1

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.

4
  • 1
    latitude/longitude values are degrees, you can not add "meter" to them. Commented Nov 17, 2022 at 11:53
  • What have you tried because unless you add a Python code attempt your question will be closed
    – gene
    Commented Nov 17, 2022 at 12:19
  • 2
    No distances in Web Mercator are correct. Being 40% off is expected at high latitudes (unless it's very high latitudes, when the error is infinite).
    – Vince
    Commented Nov 17, 2022 at 13:32
  • 1
    You're on the right track, but you picked a lousy projection for this particular purpose. UTM would be better, but will still only be reasonably accurate for distances that stay within the selected zone. Furthermore, you will have to pick the right UTM zone for your data - there is no universal zone that works for the entire globe.
    – Llaves
    Commented Nov 17, 2022 at 16:11

1 Answer 1

1

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.