13

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?

  • 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. – user23037 Oct 18 '13 at 1:32
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
  • 1
    +1, would have +10 if I could. Saved my team hours of programming. – Adam Matan Dec 1 '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 7: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 4:05
4

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
  • Note that the length for this particular example will be meaningless, as it is a geographic coordinate system (WGS84). – Mike T Aug 5 '14 at 0:38
2

I'd use ogr2ogr (http://www.gdal.org/ogr/index.html) to do it directly but if you really must use python then there are python bindings (http://pypi.python.org/pypi/GDAL/) to let you do it.

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
  • 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 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 at 15:05

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.