6 Replies Latest reply on Sep 13, 2013 4:24 PM by curtvprice

    Bug?: Creating polygons with holes using Arcpy

    lukewrogers
      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
          Dan_Patterson
          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
            lukewrogers
            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
              Dan_Patterson
              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
                dustfurn
                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
                  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
                    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
                    


                    another bugfix