I like using Plotly to create interactive plots and maps straight from Python. I find it more intuitive than Bokeh’s API. When combined with Plotly’s sister library Dash, the pair can be very powerful in creating both interactive and dynamic web based tools and visualizations.
While Bokeh has some native Datashader integration, including in the mapping components, the Plotly Mapbox mapping backend relies on Mapbox’s own image overlay capability to display arbitrary Datashader outputs. Plotly has very helpful documentation and includes this use case as an example.
That example code is shown here along with the output map which looks very cool.
import pandas as pd
df = pd.read_csv('https://raw.githubusercontent.com/plotly/datasets/master/uber-rides-data1.csv')
dff = df.query('Lat < 40.82').query('Lat > 40.70').query('Lon > -74.02').query('Lon < -73.91')
import datashader as ds
cvs = ds.Canvas(plot_width=1000, plot_height=1000)
agg = cvs.points(dff, x='Lon', y='Lat')
# agg is an xarray object, see http://xarray.pydata.org/en/stable/ for more details
coords_lat, coords_lon = agg.coords['Lat'].values, agg.coords['Lon'].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]]]
from colorcet import fire
import datashader.transfer_functions as tf
img = tf.shade(agg, cmap=fire)[::-1].to_pil()
import plotly.express as px
# Trick to create rapidly a figure with mapbox axes
fig = px.scatter_mapbox(dff[:1], lat='Lat', lon='Lon', zoom=12)
# Add the datashader image as a mapbox layer image
fig.update_layout(mapbox_style="carto-darkmatter",
mapbox_layers = [
{
"sourcetype": "image",
"source": img,
"coordinates": coordinates
}]
)
fig.show()
Although it is not particularly obvious on this scale looking at Manhattan, there is actually a projection problem if you are providing a raster that was “datashaded” in Lat/Lon space to the Mapbox image overlay since Mapbox only supports Web Mercator projections.
Mapbox supports the popular Web Mercator projection, and does not support any other projections. Web Mercator is a nearly conformal projection that is adopted by the vast majority of web maps and its use allows you to combine Mapbox’s maps with other layers in the same projection.
Commonly this projection is referred to as EPSG:900913 or EPSG:3857.
The problem becomes much more apparent when trying to plot continuous fields on larger scales. Here is some code that demonstrates the projection problem with this original approach.
import xarray as xr
import colorcet
import datashader as ds
import datashader.transfer_functions as tf
import plotly.graph_objects as go
cvs = ds.Canvas(plot_height=2000, plot_width=2000)
def to_datashader_raster(data, y="lat", x="lon"):
z = xr.DataArray(data.values, coords=[("y", data[y].values), ("x", data[x].values)])
return cvs.raster(z)
def to_shade(agg, cmap=colorcet.bmy):
return tf.shade(agg, cmap=cmap)
def to_pil(im):
return im[::-1].to_pil()
up = xr.open_dataset("Ice/201904_ice.nc")
upi = data.isel(time=0).sel(lon=slice(-170,-73))
z = upi["siconc"]
agg = to_datashader_raster(z, x="lon", y="lat")
im = to_shade(agg)
img = to_pil(im)
coords_lat, coords_lon = (z.lat.values, z.lon.values)
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]],
]
fig = go.Figure(go.Scattermapbox())
fig.update_layout(
mapbox_style="open-street-map",
mapbox_layers=[
{
"sourcetype": "image",
"source": img,
"coordinates": coordinates,
},
],
margin={"r": 0, "t": 0, "l": 0, "b": 0},
)
fig.show()
It is clear from the coastline misalignment that there is a projection problem, yet the corners of the plotted ice concentration data still correspond to the “coordinates” input parameter. One solution to this problem is to project the Lat/Lon data into Web Mercator before passing along to Datashader. Even for regular evenly spaced grids this necessarily forces the data to a 2D meshgrid in the projected space, and as a result the Datashader quadmesh transform must be used.
from pyproj import Transformer
import numpy as np
gcs_to_3857 = Transformer.from_crs(4326, 3857, always_xy=True)
def add_psuedomercator(data):
if len(data.lon.shape) < 2:
lon, lat = np.meshgrid(data.lon, data.lat)
else:
lon, lat = data.lon.values, data.lat.values
shp = lon.shape
x, y = gcs_to_3857.transform(list(lon.ravel()), list(lat.ravel()))
data["lon2"] = xr.DataArray(np.reshape(x, shp), dims=("lat", "lon"))
data["lat2"] = xr.DataArray(np.reshape(y, shp), dims=("lat", "lon"))
return up
def to_datashader_quadmesh(data, y="lat", x="lon"):
z = xr.DataArray(
data.values,
dims=["y", "x"],
coords={"Qy": (["y", "x"], data[y].values), "Qx": (["y", "x"], data[x].values)},
name="z",
)
return cvs.quadmesh(z, x="Qx", y="Qy")
By combining this projection step to the original methodology, the following example does correctly overlay the raster image layer onto the Web Mercator Plotly/Mapbox plot.
import xarray as xr
import colorcet
import datashader as ds
import datashader.transfer_functions as tf
import numpy as np
import plotly.graph_objects as go
from pyproj import Transformer
import numpy as np
cvs = ds.Canvas(plot_height=2000, plot_width=2000)
gcs_to_3857 = Transformer.from_crs(4326, 3857, always_xy=True)
def to_datashader_quadmesh(data, y="lat", x="lon"):
z = xr.DataArray(
data.values,
dims=["y", "x"],
coords={"Qy": (["y", "x"], data[y].values), "Qx": (["y", "x"], data[x].values)},
name="z",
)
return cvs.quadmesh(z, x="Qx", y="Qy")
def to_shade(agg, cmap=colorcet.bmy):
return tf.shade(agg, cmap=cmap)
def to_pil(im):
return im[::-1].to_pil()
def add_psuedomercator(data):
if len(data.lon.shape) < 2:
lon, lat = np.meshgrid(data.lon, data.lat)
else:
lon, lat = data.lon.values, data.lat.values
shp = lon.shape
x, y = gcs_to_3857.transform(list(lon.ravel()), list(lat.ravel()))
try:
data["lon2"] = xr.DataArray(np.reshape(x, shp), dims=("lat", "lon"))
data["lat2"] = xr.DataArray(np.reshape(y, shp), dims=("lat", "lon"))
except:
data["lon2"] = xr.DataArray(np.reshape(x, shp), dims=("latitude", "longitude"))
data["lat2"] = xr.DataArray(np.reshape(y, shp), dims=("latitude", "longitude"))
return up
# Read gridded NetCDF data using Xarray
up = xr.open_dataset("Ice/201904_ice.nc")
upi = data.isel(time=0).sel(lon=slice(-170,-73))
z = upi["siconc"]
# Add Web Mercator Projection to DataArray
z = add_psuedomercator(z)
# Create a Datashader quadmesh DataArray using the projected coords
agg = to_datashader_quadmesh(z, x="lon2", y="lat2")
im = to_shade(agg)
img = to_pil(im)
# TODO Maybe use the corners of the image coordinates rather than the grid point
# coordinates of the corner points?
coords_lat, coords_lon = (z.lat.values, z.lon.values)
if len(coords_lat.shape) > 1:
coords_lat = coords_lat[:, 0]
coords_lon = coords_lon[0, :]
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]],
]
fig = go.Figure(go.Scattermapbox())
fig.update_layout(
mapbox_style="open-street-map",
mapbox_layers=[
{
"sourcetype": "image",
"source": img,
"coordinates": coordinates,
},
],
margin={"r": 0, "t": 0, "l": 0, "b": 0},
)
fig.show()
Useful resources: