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!
# 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()
>>> from arcpy import Point >>> pnt = Point() >>> print str(pnt) 0 0 NaN NaN >>>
# # 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
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?
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)