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

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

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 …

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