3

I have 2 geodataframes:

  1. An original shapefile of water features read in as a geodataframe. (playas)
  2. A gdf as a result of a convex_hull operation with only polygon geometries. (polygons)

There should be intersection between these two layers, but the following code returns an empty geodataframe:

aff_playas = playas[playas.intersects(polygons)]

I am trying to think of what the problem could be with these two layers. I believe it is one of the following:

  • A projection inconsistency with the created hulls_gdf
  • function failure due to how detailed, varied the playas layer is.

The two layers plotted below have similar extent and are verified to intersect in Qgis, ArcGIS. The third is the output from the intersects().

hulls, plotted

Playas, plotted

Intersection visible in Qgis

7
  • 1
    Good idea to check the projections of the both layers. Add the result to your question.
    – user30184
    Apr 23, 2018 at 18:42
  • I added some more information about my code in the question, and I checked the crs for both layers and the seem to be consistent.
    – Stephen
    Apr 23, 2018 at 19:17
  • Doesn't "Intersects" return a boolean value? What do you expect to get as a sum of Trues and Falses?
    – user30184
    Apr 24, 2018 at 6:14
  • Good point, the following code also returns no results: aff_playas = playas[playas.intersects(polygons)]
    – Stephen
    Apr 25, 2018 at 13:59
  • 1
    Ah, I guess you're using shapely's intersects method. You could try int_gdf = gpd.sjoin(hulls_gdf, water_gdf, op='intersects').
    – Jon
    Apr 25, 2018 at 14:12

3 Answers 3

3

I have solved my problem. For some reason intersects() was not working between my two layers, even though they are shown to intersect. Here is my workaround:

sj = gpd.sjoin(playas, polygons, how='inner',op='intersects')
aff_playas = playas.geom_almost_equals(sj)
aff_playas = aff_playas[~aff_playas.index.duplicated()]
playa_subset = playas[aff_playas]

Since I couldn't get intersects() to work, I used spatial join and the geom_almost_equals function to subset the original layer to the same intended affect.

1

geopandas's .intersects works by comparing the two geodataframes based on the index values. It's an element wise comparison. If you want to compare 1 shape (eg: convex hull) against all the features, you can use your code as is, but you need to ensure that polygons is a shapely object (for example: union all the geoms to one multipolygon shapely object)

1
  • Please give a code sample and describe your solution more explicitly. Oct 5, 2020 at 16:20
1

Below code work for me.

sj = gpd.sjoin(filter_street, buffer_project_dataframe, how='inner',op='intersects')

But make sure that, buffer_project_dataframe should be a geodataframe, geosereis is not working

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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