时间戳

首页 > 解决方案 > Python中的经纬度点属于哪个多面体？

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

问题描述

我有经度和纬度，预期的结果是无论该点在哪个多面体中，我都会得到多面体的名称或 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 来完成这项工作，但我一直在努力争取在这方面取得进展。我得到的最接近的是 inside 和 contains 函数，但我一直在努力的工作领域是成功地将多面体转换为多边形，然后利用 inside 和 contains 的力量来获得所需的输出。

shapefile 已从此处下载。

标签： pythongeospatialgeopandasshapelyfiona

解决方案

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

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

推荐阅读