Saturday 28 July 2018

geometry - Updating Z coordinates of points in polygons with ArcPy


I am trying to update Z values of polygon shape with ArcPy. This script is working without error but increasing the first (0 position vertex) vertex twice. I want to update only once. How can I avoid it?


import arcpy


feature = r"D:\acrpy\VV_test\scripts\map1\ZM_KAD_00.shp" #polygon shape
fields = ["ESRI_OID", "SHAPE@Z"]
whereclause = """{0} = 7""".format(arcpy.AddFieldDelimiters(feature, fields[0]))
z_increase = 2

with arcpy.da.UpdateCursor(feature, fields, whereclause, explode_to_points=True) as cursor:
for row in cursor:
print row[1]
cursor.updateRow([row[0],row[1] + z_increase])


enter image description here



Answer



I wanted to insert a new answer rather than editing my previous one as this comes with a new solution that should handle the problematic of updating the geometries of polygons with more than one inner ring as well.


The OP in the comments pointed out that my other answer didn't work for the case "two donuts in one polygon".


Trying myself I was susprised but today I realized this problem has already come up previously in this other question where a solution is also provided.


I'll highlight the key part that makes the code you are about to see working with all types of polygons (single-to-multipart, zero-to-many inner rings):



  • an arcpy.Polygon() is made up of one or more arcpy.Array()s

  • each arcpy.Array() is constructed with at least three arcpy.Point()s (although they require four per definition, as pointed out by @Vince in the comments; yes, in ArcGIS a "Point" means a "Vertex", quite confusing...)


  • a single arcpy.Polygon() with more than one arcpy.Array() inside is a multipart polygon

  • a single arcpy.Polygon() with one or more interior rings (alias "donuts"), is NOT a multipart polygon but a singlepart polygon with interior ring(s) (quoting this)


  • a single part arcpy.Polygon() with one or more interior rings is composed by a single arcpy.Array() with a None value where the inner ring(s) begin(s), for example:


    [[x1, y1], [x2, y2], [x3, y3], [x4, y4], None, # start of the inner ring coordinatres [x5, y5], [x6, y6], [x7, y7], [x8, y8]]




  • Last but not least, arcpy.cursors list per each ring (not part, here I mean ring), a starting vertex and an end vertex which are equivalent in terms of coorindates.





  • As per the previous assessment, one needs to sikp updating (shifting) either the first or last point otherwise the coincident point is shifted twice.




The modified working code of my previous original answer is shown below:


fields = ["OID@", "SHAPE@"]
feature = r"D:\acrpy\VV_test\scripts\map1\ZM_KAD_00.shp" #polygon shape

fields = ["ESRI_OID", "SHAPE@"]

whereclause = """{0} = 7""".format(arcpy.AddFieldDelimiters(feature, fields[0]))


z_increase = 2

desc = arcpy.Describe(feature)
sr = desc.spatialReference


with arcpy.da.UpdateCursor(feature, fields, whereclause) as cursor:
for row in cursor:
print("Feature {}:".format(row[0]))

partnum = 0
# create an array for the updated feature
arr_feat = arcpy.Array()
for part in row[1]:
# create an array for the part
arr_part = arcpy.Array()
print("Part {}:".format(partnum))
i = 0
for pnt in part:
if i == 0:

print("SKIPPING FIRST VERTEX")
arr_part.add(arcpy.Point(pnt.X, pnt.Y, pnt.Z))
i += 1
continue
if pnt:
print("OLD: {}, {}, {}".format(pnt.X, pnt.Y, pnt.Z))
new_z = pnt.Z
new_z += z_increase
i += 1
print("NEW: {}, {}, {}".format(pnt.X, pnt.Y, new_z))

arr_part.add(arcpy.Point(pnt.X, pnt.Y, new_z))
else:
print("Interior Ring:")
null_point = None
arr_part.add(null_point) # MUST DO THIS WHEN INSERTING INTERIOR RING!!!
i = 0
# add part to feature array
arr_feat.add(arr_part)
partnum += 1
polygon = arcpy.Polygon(arr_feat, sr, True)


row[1] = polygon
cursor.updateRow(row)

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