经度和纬度点在Python中属于哪个多边形？

编程问答 2022-08-12

问题描述

我有经度和纬度，并且预期的结果是无论该点位于哪个多边形中，我都可以得到该多边形的名称或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，身材匀称的大猩猩和大熊猫来完成该任务，但我为此竭尽全力。我得到的最接近的是“ inside”和“ contains”功能，但是我一直在努力的工作领域是将多面多边形成功地转换为多边形，然后利用inside和contains的功能获得所需的输出。

shapefile已从here下载。

解决方法

world.crs给出{'init': 'epsg:3857'}（Web Mercator投影），因此如果要保留点的经纬度坐标系，则应首先在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

您可以在方法内整齐地检查结果：

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
fionageopandasgeospatialpythonshapely