Actualmente estoy usando Python / Numpy para tratar con datos geográficos / GPS (¡me encanta!), Y me enfrento a la tarea recurrente para calcular distancias entre puntos geográficos definidos por un par de coordenadaspn = [lon, lat].

Tengo una función que uso así:dist = geodistance(p1, p2) que es análoga a la distancia euclidiana en el álgebra lineal (resta / diferencia del vector), pero ocurre en el espacio geodésico (esférico) en lugar de en el espacio euclidiano rectangular.

Programáticamente, la distancia euclidiana está dada por

dist = ((p2[0] - p1[0])**2 + (p2[1] - p1[1])**2)**0.5

Matemáticamente, esto es equivalente a la oración "idiomática" (por falta de una palabra mejor)

dist = p1 - p1   # the "norm" of the vector difference, subtraction.

Actualmente, obtengo mi distancia como esta:

p1 = [-51.598354,-29.953363]
p2 = [-51.598701,-29.953045]
dist = geodistance(p1, p2)
print dist

>> 44.3904032407

Me gustaría hacer esto:

print p2 - p1  # these points now are from some fancy datatype

>> 44.3904032407

Y el objetivo final:

track = numpy.array([[-51.203018 -29.996149]
                     [-51.203018 -29.99625 ]
                     [-51.20266  -29.996229]
                     [-51.20229  -29.996309]
                     [-51.201519 -29.99416 ]], dtype=fancy)  # (**) or something like

print numpy.diff(track)

>> ndarray([[   0.        ]
            [   7.03531252]
            [  39.82663316]
            [  41.50958596]
            [ 172.49825765]])

Algo similar es: si tomas dosdatetime objetos y restarlos, la operación devuelve untimedelta objeto. Quiero restar dos coordenadas y obtener una distancia geodésica como resultado.

Me pregunto si una clase funcionaría, pero dtype (un "subtipo" de float32, por ejemplo) ayudaría mucho en la creación de matrices a partir de listas (** que es cómo leo las cosas de los archivos XML).

¡Muchas gracias!