How to calculate coordinates from bearing/distance?

I’m sure I miss the forest for the trees but I can’t find it in the docs:

I have a starting coordinate A N50° E8°, a bearing 310° and a distance 23 NM. I want to calculate the coordinate B where bearing and distance are pointing to.

``````> A
# lon lat
8  50

> B
# lon lat
?   ?
``````

How do I accomplish this in GMT? Your nudge in the right direction is greatly appreciated!

Refering to https://stackoverflow.com/a/33002517 there is a simple python function can do it:

``````import math

def getEndpoint(lat1,lon1,bearing,d):
R = 6371                     #Radius of the Earth
brng = math.radians(bearing) #convert degrees to radians
d = d*1.852                  #convert nautical miles to km
lat1 = math.radians(lat1)    #Current lat point converted to radians
lon1 = math.radians(lon1)    #Current long point converted to radians
lat2 = math.asin( math.sin(lat1)*math.cos(d/R) + math.cos(lat1)*math.sin(d/R)*math.cos(brng))
lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1),math.cos(d/R)-math.sin(lat1)*math.sin(lat2))
lat2 = math.degrees(lat2)
lon2 = math.degrees(lon2)
return (lat2, lon2)

lat2, lon2 = getEndpoint(lat1=50,lon1=8,bearing=310,d=23)``````

@saeedsltm thank you for bringing up that function – I’m currently writing a shell script so I try to adapt it for bash. By the looks of it there might be some issues with the equator and I haven’t tested it for the poles. That’s why I was hoping for GMT offering a solution out of the box …

The closest that I know is:

\$ Dkm=\$(echo 23*1.852 | bc -l);gmt project -C8/50 -A310 -L0/\$Dkm -G\$Dkm -Q
8 50 0
7.54112215577 50.24533433 42.596

and you could use some awk to get the 7.54112215577 50.24533433

other possibility is:

\$ echo 0 23 | gmt mapproject -R7/9/49/51 -Joa8/50/\$((310+90))/1:1 -Fn -C -I
7.54112559204 50.2453325089

I’m not sure what is more accurate / convenient.

I normally would use (lat lon order):

\$ echo 50 8 310 23 | geod +ellps=WGS84 +units=kmi -f %.6f
50.245257 7.542544 129.648938

but is not a gmt solution.

I’ve been here before but don’t remember how it ended. `project` is spherical only. Probably the best is to stick to `proj4` program `geod`. It gives the same as my Matlab version that uses the Vincenty algorithm

7.542543557083350 50.245256893178066

echo 8 50 | gmt vector -Tt310/23n -fg

But that’s spherical only too. At least it doesn’t accept -je

Yes, that is true. I could look at implmenting -je for this but not today.

Well, found the time so the geidesic version is now in master, @KristofKoch:

``````echo 8 50 | gmt vector -Tt310/23n -je
7.54254355708 50.2452568932``````
1 Like

Thank you very much @pwessel!