GeoPandas no identifica los puntos dentro del polígono

Question

Estoy tratando de intersecar algunos puntos con un polígono en GeoPandas.

Primero creo un Point GeoDataFrame a partir de un archivo de texto .

import geopandas as gp
import numpy as np
import pandas as pd
from shapely.geometry import Point, Polygon

tile_index = pd.read_csv('../tile_index.dat', names=['tile', 'x', 'y'], sep=' ')
geometry = [Point(r.x, r.y) for ix, r in tile_index.iterrows()]
tile_index = gp.GeoDataFrame(tile_index, geometry=geometry)
tile_index = tile_index[(tile_index.x.between(-10, 30)) & (tile_index.y.between(-10, 30))]
print tile_index.head()

Que produce:

     tile      x       y                                        geometry
202   226 -9.226  -0.724  POINT (-9.226000000000001 -0.7240000000000001)
203   227 -9.226   9.276                POINT (-9.226000000000001 9.276)
204   228 -9.226  19.276               POINT (-9.226000000000001 19.276)
205   229 -9.226  29.276               POINT (-9.226000000000001 29.276)
219   243  0.774  -0.724               POINT (0.774 -0.7240000000000001)

A continuación, creo un Polygon para que se crucen:

poly = Polygon([(0, 0), (0, 20), (20, 20), (20, 0)])
P = gp.GeoDataFrame(data=[1, poly]).T
P.columns = ['ID', 'geometry']
print np.any(tile_index.within(P))

Que devuelve False pero al menos 4 puntos deben estar dentro.

Si creo un punto utilizando las coordenadas del punto 237 sí se interseca.

pt = Point(10.774, 9.276)
p = gp.GeoDataFrame([[1, pt]])
p.columns = ['ID', 'geometry']
print p.within(P)

Que devuelve True .

Si lo hago P.contains(tile_index.geometry) entonces todos son False pero si lo hago np.any([P.contains(row.geometry)[0] for ix, row in tile_index.iterrows()]) entonces esto devuelve True .

Answer 1

geopandas alinea los fotogramas en el índice y, a continuación, realiza within o contains por elementos. Hay una descripción más completa en este número cuya parte más relevante es el ejemplo que se cita aquí:

>>> p1 = shapely.geometry.Polygon([(0,0), (1,0), (1,1), (0,1)])

>>> p2 = shapely.geometry.Polygon([(1,1), (2,1), (2,2), (1,2)])  
>>> p1.intersects(p2)
True
>>> s1 = geopandas.GeoSeries([p1])
>>> s2 = geopandas.GeoSeries([p2])
>>> s1.intersects(s2)
0    True
dtype: bool

# giving it another index, s1 and s2 will now first align (and in this case return False
# for each index, as there are no polygons that intersect and have matching index)
>>> s2 = geopandas.GeoSeries([p2], index=[1])
>>> s1.intersects(s2)
0    False
1    False
dtype: bool

Este comportamiento no es evidente, y sería estupendo documentarlo más claramente en geopandas . Dicho esto, dependiendo de su aplicación, un unión espacial puede ser la operación que busca.

Answer 2

La explicación de @jdmcbr me ayudó a resolver el problema, pero lo hice de una manera diferente porque mis datos de entrada eran líneas, no puntos. Una unión espacial no habría sido suficiente porque quería saber qué líneas estaban totalmente contenidas por el polígono, no sólo cuáles se cruzaban.

He utilizado una función lambda, aplicando sobre las líneas GeoDataFrame. Esto probablemente sólo funcionará si usted tiene un solo polígono (o multipolygon) en su geodataframe, para conseguir esto usted puede utilizar GeoDataFrame.dissolve()

# lines_gdf = GeoDataFrame containing lines
# poly_gdf = GeoDataFrame containing a polygon or multipolygon

lines_gdf["within"] = lines_gdf["geometry"].apply(lambda geom: poly_gdf["geometry"].contains(geom))
selected_gdf = lines_gdf[lines_gdf["within"]]