topic Re: Bug?: Creating polygons with holes using Arcpy in Python Questions

Bug?: Creating polygons with holes using Arcpy

Luke__NWGS_Rogers — Sat, 11 Dec 2021 20:49:12 GMT

I believe I have found what is a bug with the Arcpy site package and how it handles creating polygon geometries...

We are trying to generate polygon geometries using the Arcpy site package and cannot generate polygons with multiple holes. One or two holes in a polygon outer ring seem to work fine but any more than that and the polygon generation fails. In looking at the point arrays that come from an existing multiple hole polygon (using polygon.GetPart()) there is a "None" type object in the point array that indicates that subsequent points are part of a hole. However the arcpy.Array object does not allow the inserting of "None" types or null Points so it appears that holes cannot be created directly and must be inferred by the Arcpy package to create the holes.

Have a look at the attached code and see if you can make it work. The documentation has nothing about this anywhere I can find (other than this) so I have been working backwards to figure it out. The example coordinate list in the attached python file has clockwise outer rings and counter clockwise holes yet still it does not work. The output should look like this (2 polygons):

but instead looks like this:

Here is the code (in case attachment doesn't work):

import arcpy

# Globals
outputFeatureClass = r"D:\temp\geometry\example.mdb\example"
coordList = [[[[0,0],[0,10],[10,10],[10,0],[0,0]],
            [[1,1],[2,1],[2,2],[1,2],[1,1]],
            [[8,8],[9,8],[9,9],[8,9],[8,8]],
            [[8,1],[9,1],[9,2],[8,2],[8,1]],
            [[1,8],[2,8],[2,9],[1,9],[1,8]]],
            [[[10,10], [11,11], [12,10], [10,10]]]]

def main():
    array = arcpy.Array()
    point = arcpy.Point()
    # Create a list to store the features
    features = []
    # Read the coordinates
    for feature in coordList:
        for part in feature:
            for coordPair in part:
                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__':
    main()

Please help!

Re: Bug?: Creating polygons with holes using Arcpy

DanPatterson_Retired — Sat, 11 Dec 2021 20:49:15 GMT

Luke...a quick test...

# Globals
outputFeatureClass = r"c:\temp\donutdemo.shp"
coordList = [[[[0,0],[0,10],[10,10],[10,0],[0,0]],
              [[1,1],[2,1],[2,2],[1,2],[1,1]],
              [[8,8],[9,8],[9,9],[8,9],[8,8]],
              [[8,1],[9,1],[9,2],[8,2],[8,1]],
              [[1,8],[2,8],[2,9],[1,9],[1,8]]],
             [[[10,10], [11,11], [12,10], [10,10]]]]
def main():
  array = arcpy.Array()
  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:
        point.X = coordPair[0]
        point.Y = coordPair[1]
        array.add(point)
      null_point = arcpy.Point()
      array.add(null_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
  main()

Re: Bug?: Creating polygons with holes using Arcpy

Luke__NWGS_Rogers — Thu, 09 Dec 2010 17:03:56 GMT

Hi Dan-

This does indeed work which is interesting since calling arcpy.Point() creates a point object with 0.0 for X and Y which I would think would be valid coordinates so I didn't even try the solution you discovered. A little qwerky but it works!

Thanks Dan!

Re: Bug?: Creating polygons with holes using Arcpy

DanPatterson_Retired — Sat, 11 Dec 2021 20:49:18 GMT

Luke
that appears to be the case, which appears to go against the documentation which suggests that X and Y would be assigned a value of none

>>> from arcpy import Point
>>> pnt = Point()
>>> print str(pnt)
0 0 NaN NaN
>>>

| Methods inherited from arcpy.arcobjects.mixins.PointMixin:
|
| __init__(self, X=None, Y=None, Z=None, M=None, ID=None)
| Point({X}, {Y}, {Z}, {M}, {ID})

which suggests that None is being translated to 0...perhaps ESRI could provide some rationale which I presume has something to do with None not being able to be represented as a numeric value for coordinates...go figure

Re: Bug?: Creating polygons with holes using Arcpy

DustinReagan — Thu, 06 Oct 2011 20:03:06 GMT

I tried this method to create a polygon with holes, and nearly succeeded. Take a look at this thread for details. Any help explaining why the method described here fails in my case would be appreciated.

Thanks,
Dustin

Re: Bug?: Creating polygons with holes using Arcpy

curtvprice — Fri, 13 Sep 2013 14:02:17 GMT

I've been trying this method too, including a arcpy.Point(None, None) in the middle of the point array. I even got a tip from the arcpy team (through the help feedback button) to use Dan's suggestion. But it does not work for me.

I end up with a 0,0 point added to my ring when the feature is created (ie when I write the polygon to a shape).

[ATTACH=CONFIG]27437[/ATTACH]

I have an open incident going with Esri on this and hopefully will have an answer soon as we go back and forth. So far my analyst is running into the same problems. At this point we're chalking it up to missing documentation. There's GOT to be a way to do this. Someone on stack exchange apparently figured it out but did not (!) post their code.

[#NIM094838 Include a sample in the help on how to build a polygon with a hole using python. ]
[#NIM093332 Please provide some additional documentation and examples for creating donut polygons using arcpy geometry objects

Creating polygons with holes using ArcPy

curtvprice — Sat, 11 Dec 2021 20:49:21 GMT

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(
)
del Rows

another bugfix

Re: Bug?: Creating polygons with holes using Arcpy

FridjofSchmidt — Wed, 17 Nov 2021 20:50:12 GMT

Thanks Curtis, the code works great for the given examples, and also for polygons with multiple interior rings (holes). However, I can't figure out how to create a multipart polygon with a hole in one of its parts (the union of Feat0 and Feat2 - tested with ArcGIS Desktop 10.8.1). The returned polygon will only contain the first part:

Feat3 = [ [[3.0, 8.0], [1.0, 8.0], [2.0, 10.0], [3.0, 8.0]], [[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]] ]

When I try to print the returned polygon with the following method:

def printpoly(polygon): print "Area: {0}".format(polygon.area) print "Length: {0}".format(polygon.length) for part in polygon: for pnt in part: print pnt

I get this result:

Area: 2.0 Length: 6.472135955 3,0001220703125 8,0001220703125 NaN NaN 1,0001220703125 8,0001220703125 NaN NaN 2,0001220703125 10,0001220703125 NaN NaN 3,0001220703125 8,0001220703125 NaN NaN

Am I missing something or is there a bug in Esri's Polygon constructor?

Re: Bug?: Creating polygons with holes using Arcpy

lsturtevant — Thu, 24 Aug 2023 14:04:00 GMT

I realize this is a very old post, but after grappling with this issue for a while I think I might have a workaround using the symmetricDifference function for a Polygon—ArcGIS Pro | Documentation object. My specific workflow was to compute M and Z values on an existing Polygon ZM object and the only way I could figure out how to do this was to extract the exterior and interior rings into separate arrays, reassemble the parts and then run the symmetricDifference with the exterior parts and the interior parts. I've modified @curtvprice code and posted it here as well as an attachment. I believe this also helps with the issue @FridjofSchmidt was having.

#
# Write a polygon feature class with interior rings (donut holes)
#

import os
import arcpy
from arcpy import env


def makepoly(coord_list, SR):
    exteriorRing = arcpy.Array()
    interiorRing = arcpy.Array()
    for part in coord_list:
        print(part)
        exteriorSubArray = arcpy.Array()
        interiorSubArray = arcpy.Array()
        isExterior = True
        for pnt in part:
            if pnt:
                px, py = pnt
                pnt = arcpy.Point(px, py)
            if not pnt:
                isExterior = False
                if interiorSubArray.count:
                    interiorRing.add(interiorSubArray)
                interiorSubArray.removeAll()
            if pnt and isExterior:
                exteriorSubArray.add(pnt)
            if pnt and not isExterior:
                interiorSubArray.add(pnt)
        if interiorSubArray.count:
            interiorRing.add(interiorSubArray)
        exteriorRing.add(exteriorSubArray)

    exteriorPolygon = arcpy.Polygon(exteriorRing, SR, True, True)
    if interiorRing.count:
        interiorPolygon = arcpy.Polygon(interiorRing, SR, True, True)
        finalPolygon = exteriorPolygon.symmetricDifference(interiorPolygon)
    else:
        finalPolygon = exteriorPolygon
    return finalPolygon


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]]
]

Feat3 = [
    [[6.0, 11.0],
     [4.0, 11.0],
     [5.0, 13.0],
     [6.0, 11.0]],

    [[12.0, 14.0],
     [12.0, 11.0],
     [9.0, 11.0],
     [9.0, 14.0],
     [12.0, 14.0],
     None,
     [10.0, 13.0],
     [10.0, 12.0],
     [11.0, 12.0],
     [11.0, 13.0],
     [10.0, 13.0]]
]
# create the empty feature class

# with real data, provide a SR code, name or dataset for SR
SR = None
output_polygon = r'C:\ArcGIS\Default.gdb\TestPoly'

# create the polygons and write them
polyList = []
for f in [Feat0, Feat1, Feat2, Feat3]:
    print(f"coords: {f}")
    p = makepoly(f, SR)
    polyList.append(p)
    print(f"feature: {p}")

arcpy.management.CopyFeatures(polyList, output_polygon)