Sunday, 21 July 2019

arcpy - How to create multipart line and polygon with interior ring in arcgisscripting?


My objective is to swap X,Y coordinates of feature classes using Python and arcgisscripting.


Based on help examples, I managed to write script which reads and writes all geometry types.


The problem is with polyline geometry: for multipart lines, all parts are connected and they should be disjoint. It seems that connection between last vertex of one part is connected with first vertex of next part.
It's strange especially as script deals well with multipart polygons and polygons with void.


Am I missing something with multipart polyline geometry?


import arcgisscripting

gp = arcgisscripting.create(9.3)

gp.overwriteoutput = True
gp.workspace = r"D:\Work\SwapXY.gdb"
inFC = r"D:\Work\SwapXY.gdb\Line"

desc = gp.Describe(inFC)
shapeField = desc.ShapeFieldName #SHAPE field name
shapeType = desc.ShapeType


# Create FC for inserting swaped geometries (if not exist)
...

try:
rows = gp.SearchCursor(inFC) #open SearchCursor on input FC
row = rows.next()
outRows = gp.InsertCursor(outFC) #open InsertCursor on output FC

#create Array object which will contain features vertices (other than points)
vertexArray = gp.CreateObject("Array")


while row:
feature = row.getValue(shapeField) #get the SHAPE field into variable

##For point geometry there is another way of reading coordinates than for polyline/polygon
if shapeType.upper() == "POINT" or shapeType.upper() == "MULTIPOINT":
#this part works right

else:
#feature can have multiple parts - first goes iteration through parts

partNum = 0
partCount = feature.PartCount
while partNum < partCount:
part = feature.GetPart(partNum) #the output is Array of points
pnt = part.next() #take first point from Array of points

#iterate through all points in array
while pnt:
#for each geometry create new POINT obj and assign swapped Y, X. Then add vertices to ARRAY
vertex = gp.CreateObject("Point")

vertex.X = pnt.Y
vertex.Y = pnt.X
vertexArray.add(vertex)
pnt = part.next()

#If pnt is null, either the part is finished or there is an interior ring
if not pnt:
pnt = part.next()
if pnt:
print "Interior:"


partNum += 1

newFeature = outRows.newRow()
newFeature.shape = vertexArray #assign ARRAY filled with points to shape field
outRows.insertRow(newFeature)

vertexArray.RemoveAll() #clear ARRAY before processing new geometry

row = rows.next()

except:
print gp.GetMessages()



In fact script doesn't work for polygons with more than one interior ring. So the problem is wider than I thought.
In polygon case, I've found that null point between exterior and interior rings should be inserted. But when gp.CreateObject("Point") is inserted, a point with coordinates (0,0) is added.



Answer



As adding null point as it was suggested here gave a point with (0,0) cordinates, I searched existing scripts in ESRI repository and I've found another way of separating polygon rings as well as polyline parts.
The solution is based on creating two arrays: featureVertexArray for storing whole geometry and partVertexArray for storing separate rings:


try:

rows = gp.SearchCursor(inFC) #open SearchCursor on input FC
row = rows.next()
outRows = gp.InsertCursor(outFC) #open InsertCursor on output FC

#create Array object which will contain features vertices
featureVertexArray = gp.CreateObject("Array")
partVertexArray = gp.CreateObject("Array")

while row:
feature = row.getValue(shapeField) #get the SHAPE field into variable

vertex = gp.CreateObject("Point") #empty Point object to store geometry

##For point geometry there is another way of reading coordinates than for polyline/polygon
if shapeType.upper() == "POINT" or shapeType.upper() == "MULTIPOINT":
#this part works right

else:
#feature can have multiple parts - first goes iteration through parts
partNum = 0
partCount = feature.PartCount

while partNum < partCount:
part = feature.GetPart(partNum) #the output is Array of points
pnt = part.next() #take first point from Array of points

#iterate through all points in array
while pnt:
#for each geometry create new POINT obj and assign swapped Y, X. Then add vertices to ARRAY
vertex.X, vertex.Y = pnt.Y, pnt.X
partVertexArray.add(vertex)
pnt = part.next()


#If pnt is null, either the part is finished or there is an interior ring
if not pnt:
pnt = part.next()
featureVertexArray.add(partVertexArray)
partVertexArray.removeAll()
#if pnt:
#print "Interior:"

featureVertexArray.add(partVertexArray)

partNum += 1 #increment part number to run loop

newFeature = outRows.newRow() #create new row in InsCur
newFeature.shape = featureVertexArray #assign ARRAY filled with points to shape field
outRows.insertRow(newFeature) #insert new row

featureVertexArray.removeAll() #clear ARRAYS before processing new geometry
partVertexArray.removeAll()

row = rows.next()

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...