Procurando uma maneira pitônica de calcular o comprimento de uma cadeia de linhas WKT

13

Fiquei bastante insatisfeito com o cálculo do comprimento das cadeias de linhas no WGS84 em milhas . Isso me fez pensar se existe uma maneira Python mais conveniente de calcular o comprimento de uma cadeia de linhas WKT de acordo com um determinado SRID.

Eu tenho em mente algo como:

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

Estou procurando uma resposta precisa, não sin\cosaproximações.

Alguma ideia?

Adam Matan
fonte
tomkralidis, este é um site de GIS. sua resposta ignora que esta é uma distância entre coordenadas geoespaciais (consulte SRID). o bem torneado por si só não pode calcular distâncias geoespaciais, pois não tem conhecimento de projeção de mapa.

Respostas:

18

O módulo geopy fornece a fórmula Vincenty , que fornece distâncias elipsóides precisas. Junte isso ao wktcarregamento no Shapely e você terá um código razoavelmente simples:

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
fonte
1
+1, teria +10 se eu pudesse. Economizei horas de programação para minha equipe.
Adam Matan
Essa abordagem é diferente da resposta @tomkralidis se as coordenadas de entrada já estiverem no WGS-84?
LarsVegas
1
@LarsVegas sim, Shapely lida apenas com coordenadas planares - portanto, medirá as distâncias com precisão no espaço projetado, mas não geográfico (por exemplo, WGS-1984).
ACS
4

Você também pode usar a propriedade length da Shapely , ou seja:

from shapely.wkt import loads

l=loads('LINESTRING(3.0 4.0, 3.1 4.1)')
print l.length
tomkralidis
fonte
Observe que o comprimento deste exemplo em particular não terá sentido, pois é um sistema de coordenadas geográficas (WGS84).
Mike T
2

Tarde para a festa, mas com uma contribuição esperançosamente útil. Com base na resposta do scw usando geopy, escrevi uma pequena função que faz o cálculo para um objeto LineString bem torneado com arbitrariamente muitas coordenadas. Ele usa um pairsiterador do Stackoverflow.

Recurso principal: os documentos são muito mais longos que os trechos.

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
fonte
1
Isso está errado: geopy.distance.distance aceita coordenadas em (y, x), mas uma cadeia de linhas bem torneada é "uma sequência ordenada de 2 ou mais (x, y [, z])", portanto a função auxiliar do geopy lonlat () deve ser usada .
Martin Burch
@ MartinBurch: ai, você está certo. A coisa desagradável não é mesmo o [, z], mas a troca de argumento (y, x)para (x, y)o que é necessário. Obrigado por detectá-lo. Você pode observar se essa edição parece menos cheia de bugs?
Ojdo