Datashader integration for polygons in plotly mapbox

Question

I'm using plotly's Scattermapbox to overlay a map with a shaded image of polygons created by datashader's shade function (based on https://plotly.com/python/datashader/), but the projections do not seem to align, see picture below. Any suggestions how I can overcome this problem using plotly's Scattermapbox and datashader?

Reproducible example:

import geopandas as gpd
import plotly.graph_objects as go
import spatialpandas as spd
import datashader as ds
from colorcet import fire
import datashader.transfer_functions as tf

# load data
world = gpd.read_file(
    gpd.datasets.get_path('naturalearth_lowres')
)
# world = world.to_crs(epsg=3857)
# create spatialpandas DataFrame
df_world = spd.GeoDataFrame(world)
# create datashader canvas and aggregate
cvs = ds.Canvas(plot_width=1000, plot_height=1000)
agg = cvs.polygons(df_world, geometry='geometry', agg=ds.mean('pop_est'))
# create shaded image
tf.shade(agg, cmap=fire)

shaded image

# create shaded image and convert to Python image
img = tf.shade(agg, cmap=fire)[::-1].to_pil()

coords_lat, coords_lon = agg.coords["y"].values, agg.coords["x"].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]],
]

fig = go.Figure(go.Scattermapbox())
fig.update_layout(
    mapbox_style="open-street-map",
    mapbox_layers=[
        {
            "sourcetype": "image",
            "source": img,
            "coordinates": coordinates,
        }
    ]
)
fig.show()

overlayed map

I read that Scattermapbox only supports Mercator projection which I found confusing as the examples in plotly's documentation seem to be in long/lat format, but I tried converting the coordinates of the GeoDataFrame to epsg 3857, see

# world = world.to_crs(epsg=3857)

The results is that the shaded image becomes invisible. Any help would be highly appreciated.

Answer 1

Have you tried with epsg:4326? In my case, I use this one and the geometries are placed correctly.

On the other hand, with geopandas to convert the geometry column of the dataframe you have to use the parameter "inplace=True".

Answer 2

We have discovered the solution to this issue: Below is each step / function code and description:

Imports for reference :

import datashader as ds
import datashader.transfer_functions as tf
import numpy as np
import pandas as pd
import plotly.graph_objects as go
import rasterio
import shapely.geometry
import xarray as xr

_helper_add_pseudomercator_optimized : Creates array from the meshgrid with the proper mercator coordinates from the original raster with epsg:4326.

def _helper_add_pseudomercator_optimized(raster):
    """Adds mercator coordinates epsg:3857 from a raster with epsg:4326.

    Originally defined as `add_psuedomercator_adam_manuel_optimized`

    Args:
        raster: xr.DataArray: `xr.DataArray` to transform coordinates

    Returns:
        `xr.DataArray` with coordinates (x, y) transformed from epsg:4326 to epsg:3857
    """
    # Transformer that converts coordinates from epsg 4326 to 3857
    gcs_to_3857 = Transformer.from_crs(4326, 3857, always_xy=True)

    x_vals = list(raster.x.values.squeeze())  # x values from the raster dimension x
    y_vals = list(raster.y.values.squeeze())  # x values from the raster dimension x

    # Allows transformation of non-square coordinates
    y_dummy_vals = [raster.y.values[0] for v in raster.x.values]  # dummy values
    x_dummy_vals = [raster.x.values[0] for v in raster.y.values]  # dummy values

    x, _ = gcs_to_3857.transform(x_vals, y_dummy_vals)  # Obtain x output here only
    _, y = gcs_to_3857.transform(x_dummy_vals, y_vals)  # Obtain y output here only\

    # Create meshgrid with the x and y mercator converted coordinates
    lon, lat = np.meshgrid(x, y)

    # Add meshgrid to raster -> raster now has mercator coordinates for every point
    raster["x_mercator"] = xr.DataArray(lon, dims=("y", "x"))
    raster["y_mercator"] = xr.DataArray(lat, dims=("y", "x"))

    return raster

def _helper_affine_transform(raster):
    """Create new affine from a raster. Used to get new affine from the transformed affine.

    Args:
        raster: xr.DataArray: `xr.DataArray` to get the original affine and then transform

    Returns:
        New affine transform for a coarsened array
    """
    res = (raster.x[-1].values - raster.x[0].values) / raster.x.shape[0]
    scale = Affine.scale(res, -res)

    transform = (
        Affine.translation(raster.x[0].values - res / 2, raster.y[0].values - res / 2)
        * scale
    )

    return transform

def _helper_to_datashader_quadmesh(raster, y="lat", x="lon"):
    """Create lower level quadmesh with data based on flood raster. Map Flooding
    to lower level map.

    Args:
        raster: xr.DataArray: `xr.DataArray` raster of flooded regions

    Returns:
        `datashader.Canvas` based on quadmesh from original flood raster
    """
    cvs = ds.Canvas(plot_height=5000, plot_width=5000)

    z = xr.DataArray(
        raster.values.squeeze(),
        dims=["y", "x"],
        coords={
            "Qy": (["y", "x"], raster[y].values),
            "Qx": (["y", "x"], raster[x].values),
        },
        name="z",
    )

    return cvs.quadmesh(z, x="Qx", y="Qy")

def _helper_img_coordinates(raster):
    """Get coordinates of the corners of the baseline raster.

    Args:
        raster: xr.DataArray: `xr.DataArray` to get corner coordinates from

    Returns:
        coordinates of where to plot the flooded raster on the map
    """
    coords_lat, coords_lon = (raster.y.values, raster.x.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]],
    ]

    return coordinates

All operations together for the below sequence :

# Add mercator coordinates to the raster
raster = _helper_add_pseudomercator_optimized(raster)

# Create quadmesh from the burned raster
agg_mesh = _helper_to_datashader_quadmesh(raster, x="x_mercator", y="y_mercator")

# Don't plot values where the flooding is zero
agg_mesh = agg_mesh.where(agg_mesh < 0)

# Convert to datashader shade
im = tf.shade(agg_mesh, Theme.color_scale)

# Convert to image
img = im.to_pil()

# Get coordinates to plot raster on map
coordinates = _helper_img_coordinates(baseline_raster)

Then this image produced by datashader can be added to a plotly plot using the plotly objects layer, and providing this layer to the figure

layer = go.layout.mapbox.Layer(
        below="water",
        coordinates=coordinates,
        sourcetype="image",
        source=img,
    )