1

I want to check if points layer that I have intersect with polygons layer that I have, as Boolean column in the points dataframe.
I have two GeoPandas dataframes, the first one is many points and looks like this:

>>>ID  geometry
0  12  POINT (5.0279 7.4547)
1  45  POINT (6.6539 12.139)
...

and the second dataframe is one layer of many different polygons that looks something like this:

>>>name     code   geometry
0  Desert   12     POLYGON ((5.52013 13.8902, 5.5265 13.892,...)
1  Water    24     POLYGON ((5.53756 13.88472, 5.5291 13.8791,...)
...

I'm trying to check if there is any intersection between the points layer and the regions layer. For that I have determined the crs and then used intersects :

regions=regions.to_crs({'init': 'epsg:4326'})
points=points.set_crs({'init': 'epsg:4326'})
inter=points.geometry.intersects(regions.geometry)

A The script runs with the following warnings:

FutureWarning: '+init=:' syntax is deprecated. ':' is the preferred initialization method. When making the change, be mindful of axis order changes: https://pyproj4.github.io/pyproj/stable/gotchas.html#axis-order-changes-in-proj-6 return _prepare_from_string(" ".join(pjargs)) /opt/conda/lib/python3.8/site-packages/geopandas/base.py:39: UserWarning: The indices of the two GeoSeries are different.
warn("The indices of the two GeoSeries are different.")

Then when I check the results the only value is False, like all th epoints do not intersect:

inter.unique().tolist()
>>>[False]

*I saw on QGIS that the there are points that intersect and there are points that do not so there is no way this result is true

*I have checked the dtypes - each one of my geodataframes has one column that is geometry and called geometry.

My end goal: to add new column in the points geodataframe that will tell if it intersects the regions or not.

4

Intersects is row-wise operation which aligns both GeoSeries and check is a geometry A in row 0 intersects geometry B in row 0.

For check if there is any intersection, use spatial index.

import numpy as np

inp, res = regions.sindex.query_bulk(points.geometry, predicate='intersects')
points['intersects'] = np.isin(np.arange(0, len(points)), inp)

Regarding the CRS warning, just remove init part and pass it as a string.

regions=regions.to_crs('epsg:4326')
points=points.set_crs('epsg:4326')

query_bulk requires geopandas 0.8 and is not yet documented, but wraps pygeos.STRtree.query_bulk, so you can check that docs - https://pygeos.readthedocs.io/en/latest/strtree.html#pygeos.strtree.STRtree.query_bulk.

| improve this answer | |
  • thank you for your answer, the problem is that even though the geometry columns are type geometry, when I try this I get error "TypeError: Got geometry of type <class 'str'>, geometry must be a shapely geometry." – Reut Sep 30 at 8:33
  • That means that your geometry column contains strings instead of shapely.geometry. You will have to convert those incorrect ones. See stackoverflow.com/questions/56433138/… – martinfleis Sep 30 at 8:35
  • 1
    I ffound the problem- I had also to put ".geometry" after the points dataframe in the beackets. thank you for your answer! – Reut Sep 30 at 10:37
  • 1
    @hhh The code above will work even without pygeos with geopandas 0.8 – martinfleis Nov 3 at 16:53
  • 1
    @hhh only optionally - geopandas.readthedocs.io/en/latest/getting_started/… – martinfleis Nov 3 at 19:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.