Buscando una forma pitónica para calcular la longitud de una cadena lineal WKT

13

Estaba bastante insatisfecho con el cálculo de la longitud de las cadenas lineales en WGS84 en millas . Me mantuvo preguntándome si hay una manera más conveniente y pitónica de calcular la longitud de una cadena lineal WKT de acuerdo con un SRID dado.

Tengo en mente algo como:

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

Estoy buscando una respuesta precisa, no sin\cosaproximaciones.

¿Algunas ideas?

Adam Matan
fuente
tomkralidis, este es un sitio web de SIG. su respuesta ignora que esta es una distancia entre coordenadas geoespaciales (busque SRID). bien proporcionado por sí mismo no puede calcular distancias geoespaciales ya que no tiene conocimiento de la proyección del mapa.

Respuestas:

18

El módulo geopy proporciona la fórmula de Vincenty , que proporciona distancias elipsoidales precisas. wktCombine esto con la carga en Shapely, y tendrá un código razonablemente simple:

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
fuente
1
+1, tendría +10 si pudiera. Le ahorré a mi equipo horas de programación.
Adam Matan
¿Es este enfoque diferente de la respuesta @tomkralidis si las coordenadas de entrada ya están en WGS-84?
LarsVegas
1
@LarsVegas sí, Shapely solo maneja coordenadas planas, por lo que medirá distancias con precisión en el espacio proyectado, pero no geográfico (por ejemplo, WGS-1984).
scw
4

También puede usar la propiedad de longitud de Shapely , es decir:

from shapely.wkt import loads

l=loads('LINESTRING(3.0 4.0, 3.1 4.1)')
print l.length
tomkralidis
fuente
Tenga en cuenta que la longitud de este ejemplo en particular no tendrá sentido, ya que es un sistema de coordenadas geográficas (WGS84).
Mike T
2

Tarde a la fiesta, pero con una contribución útil. Sobre la base de la respuesta de scw usando geopy, escribí una pequeña función que hace el cálculo para un objeto LineString bien formado con arbitrariamente muchas coordenadas. Utiliza un pairsiterador de Stackoverflow.

Característica principal: las cadenas de documentos son mucho más largas que los fragmentos.

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:
        /programming/1257413/1257446#1257446
    """
    i = iter(lst)
    prev = i.next()
    for item in i:
        yield prev, item
        prev = item
ojdo
fuente
1
Esto está mal: geopy.distance.distance acepta coordenadas en (y, x) pero una cadena lineal bien formada es "una secuencia ordenada de 2 o más (x, y [, z])", por lo que debe usarse la función auxiliar de geopy lonlat () .
Martin Burch
@MartinBurch: ay, tienes razón. Lo desagradable ni siquiera es el [, z], pero el argumento de intercambio (y, x)a (x, y)que es necesario. Gracias por verlo. ¿Puedes ver si esta edición parece menos cargada de errores?
ojdo