Updating Z coordinates of points in polygons with ArcPy

Question

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])

Answer 1

Your case is very interesting to me as it reveals that explode_to_points=True is indeed "deconstruct(ing) a feature into its individual points or vertices." (as stated in the arcpy.da.UpdateCursor help).

HOWEVER

as in a Polygon, the first and last vertex of a feature (and all its parts) are coincidents, it seems that updating the geometry of the first vertex automatically updates that of the coincident (last) vertex. Thus, when you reach the last one which has already been shifted, it gets shifted again.

Let's visualyze this:

in ArcMap, I have a simple one feature, one part polygon. If I start an edit session and see its vertices, I can see it has 4 vertices.

But when I build a cursor and print all its same vertices, I can see there's one more at the end, which is coincident with the first one.

This was also already discussed in this post I just discovered.

Even if you avoid using the explode_to_points=True thing and go through each point in each part explicitly, you'll get the same result unfortunately.

Here you can see the same problem (code adapted from here):

This seems to me a real problem that should be handled differently by ArcPy.

SOLUTION

I guess we're left to avoid touching the first vertex of each part. Something like the same example above combined with this SO answer should work:

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))
            for i, pnt in enumerate(part):
                if i == 0:
                    arr_part.add(arcpy.Point(pnt.X, pnt.Y, pnt.Z))
                    continue
                if pnt:
                    print("OLD: {}, {}, {}".format(pnt.X, pnt.Y, pnt.Z))
                    new_z = pnt.Z
                    new_z += 2
                    print("NEW: {}, {}, {}".format(pnt.X, pnt.Y, new_z))
                    arr_part.add(arcpy.Point(pnt.X, pnt.Y, new_z))
                else:
                    print("Interior Ring:")
            # 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)

After running this on a single feature shapefile with two parts each with one donut I was able to see the updated Z correctly.

Answer 2

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)