Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
add comment

2 Answers

up vote 3 down vote accepted

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.

share|improve this answer
    
You are indeed correct. I agree with your correspondent in the above link that arcpy.Point() giving a 0,0 point is strange. That threw me. Also, because I have such a massive number of points in my rings, I had both inner and outer rings in Arrays. I just wanted to keep them that way rather than deal with individual point objects. So after finding a one-character but game-changing error in my code, I found that the technique I tried before of putting both arrays into a container array (outer ring first) and using that as input ended up working fine. Thanks for the help. –  celticflute Sep 24 '11 at 3:06
    
I've been trying writing arcpy.Point(0,0) in the middle of the array and it does not work for me even thought I've been getting advice to use it. I end up with a 0,0 point in the middle of the array, spiking my polygons with a real 0.0 point included. Would love to see a code example, celticflute! –  Curtis Price Sep 13 '13 at 14:54
add comment

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, cprice@usgs.gov

    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
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.