need help creating multipart polygon output using arcpy

[erichfisher]

Hi all,

I am working on a script that transforms the vertex coordinates of input multi-part polygon shapefiles and outputs these transformed coordinate values into new multi-part polygons.  When I run the script, I am able to generate an identical shapefile to the original (but with transformed coordinate values), but none of the individual features remain.  I've tried to debug the code myself by observing the format of my main array that is fed into the arcpy.Polygon object and compared these results against arrays from working script examples provided by ESRI.  I don't see any difference, but clearly something is not right.  I was hoping someone could have a look below and recommend a solution.  Here's the code:

import arcpy
import math

infc = "C:\\testing\\fortesting.shp"

def transform(pntX,pntY):
    global Calc_north
    global Calc_east
    global Calc_elev

    start_north = -3770671.323
    start_east = -72160.677
    end_north = -3770648.576
    end_east = -72151.373

    delta_n = start_north - end_north
    delta_e = end_east - start_east
    angle = math.degrees(math.atan(delta_n/delta_e))

    deltaY = pntX*(math.sin(math.radians(angle)))
    deltaX = pntX*(math.cos(math.radians(angle)))

    Calc_north = start_north - deltaY
    Calc_east = deltaX + start_east
    Calc_elev = pntY

array = arcpy.Array()

# Enter for loop for each feature
for row in arcpy.da.SearchCursor(infc, ["OID@", "SHAPE@"]):
    # Print the current multipoint's ID
    print("Feature {}:".format(row[0]))
    # Step through each part of the feature
    for part in row[1]:
        # Step through each vertex in the part
        sub_array = arcpy.Array()
        for pnt in part:
            # Print x,y coordinates of current point
            transform(pnt.X,pnt.Y)
            #print pnt.X, pnt.Y, Calc_east, Calc_north, Calc_elev
            new_pnt = arcpy.Point(Calc_east, Calc_north)
            sub_array.add(new_pnt)
        sub_array.add(sub_array.getObject(0))
        array.add(sub_array)
        del sub_array
multi_Part_Polygon = arcpy.Polygon(array, arcpy.SpatialReference(102484), True)
arcpy.CopyFeatures_management(multi_Part_Polygon,"C:\\testing\\test18.shp")

[DarrenWiens2]

It's hard to tell where your indentation is, but I think you need to indent the line where you create the polygon object by one level (i.e. make one polygon for each row in your cursor). Add that polygon to a new array, and write that with CopyFeatures.

Are you trying to make the same number of multipart features you had originally, or one huge multipart feature containing all parts?

[IanMurray]

Hi Erich,

For readability of your script on the forum, please read the below post on posting code blocks.

https://community.esri.com/people/curtvprice/blog/2014/09/25/posting-code-blocks-in-the-new-geonet?s...‌

I think Darren is right, you are creating a single Polygon object which each of your multipart features is being written to.  You need to create a Polygon object for each polygon feature then use that for the input.

See the example below where they create a list of Polygon objects that are used as the input for the copy features.

import arcpy

# A list of features and coordinate pairs
feature_info = [[[1, 2], [2, 4], [3, 7]],
                [[6, 8], [5, 7], [7, 2], [9, 5]]]

# A list that will hold each of the Polygon objects
features = []

for feature in feature_info:
    # Create a Polygon object based on the array of points
    # Append to the list of Polygon objects
    features.append(
        arcpy.Polygon(
            arcpy.Array([arcpy.Point(*coords) for coords in feature])))

# Persist a copy of the Polyline objects using CopyFeatures
arcpy.CopyFeatures_management(features, "c:/geometry/polygons.shp")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[erichfisher]

import arcpy
import math

 

infc = "C:\\testing\\fortesting.shp"

 

def transform(pntX,pntY):
    global Calc_north
    global Calc_east
    global Calc_elev

 

    start_north = -3770671.323
    start_east = -72160.677
    end_north = -3770648.576
    end_east = -72151.373

 

    delta_n = start_north - end_north
    delta_e = end_east - start_east
    angle = math.degrees(math.atan(delta_n/delta_e))

 

    deltaY = pntX*(math.sin(math.radians(angle)))
    deltaX = pntX*(math.cos(math.radians(angle)))

 

    Calc_north = start_north - deltaY
    Calc_east = deltaX + start_east
    Calc_elev = pntY

 

array = arcpy.Array()

 

# Enter for loop for each feature
for row in arcpy.da.SearchCursor(infc, ["OID@", "SHAPE@"]):
    # Print the current multipoint's ID
    print("Feature {}:".format(row[0]))
    # Step through each part of the feature
    for part in row[1]:
        # Step through each vertex in the part
        sub_array = arcpy.Array()
        for pnt in part:
            # Print x,y coordinates of current point
            transform(pnt.X,pnt.Y)
            #print pnt.X, pnt.Y, Calc_east, Calc_north, Calc_elev
            new_pnt = arcpy.Point(Calc_east, Calc_north)
            sub_array.add(new_pnt)
        sub_array.add(sub_array.getObject(0))
        array.add(sub_array)
        del sub_array
multi_Part_Polygon = arcpy.Polygon(array, arcpy.SpatialReference(102484), True)
arcpy.CopyFeatures_management(multi_Part_Polygon,"C:\\testing\\test18.shp")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Hi everyone, thanks for the quick responses.  I did not know about the syntax highlighter.  Thanks for pointing that out. I have re-pasted the code again and I hope this is easier to read. 

See, I thought it was the creation of the arcpy.Polygon object too and I have tested the script by indenting the code by one level (so it is part of the "for row in arcpy.da.Searchcursor") and two levels (part of the "for part in row[1]") and neither works. 

I was wondering if there was something wrong with my input data?  the example that you provide Ian, for example, was one I used for reference and other like it all rely on list inputs.  My values are drawn from an existing shapefile, transformed in the function, and then returned as numeric values.  Would that make any difference?

Thanks for all the help,

Erich

[DarrenWiens2]

You can just update the geometry within an UpdateCursor.

[IanMurray · Accepted solution]

I think you missed the point.  You need to create a Polygon Object for each row, not for the entire feature class.  Right now you are creating only a single polygon object, hence the single polygon output you are getting of all features together.  You need to create a list of Polygon objects as the input for Copy Features so you will have a record for each original feature.  This is why in the example, features was the input for Copy Features, not a single Polygon Object, since there were two Polygons created in the example script for each list of coordinates.

I believe if you adapted it to something like this it would take care of things, sorry about the gratuitous whitespace.

features = []

# Enter for loop for each feature
for row in arcpy.da.SearchCursor(infc, ["OID@", "SHAPE@"]):
    # Print the current multipoint's ID
    print("Feature {}:".format(row[0]))

    polygon = arcpy.Polygon()

    array = arcpy.Array()
    # Step through each part of the feature
    for part in row[1]:
        # Step through each vertex in the part
        sub_array = arcpy.Array()
        for pnt in part:
            # Print x,y coordinates of current point
            transform(pnt.X,pnt.Y)
            #print pnt.X, pnt.Y, Calc_east, Calc_north, Calc_elev
            new_pnt = arcpy.Point(Calc_east, Calc_north)
            sub_array.add(new_pnt)
        sub_array.add(sub_array.getObject(0))
        array.add(sub_array)
        del sub_array

  polygon.add(array)

  features.append(polygon)

  del array

  del polygon

#multi_Part_Polygon = arcpy.Polygon(array, arcpy.SpatialReference(102484), True)
arcpy.CopyFeatures_management(features,"C:\\testing\\test18.shp")
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[erichfisher]

Thanks so much for the explanation. I was able to generate the individual parts in my output now.