Datashader如何在Plotly Mapbox中进行多边形绘制？

Question

我正在使用plotly的Scattermapbox将地图与由datashader的shade函数创建的多边形阴影图像叠加（基于https://plotly.com/python/datashader/），但是投影似乎不对齐，如下图所示。有什么建议可以使用plotly的Scattermapbox和datashader解决这个问题吗？

可重现的示例：

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)

带阴影的图片

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

覆盖地图

我看到Scattermapbox只支持Mercator投影，这让我感到困惑，因为plotly文档中的示例似乎是以经纬度格式给出的。但我尝试将GeoDataFrame坐标转换为epsg 3857，具体操作请参考

# world = world.to_crs(epsg=3857)

结果是阴影图像变得不可见。非常感谢您的帮助。

- Manuel Diehn

---

非常有趣 - Adam

我已经从多个方面进行了调查...我怀疑坐标可能有误，但微调几乎没有什么影响。我认为由Datashader创建的xarray可能是错误的 - 但是如果我重新创建一个pandas数据帧并执行密度mapbox，则一切正常。所以现在陷入了死胡同... - Rob Raymond

我建议使用基于HoloViews的指南来使用Datashader与plotly，网址为https://dash.plotly.com/holoviews。但是，如果这些指南不包括多边形或者您只想手动完成，您可以使用datashader.utils.lnglat_to_meters函数将原始坐标转换为米（请参见https://dev59.com/qazla4cB1Zd3GeqPE_K1#51385389）。 - James A. Bednar

@JamesA.Bednar 感谢您的建议。我想使用plotly的图形对象库中的Scattermapbox，原因是可以在“update_layout”调用中使用选项“'below': 'water'”添加层。我没有找到类似的东西在holoviews中。 - Manuel Diehn

我不确定Scattermapbox具体支持什么，但如果您在HoloViews中有绘图a、b和c，只需执行a*b*c即可显示c叠加在b上，再叠加在a上。 - James A. Bednar

Answer 1

我们已经发现了这个问题的解决方法：以下是每个步骤/函数代码和描述：

引用以供参考：

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：使用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

以下序列的所有操作：

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

然后，使用datashader生成的此图像可以通过plotly对象层添加到plotly绘图中，并将该层提供给图形。

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

- Adam

Answer 2

你尝试过使用 epsg:4326 吗？在我的情况下，我使用这个并且几何图形被正确地放置了。

另一方面，使用 geopandas 转换数据框的几何列时，你必须使用参数 "inplace=True"。

- Francisco Javier López Andreu

---

请使用注释来提问。 - Tobi208