绘图地图框中多边形的数据着色器集成

Question

我正在使用 plotly 的 Scattermapbox 用由 datashader 的阴影功能（基于 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)

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

我读到 Scattermapbox 仅支持墨卡托投影，我发现这很令人困惑，因为 plotly 文档中的示例似乎是 long/lat 格式，但我尝试将 GeoDataFrame 的坐标转换为 epsg 3857，请参阅

# world = world.to_crs(epsg=3857)

结果是阴影图像变得不可见。任何帮助将不胜感激。

Answer 1

您是否尝试过使用 epsg:4326？就我而言，我使用了这个，并且几何图形放置正确。

另一方面，使用 geopandas 来转换数据框的几何列，您必须使用参数“inplace=True”。

Answer 2

我们已经找到了解决这个问题的方法：下面是每个步骤/功能代码和描述：

进口参考：

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)

然后可以使用 plotly objects 层将 datashader 生成的图像添加到 plotly plot 中，并将该层提供给图形

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