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: https://github.com/jcfbeardsley/tastemp/tree/datashader
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('temp.nc')
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)
gdf.set_index('id',inplace=True)
# 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
else:
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))
else:
import plotly.express as px
fig = px.scatter_mapbox(df.dropna(), lat='latitude', lon='longitude')
# Add the datashader image as a mapbox layer image
fig.update_layout(mapbox_style="carto-darkmatter",
mapbox_layers = [
{
"sourcetype": "image",
"source": img,
"below": 'traces',
"coordinates": coordinates,
"opacity": 0.5
}]
)
fig.write_html('datashader_demo.html', auto_open=True)