Aggregate Polygons seems to be creating an invalid geometry

463
4
03-03-2022 10:45 AM
williamwinner
New Contributor III

I'm running a script that takes a raster file and creates depth areas from it.  The process involves a number of tools including Aggregate Polygons.

It is this tool that seems to be causing me problems.  After I run this tool, I'm running a custom script to reduce the number of vertices but when the script runs, one of the holes in a feature seems to be completely missing.  (I've run this on dozens of rasters dozens of times and only had this one feature out of thousands have a problem.)

When I process the file, the end result is always missing a very large hole.  So, I've investigated by exporting all of the feature vertices to points.  There are 260,381 vertices in the feature.  But, when I run the following:

 

cursor = arcpy.da.SearchCursor(inFile, ["SHAPE@"])

for feat in cursor:
    arcpy.AddMessage("part count=" + str(feat[0].partCount))
    parts = feat[0].getPart()

    for part in parts:
        arcpy.AddMessage("  -" + str(part.count))

 

it returns:

part count=1
-205085

I converted the feature to lines and then exported just the missing hole and it has 55,313 vertices.  There should be 18 holes in the feature.  Without that hole, I would expect to find 205,085 points in the feature (205,068 + 17 null points indicating the start of an inner ring).

So, why is the shape geometry not finding all of the points?  When I export it to it's vertices, it finds them all but when I run the above code, it misses an entire hole.  And repair geometry doesn't fix it.  After, it still shows one part with 205,085 vertices.  But, when I open it in ArcMap, it shows the hole fine.

I'm at a loss to figure out how to see this particular hole.  I can't fix the problem if I can't figure out how to find it.

Any suggestions?  I've tried everything I can think of to get the data out.  I've even tried exporting to JSON and then back from JSON but it has a problem parsing the JSON file.

 

UPDATE:

I just tried to take the raw input from a SearchCursor and output it to an InsertCursor and the hole also disappears.  I've also submitted a Bug Report because it seems like the Aggregate Polygons tool has broken this feature.  Here's the code I just ran:

 

import os

inFile = arcpy.GetParameterAsText(0)
outFile = arcpy.GetParameterAsText(1)
path, name = os.path.split(outFile)
arcpy.AddMessage(path + "\\" + name)
arcpy.CreateFeatureclass_management(path, name, "POLYGON", spatial_reference=arcpy.Describe(inFile).spatialReference)

cursor = arcpy.da.SearchCursor(inFile, ["SHAPE@"])
outputCursor = arcpy.da.InsertCursor(outFile, ["SHAPE@"])
for feat in cursor:
    arcpy.AddMessage("part count=" + str(feat[0].partCount))
    parts = feat[0].getPart()

    outputCursor.insertRow([arcpy.Polygon(parts)])

    arcpy.AddMessage(str(parts.count))
    for part in parts:
        arcpy.AddMessage("  -" + str(part.count))

del cursor, outputCursor

 

The feature went from:

williamwinner_0-1646337675340.png

to:

williamwinner_1-1646337698750.png

 

0 Kudos
4 Replies
DanPatterson
MVP Esteemed Contributor

Check Geometry? any reports from it?

Check Geometry (Data Management)—ArcGIS Pro | Documentation


... sort of retired...
0 Kudos
williamwinner
New Contributor III

No, it came back without any outputs and nothing in the table.

0 Kudos
DanPatterson
MVP Esteemed Contributor

bad_sp.png

This is the shapefile in a geodatabase, 260381 vertices, 1 part, 18 holes and 0 curves

bad_sp2.png

of course, the channel isn't closed at either end, perhaps due to a conversion or clipping artifact


... sort of retired...
0 Kudos
williamwinner
New Contributor III

I've worked on it in both a GDB and the in_memory workspace.  My larger script uses the in_memory workspace to create temporary feature classes.  But, during de-bugging, I've been running through the script step by step and storing everything in a GDB.  I have the same issues whether in the GDB or not.  I haven't tried it in ArcPro yet as we're only starting to switch over and the script needs to work in ArcMap.

0 Kudos