3

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.

0

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.

  • Indeed the problem seems to be with the shape of the array. I removed the list() call but I still get the error. The issue is that cartopy needs a shapely LineString to compute geometry_length of, but while printing the LineString, it seems to be looking like this: LINESTRING (49.009722 2.547778, 49.015556 2.573611, 49.021389 2.599167, 49.039167 2.676389, 49.048056 2.715,.... . So it seems like a deadlock. – VictorX Aug 1 at 9:05
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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