Common Python geometry calculations in the Field Calculator

[wilfredwaters]

Here it has Common VBA geometry calculations in the Field Calculator, provide as a number of code blocks. I have pasted these below and am wondering if someone can make a set of equivalent code blocks in Python, or can point me where someone has already made them?

Common VBA geometry calculations in the Field Calculator
Below are some code samples you can use to perform geometric calculations with VBA statements, rather than having the Calculate Geometry command do them for you.
To use the samples, open the Field Calculator, check Advanced, and type the VBA statement in the first text box. Then, type the variable name in the text box directly under the field name.

Area

Dim dblArea as double
Dim pArea as IArea
Set pArea = [shape]
dblArea = pArea.area

Perimeter

Dim dblPerimeter as double
Dim pCurve as ICurve
Set pCurve = [shape]
dblPerimeter = pCurve.Length

Length

Dim dblLength as double
Dim pCurve as ICurve
Set pCurve = [shape]
dblLength = pCurve.Length

X-coordinate of a point

Dim dblX As Double
Dim pPoint As IPoint
Set pPoint = [Shape]
dblX = pPoint.X

Y-coordinate of a point

Dim dblY As Double
Dim pPoint As IPoint
Set pPoint = [Shape]
dblY = pPoint.Y

X-coordinate of a polygon centroid

Dim dblX As Double
Dim pArea As IArea
Set pArea = [Shape]
dblX = pArea.Centroid.X

Y-coordinate of a polygon centroid

Dim dblY As Double
Dim pArea As IArea
Set pArea = [Shape]
dblY = pArea.Centroid.Y

To add the z-value of the start point of a polyline

Dim dblZ As Double
Dim pLine As IPolyline
Dim pPoint as IPoint
Set pLine = [Shape]
Set pPoint = pLine.FromPoint
dblZ= pPoint.Z

[wilfredwaters]

area:

!shape.area@hectares!

(or whatever other unit is required)

[DanPatterson_Retired]

did you rule out those contained in the help files?
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Calculate_Field_examples/005s0000002m0...

[LoganPugh]

EasyCalculate is a free set of 100+ field calculator scripts, definitely check them out: http://www.ian-ko.com/free/free_arcgis.htm

[wilfredwaters]

EasyCalculate is a free set of 100+ field calculator scripts, definitely check them out: http://www.ian-ko.com/free/free_arcgis.htm

Wow this looks stupendiforously useful, however I note it only mentions up to Arc 9.x in the documentation, not Arc 10. Will there be compatibility issues?

Thanks,

Wil

[wilfredwaters]

Ianko responded to my email with this link for arc 10:

http://www.ian-ko.com/free/EC10/EC10_main.htm

[RichardFairhurst]

ArcGIS 10 supports all of the common field calculator geometry expressions in Python as simple expressions (no advanced calculations necessary, unlike VBA).  Here they are:

Area

!shape.area!

Perimeter

!shape.length! (for polygons length returns the perimeter length)

Length

!shape.length! (for polylines length returns the line length)

X-coordinate of a point

!shape.firstpoint.X!

Y-coordinate of a point

!shape.firstpoint.Y!

X-coordinate of a polygon centroid

!shape.centroid.X! (works for polyline too, but is not necessarily the midpoint, see midpoint/truecentroid below)

Y-coordinate of a polygon centroid

!shape.centroid.Y! (works for polyline too, but is not necessarily the midpoint, see midpoint/truecentroid below)

To add the z-value of the start point of a polyline

!shape.firstpoint.Z!

To add the z-value of the end point of a polyline

!shape.lastpoint.Z!

To add the z-value of the mid point of a polyline

!shape.truecentroid.Z!

All of the above calculations can get the X, Y, Z and M values of a geometry (if Z and M are enabled) by substituting the appropriate letter at the end of any of the expressions that return a coordinate.

Using the expression pattern below units can be converted in one expression to any of the following units for area:  ACRES | ARES | HECTARES | SQUARECENTIMETERS | SQUAREDECIMETERS | SQUAREINCHES | SQUAREFEET | SQUAREKILOMETERS | SQUAREMETERS | SQUAREMILES | SQUAREMILLIMETERS | SQUAREYARDS | SQUAREMAPUNITS | UNKNOWN (do not use with geographic coordinate systems that use angular units as the results will be questionable)

!shape.area@squarekilometers!

