polygon with holes + arcpy.AsShape is broken?

DustinReagan · ‎10-06-2011

Hi,

I was searching for a way to programmatically create a polygon with a hole, as there isn't any obvious way of doing this in the documentation, I started fumbling around on my own. I came up with this test script:

gjPolygonWithHole  = {
    "type": "Polygon",
    "coordinates": [
                    [[-88.684979400000032, 38.154361399999985], [-88.666690199999962, 38.154329800000042], [-88.67131005539963, 38.154337782166039], [-88.675899944207302, 38.154345712555887], [-88.675899944207302, 38.150691003182175], [-88.67131005539963, 38.150691003182175], [-88.67131005539963, 38.147037172952786], [-88.675899944207302, 38.147037172952786], [-88.685066075460313, 38.147037172952786], [-88.684979400000032, 38.154361399999985]],
                    [[-88.680463013034853, 38.152519617763254], [-88.680463013034853, 38.15068957953018], [-88.678182297717754, 38.15068957953018], [-88.678182297717754, 38.152519617763254], [-88.680463013034853, 38.152519617763254]]]
    }

def main():

    outputShapeDir = r"c:\temp"
    outputShapeFile = r"test.shp"
    spatialRef = r'X:\Nadel&Gussman\McLeansboro\Sections.prj'
    newFC = arcpy.CreateFeatureclass_management( outputShapeDir, outputShapeFile, "POLYGON", "", "", "", spatialRef )
    cur = arcpy.InsertCursor(newFC)
    
    array = arcpy.Array()
    point = arcpy.Point()
    # Create a list to store the features
    features = []
    # Read the coordinates
    for part in gjPolygonWithHole["coordinates"]:
        for coordPair in part:
            point.X = coordPair[0]
            point.Y = coordPair[1]
            array.add(point)
        
    feat = cur.newRow()
    feat.shape = array
    cur.insertRow( feat )

if __name__ == '__main__':
    import arcpy
    arcpy.env.overwriteOutput = True
    main()

Amazingly, it works (I'm not really sure how, since I don't do anything to separate the interior from the exterior, other than suddenly switch winding order):

Unfortunately, this doesn't work in all polygon-with-a-hole cases that I am working with.

So...I noticed the seemingly useful "AsShape" functionality here. Specifically, (notice the AsShape example on doing just this at the bottom of the page i linked).

I put together a little sample script to test this:

gjPolygonWithHole  = {
    "type": "Polygon",
    "coordinates": [
                    [[-88.684979400000032, 38.154361399999985], [-88.666690199999962, 38.154329800000042], [-88.67131005539963, 38.154337782166039], [-88.675899944207302, 38.154345712555887], [-88.675899944207302, 38.150691003182175], [-88.67131005539963, 38.150691003182175], [-88.67131005539963, 38.147037172952786], [-88.675899944207302, 38.147037172952786], [-88.685066075460313, 38.147037172952786], [-88.684979400000032, 38.154361399999985]],
                    [[-88.680463013034853, 38.152519617763254], [-88.680463013034853, 38.15068957953018], [-88.678182297717754, 38.15068957953018], [-88.678182297717754, 38.152519617763254], [-88.680463013034853, 38.152519617763254]]]
    }

def main():

    outputShapeDir = r"c:\temp"
    outputShapeFile = r"test.shp"
    spatialRef = r'X:\Nadel&Gussman\McLeansboro\Sections.prj'
    newFC = arcpy.CreateFeatureclass_management( outputShapeDir, outputShapeFile, "POLYGON", "", "", "", spatialRef )
    cur = arcpy.InsertCursor(newFC)
    
    feat = cur.newRow()
    feat.shape = arcpy.AsShape( gjPolygonWithHole )
    cur.insertRow( feat )

if __name__ == '__main__':
    import arcpy
    arcpy.env.overwriteOutput = True
    main()

FAILURE:

The hole is completely ignored. Any clue as to why?

Regarding holes, I found this thread, and adapted their method to my test script + data. It almost works:

Here's my example code:

gjPolygonWithHole  = {
    "type": "Polygon",
    "coordinates": [
                    [[-88.684979400000032, 38.154361399999985], [-88.666690199999962, 38.154329800000042], [-88.67131005539963, 38.154337782166039], [-88.675899944207302, 38.154345712555887], [-88.675899944207302, 38.150691003182175], [-88.67131005539963, 38.150691003182175], [-88.67131005539963, 38.147037172952786], [-88.675899944207302, 38.147037172952786], [-88.685066075460313, 38.147037172952786], [-88.684979400000032, 38.154361399999985]],
                    [[-88.680463013034853, 38.152519617763254], [-88.680463013034853, 38.15068957953018], [-88.678182297717754, 38.15068957953018], [-88.678182297717754, 38.152519617763254], [-88.680463013034853, 38.152519617763254]]]
    }

def main():

    outputShapeDir = r"c:\temp"
    outputShapeFile = r"test.shp"
    spatialRef = r'X:\Nadel&Gussman\McLeansboro\Sections.prj'
    newFC = arcpy.CreateFeatureclass_management( outputShapeDir, outputShapeFile, "POLYGON", "", "", "", spatialRef )
    cur = arcpy.InsertCursor(newFC)
    
    array = arcpy.Array()
    point = arcpy.Point()
    # Create a list to store the features
    features = []
    # Read the coordinates
    for part in gjPolygonWithHole["coordinates"]:
        for coordPair in part:
            point.X = coordPair[0]
            point.Y = coordPair[1]
            array.add(point)
        null_point = arcpy.Point()
        array.add(null_point)
    feat = cur.newRow()
    feat.shape = arcpy.Polygon( array )
    cur.insertRow( feat )

