Determining Lengths of Polyline Features Regardless of Spatial Reference using ArcPy

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