2

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

  • 1
    Good idea to check the projections of the both layers. Add the result to your question. – user30184 Apr 23 '18 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 '18 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 '18 at 6:14
  • Good point, the following code also returns no results: aff_playas = playas[playas.intersects(polygons)] – Stephen Apr 25 '18 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 '18 at 14:12
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.

| improve this answer | |
0

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)

| improve this answer | |

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.