33

I have a GeoDataFrame of polygons (~30) and a GeoDataFrame of Points (~10k)

I'm looking to create 30 new columns (with appropriate polygon names) in my GeoDataFrame of Points with a simple boolean True/False if the point is present in the polygon.

As an example, the GeoDataFrame of Polygons is this:

id  geometry
foo POLYGON ((-0.18353,51.51022, -0.18421,51.50767, -0.18253,51.50744, -0.1794,51.50914))
bar POLYGON ((-0.17003,51.50739, -0.16904,51.50604, -0.16488,51.50615, -0.1613,51.5091))

The GeoDataFrame of Points is like this:

counter     points
   1     ((-0.17987,51.50974))
   2     ((-0.16507,51.50925))

Expected output:

counter          points        foo    bar
   1    ((-0.17987,51.50974))  False  False
   1    ((-0.16507,51.50925))  False  False

I can do this manually by:

foo = df_poly.loc[df_poly.id=='foo']
df_points['foo'] = df_points['points'].map(lambda x: True if foo.contains(x).any()==True else False

But given that I have 30 polygons, I was wondering if there is a better way. Appreciate any help!

Kvothe
  • 1,341
  • 7
  • 20
  • 33

1 Answers1

45

Not really clear what kind of data structures you actually have. Also, all your expected results are False, so that's kind of hard to check. Assuming GeoSeries and GeoDataFrames, I would do this:

from shapely.geometry import Point, Polygon
import geopandas

polys = geopandas.GeoSeries({
    'foo': Polygon([(5, 5), (5, 13), (13, 13), (13, 5)]),
    'bar': Polygon([(10, 10), (10, 15), (15, 15), (15, 10)]),
})

_pnts = [Point(3, 3), Point(8, 8), Point(11, 11)]
pnts = geopandas.GeoDataFrame(geometry=_pnts, index=['A', 'B', 'C'])
pnts = pnts.assign(**{key: pnts.within(geom) for key, geom in polys.items()})

print(pnts)

And that gives me:

        geometry    bar    foo
A    POINT (3 3)  False  False
B    POINT (8 8)  False   True
C  POINT (11 11)   True   True
Paul H
  • 65,268
  • 20
  • 159
  • 136
  • Nice answer! I wonder, how would you do this so that instead that each key is a column, there's a single column and the content is the key rather than a True/False? – WGS Jan 18 '18 at 07:33
  • @Manhattan points can be in more than one polygon, so that model breaks down if you want to keep your (geo)dataframe tidy http://vita.had.co.nz/papers/tidy-data.pdf. That said, you could simply `melt` this dataframe and query out the false values. – Paul H Jan 18 '18 at 15:25
  • Aha, nice use of `melt`. Thanks! – WGS Jan 22 '18 at 09:21
  • It fails for Point(35.68333333, 139.75) & Polygon([(90,0),(0,180),(0,-180)]) – jolly Feb 20 '18 at 02:40
  • 3
    @pksingh an important nuance of `shapely` is that it's for geometric, not geographic operations. In other words, `Polygon([(90,0),(0,180),(0,-180)]` is a triangle, not a hemisphere – Paul H Feb 20 '18 at 02:45
  • Thanks @PaulH , is there a direct way for the same? – jolly Feb 20 '18 at 04:50
  • @jolly you need to project your lat/lon into a UTM or other projected coordinate system. BTW, what you mean by "fails"? Do you mean, raises an error or returns an unexpected result? In planar coordinates your point is outside of your triangle. – Paul H Feb 20 '18 at 06:28
  • Thanks @PaulH , Fail is exactly what you said. Point is indeed outside 2d triangle. But is there a way to check if it is inside 3d triangle – jolly Feb 20 '18 at 18:56
  • @jolly `shapely` is for planar geometries. I don't work with volumes so I can't help you there. – Paul H Feb 20 '18 at 20:54
  • For this nice example given by Paul, I get an error: AttributeError: 'GeoSeries' object has no attribute 'items' File "", line 11, in pnts = pnts.assign(**{key: pnts.within(geom) for key, geom in polys.items()}) File "/usr/local/lib/python2.7/dist-packages/pandas/core/generic.py", line 3077, in __getattr__ return object.__getattribute__(self, name) AttributeError: 'GeoSeries' object has no attribute 'items' – Solomon Vimal May 07 '19 at 18:03
  • @SolomonVimal which versions of geopandas and pandas do you have? This still works for me on geopandas v0.4.0 and pandas v0.24.1 – Paul H May 07 '19 at 18:43
  • I have geopandas==0.3.0+48.ga577b12 & geopandas-osm==0.0.1 & pandas==0.20.3 Not sure why it isn't working for me. – Solomon Vimal May 07 '19 at 18:48
  • the error is coming from pandas, so upgrade that first – Paul H May 07 '19 at 18:49
  • I can try to upgrade mine to your version. – Solomon Vimal May 07 '19 at 18:49
  • Now I have geopandas==0.4.0 and geopandas-osm==0.0.1 and pandas==0.24.1, and it works! Thanks! – Solomon Vimal May 07 '19 at 18:51
  • Maybe this could alternatively be done with `sjoin` and then `pivot`? https://geopandas.org/en/stable/gallery/spatial_joins.html#Joins – Richard DiSalvo Apr 12 '23 at 20:36