What is unit of Shapely length attribute?

Question

I'm doing a very simple calculation of the length of a Polyline using shapely:

from shapely.geometry import LineString
... 
xy_list = [map(float, e) for e in xy_intm]
line = LineString(xy_list)
s = '%s,%s,%s' % (fr, to, line.length)

My coordinates are in WGS84. I can't seem to find any information about Shapely's length attribute. What is the unit of the length attribute? Is there an easy way to convert to km or meters?

Answer 1

Coordinate Systems

[...] Shapely does not support coordinate system transformations. All operations on two or more features presume that the features exist in the same Cartesian plane.

Source: http://toblerity.org/shapely/manual.html#coordinate-systems

shapely is completely agnostic in reference to SRS. Therefore, the length attribute is expressed in the same unit of coordinates of your linestring, i.e. degrees. In fact:

>>> from shapely.geometry import LineString
>>> line = LineString([(0, 0), (1, 1)])
>>> line.length
1.4142135623730951

Instead, if you want to express length in meters, you have to transform your geometries from WGS84 to a projected SRS using pyproj (or, better, execute geodesic distance calculation, see Gene's answer). In detail, since version 1.2.18 (shapely.__version__), shapely supports the geometry transform functions ( http://toblerity.org/shapely/shapely.html#module-shapely.ops) that we can use it in conjunction with pyproj. Here's a quick example:

from shapely.geometry import LineString
from shapely.ops import transform
from functools import partial
import pyproj
line1 = LineString([(15.799406, 40.636069), (15.810173,40.640246)])
print(str(line1.length) + " degrees")
0.0115488362184 degrees
Geometry transform function based on pyproj.transform
project = partial(
    pyproj.transform,
    pyproj.Proj('EPSG:4326'),
    pyproj.Proj('EPSG:32633'))
line2 = transform(project, line1)
print(str(line2.length) + " meters")
1021.77585965 meters

Answer 2

As alfaciano says in shapely, the distance is the Euclidean Distance or Linear distance between two points on a plane and not the Great-circle distance between two points on a sphere.

from shapely.geometry import Point
import math


point1 = Point(50.67,4.62)
point2 = Point(51.67, 4.64)

# Euclidean Distance
def Euclidean_distance(point1,point2):
     return math.sqrt((point2.x()-point1.x())**2 + (point2.y()-point1.y())**2)

print Euclidean_distance(point1,point2)
1.00019998 # distance in degrees (coordinates of the points in degrees)

# with Shapely
print point1.distance(point2)
1.0001999800039989 #distance in degrees (coordinates of the points in degrees)

For the great-circle distance, you need to use algorithms as the law of cosines or the Haversine formula (look at Why is law of cosines more preferable than haversine when calculating distance between two latitude-longitude points?) or use the module pyproj that performs geodetic calculations.

# law of cosines
distance = math.acos(math.sin(math.radians(point1.y))*math.sin(math.radians(point2.y))+math.cos(math.radians(point1.y))*math.cos(math.radians(point2.y))*math.cos(math.radians(point2.x)-math.radians(point1.x)))*6371
print "{0:8.4f}".format(distance)
110.8544 # in km
# Haversine formula
dLat = math.radians(point2.y) - math.radians(point1.y)
dLon = math.radians(point2.x) - math.radians(point1.x)
a = math.sin(dLat/2) * math.sin(dLat/2) + math.cos(math.radians(point1.y)) * math.cos(math.radians(point2.y)) * math.sin(dLon/2) * math.sin(dLon/2)
distance = 6371 * 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
print "{0:8.4f}".format(distance)distance
110.8544 #in km

# with pyproj
import pyproj
geod = pyproj.Geod(ellps='WGS84')
angle1,angle2,distance = geod.inv(point1.x, point1.y, point2.x, point2.y)
print "{0:8.4f}".format(distance/1000)
110.9807 #in km

You can test the result at Longitude Latitude Distance Calculator

Answer 3

You can now use shapely w/ pyproj to get geodesic length in meters:

From: https://pyproj4.github.io/pyproj/stable/api/geod.html#pyproj.Geod.geometry_length

from pyproj import Geod
from shapely.geometry import Point, LineString
line_string = LineString([Point(1, 2), Point(3, 4)])
geod = Geod(ellps="WGS84")
f"{geod.geometry_length(line_string):.3f}"
'313588.397'

Answer 4

Use GeographicLib, which is written in many programming languages. The Python implementation is pure Python and does not have any dependencies, so it is easy to install: pip install geographiclib.

Here is an example:

from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84
parameters: lat1, lon1, lat2, lon2
g = geod.Inverse(4.62, 50.67, 4.64, 51.67)
print("The distance is {:.3f} m.".format(g['s12']))
The distance is 110980.677 m.

Or with shapely points:

from shapely.geometry import Point
order is: lon, lat
p1 = Point(50.67, 4.62)
p2 = Point(51.67, 4.64)
g = geod.Inverse(p1.y, p1.x, p2.y, p2.x)
print("The distance is {:.3f} m.".format(g['s12']))
The distance is 110980.677 m.

Answer 5

Update of the most popular answer using the latest version of pyproj (presently 3.3.1), and with some extra tips for different co-ordinate systems.

As mentioned shapely is agnostic to coordinate system, so use pyproj to convert to a co-ordinate system in metres.

from shapely.ops import transform
from shapely.geometry import LineString
from pyproj import Transformer
transformer = Transformer.from_crs("epsg:4326", "epsg:32633")
line1 = LineString([(15.799406, 40.636069), (15.810173,40.640246)])
line2 = transform(transformer.transform, line1)
print(str(line2.length) + " meters")
1400.1203731431115 meters

If like me you are using co-ordinates in (lng, lat) or (x,y) form, then you'll need to add a parameter to the transformer creation method:

transformer = Transformer.from_crs("epsg:4326", "epsg:32633", always_xy=True)

In this case I'm converting to the co-ordinate system used in the answer I linked, but your co-ordinates might be in a different part of the world. You can use https://epsg.io/ to search for your particular country or co-ordinate system. In my case I had co-ordinates in Mexico, so I used "epsg:4484":

transformer = Transformer.from_crs("epsg:4326", "epsg:4484", always_xy=True)