pythongisgeospatialcartopycartography

How to transform X and Y coord from pygbif in lat, lon for Cartopy

enter image description here

The points plotted on the final map (especially the Mercator projection) do not match the continental boundaries. This misalignment suggests that there might be an error in how I'm transforming the coordinates. enter image description here

I'm using Python to plot occurrences of a species on a world map, but the points do not align correctly with the continents. The data is fetched using the pygbif library, decoded with mapbox_vector_tile x and y coords between 0...512.

And then visualized using matplotlib and cartopy. However, the plotted points appear misaligned with the continental boundaries.

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)
Solution

  • Your data is provided in epsg 3857 not in epsg 4326!

    So you have to convert the coordinates to x/y in epsg 3857 before re-projecting them to lon/lat in epsg 4326...

    Changing the following lines will give you the proper coordinates:

    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

    Here's a plot with eomaps to show that it worked:

    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))

    enter image description here