ArcGIS 10.0/arcpy: How to create a polygon geometry from inner and outer ring point arrays?

Question

I have two arrays, each containing a closed ring of points, one for an outer ring and the other for an inner ring, the latter running counterclockwise, of course. I am using the arcpy.Polygon() factory function to build the polygon. My question: how do I input the points to Polygon()?

Documentation gives the general syntax as Polygon(inputs, spatial reference, has_Z, has_M). The inputs part is what is sketchy. The inputs supposedly can be point objects or arrays. My problem is that I have a polygon with a doughnut hole. I can build a simple polygon with a single point array as input ... no problem; however, in what way do I add the two sets of points? In the "Reading Geometries" help page it says that parts (in this case rings) are separated by a null point object. For polygons with holes, the outer ring goes first. So, I tried adding the outer array then the inner array to a container array. That didn't work: no hole. I tried the same thing with a null point object between them. No joy: same result. Then, I tried putting both arrays into a Python list as the input. That threw an error. The sample scripts in the help files never deal with writing polygons with holes. So, I really don't know what format to use for the inputs. Does anyone have any insight into this one? I would surely appreciate the help. Thanks, all.

Answer 1

Since celticflute did not post his code, I had to figure it out again.

Finally got a working sample. This seems to work pretty well.

The key thing is that you do not write polygon parts and holes differently!

Instead, you construct a polygon from nested arrays which represent parts made of up rings. When you write the geometry by writing it to a shape field or use it in a tool, behind the scenes arc objects planarizes the rings within in each polygon part and determines what's a hole and what isn't. And sorts the coordinates clockwise, and burns the xys into the coordinate system and domain of the geodatabase feature class.

#
# Write a polygon feature class
#

import os
import arcpy
from arcpy import env

def makepoly(coord_list, SR=None):
    """Convert a Python list of coordinates to an ArcPy polygon feature

    Author: Curtis Price, USGS, [email protected]

    Examples, from Desktop Help 10.x: Reading Geometries

    Feat0 = [
            [[3.0, 8.0],
            [1.0, 8.0],
            [2.0, 10.0],
            [3.0, 8.0]]
            ]
    Feat1 = [
            [[5.0, 3.0],
            [3.0, 3.0],
            [3.0, 5.0],
            [5.0, 3.0]],

            [[7.0, 5.0],
            [5.0, 5.0],
            [5.0, 7.0],
            [7.0, 5.0]],
            ]

            # this feature has an interior ring (donut)

            Feat2 = [
            [[9.0, 11.0],
            [9.0, 8.0],
            [6.0, 8.0],
            [6.0, 11.0],
            [9.0, 11.0],
            None,
            [7.0, 10.0],
            [7.0, 9.0],
            [8.0, 9.0],
            [8.0, 10.0],
            [7.0, 10.0]]
            ]
"""

    parts = arcpy.Array()
    rings = arcpy.Array()
    ring = arcpy.Array()
    for part in coord_list:
        for pnt in part:
            if pnt:
                ring.add(arcpy.Point(pnt[0], pnt[1]))
            else:
                # null point - we are at the start of a new ring
                rings.add(ring)
                ring.removeAll()
        # we have our last ring, add it
        rings.add(ring)
        ring.removeAll()
        # if we only have one ring: remove nesting
        if len(rings) == 1:
            rings = rings.getObject(0)
        parts.add(rings)
        rings.removeAll()
    # if single-part (only one part) remove nesting
    if len(parts) == 1:
        parts = parts.getObject(0)
    return arcpy.Polygon(parts, SR)

# test data from:
# Desktop Help 10.0: Reading Geometries

Feat0 = [
        [[3.0, 8.0],
        [1.0, 8.0],
        [2.0, 10.0],
        [3.0, 8.0]]
        ]
Feat1 = [
        [[5.0, 3.0],
        [3.0, 3.0],
        [3.0, 5.0],
        [5.0, 3.0]],

        [[7.0, 5.0],
        [5.0, 5.0],
        [5.0, 7.0],
        [7.0, 5.0]],
        ]

# this last feature has an interior ring (donut)

Feat2 = [
        [[9.0, 11.0],
        [9.0, 8.0],
        [6.0, 8.0],
        [6.0, 11.0],
        [9.0, 11.0],
        None,
        [7.0, 10.0],
        [7.0, 9.0],
        [8.0, 9.0],
        [8.0, 10.0],
        [7.0, 10.0]]
        ]

# test code

# create the empty feature class

# with real data, provide a SR code, name or dataset for SR
# SR = arcpy.SpatialReference(4326)
SR = None
env.workspace = env.scratchGDB
Data = arcpy.CreateScratchName("","","featureclass",env.workspace)
print "writing: " + Data
print
arcpy.CreateFeatureclass_management(os.path.dirname(Data),
                                    os.path.basename(Data),"Polygon",
                                    spatial_reference=SR)

# create the polygons and write them

Rows = arcpy.da.InsertCursor(Data, "SHAPE@")
for f in [Feat0, Feat1, Feat2]:
    print "coords: " +  repr(f)
    p = makepoly(f)
    print "feature: " + repr(p)
    Rows.insertRow([p])
del Rows

Answer 2

Creating donut polygons was discussed on in this link with a python code example (for use in ArcMap). The trick being to add a null point between the list of coordinates between the outer and inner ring. Note, the outer ring points are clockwise and the inner ring points are counterclockwise.