if __name__ == '__main__':
    import arcpy
    arcpy.env.overwriteOutput = True
    main()

Any ideas why this doesn't work? Why is creating geometry such a hassle with arcpy? Why do some methods work and some do not?

Most likely I'm missing something obvious...
Thanks,
Dustin

KimOllivier · ‎10-06-2011

The key to understanding "donuts" is that there is a [null point] separating the nested rings list of points in each part. Not between parts.

There are some conference technical workshop presentations that set this out very well.

A polygon array in a python pseudocode list would look like this: (of course it has to be in an Arcpy array in internal binary, you cannot just pour this into the array, maybe you can with AsShape, but assembling the list as a JSON string is just as hard). Your example list does not distinguish between parts and donuts.

[[pt,pt,pt,,pt,pt,pt],[pt,pt,pt,pt]]

This is a polygon containing two separate parts, with the first part containing a hole.
You don't have to reverse the directions, ArcGIS does that for you automatically to generate a valid shape. You can have islands inside donuts and so on. Lakes can be like this, an island in the lake that has a small lake in the island.

When you assemble the part array of points and you get to the donut, just add a null point and carry on. I think of it as a 'pen down - pen up' operation that we used to program for pen plotters.

It appears that your error is to add a null point at the end of a part, rather than making the hole within the same part. Likewise adding a null point between parts would not make a donut, just a null part which would be erases by the validator.

Earlier when you reversed the coordinates, it seems that ArcObjects detected that and added in a null for you to make a donut, but I wouldn't rely on that.

# donut example
# with explicit null pair to flag a donut
# this example only has one feature and one part in the feature
# coordinate pairs are shown as tuples for clarity
# Kim Ollivier
coordList = [[
            [(-88.684979400000032, 38.154361399999985), (-88.666690199999962, 38.154329800000042),
            (-88.67131005539963, 38.154337782166039), (-88.675899944207302, 38.154345712555887),
            (-88.675899944207302, 38.150691003182175), (-88.67131005539963, 38.150691003182175),
            (-88.67131005539963, 38.147037172952786), (-88.675899944207302, 38.147037172952786),
            (-88.685066075460313, 38.147037172952786), (-88.684979400000032, 38.154361399999985),
            (),
            (-88.680463013034853, 38.152519617763254), (-88.680463013034853, 38.15068957953018),
            (-88.678182297717754, 38.15068957953018), (-88.678182297717754, 38.152519617763254),
            (-88.680463013034853, 38.152519617763254)]
            ]]
def main():
  array = arcpy.Array()
  point = arcpy.Point()
  null_point = arcpy.Point()
  # Create a list to store the features
  features = []
  # Read the coordinates
  for feature in coordList:
    print "feature", feature
    for part in feature:
      for coordPair in part:
        print coordPair
        if len(coordPair) == 0:
            array.add(null_point)
        else:
            point.X = coordPair[0]
            point.Y = coordPair[1]
            array.add(point)
      # Create the polygon object
    polygon = arcpy.Polygon(array)
    # Clear the array for the next feature
    array.removeAll()
    # Append to the feature list
    features.append(polygon)

    # Copy the features to an output feature class
    arcpy.CopyFeatures_management(features, outputFeatureClass)

if __name__ == '__main__':
    import arcpy
    arcpy.env.overwriteOutput = True
    outputShapeDir = "c:/temp"
    outputShapeFile = "test.shp"
    if not arcpy.Exists(outputShapeDir+"/"+outputShapeFile):
        spatialRef = arcpy.SpatialReference(4326) #  or "GCS_WGS_1984"
        outputFeatureClass = arcpy.CreateFeatureclass_management( outputShapeDir, outputShapeFile, "POLYGON", "", "", "", spatialRef )
    else:
        outputFeatureClass = outputShapeDir+"/"+outputShapeFile
    main()

To go back the other way, see my resource example tool "Fill Donut" on the resources page.
http://resources.arcgis.com/gallery/file/Geoprocessing-Model-and-Script-Tool-Gallery/details?entryID...
I know that Esri has added a tool to do this and you can manually edit donuts, but this is a script to show how to program them.

DustinReagan · ‎10-07-2011

Thank you very much for the reply.

Unfortunately your script also produces incorrect (but close) results:

Additionally, your method appears to be functionally equivalent to the last method that I posted.

You say:

It appears that your error is to add a null point at the end of a part, rather than making the hole within the same part. Likewise adding a null point between parts would not make a donut, just a null part which would be erases by the validator.

If you'll look closely at my code, you'll see that I only ever add a single part to "array", putting a null Point between the points that describe exterior and the points that describe the interior (all in the same part).

array = arcpy.Array()
    point = arcpy.Point()
    # Create a list to store the features
    features = []
    # Read the coordinates
    for part in gjPolygonWithHole["coordinates"]:
        for coordPair in part:
            point.X = coordPair[0]
            point.Y = coordPair[1]
            array.add(point)
        null_point = arcpy.Point()
        array.add(null_point)
    feat = cur.newRow()
    feat.shape = arcpy.Polygon( array )
    cur.insertRow( feat )

Unless I am mistaken, adding multiple parts to a polygon would involve something like:

point = arcpy.Point()
polygonArray = arcpy.Array()
partArray = arcpy.Array()
for part in coordinateList:
    for coordPair in part:
        point.X = coordPair[0]
        point.Y = coordPair[1]
        partArray.add(point)
    #add the null point
    partArray.add( arcpy.Point() )
    polygonArray.add(partArray)
    partArray.removeAll()
newMultiPartPolygon = arcpy.Polygon( polygonArray )
polygonArray.removeAll()

Thanks for your time,
Dustin