Geopandas 叠加层(交叉点)Returns 零行

Geopandas Overlay (Intersection) Returns Zero Rows

我创建了一个带有 Points geometry 列的单行 GeoDataFrame,如下所示:

1。创建df

df = pd.DataFrame([[51.502687, -3.538329, 2242, 1, 47]],
                   columns = ["lat","long","num_of_trucks","performance","num_of_routes"]
                 )

df 

2。从 df

创建 gdf
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lat"],df["long"]),
                                                                 crs={"init": "epsg:4326"})
gdf

3。重新投影到使用米的 CRS

gdf.to_crs(epsg=3395,inplace=True)

4。计算 Point 周围的 10KM 缓冲区半径,并将 Point 的原始 geometry 列替换为新的缓冲区半径 Polygon

gdf["geometry"] = gdf.geometry.buffer(10000)

上述过程生成了我想要的单行 GeoDataFrame gdf,其中包含一个包含单个多边形的新 geometry 列。

我有第二个 GeoDataFrame,名为 map_selection。第二个 GeoDataFrame 由这个 shapefile 中的几行组成:https://geoportal.statistics.gov.uk/datasets/local-authority-districts-may-2020-boundaries-uk-buc。 此 GeoDataFrame 的每一行代表一个英国地方当局区(作为其 geometry 列中的 Polygon/Multigon)。请注意,此 shapefile 最初使用 CRS = EPSG 27700,但我在加载它时将其更改为 CRS = ESPG 4326,因为我需要在 Geoviews 中绘制多边形。

考虑下图:

新生成的gdfgeometry列中包含的Polygon就是红圈。 每个蓝色多边形代表 map_selection geometry 列中包含的地方当局区域 Polygons/Multigons。

我正在尝试计算 map_selection 中每个蓝色多边形的面积,它们加在一起构成 gdf 中的单个红色多边形。请注意,红色多边形不会与所有蓝色多边形相交,但肯定会与一些蓝色多边形相交。 我的最终目标是根据适当的蓝色多边形的数量来估计红色多边形的数量。

我认为下面的内容会让我得到多边形之间的交集,但是它 returns 零行:

5。重新投影 gdf CRS 以匹配 map_selection:

gdf.to_crs(epsg=4326,inplace=True)

6. Return gdfmap_selection GeoDataFrames 之间的交集:

gpd.overlay(map_selection,gdf,how="intersection")

输出为零行,以下项的输出也是如此:

gpd.overlay(gdf,map_selection,how="intersection")

如有任何想法或建议,我将不胜感激。

谢谢

您很可能有不同的预测。您的管理边界在 EPSG:27700.

为确保您 re-projecting 到其他数据帧的 CRS,请始终直接传递该 CRS。

gdf.to_crs(map_selection.crs, inplace=True)

所以我按照@martinfleis 的建议绘制了我的 2 个 GeoDataFrames,这表明我的单个多边形的纬度和经度值是倒置的。奇怪的是,当您查看 GeoDataFrame 的行时,它们并没有以这种方式出现,只有在绘制时才会出现,这就是我没有选择它的原因。

原来我的代码 Step 2 应该是:

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["long"],df["lat"]),
                                                                 crs={"init": "epsg:4326"})

而不是

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["lat"],df["long"]),
                                                                 crs={"init": "epsg:4326"})

重申一下 - 如果您在 运行 之后查看 gdf GeoDataFrame,两种情况下的纬度和经度值将相同。只有当我在 Geoviews 中绘制这两个 GeoDataFrames 的叠加层时,纬度和经度的反转才变得明显。

我的覆盖现在成功 returns 预期的交叉行。