如何将pygbif的X/Y坐标转为Cartopy经纬度？坐标匹配异常

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

我正在使用Python在世界地图上绘制某物种的分布点，但这些点与大陆无法正确对齐。数据通过pygbif库获取，经mapbox_vector_tile解码后得到0...512范围内的X、Y坐标。随后我使用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)

# 将x、y转换为墨卡托WGS84范围：-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散点图
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()

# 墨卡托散点图
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)

你当前的坐标转换逻辑错误，因为Mapbox矢量瓦片的0-512像素坐标不能直接线性映射到经纬度范围，它们对应的是Web Mercator（EPSG:3857）投影下的瓦片像素坐标，需要先转换为EPSG:3857的空间坐标，再转成WGS84经纬度（EPSG:4326）。

修正步骤： #

1. 先将0-512的瓦片像素坐标转换为EPSG:3857空间坐标，EPSG:3857的世界范围是[-20037508.34, 20037508.34]（x和y方向）。
2. 再将EPSG:3857坐标转换为WGS84经纬度。

修改后的代码： #

需要先安装pyproj库（用于坐标转换）：

pip install pyproj

替换原代码中的坐标转换部分，完整代码如下：

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
from pyproj import Transformer

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)

# 1. 将0-512的瓦片像素坐标转换为EPSG:3857空间坐标
web_mercator_min = -20037508.342789244
web_mercator_max = 20037508.342789244
df['x_3857'] = (df['x'] / extent) * (web_mercator_max - web_mercator_min) + web_mercator_min
df['y_3857'] = (df['y'] / extent) * (web_mercator_max - web_mercator_min) + web_mercator_min

# 2. 将EPSG:3857转换为WGS84经纬度（EPSG:4326）
transformer = Transformer.from_crs("EPSG:3857", "EPSG:4326", always_xy=True)
df['lon'], df['lat'] = transformer.transform(df['x_3857'], df['y_3857'])

# XY散点图
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坐标')
plt.savefig('xy.png')
plt.show()

# 墨卡托投影散点图
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('墨卡托坐标')
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('正射投影坐标')
plt.savefig('ortho.png')
plt.show()

display(df)

关键说明： #

• 你之前直接将像素坐标线性映射到经纬度范围，忽略了Web Mercator投影的非线性特性（尤其是纬度方向的拉伸），导致点的位置错位。
• 使用pyproj进行严格的坐标转换，能保证坐标系统的正确性，让分布点和大陆边界准确对齐。

内容的提问来源于stack exchange，提问作者Maxim