I have the following ArcPy script to calculate total length of Polylines in a Feature Class/Shapefile:

arcpy.Project_management('lines.shp', 'lines_projected.shp', arcpy.SpatialReference(TODO EPSG CODE))
arcpy.MakeFeatureLayer_management('lines_projected.shp', 'lines_layer')
geometries = arcpy.CopyFeatures_management('lines_layer', arcpy.Geometry())
print 'Number of lines: ' + str(len(geometries))
total_length = sum([g.length for g in geometries])
print 'Total length: ' + str(total_length/1000) + 'km'

What is the methodology to find appropriate Coordinate Reference System assuming the data associated with well known region (For example "Berlin" or "France")?

For example why a simple approach of choosing the appropriate UTM zone out of the 60 won't work

  • 1
    Berlin and France are very different things and the sum of the length of lines could vary greatly depending on the scale at which they were captured, in addition to varying by coordinate reference. This question is quite broad, certainly broader than could be answered in a paragraph or three, which makes it a poor fit for our "Focused question / Best answer" model. – Vince Mar 2 '18 at 16:45
  • According to this: qgistutorials.com/en/docs/calculating_line_lengths.html "Since we are interested in calculating length, let’s select an equidistance projection" & "North_America_Equidistant_Conic EPSG:102010 as the CRS" it seems there is a subset of EPSG that are appropriate for this task. Of course accuracy will be compromised, but the question is how much? – michael Mar 2 '18 at 17:56
  • 2
    You misunderstand the meaning of "equidistant projection." Except between special points or along special lines, so-called equidistant projections are among the worst possible choices for computing general distances within any region. Perhaps the best general technique for selecting a good projection and assessing its accuracy is to study the Tissot Indicatrices of candidate projections: see gis.stackexchange.com/questions/5068, gis.stackexchange.com/questions/185100, gis.stackexchange.com/questions/50692, and gis.stackexchange.com/questions/32115. – whuber Mar 2 '18 at 18:45
  • @whuber Thanks, by the why a simple approach of choosing the appropriate UTM zone out of the 60 won't work? – michael Mar 2 '18 at 20:30
  • 2
    UTM certainly will work (and I don't see anything in this thread or its comments suggesting otherwise), provided (a) your accuracy requirements are consistent with UTM and (b) your region doesn't straddle two or more UTM zones. – whuber Mar 2 '18 at 20:37

If using ArcGIS Pro: Use getLength on a geometry instead which would let you specify the measurementType of geodesic. Look under AngleAndDistanceTo for description of measurement types. Also available in ArcGIS Desktop.

Or use the Add Geometry Attributes tool with the length_geodesic option. It's also available in ArcGIS Desktop.

Disclosure: I work for Esri.

| improve this answer | |
  • Thanks. I also want to understand the underlying process since probably the methods doing my code but the missing piece in the puzzely is the spatial reference. – michael Mar 2 '18 at 18:38
  • By the way seems like in the end all in all this will be the most practical solution like many times ('unfortunately'): Just use the black box ArcPy API without understand the underlying stuff. – michael Mar 2 '18 at 20:34

Theoretically (rather than practically), to select a projection...find the centroid of your point, get of list of all the projections, subset that list to ones where their 'bounding box' (AKA projected bounds) contains your centroid. That'll cut it down from thousands to hundreds. Then do a compare again, comparing the bounding box of your file to the projected bounds of the projection. Do something like correlation between the numbers to compare; use the one that fits best.

| improve this answer | |

Not the answer you're looking for? Browse other questions tagged or ask your own question.