Python中的经纬度点属于哪个多面体？

Question

我有经度和纬度，预期结果是无论点在哪个多边形中，我都会得到多边形的名称或 ID。

import geopandas as gpd

world = gpd.read_file('/Control_Areas.shp')
world.plot()

输出

   0    MULTIPOLYGON (((-9837042.000 6137048.000, -983...
   1    MULTIPOLYGON (((-11583146.000 5695095.000, -11...
   2    MULTIPOLYGON (((-8542840.287 4154568.013, -854...
   3    MULTIPOLYGON (((-10822667.912 2996855.452, -10...
   4    MULTIPOLYGON (((-13050304.061 3865631.027, -13.

之前的尝试： 我已经尝试过 fiona、shapely 和 geopandas 来完成这项工作，但我一直在努力取得进展。我得到的最接近的是 within 和 contains 函数，但我一直在努力的工作领域是成功地将多面体转换为多边形，然后利用 within 和 contains 的力量来获得所需的输出。

形状文件已从 here 下载。

Answer 1

world.crs 给出 {'init': 'epsg:3857'}（Web 墨卡托投影），因此如果您想保留点的 latitude-longitude 坐标系，您应该首先在 WGS84 投影中重新投影 GeoDataFrame。

world = world.to_crs("EPSG:4326")

然后你可以使用GeoPandas的intersects方法来找到包含你的点的多边形的索引。

例如纽约市：

from shapely.geometry import Point
NY_pnt = Point(40.712784, -74.005941)
world[["ID","NAME"]][world.intersects(NY_pnt)]

这导致：

ID  NAME
20  13501   NEW YORK INDEPENDENT SYSTEM OPERATOR

您可以使用 shapely within 方法检查结果：

NY_pnt.within(world["geometry"][20])

如果有多个点，可以创建一个GeoDataFrame并使用sjoin方法：

NY_pnt = Point(40.712784, -74.005941)
LA_pnt = Point(34.052235, -118.243683)
points_df = gpd.GeoDataFrame({'geometry': [NY_pnt, LA_pnt]}, crs='EPSG:4326')
results = gpd.sjoin(points_df, world, op='within')
results[['ID', 'NAME']]

输出：

ID  NAME
0   13501   NEW YORK INDEPENDENT SYSTEM OPERATOR
1   11208   LOS ANGELES DEPARTMENT OF WATER AND POWER

Python中的经纬度点属于哪个多面体？

Which multipolygon does the point of longtitude and latitude belong to in Python?

python

geospatial

shapely

fiona

geopandas