Sunday, 15 February 2015

arcpy - Calculate area within Python script in ArcMap


I am trying to calculate the area of a polygon within my Python script. I create a new polygon from merging two together, and I'd like to add the area of the resulting polygon to a field in the output file. The polygon is stored in a regular shapefile, and is projected. Area preferably in map units.


I would have thought this to be a quite common and simple task, but despite a lot of Googleing I have been unable to find a working solutions so far.


I was planning on using arcpy.updateCursor to insert the value once it's calculated (there's only one feature in the FC at this stage), so easiest is if it can be returned as a variable. Any alternative solution that accomplishes the same task (getting the area value into correct field) will also work.



I have also tried the Field calculator from Python. Modified from the help pages I thought the following would work, but so far no luck.


arcpy.AddField_management(tempPgs, "Shape_area", 'DOUBLE')
exp = "float(!SHAPE.AREA!.split())"
arcpy.CalculateField_management(tempPgs, "Shape_area", exp)

Running ArcGIS Basic 10.1 SP1 with Python 2.7 on Windows 7.


Relevant parts of my current code looks like this:


#/.../
arcpy.Copy_management(inpgs, outpgs)
arcpy.AddField_management(outpgs, 'Shape_area', 'LONG')

fields = AM.FieldLst(outpgs)

#/.../

# Identify and search for shapes smaller than minimum area
where1 = '"' + 'Shape_Area' + '" < ' + str(msz)
polyrows = arcpy.SearchCursor(inpgs, where1)

for prow in polyrows:
grd1 = prow.GridID # GridID on the current polygon

grd2 = nDD.get(grd1) # GridID on the polygon downstream

# Update features
if grd2
geometry1 = prow.Shape
geometry2 = geometryDictionary[grd2]

# Update temporary features
arcpy.Merge_management([geometry1, geometry2], tempMerged)
arcpy.Dissolve_management(tempMerged, tempPgs)


fds = AM.FieldLst(tempPgs)

for field in fields[2:]:
arcpy.AddField_management(tempPgs, field, 'DOUBLE')

for fd in fds[2:]:
arcpy.DeleteField_management(tempPgs, fd)

exp = "float(!SHAPE.AREA!.split())"

arcpy.CalculateField_management(tempPgs, "Shape_area", exp)

# Append them to output FC
try:
arcpy.Append_management(tempPgs, outpgs, "TEST")
except arcgisscripting.ExecuteError:
arcpy.Append_management(tempPgs, outpgs, "NO_TEST")

elif ...


else ...

Answer



There are three different ways to find and store polygon area into a feature class with arcpy: 1) field calculator, 2) "classic" arcpy cursors, and 3) arcpy.da cursors. Some of this is borrowed from my previous answer about using SearchCursor.




1. Field calculator




  • When using field calculator, there are three different expression types that use different expression parsers. This is specified in the third parameter of the Calculate Field geoprocessing tool. When accessing the properties of the Geometry object using like in !shape.area!, you should use the Python 9.3 parser.





  • The expression you had before performed a split() command on the result of !SHAPE.AREA!. This returns a Python list object, which cannot be cast into a float object.




  • In your expression, you can specify the unit of the returned area by using the @SQUAREKILOMETERS flag, replacing SQUAREKILOMETERS with the units on the Calculate Field help page.




Here's the Python code I would use for this method:


tempPgs = "LayerName"
arcpy.AddField_management(tempPgs, "Shape_area", "DOUBLE")
exp = "!SHAPE.AREA@SQUAREKILOMETERS!"

arcpy.CalculateField_management(tempPgs, "Shape_area", exp, "PYTHON_9.3")



2. Arc 10.0 - "Classic" cursors




  • When using classic cursors (i.e. arcpy.UpdateCursor) the cursor object is an iterable containing row objects. You need to use the getValue and setValue methods to get the geometry from the row (as a geometry object and set the area value to the row as a float.




  • Your output row is stored in a temporary scratch space until you call the updateRow method on the cursor. This saves the new data into the actual dataset.





Here's the Python code I would use for this method:


tempPgs = "LayerName"
arcpy.AddField_management(tempPgs, "Shape_area", "DOUBLE")
geometryField = arcpy.Describe(tempPgs).shapeFieldName #Get name of geometry field
cursor = arcpy.UpdateCursor(tempPgs)
for row in cursor:
AreaValue = row.getValue(geometryField).area #Read area value as double
row.setValue("Shape_area",AreaValue) #Write area value to field

cursor.updateRow(row)
del row, cursor #Clean up cursor objects



3. Arc 10.1 - arcpy.da cursors




  • When using the new cursors in the data access module (i.e. arcpy.da.UpdateCursor) you need to pass in a list of field names as the second parameter in the cursor constructor. This requires some more work up-front, but the resulting row objects are Python lists, which makes it easier to read and write data when iterating through cursor rows. arcpy.da.UpdateCursor also has better performance than arcpy.UpdateCursor, partially because it skips unimportant fields, especially geometry.





  • When reading geometry, you can choose one of a number of geometry tokens, e.g. SHAPE@TRUECENTROID, SHAPE@AREA, or SHAPE@. Using a "simpler" token greatly improves performance compared to SHAPE@, which contains all of the geometry information. The full list of tokens is at the arcpy.da.UpdateCursor help page.




  • As before, your output row is stored in a temporary scratch space until you call the updateRow method on the cursor. This saves the new data into the actual dataset.




Here's the Python code I would use for this method:


tempPgs = "LayerName"
arcpy.AddField_management(tempPgs, "Shape_area", "DOUBLE")
CursorFieldNames = ["SHAPE@AREA","Shape_area"]

cursor = arcpy.da.UpdateCursor(tempPgs,CursorFieldNames)
for row in cursor:
AreaValue = row[0].area #Read area value as double
row[1] = AreaValue #Write area value to field
cursor.updateRow(row)
del row, cursor #Clean up cursor objects

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