Edit 2017:
After all the below work, I've come to learn that there's a much easier way to determining feature lengths regardless of spatial reference with use of the field calcutor.
This blog explains how.
Here's an example:
arcpy.CalculateField_management ("InputLayer", "LenMiles", "!shape.length@miles!", "PYTHON_9.3")
Original Blog:
When working with GIS, knowing lengths of polyline features can be fundamental to accomplishing your goals. However, something so basic isn't necessarily always simple. I wanted to create a python function that would take any feature class and return each feature's length in a standard unit (meters). I wanted the code to work for any feature class, regardless of projection or how large or small an area of the globe the feature class's features covered.
Determining lengths of features in a projected feature class is easy. ArcPy allows you to do so by the feature class's
spatial reference object. This object has the metersPerUnit method, which returns the conversion factor for the feature class's units to meters. Once this is known, a
cursor can be used to obtain the length of each feature. This distance can be multiplied by the conversion factor for the standardized distance.
Geographic coordinate systems do not have a measure unit, so determining each feature's length is not as straight forward a task. Initially I couldn't find a solution, and so asked
this question on GIS StackExchange. The answer presented the solution of setting the feature class's coordinate system within a cursor. Picking which coordinate system also wasn't easy. The answer suggested WGS Web Mercator. However, one of the commenters warned that WGS Web Mercator isn't actually in meters. He suggested using the UTM zone the data fell in as the projection. This was a good suggestion, except for the fact that I wanted my code to be robust enough to work regardless of where it was geographically. I had to find a way to determine which UTM zone each feature fell within. With a little poking around I came across
Calculate UTM Zone. This tool appears to be designed for use with data driven pages, but I found it suited my needs relatively well. The output, which populates a field with WKT, wasn't perfect however. I asked
this question for suggestions on how I could determine each feature's UTM zone from its output consistently, but didn't find an answer. Initially I came up with my own string manipulation method, but later discovered the loadFromString () method of the spatial reference object, that lets you update a spatial reference from well known text. From here, it was simply a matter of projecting each feature into its corresponding UTM zone and creating the dictionary.
The code:
def FeatureLengths (inFc):
"""
input - inFc: the input feature class
returns dictionary -
key: feature OID
value: distance in meters
"""
#get spatial reference type
sr = arcpy.Describe (inFc).spatialReference
srType = sr.type
#Projected coordinate systems
if srType == "Projected":
#meters per unit
metersPerUnit = sr.metersPerUnit
#out Di
return dict ([(oid, dist * metersPerUnit) for
oid, dist in
arcpy.da.SearchCursor
(inFc, ["OID@", "SHAPE@LENGTH"])])
#Geographic coordinate systems
#add UTM field to feature class
flds = [f.name for f in arcpy.ListFields (inFc)]
utmFld = "UTM"
i = 0
#Make sure field name is not in use
while utmFld in flds:
utmFld = "UTM{0}".format (i)
i += 1
arcpy.AddField_management (inFc,
utmFld,
"TEXT",
field_length = 1000)
#calculate UTM zone
arcpy.CalculateUTMZone_cartography (inFc, utmFld)
#add sql field delimiters
delimUtmFld = arcpy.AddFieldDelimiters (inFc, utmFld)
#output dictionary that will hold oids and distances in meters
outDi = {}
#list to hold utm zones
utms = []
#get unique layer name
layer = "lyr"
while arcpy.Exists (layer):
i += 1
layer = "lyr" + str (i)
#access table
with arcpy.da.SearchCursor (inFc, [utmFld]) as cursor:
for utmWkt, in cursor:
#check if utm has been iterated
if utmWkt in utms:
continue
utms += [utmWkt]
#create new spatial reference
sr = arcpy.SpatialReference ()
#update spatial reference with wkt
sr.loadFromString (utmWkt)
#create sql to filter for utm zone
sql = "{0} = '{1}'".format (delimUtmFld, utmWkt)
#create feature layer per utm zone
arcpy.MakeFeatureLayer_management (inFc, layer, sql)
#iterate layer with projection
with arcpy.da.SearchCursor (layer,
["OID@", "SHAPE@LENGTH"],
spatial_reference
= sr) as lyrCursor:
for oid, shpLen in lyrCursor:
#update dictionary with length
if shpLen:
outDi [oid] = shpLen
#if shape error, return zero
else:
outDi [oid] = 0
del lyrCursor
#Delete layer
arcpy.Delete_management (layer)
del cursor
#delete utm zone field
arcpy.DeleteField_management (inFc, utmFld)
#return out dictionary
return outDi
