Calculating geodesic distance along a path (lat/lon points) at once?

Question

Using Python 3.7 with Jupyter Notebook on a Win10 x64.

I have a list of lat-long tuples representing a path and I want to calculate the total length along that path in meters (or km, miles etc). I would like to avoid calculating each segment distance and then adding them all together as I have to do this for 7 million paths. So time efficiency is key. Adding all segments after computing individual distances takes 7ms per path. I want to make it at least 1ms fast.

I need to calculate distances using the WGS84 ellipsoid, so spherical (Haversine) is not sufficient. I think I could work with 1m accurracy. Points are distributed randomly along paths. Some can have 20 kilometers between them, and some have below 1 km.

Here is the path with points (decimal lat/lon):

[(49.009722, 2.547778), (49.015556, 2.573611), (49.021389, 2.599167), (49.039167, 2.676389), (49.048056, 2.715), (49.044444, 2.835), (49.041667, 2.928333), (49.042778, 2.942222), (49.051667, 3.066667), (49.061389, 3.205), (49.072222, 3.357222), (49.085, 3.536944), (49.086111, 3.550833), (49.097778, 3.729444), (49.113056, 3.963056), (49.130833, 4.238056), (49.138889, 4.361667), (49.1925, 4.564444), (49.306667, 4.995556), (49.333611, 5.096944), (49.395, 5.329167), (49.490556, 5.690833), (49.514444, 5.781111), (49.53, 5.845833), (49.599444, 6.127778), (49.637222, 6.281667), (49.673333, 6.440278), (50.0475, 8.078333), (50.053611, 8.637222), (50.056667, 8.800278), (50.063056, 9.19), (50.066944, 9.486389), (50.07, 9.783056), (50.072778, 10.098611), (50.073333, 10.242778), (50.075278, 10.728889), (50.046667, 10.863333), (50.0325, 10.930278), (49.981111, 11.172222), (49.969722, 11.225833), (49.961111, 11.491389), (49.959444, 11.547222), (49.957222, 11.617222), (49.946111, 11.9325), (49.937222, 12.343889), (49.933333, 12.47), (49.9325, 12.498056), (49.928611, 12.624167), (49.924167, 12.764444), (49.919444, 12.918611), (49.910833, 13.199167), (49.909444, 13.241111), (49.907778, 13.283611), (49.900556, 13.481944), (50.077222, 13.840278), (50.124167, 13.995556), (50.182778, 14.189722), (50.220278, 14.315), (50.268889, 14.478056), (50.211389, 14.403611), (50.166389, 14.345), (50.133611, 14.3025), (50.100833, 14.26)]

I have discovered cartopy - a package providing geodesic calculations that uses shapely and proj. However the docs on cartpopy do not give a relevent example and I am stuck at this point. Basically the geometry_length gives the length of a shapely object in one go, so I am doing it in the following way:

#defining the geoid on which to make calculations
myGeod = geodesic.Geodesic(6378137.0,1 / 298.257223563)

#making my list of latlon (in decimal degrees) into a shapely 
shapelyObject = LineString(list(latlon_dd))

#applying the method on the shapelyObject given the defined ellipsoid
myGeod.geometry_length(shapelyObject)

I want to compute the length in meters, which should be around 917,315.3 meters. Instead, I get this ValueError:

ValueError                                Traceback (most recent call last)
<ipython-input-243-7c75042775e3> in <module>
      6 
      7 #applying the method on the shapelyObject given the defined ellipsoid
----> 8 myGeod.geometry_length(shapelyObject)

lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.geometry_length()

lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.geometry_length()

lib\cartopy\geodesic\_geodesic.pyx in cartopy.geodesic._geodesic.Geodesic.inverse()

ValueError: Expecting input points to be (N, 2), got (1, 63)

I can also use a different package, as long as it's fast. GeoPy and GeographicLib did not give better results.

Answer 1

ValueError: Expecting input points to be (N, 2), got (1, 63)

This is a problem with the 'shape' of your array. It suggests the code is expecting something like this.. a list of coordinates, where each coordinate is an array (or tuple) of length 2.

[[x1,y1],[x2,y2], ... ,[xN,yN]]

but somewhere, your code has 'flattened' the array, so it looks like this

[x1,y1,x2,y2 ... ,xN,yN]

where N is 63.

I suspect it's the list() call in this line

shapelyObject = LineString(list(latlon_dd))

Try removing the list() and see if that works?

My numpy/python is rusty now, but hopefully that will help you, or prompt someone else to suggest a working example. Have a look at numpy reshape too as this may also be useful in situations like this.

Answer 2

Going to answer my question based on the reply found here. Apparently it is a bug, and the current workaround is to use:

myGeod.geometry_length(np.array(shapelyObject.coords))

instead of

myGeod.geometry_length(shapelyObject)

Will update when a final solution is available.