¿Trabajar con coordenadas geográficas en shapely?

Question

Tengo una pregunta muy básica sobre el trabajo con el shapely en Python. Tengo unas coordenadas geográficas recibidas de un GPS. Me gustaría utilizar shapely para calcular la gran distancia cíclica en metros entre dos puntos. Empecé con:

>>> from shapely.geometry import Point
>>> p1 = Point(43.374880, -78.119956)
>>> p2 = Point(43.374868, -78.119666)

Creo que esto me da dos puntos en el sistema de coordenadas cartesianas, lo que no va a ser muy útil. Estas coordenadas provienen de un GPS que presumiblemente está utilizando el WGS84 CRS. He visto algún ejemplo usando pyproj y shapely.ops.transform para transformar entre sistemas de coordenadas, pero todos ellos parecen implicar puntos de datos que ya están en algún sistema de coordenadas cartográfico válido con un identificador EPSG. No estoy seguro de qué hacer con puntos ingenuos como los anteriores que no tienen un CRS geográfico asociado.

Cómo asociar estos puntos con el SIR WGS84 de manera que p1.distance(p2) ¿dará la distancia en metros? ¿Estoy utilizando las herramientas equivocadas?

Answer 1

El problema que se intenta resolver es el Problema geodésico inverso . Afortunadamente la Python geographiclib puede hacer esto por ti:

from geographiclib.geodesic import Geodesic

p1_lat, p1_lon = 43.374880, -78.119956
p2_lat, p2_lon = 43.374868, -78.119666
geod = Geodesic.WGS84

# note the Inverse method expects latitude then longitude
g = geod.Inverse(p1_lat, p1_lon, p2_lat, p2_lon)

print("Distance is {:.2f}m".format(g['s12']))

Dando el resultado:

Distance is 23.54m

Dicho esto, dado que es probable que los puntos estén a sólo 10 metros de distancia, trabajar en unidades proyectadas será mucho más fácil.

Answer 2

En realidad, creo que lo he descubierto. Los dos puntos anteriores se encuentran a orillas del lago Ontario, por lo que puedo convertir el (lat, lon) en una referencia de cuadrícula utilizando EPSG:32117 ( NAD83/Nueva York Oeste ). Terminé con:

>>> from shapely.geometry import Point
>>> from pyproj import Proj
>>> nys = Proj(init='EPSG:32117')
>>> p1 = Point(43.374880, -78.119956)
>>> p2 = Point(43.374868, -78.119666)
>>> p1_proj = nys(p1.y, p1.x)
>>> p2_proj = nys(p2.y, p2.x)
>>> d = Point(p1_proj).distance(Point(p2_proj))
>>> d
23.539335485921658

Lo que parece razonable para mis necesidades (los puntos que estoy utilizando nunca van a estar a más de decenas de metros de distancia).

Se trata de una distancia lineal en lugar de una distancia ortodrómica, y requiere conocer la CRS adecuada para las coordenadas dadas. Si hay una solución más general, sería genial.