如何将 pygbif 中的 X 和 Y 坐标转换为 Cartopy 的经纬度

Question

最终地图上绘制的点（尤其是墨卡托投影）与大陆边界不匹配。这种未对准表明我转换坐标的方式可能存在错误。

我正在使用 Python 在世界地图上绘制物种的出现情况，但这些点与大陆没有正确对齐。使用

pygbif

库获取数据，并使用

mapbox_vector_tile

x 和 y 坐标在 0...512 之间进行解码。

然后使用

matplotlib

和

cartopy

进行可视化。然而，绘制的点似乎与大陆边界不一致。

from pygbif import maps
import mapbox_vector_tile

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import pandas as pd
import cartopy.feature as cfeature

x = maps.map(
    taxonKey = 3034824,
    srs='EPSG:3857', # Web Mercator
    format = ".mvt"
)

content = mapbox_vector_tile.decode(x.response.content)
extent = content['occurrence']['extent'] # 512
print('extent:', extent)

data = content['occurrence']['features']

records = [
    {
        'x': item['geometry']['coordinates'][0],
        'y': item['geometry']['coordinates'][1],
        'total': item['properties']['total']
    } for item in data
]

df = pd.DataFrame(records)

# conver x and y to mercator WGS84 Bounds: -180.0, -85.06, 180.0, 85.06
df['lon'] = (df['x'] - extent / 2) / (extent / 2) * 180
df['lat'] = (df['y'] - extent / 2) / (extent / 2) * 85.06
#df['lat'] = (df['y'] - extent / 2) / (extent / 2) * 90

# XY scatter plot
fig, ax = plt.subplots(figsize=(5, 5))
plt.scatter(df['x'], df['y'], c=df['total'], cmap=plt.cm.jet, s=0.1)
plt.xlim(0, extent)
plt.ylim(0, extent)
plt.title('xy coordinates')
plt.savefig('xy.png')
plt.show()

# Mercator scatter plot
fig, ax = plt.subplots(figsize=(5, 5))
plt.scatter(df['lon'], df['lat'], c=df['total'], cmap=plt.cm.jet, s=0.1)
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.title('mercator coordinates')
plt.savefig('mercator.png')
plt.show()

fig, ax = plt.subplots(
    figsize=(5, 5), 
    subplot_kw={'projection': ccrs.Orthographic(central_latitude=38, central_longitude=55)}
)
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.scatter(
    df['lon'], df['lat'], 
    c=df['total'], cmap=plt.cm.jet, s=1,
    transform=ccrs.PlateCarree()
)
ax.set_global()
ax.spines['geo'].set_visible(False)
plt.title('ortho coordinates')
plt.savefig('ortho.png')
plt.show()

display(df)

Answer 1

您的数据是在 epsg 3857 中提供的，而不是在 epsg 4326 中提供的！

因此，您必须将坐标转换为 epsg 3857 中的 x/y，然后再将它们重新投影到 epsg 4326 中的 lon/lat...

更改以下几行将为您提供正确的坐标：

import numpy as np

df['x_3857'] = (df['x'] - extent / 2) / (extent / 2) * np.pi * 6378137
df['y_3857'] = (df['y'] - extent / 2) / (extent / 2) * np.pi * 6378137

---

这是一个带有 eomaps 的图，表明它有效：

from eomaps import Maps

m = Maps()
m.add_feature.preset("coastline", "ocean", "land")
m.set_data(df, x="x_3857", y="y_3857", parameter="total", crs=3857)
m.set_shape.rectangles()
m.plot_map(cmap="RdYlBu_r", set_extent=False)
m.add_colorbar(hist_kwargs=dict(log=True))

Answer 2

感谢 raphael，这里是 EOmaps lib 的最终代码和结果。

from pygbif import maps
import mapbox_vector_tile

import pandas as pd
from eomaps import Maps

# Получение данных
x = maps.map(
    taxonKey=3034824,
    srs='EPSG:3857',  # Web Mercator
    format=".mvt",
    year=2022
)

content = mapbox_vector_tile.decode(x.response.content)
data = content['occurrence']['features']

records = [
    {
        'x': item['geometry']['coordinates'][0],
        'y': item['geometry']['coordinates'][1],
        'total': item['properties']['total']
    } for item in data
]
df = pd.DataFrame(records)

extent = content['occurrence']['extent']  # 512
df['x_3857'] = (df['x'] - extent / 2) / (extent / 2) * np.pi * 6378137
df['y_3857'] = (df['y'] - extent / 2) / (extent / 2) * np.pi * 6378137

# Преобразование координат Web Mercator в широту и долготу
def web_mercator_to_latlon(x, y):
    lon = (x / 6378137.0) * 180.0 / np.pi
    lat = (y / 6378137.0) * 180.0 / np.pi
    lat = 180.0 / np.pi * (2.0 * np.arctan(np.exp(lat * np.pi / 180.0)) - np.pi / 2.0)
    return lon, lat

df['lon'], df['lat'] = zip(*df.apply(lambda row: web_mercator_to_latlon(row['x_3857'], row['y_3857']), axis=1))


m = Maps(crs=Maps.CRS.Orthographic(37.6176, 55.7558), figsize=(10, 10))

m.add_feature.preset.ocean(facecolor="whitesmoke")
m.add_feature.preset.land(facecolor="gainsboro")
m.add_feature.preset.coastline(lw=0.1, edgecolor="black")

m.set_data(df, x="x_3857", y="y_3857", parameter="total", crs=3857)
m.set_shape.rectangles()
m.plot_map(cmap="RdYlBu_r", set_extent=False)

m.add_title("Распространение Борщевика в 2024 году", fontsize=26, fontweight="bold")

m.savefig("globe.png", dpi=300, bbox_inches="tight")
m.show()