Thursday, 6 June 2013

Donut Polygons using Python

The thread listed here from the ESRI forum contains a discussion on creating donut polygons using Python and arcpy.  It also demonstrates how to use "def" 's and the if __name__ == '__main__':  to produce a standalone script to run for testing purposes.  Also note that the "raw" formatting of the output shapefile name.  This is equivalent to using \\ or / slashes in the output filename.  Load the shapefile into ArcMap to see the results.  Note, that the shapefile is defined as a list of list of lists of coordinate values.  Arc* defines a polygon in clockwise order and holes (holes) in counter-clockwise order and multipart features (whether they are discrete parts or multiple holes) as a list of lists.  In this example the outer bounding box is defined as
[[0,0],[0,10],[10,10],[10,0],[0,0]]  which is a list of coordinates... the remaining holes etc follow this logic, so build slowly if coding by hand.


# 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()

No comments:

Post a Comment