[GIS] Difference in destination location between pyproj and geopy

geodesygeopypyprojpythonvincenty-formulae

I am looking at ways to calculate (in Python) the destination location from a point given bearing and range.

Based on the results comparison from the 2 libraries in subject(geopy and pyproj), I noticed that there is an increasing difference in the final output. For instance, at 100 km is roughly of the order of decimeters. This is a minimal example of what I mean:

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

I got these results:

        | 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

My main question is whether I am doing something wrong (both settings should be using WGS84).

If not, is the difference due to different formulas in use (Vincenty for geopy vs. Karney for pyproj)? E.g., the round-off error cited here.

Best Answer

It looks like you've done everything correctly. You can evaluate the errors from each method by performing the inverse calculations to find the distance given the origin and destination coordinates, then evaluate the residuals of distances. This is a round-trip exercise.

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

Shows:

> inverse distance residule (m)
geopy  0.1434377
prj    0.0000000

So you can see that Vincenty's method determines an inverse distance that is over a decimetre different for the same coordinates. Karney's method has errors within machine precision, which is less than 15 nanometers. In this example the error is 0.1455 nm, which is around the diameter of a hydrogen atom.

The problem is probably with geopy's destination method. Let's compare a second implementation of Vincenty's method with PostGIS versions 2.1, shown here. (PostGIS version 2.2 with Proj 4.9 and later use Karney's methods). The round-trip distance residuals from PostGIS 2.1 is always less than 1 cm. For this example it is 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
Related Question