Arcpy (2.6) modify multiparts polygon vertice xy vlaue

2235
17
Jump to solution
08-30-2016 02:14 AM
tinyMee
New Contributor

Dear all,

I am using python2.6 (together with arcgis 10.0) and try to modify the vertice xy values for multiple parts polygon using the following scripts, no error thrown but the vertice xy value not changed. Please help to advice what is wrong in my script:

def modify(f, n):
     s = '%.12f' % f
     i, p, d = s.partition('.')
     return '.'.join([i, (d+'0'*n)[:n]])

import arcpy
import logging
from arcpy import env

env.workspace = r"C:\modifyTest"
infc = "testing.mdb/donutTesting1"
desc = arcpy.Describe(infc)
shapefieldname = desc.ShapeFieldName
my_field = "SHAPE"
rows = arcpy.UpdateCursor(infc) 
for row in rows:
     feat = row.getValue(shapefieldname) 
     arr_pol = arcpy.Array()
     partnum = 0
     if feat.isMultipart == True: 
          for part in feat:
               parts = []
               arr_parts = arcpy.Array()
 
               while partnum < feat.partCount:
                    part = feat.getPart(partnum)
                    pnt = part.next()
                    while pnt: 
                         X_value = modify( pnt.X , 2)
                         Y_value = modify( pnt.Y , 2)
                         pnt.X= X_value
                         pnt.Y= Y_value
                         #print pnt.X, pnt.Y
                         parts.append([X_value, Y_value])
                         XYPoint = arcpy.Point(X_value, Y_value)
                         arr_parts.add(XYPoint) 
                         arr_pol.add(arr_parts)
 
                         pnt = part.next()
                         if pnt is None:
                              pnt = part.next()
                              if pnt:
                                   X_value = modify( pnt.X , 2)
                                   Y_value = modify( pnt.Y , 2)
                                   pnt.X= X_value
                                   pnt.Y= Y_value
                                   parts.append([X_value, Y_value])
                                   XYPoint = arcpy.Point(X_value, Y_value)
                                   arr_parts.add(XYPoint) 
                                   arr_pol.add(arr_parts) 
                    partnum += 1
              print parts
      polygon = arcpy.Polygon(arr_pol)
      row.setValue(my_field,polygon)
      rows.updateRow(row)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
Tags (2)
0 Kudos
17 Replies
DanPatterson_Retired
MVP Emeritus

how are you assessing that it has not changed?  

