GeoPandas – Troubleshooting No Intersection Found Issue

geopandasintersectionpython

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.

Best Answer

Intersects is row-wise operation which aligns both GeoSeries and check is a geometry A in row 0 intersects geometry B in row 0. See the documentation on that: https://geopandas.readthedocs.io/en/stable/docs/reference/api/geopandas.GeoSeries.intersects.html#geopandas.GeoSeries.intersects

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. You can check the docs - https://geopandas.readthedocs.io/en/stable/docs/reference/api/geopandas.sindex.SpatialIndex.query_bulk.html#geopandas.sindex.SpatialIndex.query_bulk

Related Question