Différence d'emplacement de destination entre pyproj et geopy

Je cherche des moyens de calculer (en Python) l'emplacement de destination à partir d'un point et d'un angle donnés.

Sur la base de la comparaison des résultats des 2 bibliothèques dans le sujet ( geopyet pyproj), j'ai remarqué qu'il y a une différence croissante dans la sortie finale. Par exemple, à 100 km est à peu près de l'ordre du décimètre. Ceci est un exemple minimal de ce que je veux dire:

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))

J'ai obtenu ces résultats:

        | 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

Ma principale question est de savoir si je fais quelque chose de mal (les deux paramètres doivent être utilisés WGS84).

Sinon, la différence est-elle due aux différentes formules utilisées (Vincenty pour geopyvs Karney pour pyproj)? Par exemple, l'erreur d'arrondi citée ici .

Il semble que vous ayez tout fait correctement. Vous pouvez évaluer les erreurs de chaque méthode en effectuant les calculs inverses pour trouver la distance en fonction des coordonnées d'origine et de destination, puis évaluer les résidus des distances. Il s'agit d'un exercice aller-retour.

# 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))

Montre:

> inverse distance residule (m)
geopy  0.1434377
prj    0.0000000

Vous pouvez donc voir que la méthode de Vincenty détermine une distance inverse qui est sur un décimètre différent pour les mêmes coordonnées. La méthode de Karney présente des erreurs dans la précision de la machine, qui est inférieure à 15 nanomètres. Dans cet exemple, l'erreur est de 0,1445 nm, soit environ le diamètre d'un atome d'hydrogène.

Le problème vient probablement de la méthode de destination du géopy . Comparons une deuxième implémentation de la méthode de Vincenty avec les versions PostGIS 2.1, illustrées ici . (PostGIS version 2.2 avec Proj 4.9 et versions ultérieures utilisent les méthodes de Karney ). La distance aller-retour résiduelle de PostGIS 2.1 est toujours inférieure à 1 cm. Pour cet exemple, c'est 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

La demande d'extraction de geopy résout votre problème avec geopy. Vous devrez d'abord installer le paquetage géographique Python avec

pip install geographiclib