Shapely linestring .length units

user_7732337 0
Viewed 4584

I'm trying to find the length (in metres) of linestrings in Shapely, but can't seem to achieve the expected result. It's almost guaranteed I'm making some mistake in terms of coordinate systems, but I haven't been able to figure it out.

Here's some simplified code for a single line:

from shapely.geometry import LineString
line = LineString([(12875996.563923, -3940011.116702), (12872802.929335, -3937989.118438)])
line.crs = 'EPSG:3857'
line.length

The output is 3779.92m. However, according to my measurements, it should be ~3159m.

That is based on checking the location of the points (the center of two roundabouts) at the following site, and then measuring between them in Google Earth.

Point 1: http://epsg.io/map#srs=3857&x=12875996.563923&y=-3940011.116702&z=14&layer=streets
Point 2: http://epsg.io/map#srs=3857&x=12872802.929468&y=-3937989.118538&z=17&layer=streets

shapely coordinate-systems spatial python
1 Answers

To show that, line.length from shapely geometry is just a Euclidean distance, one can compute like this:

import math
x1, y1 = (12875996.563923, -3940011.116702)
x2, y2 = (12872802.929335, -3937989.118438)
math.hypot(x2-x1,  y2-y1)  # 3779.91783790157

To get correct distance with pyproj module, follow these steps:

import pyproj  #v2.4.1
x1, y1 = (12875996.563923, -3940011.116702)
x2, y2 = (12872802.929335, -3937989.118438)
lon1, lat1 = pyproj.Proj("epsg:3857")(x1, y1, inverse=True)
lon2, lat2 = pyproj.Proj("epsg:3857")(x2, y2, inverse=True)
_,_,dist_km = pyproj.Geod(ellps='WGS84').inv(lon1, lat1, lon2, lat2)
dist_km   # 3157.214113925091
user_21774130
Related