К какому мультиполигону принадлежат точки долготы и широты в Python?

У меня есть долгота и широта, и ожидаемый результат состоит в том, что в каком бы мультиполигоне ни находилась точка, я получаю имя или идентификатор мультиполигона.

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, чтобы сделать это, но мне ужасно тяжело, чтобы добиться в этом прогресса. Самое близкое, что я получил, - это функция внутри и содержит, но область работы, с которой я боролся, также успешно преобразует мультиполигон в многоугольник, а затем использует силу внутри и содержит для получения желаемого результата.

Шейп-файл был загружен с сайта здесь.


python geopandas geospatial shapely fiona
person Aizaz Ali    schedule 03.09.2020    source источник

Ответы (1)


arrow_upward
2
arrow_downward

world.crs дает {'init': 'epsg:3857'} (проекция Веб-Меркатора), поэтому вам следует сначала перепроецировать ваш GeoDataFrame в проекции WGS84, если вы хотите сохранить систему координат широты-долготы вашей точки.

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

Затем вы можете использовать метод пересечений GeoPandas, чтобы найти индексы полигонов, содержащих вашу точку.

Например, для города Нью-Йорк:

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
person J-B    schedule 03.09.2020