Python: как отфильтровать точки только в определенной границе с помощью геопанд? - PullRequest
Новая
       61

Python: как отфильтровать точки только в определенной границе с помощью геопанд?

0 голосов
/ 19 апреля 2019

У меня есть шейп-файл с некоторыми точками. 

import geopandas as gpd
from geopandas.tools import sjoin
import osmnx as ox

myShape = gpd.read_file('myShape.shp')
myShape.head(3)

geometry
0   POINT (-72.09513801999077 18.6526410972226)
1   POINT (-72.21044508038457 19.61039786418674)
2   POINT (-72.27903674995586 18.52939294725028)

Затем я извлек город из открытой карты улиц и ее границ.Я хочу сохранить только те точки, которые находятся внутри города.

Вот что я делаю: 

city = ox.gdf_from_place('Paris', which_result=2)
city = gpd.GeoDataFrame(city)

myShapeTmp = myShape.copy()
for i in myShape.index:
    if (myShape['geometry'][i] in gdf['geometry']) == False:
        myShapeTmp = myShapeTmp.drop([i], axis=0)

Но это занимает вечность.Есть ли другой способ сделать это?

Ответы [ 2 ]

1 голос
/ 19 апреля 2019

Я решил таким образом 

city = ox.gdf_from_place('Paris', which_result = 2)
city = gpd.GeoDataFrame(city)
boundary = city.ix[0].geometry
myShapeTmp = school[myShape.geometry.within(boundary)]
0 голосов
/ 19 апреля 2019

Другой способ получить points in polygon может быть таким:

Давайте создадим файл port_au_princ.csv со следующим содержимым: 

index,geometry
0,POINT(-72.09513801999077 18.6526410972226)
1,POINT(-72.21044508038457 19.61039786418674)
2,POINT(-72.27903674995586 18.52939294725028)

Вот код, которыйвыполняет все необходимые операции. 

import pandas as pd
import geopandas as gpd
from geopandas.tools import sjoin
import osmnx as ox
from shapely.wkt import loads

# get polygon of the city
city_pgon = ox.gdf_from_place('Port-au-Prince', which_result = 2)
# create a geoDataFrame out of it, set appropriate CRS
city_pgon = gpd.GeoDataFrame(city_pgon, crs={'init': 'epsg:4326'})
# read the file
pap = pd.read_csv('port_au_princ.csv', sep=',')
geometry = [loads(pt) for pt in pap['geometry']]  # prep geometry
# create geoDataFrame from `pap` dataframe
geo_pap = gpd.GeoDataFrame(pap, geometry=geometry, crs={'init': 'epsg:4326'})
pts_in_poly = gpd.sjoin(geo_pap, city_pgon, op='within', how='inner') # do spatial join

ax1 = city_pgon.plot()   # plot the polygon
geo_pap.plot(ax=ax1, color='red', zorder=5)         # plot red points
pts_in_poly.plot(ax=ax1, color='yellow', zorder=6)  # plot yellow point

Сюжет будет:

point-in-poly

Распечатайте некоторую информацию pts_in_poly сэтот код: 

print(pts_in_poly.index, pts_in_poly.geometry)

выход, 

Int64Index([2], dtype='int64') 2    POINT (-72.27903674995586 18.52939294725028)
Name: geometry, dtype: object
...