Using the expression pattern below units can be converted in one expression to any of the following units for length:  CENTIMETERS | DECIMALDEGREES | DECIMETERS | FEET | INCHES | KILOMETERS | METERS | MILES | MILLIMETERS | NAUTICALMILES | POINTS | UNKNOWN | YARDS (geodesic conversions occur for geographic coordinate systems that use angular units)

!shape.length@kilometers!

See the Working with geometry in Python help for calculating other geometry properties, such as extent, hullRectangle, isMultipart, partCount, pointCount, labelPoint, and type.

[AnnStark]

Thank you rfairhur24.   

This is getting me in the right direction , however I'm still stumped as to how I can calculate the lat/long of a point.   

!shape.length@DECIMALDEGREES!  returns zero.
!shape.firstpoint.X! returns the location in my State plane coordinates.
!shape.firstpoint.X@DECIMALDEGREES! doesn't work .  (just taking a stab at things)

Any one know how to do this?
A couple clicks in the desktop software does it easily, but I need to incorporate this functionality into a python script that will do this each night with new points.

Ann

[BryanTaylor]

I'm having the same problem Ann (except I think mine is worse.  I've lost access to field calculator and geometry calculator.).  The only thing I've found as a possible answer is that you have to reproject your stateplane layer to a geographic coordinate system like wgs84 and then add the x/y table values for lat/long.  Depending on how you want to use the data (I want lat/longs but within a stateplane projected layer), you will have to extract the field data and write it back into your stateplane layer record by record, or reproject the geographic layer back to stateplane.  It's a pain in the rear for sure.  I can't believe it's this hard to do something so basic.  And I hope nobody responds with "just use the x/y data tool",  or "use the geometry calculator", because neither of them work properly, nor do they throw errors or give any clue as to why they don't work properly.  I hope my 2 cents might help.  I'll keep checking back to see if a more viable solution arrives that I can use also.....I'd hate to have to write my own tool, when ESRI has one they need to fix.

[RichardFairhurst]

The only way I have found to avoid using the Project tool to do a conversion from a spatial projection that uses linear units (like State Plane coordinates) to angular units (like Decimal Degrees) is using the Geoprocessing Environment settings to set up any transformation option and the spatial reference parameter of a Python UpdateCursor.  If you have ArcGIS Desktop 10.1 or later you should only use a Data Access UpdateCursor (da UpdateCursor), not the original Python UpdateCursor syntax.

A note of warning, if you use this approach with an XY Event layer that has its original coordinates in linear units, like State Plane coordinates, you need to first duplicate those coordinates into a pair of new fields that you want to hold your Latitude (Y coordinate) and Longitude (X coordinate) and use those new fields with the XY Event Layer tool, because if you include the SHAPE@X or SHAPE@Y feilds in the UpdateCursor field list and use the updateRow method, the fields used to set the coordinates of the XY Event Layer will convert and be overwritten with new coordinate values that match the new spatial reference and transformation.

The script below can create a new Latitude and Longitude field in an XY Event Layer table that has only State Plane coordinates, make a copy of the original state plane coordinates into the new Latitude and Longitude fields, and use an XY Event layer and a da UpdateCursor opened with a different spatial reference to convert the values in the Latitude and Longitude fields into the coordinates matching the new spatial reference.  It can run the entire script to process over 131K records in an XY Event Layer in approximately 36 seconds, which is about 5 to 10 times faster than methods that used to use the Field Calculator.

# Import arcpy module
from time import strftime 

print "Start script: " + strftime("%Y-%m-%d %H:%M:%S")

import arcpy

# Set or comment out Geographic Transformation
arcpy.env.geographicTransformations = "NAD_1983_To_WGS_1984_5"

# Use a table with State Plane NAD 83 coordinates
fc = r"C:\Users\RFAIRHUR\Documents\ArcGIS\DEFAULT.gdb\CL_INTERSECTIONS_PAIRS"

print "Finished Setup: " + strftime("%Y-%m-%d %H:%M:%S")

