Diferencia en la ubicación de destino entre pyproj y geopy

8

Estoy buscando formas de calcular (en Python) la ubicación de destino desde un punto dado rumbo y rango.

Basado en la comparación de resultados de las 2 bibliotecas en el tema ( geopyy pyproj), noté que hay una diferencia creciente en el resultado final. Por ejemplo, a 100 km es aproximadamente del orden de decímetros. Este es un ejemplo mínimo de lo que quiero decir:

from __future__ import absolute_import, division, print_function

long_1 = -1.729722
lat_1 = 53.320556
bearing = 96.021667
distance = 124.8  #in km

# using geopy

import geopy
from geopy.distance import VincentyDistance

origin = geopy.Point(lat_1, long_1)
destination = VincentyDistance(kilometers=distance).destination(origin, bearing)
gp_lat_2 = destination.latitude
gp_long_2 = destination.longitude

# using pyproj

from pyproj import Geod
g = Geod(ellps='WGS84')
prj_long_2, prj_lat_2, prj_bearing_2 = g.fwd(long_1, lat_1, bearing, distance*1000)

# results comparison
print("        | pyproj        | geopy")
print("long_2   %.6f     %.6f" % (prj_long_2, gp_long_2))
print("lat_2    %.6f     %.6f" % (prj_lat_2, gp_lat_2))

print("> DELTA pyproj, geopy")
print("delta long_2   %.7f" % (prj_long_2 - gp_long_2))
print("delta lat_2    %.7f" % (prj_lat_2 - gp_lat_2))

Obtuve estos resultados:

        | pyproj        | geopy
long_2   0.127201         0.127199
lat_2    53.188432        53.188432

> DELTA pyproj, geopy
delta long_2   0.0000021
delta lat_2    -0.0000002

Mi pregunta principal es si estoy haciendo algo mal (ambas configuraciones deberían estar usando WGS84).

Si no, ¿la diferencia se debe a las diferentes fórmulas en uso (Vincenty para geopyvs. Karney para pyproj)? Por ejemplo, el error de redondeo citado aquí .

gmas80
fuente

Respuestas:

6

Parece que has hecho todo correctamente. Puede evaluar los errores de cada método realizando los cálculos inversos para encontrar la distancia dadas las coordenadas de origen y destino, luego evalúe los residuos de las distancias. Este es un ejercicio de ida y vuelta.

# For Vincenty's method:
geopy_inv_dist = geopy.distance.vincenty(origin, destination).m
# For Karney's method:
prj_inv_dist = g.inv(long_1, lat_1, prj_long_2, prj_lat_2)[2]  # s12

print("> inverse distance residule (m)")
print("geopy  %.7f" % (distance * 1000 - geopy_inv_dist))
print("prj    %.7f" % (distance * 1000 - prj_inv_dist))

Muestra:

> inverse distance residule (m)
geopy  0.1434377
prj    0.0000000

Entonces puede ver que el método de Vincenty determina una distancia inversa que está sobre un decímetro diferente para las mismas coordenadas. El método de Karney tiene errores en la precisión de la máquina, que es inferior a 15 nanómetros. En este ejemplo, el error es 0.1455 nm, que es alrededor del diámetro de un átomo de hidrógeno.


El problema es probablemente con el método de destino de geopy . Comparemos una segunda implementación del método de Vincenty con las versiones 2.1 de PostGIS, que se muestran aquí . (PostGIS versión 2.2 con Proj 4.9 y posterior usa los métodos de Karney ). Los residuos de distancia de ida y vuelta de PostGIS 2.1 siempre son inferiores a 1 cm. Para este ejemplo es 255 nm:

SELECT PostGIS_Version(),
  ST_AsText(origin) AS origin,
  ST_AsText(ST_Project(origin, distance, azimuth)) AS destination,
  ST_Distance(ST_Project(origin, distance, azimuth), origin) AS roundtrip_distance,
  distance - ST_Distance(ST_Project(origin, distance, azimuth), origin) AS postgis_residual
FROM (
  SELECT 124.8 * 1000 AS distance, radians(96.021667) AS azimuth,
    ST_SetSRID(ST_MakePoint(-1.729722, 53.320556), 4326)::geography AS origin
) AS f;
-[ RECORD 1 ]------+-----------------------------------------
postgis_version    | 2.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
origin             | POINT(-1.729722 53.320556)
destination        | POINT(0.12720134063267 53.1884316458524)
roundtrip_distance | 124799.999999745
postgis_residual   | 2.54993210546672e-007
Mike T
fuente
¡Gracias por esta respuesta extremadamente completa, extremadamente inteligente!
Jason Scheirer
@ Mike-t, gracias por la buena respuesta! Me gusta la idea de una tercera fuente. He abierto una geopyentrada: github.com/geopy/geopy/issues/140
gmas80