0

I know the same question has been answered elsewhere; however, the solution did not work for me. I already have both of my objects as shapely geometry.polygon.Polygon objects. I have two objects that intersect clearly, but are returning False when I want to output whether they intersect (i.e., it should return True).

Coordinates of Poly1 = (-109.91824, 32.715309) (-109.97012, 32.70005) (-110.09524, 32.632911) (-110.24172, 32.605446) (-110.35159, 32.584083) (-110.47061, 32.571876) (-110.54995, 32.602394) (-110.55606, 32.666481) (-110.51028, 32.727516) (-110.4523, 32.782448) (-110.38516, 32.822121) (-110.26614, 32.840431) (-110.15017, 32.852638) (-110.05251, 32.849586) (-109.96401, 32.797706) (-109.91519, 32.758034) (-109.91824, 32.715309)

Coordinates of Poly2 = (-126.0, 43.0) (-103.0, 43.0) (-126.0, 26.0) (-103.0, 26.0) (-126.0, 43.0)

Notice that Poly2 is just a large rectangle. I essentially want to get every shapefile within this rectangle. It also might be helpful to note that Poly1 was read from a GIS shapefile. Poly2 was created by me using the following code:

p1 = geometry.Point(-126,43)
p2 = geometry.Point(-103,43)
p3 = geometry.Point(-126,26)
p4 = geometry.Point(-103,26)
pointList = [p1, p2, p3, p4, p1]
poly2 = geometry.Polygon(pointList)

Does anyone have any idea why this might be returning False when it should be returning true?

It is also important to note that poly1 is created as the jth row of a geopandas dataframe, but I have confirmed that the output of that is a geometry.polygon.Polygon object:

poly1 = gdf_i.geometry.iloc[j]

Lastly: For reference, that dark red shape right below Phoenix (overlayed on the orange shape) is Poly1. You can visually see that it definitely intersects with Poly2 (the bounds of the map). enter image description here

Any help is appreciated. Thank you!

Edit: Code for intersection is as follows:

if poly2.intersects(poly1) == True:
        gdf_j = gdf_i.iloc[j]
Palle Due
  • 5,929
  • 4
  • 17
  • 32
Richard Strott
  • 41
  • 5
  • Does this answer your question? [intersection of two geopandas GeoSeries gives warning "The indices of the two GeoSeries are different." and few matches](https://stackoverflow.com/questions/72185667/intersection-of-two-geopandas-geoseries-gives-warning-the-indices-of-the-two-ge) – Michael Delgado Jul 12 '22 at 14:40
  • 1
    Please post the code you’re using to get the behavior you’re seeing. I’m assuming you’re using gdf.intersection? Use sjoin instead. – Michael Delgado Jul 12 '22 at 14:41
  • @MichaelDelgado See the edit above. The code wasn't working well here. I'll try sjoin. Thanks! – Richard Strott Jul 12 '22 at 14:46
  • 1
    Maybe you're confusing intersect and contain ? On the picture you added, the dark red shape doesn't intersect with the bounds of the Arizona state but is contained within. – Pauuuuuuul Jul 12 '22 at 14:51
  • @Pauuuuuuul intersect is a superset of contains in the shapely predicates grammar. The problem is a common misunderstanding of using binary predicates vs. sjoin in geopandas. – Michael Delgado Jul 12 '22 at 15:03
  • @Pauuuuuuul intersect should work to though. All of the other shapes (the orange, and the white) are showing up in the final output of this task. – Richard Strott Jul 12 '22 at 15:05
  • 1
    Can you post all of your code? This looks like it’s inside a for loop? What’s j? Please see the guide to [ask] and try to create a [mre] – Michael Delgado Jul 12 '22 at 15:05
  • are you sure your polygons are valid? – Ian Turton Jul 12 '22 at 15:18
  • @MichaelDelgado j is just an integer. I'm looping through each row of a geodataframe with all polygons, then assigning the geometry to a variable. I think what is important is that I checked that poly1 and poly2 are both geometry.polygon.Polygons – Richard Strott Jul 12 '22 at 15:46
  • @IanTurton How do I determine a valid polygon? – Richard Strott Jul 12 '22 at 15:47
  • @IanTurton I made sure the polygons were valid and it didn't help. – Richard Strott Jul 12 '22 at 16:02
  • @RichardStrott did you actually test with [`poly2.is_valid`](https://shapely.readthedocs.io/en/stable/manual.html#object.is_valid)? I think that's what Ian Turton is talking about – Michael Delgado Jul 12 '22 at 16:49
  • @MichaelDelgado I did the .is_valid test. Poly was not valid. So I changed it to be valid with a method I saw online. It still wasn't picking up on my intersection, so I assumed that the validation wasn't the problem. It turns out I had made it valid but still maintained the hourglass shape, which cut out a lot of shapes that I wanted to keep. – Richard Strott Jul 12 '22 at 18:30
  • @IanTurton See comment above. Thanks for the help! – Richard Strott Jul 12 '22 at 18:30
  • fwiw, this is why we ask for a *complete* [mcve]. we didn't even have the shape you were using! glad it worked – Michael Delgado Jul 12 '22 at 19:15
  • @MichaelDelgado Touche - I'll make sure to include in the future. – Richard Strott Jul 13 '22 at 01:27

1 Answers1

0

Your bounding box forms an hour glass, not a square. You should change your poly2 definition to:

p1 = geometry.Point(-126,43)
p2 = geometry.Point(-103,43)
p3 = geometry.Point(-103,26)
p4 = geometry.Point(-126,26)
pointList = [p1, p2, p3, p4, p1]
poly2 = geometry.Polygon(pointList)
Michael Delgado
  • 13,789
  • 3
  • 29
  • 54