Datashader and Mapbox - Projection issues?

Hi team,

Using the Datashader example here, I’ve managed to write a small script that will extract temperature from a given NetCDF file and create Datashader polygons for plotting on a Mapbox (see below). I’ve opted for Datashader polygons due to the volume of points and performance limitations with Choroplethmapbox.

The script runs successfully, but I’m finding issues with the projected location of the Datashader image when rendered. The original dataset is in EPSG 4326, and when plotted using Scattermapbox or Choroplethmapbox, the points are in the correct locations (i.e. the ocean relative to the coastline). However, the Datashader image is suffering distortions in part of the image, as you’ll see from the mismatch between the Datashader and Scattermapbox results below. The coordinates appear to line up okay in the south-west and north-east corners of the grid, but are distorted in the north-west, south-east and central parts (differing by approx 3 grid cells in the Y direction). The Datashader polygons should be overlapping with all those points shown from Scattermapbox.

Given that there is nothing wrong with the geojson in the supplied polygon boundaries (as evidenced by plotting Choroplethmapbox in the correct location), I’m trying to get my head around what is happening in Datashader to create this distortion. Can someone please explain?

I read in this post that all Datashader data should be in EPSG 3857 for plotting with Mapbox, but in my instance that results in an even larger discrepancy. Is this simply a case of a different projection mismatch, or is something bigger going on here?

Thanks in advance.

Full script and datasets available for download here:

import xarray as xr
import json
import geopandas as gpd
import spatialpandas as spd
from colorcet import fire
import datashader.transfer_functions as tf

'''Simple DataShader Projection Issue Example'''

# Set whether to reproject the data to web mercator (from EPSG:4326) before running datashader:
reprojectToWebMercator = False

# Set whether to generate datashader polygons. If false, generates points instead:
datashaderPolygons = True

# Set whether to overlay Mapbox Choropleth alongside Datashader result. If false, renders Mapbox Scatter alongside instead:
mapboxChoropleth = False

# Read the target dataset
nc = xr.open_dataset('')
df = nc['temp'][0, -1, 1:-1, 1:-1].to_dataframe()
df.index = ['c' + ('_r'.join(map(str, tuple(v+1 for v in x)))) for x in df.index.values]

# Read in polygon boundaries:
with open('cells.geojson') as json_file:
    cells = json.load(json_file)
gdf = gpd.read_file('cells.geojson')

if reprojectToWebMercator:
    gdf = gdf.to_crs(epsg=3857)

# Merge with data values from NetCDF DF based on index:
gdf = gdf.join(df)

# Create datashader polygons:
import datashader as ds
cvs = ds.Canvas(plot_width=2000, plot_height=2000)
if datashaderPolygons:
    agg = cvs.polygons(spd.GeoDataFrame(gdf), geometry='geometry',agg=ds.mean('temp'))
    coords_lat, coords_lon = agg.coords['y'].values, agg.coords['x'].values
    agg = cvs.points(df, x='longitude', y='latitude',agg=ds.mean('temp'))
    coords_lat, coords_lon = agg.coords['latitude'].values, agg.coords['longitude'].values

# Corners of the image, which need to be passed to mapbox
coordinates = [[coords_lon[0], coords_lat[0]],
               [coords_lon[-1], coords_lat[0]],
               [coords_lon[-1], coords_lat[-1]],
               [coords_lon[0], coords_lat[-1]]]

img = tf.shade(agg, cmap=fire)[::-1].to_pil()

if mapboxChoropleth:
    import plotly.graph_objects as go
    fig = go.Figure(go.Choroplethmapbox(geojson=cells, locations=df.cellid))
    import as px
    fig = px.scatter_mapbox(df.dropna(), lat='latitude', lon='longitude')

# Add the datashader image as a mapbox layer image
                 mapbox_layers = [
                    "sourcetype": "image",
                    "source": img,
                    "below": 'traces',
                    "coordinates": coordinates,
                    "opacity": 0.5
fig.write_html('datashader_demo.html', auto_open=True)
1 Like


Was this issue ever resolved?

I’m facing something very similar!

It is normal that you don’t see a difference with point location and you see it with a raster. The webmercator assumes a spherical Earth whereas most of the lat lon datasets are in WGS84 where your datum is an obloid Earth. Why plotly mapbox thought that it was a good idea to not offer to put the coordinates in the native projection system of mapbox and stick to lat lon coordinates does not make much sense.
Basically you can warp your dataset so that the curvature of your dataset is matching the one of a sphere.
This page has a few minor bugs but you surely will find your way

The other solution is to use something like xarray interp but based on their example, I could not really figure out what to do.