visually?      (1 mm shift isn't going to show on a map)  

Within the table?     (coordinates values in tables aren't updated automagically)

If you are getting no errors and you have fixed the change method to do proper rounding, how are you assessing whether a point has been rounded (ie shifted) or not?  and how are you then assessing whether it doesn't work for multiparts?  And of course this isn't in a geodatabase or some other container that controls coordinates.  Did you try a shapefile whose coordinates are solely controlled by it

0 Kudos
tinyMee
New Contributor

Hi Dan,

The script from stand alone, it means I run the script without open arcmap. After running the script I manually check the vertex xy value in editing session.

I will test on a shapefile as you suggested, and update the result. Thank you

0 Kudos
DanPatterson_Retired
MVP Emeritus

Checking in an edit session isn't what is going to be what is shown if you calculate the X, Y coordinates and throw them into the table using the Add XY coordinates tool.  Compare the shapefile results vs that of your original file

0 Kudos
XanderBakker
Esri Esteemed Contributor

Apart from the indentation there are some funcional errors in your script. You only process the multipart polygons, while the final polygon is created from a arr_pol, which will be empty when the feature is not a multipart. You should process all features the same way regardless of it being multi or single part (all  the features in your sample data have only 1 part). If you don't do that and you write the original polygons when it is single part, you will get this:

And yes rounding the coordinates on 2 decimales is visible in you data. I didn't your code, but the da cursor (ArcGIS 10.1 and higher) which is far more efficient. 

Code used (to produce the "problem" above):

def main():
    import arcpy

    fc = r'D:\Xander\GeoNet\Donuts\shp\donutTesting.shp'
    fc_out = r'D:\Xander\GeoNet\Donuts\shp\donutTesting_out.shp'
    sr = arcpy.Describe(fc).spatialReference

    lst_feats = []
    with arcpy.da.UpdateCursor(fc, ('SHAPE@', 'OID@')) as curs:
        for row in curs:
            polygon = row[0]
            oid = row[1]
            print oid
            if polygon.isMultipart:
                print "multipart, # of parts:", polygon.partCount
                lst_pol = []
                for i in range(polygon.partCount):
                    part = polygon.getPart(i)
                    print" - part", i
                    lst_part = []
                    for pnt in part:
                        pnt2 = roundCoords(pnt, 2)
                        print "   - pnt", pnt, pnt2
                        lst_part.append(pnt2)
                    lst_pol.append(lst_part)
                polygon2 = arcpy.Polygon(arcpy.Array(lst_pol), sr)
                lst_feats.append(polygon2)
            else:
                print "singlepart, # of parts:", polygon.partCount
                lst_feats.append(polygon)

    arcpy.CopyFeatures_management(lst_feats, fc_out)


def roundCoords(pnt, decimals):
    if pnt:
        return arcpy.Point(round(pnt.X, decimals), round(pnt.Y, decimals))
    else:
        return pnt

if __name__ == '__main__':
    main()

Output messages:

0
singlepart, # of parts: 1
1
multipart, # of parts: 1
 - part 0
   - pnt 48371,4859999996 42834,2699999996 NaN NaN 48371,49 42834,27 NaN NaN
   - pnt 48372,1289999997 42966,5089999996 NaN NaN 48372,13 42966,51 NaN NaN
   - pnt 48371,9040000001 42977,0130000003 NaN NaN 48371,9 42977,01 NaN NaN
   - pnt 48400,7130000005 42992,7489999998 NaN NaN 48400,71 42992,75 NaN NaN
   - pnt 48620,4069999997 42952,8029999994 NaN NaN 48620,41 42952,8 NaN NaN
   - pnt 48627,1799999997 42891,8430000003 NaN NaN 48627,18 42891,84 NaN NaN
   - pnt 48481,5530000003 42841,0429999996 NaN NaN 48481,55 42841,04 NaN NaN
   - pnt 48371,4859999996 42834,2699999996 NaN NaN 48371,49 42834,27 NaN NaN
   - pnt None None
   - pnt 48435,8629999999 42969,7170000002 NaN NaN 48435,86 42969,72 NaN NaN
   - pnt 48466,034 42961,6860000007 NaN NaN 48466,03 42961,69 NaN NaN
   - pnt 48468,5360000003 42963,7890000008 NaN NaN 48468,54 42963,79 NaN NaN
   - pnt 48470,3729999997 42965,7080000006 NaN NaN 48470,37 42965,71 NaN NaN
   - pnt 48435,8629999999 42969,7170000002 NaN NaN 48435,86 42969,72 NaN NaN
   - pnt None None
   - pnt 48391,9639999997 42944,5050000008 NaN NaN 48391,96 42944,51 NaN NaN
   - pnt 48408,0039999997 42942,8200000003 NaN NaN 48408 42942,82 NaN NaN
   - pnt 48416,9890000001 42950,9759999998 NaN NaN 48416,99 42950,98 NaN NaN
   - pnt 48391,4349999996 42957,2050000001 NaN NaN 48391,43 42957,21 NaN NaN
   - pnt 48391,9639999997 42944,5050000008 NaN NaN 48391,96 42944,51 NaN NaN
2
singlepart, # of parts: 1

Doing this (the right way) in ArcGIS 10.0, you could use this code:

def main():
    import arcpy

    fc = r'D:\Xander\GeoNet\Donuts\shp\donutTesting.shp'
    fc_out = r'D:\Xander\GeoNet\Donuts\shp\donutTesting_out2.shp'
    sr = arcpy.Describe(fc).spatialReference
    fld_shp = arcpy.Describe(fc).ShapeFieldName
    fld_oid = arcpy.Describe(fc).OIDFieldName

    curs = arcpy.UpdateCursor(fc)
    for row in curs:
        polygon = row.getValue(fld_shp)
        oid = row.getValue(fld_oid)
        print "oid", oid
        lst_pol = []
        for i in range(polygon.partCount):
            part = polygon.getPart(i)
            print" - part", i
            lst_part = []
            for pnt in part:
                pnt2 = roundCoords(pnt, 2)
                print "   - pnt", pnt, pnt2
                lst_part.append(pnt2)
            lst_pol.append(lst_part)
        polygon2 = arcpy.Polygon(arcpy.Array(lst_pol), sr)
        row.setValue(fld_shp, polygon2)
        curs.updateRow(row)
    del curs, row

def roundCoords(pnt, decimals):
    if pnt:
        return arcpy.Point(round(pnt.X, decimals), round(pnt.Y, decimals))
    else:
        return pnt

if __name__ == '__main__':
    main()

Zoom in to coordinates 48391.435  42957.205 Meters to see the effect. 

tinyMee
New Contributor

Thanks so much for your help Dan and Xander, I use the 2nd script (yes this is the right way) but I get runtime error when process the multiple parts polygon:

oid 0
- part 0
- pnt 48408.0001220703 42942.8201293945 NaN NaN 48408 42942.82 NaN NaN
- pnt 48391.9600830078 42944.5001220703 NaN NaN 48391.96 42944.5 NaN NaN
- pnt 48391.4301147461 42957.2000732422 NaN NaN 48391.43 42957.2 NaN NaN
- pnt 48416.9801025391 42950.9700927734 NaN NaN 48416.98 42950.97 NaN NaN
- pnt 48408.0001220703 42942.8201293945 NaN NaN 48408 42942.82 NaN NaN
oid 1
Traceback (most recent call last):
File "C:\Script\dobut.py", line 38, in <module>
main()
File "C:\Script\dobut.py", line 26, in main
polygon2 = arcpy.Polygon(arcpy.Array(lst_pol), sr)
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\arcobjects\mixins.
py", line 196, in dunu
*gp_fixargs(args, True))
File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\geoprocessing\_bas
e.py", line 474, in <lambda>
return lambda *args: val(*gp_fixargs(args))
RuntimeError: Object: CreateObject cannot create geometry from inputs

I run the python file through batch file. 

0 Kudos
DanPatterson_Retired
MVP Emeritus

you changed the file names?
did you try it with one file so named in the script?
what do you mean a batch file? the script isn't set up to be batched unless you have made modifications to it or your batch file provides the right parameters

XanderBakker
Esri Esteemed Contributor

I see that the data is different from the polygons you shared. Could you share this data to see what is different? And please clarify (as Dan already asked) what you mean with running the python file through a batch file (and why you want to do that)? 

0 Kudos
tinyMee
New Contributor

Hi Xander, I run the python file from a .bat file, the content in the .bat file is "C:\Python26\Arcgis10.0\python.exe C:\testing\donut.py" since I need add to the batch file in a scheduled job. It seems it is not the correct way?

Now I am using da cursor (ArcGIS 10.1) with the code shared by you, it works properly. Thank you.

0 Kudos