How to calculate coordinates from bearing/distance?

Q&A
1

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.

brgdist
brgdist600×600 27.6 KB

> A
# lon lat
    8  50

> B
# lon lat
    ?   ?

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

2

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

@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 …

4

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.

5

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

6

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

7

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

8

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

9

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
10

Thank you very much @pwessel!