# Process: Add Latitude Field
arcpy.AddField_management(fc, "LATITUDE", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Add Longitude Field
arcpy.AddField_management(fc, "LONGITUDE", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

print "Finish Adding Fields: " + strftime("%Y-%m-%d %H:%M:%S")

# Specify original X and Y fields and corresponding Longitude and Latitude output fields.
fields = ( 'LONGITUDE', 'X_COORD', 'LATITUDE', 'Y_COORD')

#Open a da UpdateCursor and duplicate the X and Y coordinates in the Longitude and Latitude fields.
with arcpy.da.UpdateCursor(fc, fields) as updateRows:
    for updateRow in updateRows:
        updateRow[0] = updateRow[1]
        updateRow[2] = updateRow[3]
        updateRows.updateRow(updateRow)

print "Finished UpdateCursor Longitude and Latitude edits " + strftime("%Y-%m-%d %H:%M:%S")

# Set a name for your XY Event Layer
CL_INTERSECTIONS_PAIRS_Events = "CL_INTERSECTIONS_PAIRS Events"

# Make an X/Y Event Layer with NAD 83 projection.
arcpy.MakeXYEventLayer_management(fc, "LONGITUDE", "LATITUDE", CL_INTERSECTIONS_PAIRS_Events, "PROJCS['NAD_1983_StatePlane_California_VI_FIPS_0406_Feet',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Lambert_Conformal_Conic'],PARAMETER['False_Easting',6561666.666666666],PARAMETER['False_Northing',1640416.666666667],PARAMETER['Central_Meridian',-116.25],PARAMETER['Standard_Parallel_1',32.78333333333333],PARAMETER['Standard_Parallel_2',33.88333333333333],PARAMETER['Latitude_Of_Origin',32.16666666666666],UNIT['Foot_US',0.3048006096012192]];-118608900 -91259500 3048.00609601219;-100000 10000;-100000 10000;3.28083333333333E-03;0.001;0.001;IsHighPrecision", "")

print "Finished Making XY Event Layer: " + strftime("%Y-%m-%d %H:%M:%S")

# Use the Event layer as the cursor input
fcl = CL_INTERSECTIONS_PAIRS_Events

# Specify shape coordinate fields.
fields = ("SHAPE@X", "SHAPE@Y")

# Specify the Spatial Reference of the Cursor
latLonRef = "Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"

# Create update cursor for feature class with a new Spatial Reference.
with arcpy.da.UpdateCursor(fcl, fields, ' ', latLonRef) as updateRows:
    for updateRow in updateRows:
        # Lat and Long convert just by using the updateRow method on the Shape fields
        updateRows.updateRow(updateRow)

print "Finished Converting Longitude and Latitude fields: " + strftime("%Y-%m-%d %H:%M:%S")

# The coordinates in the Longitude and Latitude fields are now converted to Decimal Degrees.
print "Finish script: " + strftime("%Y-%m-%d %H:%M:%S")

Below is a script that adds a Latitude and Longitude field to a State Plane NAD 83 Point Feature Class and using an UpdateCursor with a different spatial reference is able to fill those fields in with the Decimal Degree coordinates of a WGS 1984 spatial reference.  The underly geometry of the points is unaffected by the script.  The entire script takes approximately 25 seconds to run on over 131K points.

from time import strftime  

print "Start script: " + strftime("%Y-%m-%d %H:%M:%S")

import arcpy

# Set or comment out Geographic Transformation
arcpy.env.geographicTransformations = "NAD_1983_To_WGS_1984_5"

# Use a table with State Plane NAD 83 coordinates
fc = r"C:\Users\RFAIRHUR\Documents\ArcGIS\DEFAULT.gdb\CL_INT_PAIRS_POINTS"

print "Finished Setup: " + strftime("%Y-%m-%d %H:%M:%S")

# Process: Add Latitude Field
arcpy.AddField_management(fc, "LATITUDE", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

# Process: Add Longitude Field
arcpy.AddField_management(fc, "LONGITUDE", "DOUBLE", "", "", "", "", "NULLABLE", "NON_REQUIRED", "")

print "Finish Adding Fields: " + strftime("%Y-%m-%d %H:%M:%S")

# Specify shape coordinate fields.
fields = ("LONGITUDE", "SHAPE@X", "LATITUDE", "SHAPE@Y")

# Specify the Spatial Reference of the Cursor
latLonRef = "Coordinate Systems\Geographic Coordinate Systems\World\WGS 1984.prj"

# Create update cursor for feature class with a new Spatial Reference.
with arcpy.da.UpdateCursor(fc, fields, ' ', latLonRef) as updateRows:
    for updateRow in updateRows:
        updateRow[0] = updateRow[1]
        updateRow[2] = updateRow[3]
        updateRows.updateRow(updateRow)

print "Finished Converting Longitude and Latitude fields: " + strftime("%Y-%m-%d %H:%M:%S")

# The coordinates in the Longitude and Latitude fields are now converted to Decimal Degrees.
print "Finish script: " + strftime("%Y-%m-%d %H:%M:%S")