2

I am working on a script where I identify polygons that contain other polygons, and I've run into a very perplexing issue. In the image below, green polygons have been identified as "container" polygons, and yellow polygons are supposedly polygons that are either not containers or are not contained.

Parcels

The red highlighted polygon has holes where the smaller "polygon-itos" (see turquoise polygon) are located, but it should still be green. When checking how the script might have gone awry, I found that the polygon object was not reporting the bounding box of the container polygon correctly. I obtained maximum coordinates for each polygon using a search cursor:

polygon_dict = {}
with arcpy.da.SearchCursor(polygon_fc, ['SHAPE@', ID]) as cur:
    for row in cur:
        xmax = row[0].extent.XMax; ymax = row[0].extent.YMax
        xmin = row[0].extent.XMin; ymin = row[0].extent.YMin
        polygon_dict[row[1]] = [xmax, ymax, xmin, ymin]

polygon_dict[red_polygon] = [3072744.5390000045,
                             1254989.898000002,
                             3072732.677000001,
                             1254981.4430000037]
polygon_dict[turquoise_polygon] = [3072811.898000002,
                                   1255530.0219999999,
                                   3072794.446999997,
                                   1255488.5189999938]

turq = polygon_dict[turquoise_polygon]
red = polygon_dict[red_polygon]

When I compared the maximum coordinates of the red polygon and the turquoise polygon using turq[0] < red[0] and turq[1] < red[1] and turq[2] > red[2] and turq[3] > red[3], I expected to see that turquoise would be completely within the red, but that was not the case. Just looking at the picture indicated that the turquoise polygon should at the very least touch the red one, but even testing that theory with the polygon geometry method red_poly.touches(turquoise_poly) returned false.

Why is this happening when the turquoise polygon's bounding box is clearly completely inside the red polygon's bounding box?

  • Do you have the polygons backwards? Just looking at the y min/max values, the "turquoise" is much taller than the "red" one AKA they may be switched. – mkennedy Feb 9 '18 at 1:36
  • @mkennedy I'm sure they're not switched, but I wish they were! – jesnes Feb 9 '18 at 16:23
1

I assume the values you gave are the ones directly from the extent objects. Just to note: I don't know if you copied or typed the code in but there is a typo (ximn vs xmin). I modified the code a little.

polygon_dict = {}
red_polygon = "r"
turquoise_polygon = "t"
polygon_dict[red_polygon] = [3072744.5390000045,
                             1254989.898000002,
                             3072732.677000001,
                             1254981.4430000037]
polygon_dict[turquoise_polygon] = [3072811.898000002,
                                   1255530.0219999999,
                                   3072794.446999997,
                                   1255488.5189999938]

turq = polygon_dict[turquoise_polygon]
red = polygon_dict[red_polygon]
print (turq[0] < red[0])
print (turq[2] > red[2])
print (turq[3] > red[3])
print (turq[1] < red[1])

Assuming the list order is xmax, ymax, xmin, ymin, then those should all be true when it prints. I get False True True False. So the xmax for the bigger polygon is not greater than the xmax for the smaller one (same for ymax). So the problem is with the maximum values. Either somewhere in your dictionary creation or with the coordinates you've provided. Is it possible you have duplicate ids in the dataset, so you are overwriting the dictionary values inadvertently?

You could also try using contains (instead of touches) directly from the extent object to test if the little polygon is inside the bigger one: https://resources.arcgis.com/en/help/main/10.1/index.html#//018z00000072000000

Store the extent object directly in the dictionary if you need to instead of a list.

Also, is it possible you have multipart features? Maybe that is giving you the bigger than expected extent.

Another suggestion, you can write out the extents as a polygon in a new feature class. That would visually show you the problem areas.

  • Thanks for your time and noticing the code inconsistency (just a copy paste error). Yes, these are the values I get on the red and turq polygons. It looks like this specific case was duplicate IDs, though that's not the case for every false-negative container polygon in the dataset, so I'm going to have to do a little more digging. – jesnes Feb 9 '18 at 17:10
  • Maybe use the OID/FID as a unique identifier for each feature to test that your code is working. – dslamb Feb 9 '18 at 17:21

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.