2

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

5
  • 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
    Commented Mar 2, 2018 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
    Commented Mar 2, 2018 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
    Commented Mar 2, 2018 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
    Commented Mar 2, 2018 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
    Commented Mar 2, 2018 at 20:37

2 Answers 2

2

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.

2
  • 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
    Commented Mar 2, 2018 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
    Commented Mar 2, 2018 at 20:34
0

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.

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