绘图地图不显示几何图形

Question

我正在尝试使用 plotly 绘制构成不列颠哥伦比亚省一部分的多边形地理数据框。我已经使用 geopandas 绘制了 gdf，所以我知道几何形状没问题。当我尝试使用 plotly 语法使用几何列代替 json 文件来绘制地理数据框时，绘图会在浏览器中打开一个窗口，其中包含图例和颜色条，但框中没有地图。

pacificrange_CP_web.to_crs(pyproj.CRS.from_epsg(4326), inplace=True)

fig = px.choropleth(pacificrange_CP_web, 
                geojson=pacificrange_CP_web.geometry, 
                locations=pacificrange_CP_web.polyid, 
                color="protected")
fig.update_geos(fitbounds="locations", 
            visible=False)

fig.show()

据我了解，plotly 并不像 geopandas 那样支持任何 crs。所以在第一行我重新投影。

此外，我对位置提示感到不安。我认为这是一种识别等值线内单个 plygon 的方法？

我也尝试将 gdf ​​转换为 json 文件并将其用于几何图形并将两者链接起来，但这是另一个问题。

如果有人能指出我哪里出错了，将不胜感激。

Answer 1

• 使用了来自 GitHub 的示例几何体
• 很明显，这个几何图形的部分太多，无法使用 plotly 进行有效绘图
  • 创建了实用函数 reduce_geometry()，它具有三种减少几何形状的方法，即 MultiPolygon
  • 可以使用size、percentile 或topn。演示了 topn 在 MultiPolygon 中仅使用最大的 N 个几何图形
  • 此函数还具有获取其所做工作透明度的模式。 join() 将此信息放到 GeoDataFrame 上（在 hover_data 中使用）
  • MultiGeometry 仍然意味着悬停文本在它出现的地方有些奇怪。可以选择将 explode() 几何体转换为多边形
• 它不是 EPSG:4326，所以预计它可以与 plotly 一起使用

import geopandas as gpd
import shapely.geometry
import numpy as np
import plotly.express as px
import requests
from pathlib import Path
from zipfile import ZipFile
import urllib
import pandas as pd

# fmt: off
# download boundaries
url = "https://github.com/maxduso/pacificrange_CP_web/blob/85b3005c0d95e838f9e18e1e7923e90adfbba682/pacificrange_subset.zip?raw=true"

f = Path.cwd().joinpath(urllib.parse.urlparse(url).path.split("/")[-1])
# fmt: on

if False and f.exists():
    f.unlink()

if not f.exists():
    r = requests.get(url, stream=True, headers={"User-Agent": "XY"})
    with open(f, "wb") as fd:
        for chunk in r.iter_content(chunk_size=128):
            fd.write(chunk)
    zfile = ZipFile(f)
    zfile.extractall(f.stem)

# load downloaded boundaries
gdf2 = gpd.read_file(str(f.parent.joinpath(f.stem).joinpath(f"{f.stem}.shp")))

# utility function to reduce number of polygons in multipolygon
# one of following can be passed
#   size - minimum size of a polygon within multiploygon
#   percentile - for example 95, take 5% largest polygons
#   topn - take largest n polygons
def reduce_geometry(g, size=None, percentile=None, topn=None, info=False):
    if isinstance(g, shapely.geometry.Polygon):
        if info:
            return {"minarea": g.area, "polycount": 1, "kept": 1}
        else:
            return g
    if percentile:
        size = np.percentile([p.area for p in g.geoms], percentile)
    elif topn:
        topn = min(topn, len(g.geoms))
        size = sorted([p.area for p in g.geoms])[-topn]
    polys = [p for p in g.geoms if p.area >= size]
    infod = {"minarea": size, "polycount": len(g.geoms), "kept": len(polys)}

    if info:
        return infod
    if len(polys) == 1:
        return polys[0]
    elif len(polys) == 0:
        return g.geoms[np.argmax([p.area for p in g.geoms])]
    else:
        return shapely.geometry.MultiPolygon(polys)


# simplify geometry, take biggest n polygons in each multipolygon
# join info of this process onto data frame for transparency
TOPN = 20
gdf2 = gdf2.join(
    gdf2["geometry"].apply(reduce_geometry, topn=TOPN, info=True).apply(pd.Series)
)
gdf2["geometry"] = gdf2["geometry"].apply(reduce_geometry, topn=TOPN)

# optionally explode multipolygons into polygons (means hover text is better...)
EXPLODE=True
if EXPLODE:
    gdf2 = pd.merge(
        gdf2.drop(columns="geometry"),
        gdf2["geometry"].explode(index_parts=True).reset_index(),
        left_index=True,
        right_on="level_0",
    ).assign(
        source_polyid=lambda d: d["polyid"],
        polyid=lambda d: d.loc[:, ["polyid", "level_1"]]
        .astype(str)
        .apply("_".join, axis=1)
    )

# make geopandas data frame compatible with question code...
pacificrange_CP_web = (
    gdf2.to_crs("EPSG:4326")
    .set_index("polyid", drop=False)
)

fig = px.choropleth(
    pacificrange_CP_web,
    geojson=pacificrange_CP_web.geometry,
    locations=pacificrange_CP_web.polyid,
    hover_name="name_e",
    hover_data=["polycount","kept"],
    color="protected",
)
fig.update_geos(fitbounds="locations", visible=False).update_layout(
    margin={"l": 0, "r": 0, "t": 0, "b": 0}
)

mapbox choropleth

layout = dict(
    mapbox={
        "style": "carto-positron",
        "center": {
            "lon": sum(pacificrange_CP_web.total_bounds[[0, 2]]) / 2,
            "lat": sum(pacificrange_CP_web.total_bounds[[1, 3]]) / 2,
        },
        "zoom": 7,
    },
    margin={"l": 0, "r": 0, "t": 0, "b": 0},
)

px.choropleth_mapbox(
    pacificrange_CP_web,
    geojson=pacificrange_CP_web.geometry,
    locations="polyid",
    hover_name="name_e",
    hover_data=["polycount", "kept"],
    color="protected",
).update_layout(layout)