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:
to:
Check Geometry? any reports from it?
Check Geometry (Data Management)—ArcGIS Pro | Documentation
No, it came back without any outputs and nothing in the table.
This is the shapefile in a geodatabase, 260381 vertices, 1 part, 18 holes and 0 curves
of course, the channel isn't closed at either end, perhaps due to a conversion or clipping artifact
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.