5

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

Georgy
  • 12,464
  • 7
  • 65
  • 73
EllRob
  • 145
  • 1
  • 3
  • 11
  • 1
    EPSG:3857 is a Spherical Mercator projection coordinate system. Distortion in length is unacceptable in most situation. – swatchai May 22 '20 at 13:28
  • Shapely [doesn't support](https://shapely.readthedocs.io/en/latest/manual.html#coordinate-systems) coordinate system transformations. Assigning the `crs` attribute doesn't have any effect. Comment that line and you will see that the result is the same. I think you need to use pyproj library instead. Check [this example in Shapely docs](https://shapely.readthedocs.io/en/latest/manual.html#shapely.ops.transform) and [this entry in pyproj docs](https://pyproj4.github.io/pyproj/stable/api/geod.html#pyproj.Geod.geometry_length). – Georgy May 22 '20 at 13:41

1 Answers1

4

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
swatchai
  • 17,400
  • 3
  • 39
  • 58