绘图地图不显示几何图形

Posted 2023-02-18

【中文标题】绘图地图不显示几何图形

【英文标题】：plotly map does not display geometries

【发布时间】：2021-12-22 09:50:27

【问题描述】：

我正在尝试使用 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 文件并将其用于几何图形并将两者链接起来，但这是另一个问题。

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

【参考方案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)

【讨论】：

哇，多么深入，非常感谢 Rob！嗯，所以我从你的 fig=... 开始，这创造了一个情节，但它看起来像是从中心点发出的带有颜色的线条。我假设这是因为几何在某些时候存在问题。你认为 QGIS 中的“修复几何”可以解决这个问题吗？然后我添加了使用更新的地理数据框设置索引的代码块，并且没有添加我已经拥有的“受保护”列。结果又回到了图形应该在的地方。

我没有使用过 QGIS，所以不能这样评论。您是否将形状文件加载到pacificrange_CP_web 中？我怀疑多边形的几何/定义有问题。 pacificrange_CP_web.boundary.plot() 看起来怎么样？显示BC地区的边界？

所以 pacificrange_CP_web 只是 BC 的一个子集。但是是的，它确实绘制了使用该线时打算使用的区域的边界。不知道这里有没有在评论中插入图片的好方法。此外，我确实发现了一些几何问题并在 QGIS 中修复了它们。但是，自从我发表上述评论以来，这样做并没有改变我的结果。为了澄清，我只使用了您的代码的以下部分，因为我相信您用来创建与我的可比较数据集的其余部分。

pacificrange_CP_web = ( pacificrange_CP_comp.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, color="protected", ) fig.update_geos(fitbounds="locations", visible=False)

两件事：1）你运行我所有的代码你得到相同的结果吗？ （即排除环境问题）。 2) 你能把你的形状数据作为 zip 文件分享到云驱动器/GitHub 上吗？ （所以我可以检查一下）