12

I was quite unsatisfied with Calculating Length of Linestrings in WGS84 in Miles. It kept me wondering if there is a more convenient, Pythonic way to calculate the length of a WKT linestring according to a given SRID.

I have in mind something like:

srid="WGS84"
line="LINESTRING(3.0 4.0, 3.1 4.1)"
print length(line, srid)

I'm looking for an accurate answer, not sin\cos approximations.

Any ideas?

Community
  • 1
Adam Matan
  • 6,838
  • 7
  • 37
  • 50
  • 1
    tomkralidis, this is a GIS website. your answer ignores that the this is a distance between geospatial coordinates (look up SRID). shapely in of itself cannot compute geospatial distances as it has no knowledge of map projection. –  Oct 18 '13 at 01:32

5 Answers5

18

The geopy module provides the Vincenty formula, which provides accurate ellipsoid distances. Couple this with the wkt loading in Shapely, and you have reasonably simple code:

from geopy import distance
from shapely.wkt import loads

line_wkt="LINESTRING(3.0 4.0, 3.1 4.1)"

# a number of other elipsoids are supported
distance.VincentyDistance.ELLIPSOID = 'WGS-84'
d = distance.distance

line = loads(line_wkt)

# convert the coordinates to xy array elements, compute the distance
dist = d(line.xy[0], line.xy[1])

print dist.meters
scw
  • 16,391
  • 6
  • 64
  • 101
  • 1
    +1, would have +10 if I could. Saved my team hours of programming. – Adam Matan Dec 01 '10 at 22:03
  • Is this approach any different from @tomkralidis answer if the input coordinates are in WGS-84 already? – LarsVegas Dec 18 '13 at 07:54
  • 1
    @LarsVegas yes, Shapely only handles planar coordinates -- so it will measure distances accurately in projected space, but not geographic (e.g. WGS-1984). – scw Dec 23 '13 at 04:05
3

You could also use Shapely's length property, i.e.:

from shapely.wkt import loads

l=loads('LINESTRING(3.0 4.0, 3.1 4.1)')
print l.length
tomkralidis
  • 1,297
  • 6
  • 8
  • 1
    Note that the length for this particular example will be meaningless, as it is a geographic coordinate system (WGS84). – Mike T Aug 05 '14 at 00:38
2

Late to the party, but with a hopefully useful contribution. Building on scw's answer using geopy, I wrote a small function that does the calculation for a shapely LineString object with arbitrarily many coordinates. It uses a pairs iterator from Stackoverflow.

Main feature: the docstrings are much longer than the snippets.

def line_length(line):
    """Calculate length of a line in meters, given in geographic coordinates.
    Args:
        line: a shapely LineString object with WGS 84 coordinates
    Returns:
        Length of line in meters
    """
    # Swap shapely (lonlat) to geopy (latlon) points
    latlon = lambda lonlat: (lonlat[1], lonlat[0])
    total_length = sum(distance(latlon(a), latlon(b)).meters
                       for (a, b) in pairs(line.coords))
    return round(total_length, 0)


def pairs(lst):
    """Iterate over a list in overlapping pairs without wrap-around.

    Args:
        lst: an iterable/list

    Returns:
        Yields a pair of consecutive elements (lst[k], lst[k+1]) of lst. Last 
        call yields the last two elements.

    Example:
        lst = [4, 7, 11, 2]
        pairs(lst) yields (4, 7), (7, 11), (11, 2)

    Source:
        https://stackoverflow.com/questions/1257413/1257446#1257446
    """
    i = iter(lst)
    prev = i.next()
    for item in i:
        yield prev, item
        prev = item
ojdo
  • 343
  • 1
  • 2
  • 10
  • 1
    This is wrong: geopy.distance.distance accepts coordinates in (y, x) but a shapely linestring is "an ordered sequence of 2 or more (x, y[, z]) " so geopy's helper function lonlat() must be used. – Martin Burch Apr 10 '19 at 12:47
  • @MartinBurch: ouch, you are right. The nasty thing is not even the [, z], but the argument swap (y, x) to (x, y) that is necessary. Thanks for spotting it. Can you eyeball whether this edit looks less bug-ridden? – ojdo Apr 10 '19 at 15:05
2

I'd use ogr2ogr to do it directly but if you really must use Python then there are Python bindings to let you do it.

Glorfindel
  • 1,096
  • 2
  • 9
  • 14
Ian Turton
  • 81,417
  • 6
  • 84
  • 185
0

The solution suggested by @ojdo is very nice and it works.

Another way is to do a lat, lon conversion to UTM e.g. using https://gist.github.com/twpayne/4409500. And then the distance is immediate.

TodoranIonut
  • 1