Lat/lon after moving x meters to a specified direction

gpshaversinepyprojpythonwgs84

I have recently been working with coordinate systems in Python with pyproj package. I'm facing the following problem: I know my current position from a GPS receiver as EPSG:4326 latitude/longitude values. If I move x meters to a specified direction a, say to a direction of 15-degrees-anticlockwise-from-east, what is my new position in the same EPSG:4326 latitude/longitude values? How is this problem solved mathematically and are there any functions in pyproj that can assist me?

Edit: Here is a naive approach, where I transform to EPSG:3857 (Web Mercator projection) projected coordinate system, move 10000 units there, and transform back to EPSG:4326, which does not produce correct results:

import numpy as np
from pyproj import Transformer
from pyproj import CRS

# Coordinate systems
CRS_GEO = CRS.from_epsg(4326)  # Geodetic coordinate reference system WGS 84
CRS_WM = CRS.from_epsg(3857) # Web Mercator projected coordinate system

# Transformers
T_GEO_WM = Transformer.from_crs(CRS_GEO, CRS_WM, always_xy=True)
T_WM_GEO = Transformer.from_crs(CRS_WM, CRS_GEO, always_xy=True)

# Original places
big_ben_geo = (-0.1246547718409441, 51.50075527777225) # Geodetic coordinates of Big Ben (lon, lat)
big_ben_wm = T_GEO_WM.transform(*big_ben_geo)

# Move 10 km to a direction of 15 deg anticlockwise from east
r = 10000
a = np.radians(15)
movement = r * np.array([np.cos(a), np.sin(a)])
big_ben_wm_moved = (big_ben_wm[0] + movement[0], big_ben_wm[1] + movement[1])

# Transform back to geodetic coordinates
big_ben_wm_moved_geo = T_WM_GEO.transform(*big_ben_wm_moved)
print(big_ben_wm_moved_geo)
# Prints (-0.03788417853281929, 51.51522627382787)

Now if I check the distance between those two locations from e.g. Google Maps, the actual distance is ~6.2 km although I wanted to move 10 km.

Edit2: I think I found the correct search term: Inverse/reverse haversine.
I cannot test these right now, but I think using Python packages in these links should work for me:

Edit3: Using pyproj.Geod works perfectly for my needs, as suggested by raphael.

Best Answer

You don't need to re-invent the wheel... what you want is the Pyproj.geod module!

https://pyproj4.github.io/pyproj/stable/api/geod.html

From the docs:

performs forward and inverse geodetic, or Great Circle, computations. The forward computation (using the ‘fwd’ method) involves determining latitude, longitude and back azimuth of a terminus point given the latitude and longitude of an initial point, plus azimuth and distance. The inverse computation (using the ‘inv’ method) involves determining the forward and back azimuths and distance given the latitudes and longitudes of an initial and terminus point.

from pyproj import Geod
geod = Geod(ellps="WGS84")
geod.fwd(lons=20, lats=30, az=15, dist=1000)
Out[21]: (20.00268267840205, 30.008713584787657, -164.99865848413666)
Related Question