+ Reply to Thread
Results 1 to 7 of 7

Thread: Bug?: Creating polygons with holes using Arcpy

  1. #1
    Luke (NWGS) Rogers
    Join Date
    Jul 2010
    Posts
    2
    Points
    0
    Answers Provided
    0


    0

    Default Bug?: Creating polygons with holes using Arcpy

    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):

    Code:
    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!
    Attached Files
    Last edited by lukewrogers; 12-09-2010 at 09:55 AM.

  2. #2
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,753
    Points
    442
    Answers Provided
    41


    0

    Default Re: Bug?: Creating polygons with holes using Arcpy

    Luke...a quick test...
    Code:
    # 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()
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  3. #3
    Luke (NWGS) Rogers
    Join Date
    Jul 2010
    Posts
    2
    Points
    0
    Answers Provided
    0


    0

    Thumbs up Re: Bug?: Creating polygons with holes using Arcpy

    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!

  4. #4
    Dan Patterson

    Join Date
    Apr 2010
    Posts
    1,753
    Points
    442
    Answers Provided
    41


    0

    Default Re: Bug?: Creating polygons with holes using Arcpy

    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
    Code:
    >>> 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
    Geomatics, Carleton University, Ottawa, Canada
    http://obidangis.blogspot.ca/

  5. #5
    Dustin Reagan
    Join Date
    Jul 2011
    Posts
    19
    Points
    0
    Answers Provided
    0


    0

    Default Re: Bug?: Creating polygons with holes using Arcpy

    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

  6. #6
    Curtis Price

    Join Date
    Oct 2009
    Posts
    1,796
    Points
    874
    Answers Provided
    127


    0

    Default Re: Bug?: Creating polygons with holes using Arcpy

    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).

    Click image for larger version

Name:	bad_write_geometry.PNG
Views:	25
Size:	5.0 KB
ID:	27437

    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
    Last edited by curtvprice; 09-13-2013 at 07:05 AM.

  7. #7
    Curtis Price

    Join Date
    Oct 2009
    Posts
    1,796
    Points
    874
    Answers Provided
    127


    1

    Default Creating polygons with holes using ArcPy

    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.

    Code:
    #
    # 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
    another bugfix
    Last edited by curtvprice; 09-16-2013 at 10:36 AM. Reason: script fix

+ Reply to Thread

Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts