# calculate bearing from one MapPoint to another MapPoint

I am using Esri.ArcGISRuntime package inside a C# program we are writing.

I need to create specific shapes from a reference point and project them in different directions. This is easy to create a new point 5000m away on a 45 degree bearing. ie

``````var startPoint = new MapPoint(125.4, -33.2, SpatialReferences.Wgs84);
var endPoint = GeometryEngine.GeodesicMove(startPoint, 5000, LinearUnits.Meters, 45);
``````

Later on when I need to reverse engineer the points, I can easily recalculate the distance between these points using

``````var distance = GeometryEngine.GeodesicDistance(startPoint, endPoint, LinearUnits.Meters);
``````

However, I cannot find anyway to determine what the original angle or bearing was between these 2 points.

I have tried using Arctan2, but the MapPoint coordinates are skewed by earth curvature (I think) and therefore straight forward calculation is incorrect.

``````var degrees = Math.Atan2(endPoint.Y - startPoint.Y, endPoint.X - startPoint.X) / Math.Pi * 180;
``````

The above returns a figure of roughly 39 degrees (PS I tried changing X & Y around as well, no better).

Take 39 as first estimate, 29 as low boundary, 49 as high boundary. Try to minimise distance between known end point and new end point at given bearing using 1st equation, see how golden section works https://en.wikipedia.org/wiki/Golden_section_search

Update on original answer I don’t have C# etc, so this is solution using arcpy

``````import  arcpy, traceback, os, sys, math
from arcpy import env
env.overwriteOutput = True
tbl=r"d:/scratch/excel.dbf"
line=r"D:/Scratch/theLine.shp"
g=arcpy.Geometry()
gr=(math.sqrt(5)-1)/2
def gss(a,b,tol):
c=b-gr*(b-a)
d=a+gr*(b-a)
while abs(c-d)>tol:
fc=f(c);fd=f(d)
if fc<fd:
b=d
d=c
c=b-gr*(b-a)
else:
a=c
c=d
d=a+gr*(b-a)
return (b+a)/2

def f(x):
startPointX=125.4
startPointY=-33.2
endPointX=125.437904
endPointY=-33.168116
L=5000
D=x
arcpy.DeleteRows_management(tbl)
dFields=("X","Y","L","D")
curT = arcpy.da.InsertCursor(tbl,dFields)
dOut=(startPointX,startPointY,5000,float(x))
curT.insertRow(dOut)
arcpy.BearingDistanceToLine_management(tbl,line,"X","Y","L","METERS","D","DEGREES","GEODESIC","#","GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119522E-09;0.001;0.001;IsHighPrecision")
theLine=arcpy.CopyFeatures_management(line,g)
newEndP=theLine[0].lastPoint
newEndX=newEndP.X
newEndY=newEndP.Y
dist=abs(newEndX-endPointX)+abs(newEndY-endPointY)
print x
return dist
try:
angle=gss(39.0,49.0,1e-8)
print 'Solution found %s' %angle

except:
tb = sys.exc_info()[2]
tbinfo = traceback.format_tb(tb)[0]
pymsg = "PYTHON ERRORS:\nTraceback Info:\n" + tbinfo + "\nError Info:\n    " + \
str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"
``````

Result:

• @simo.379209 have a look at update in my answer Oct 8 '15 at 0:44
• thanks for the update. I can see a bit more now, even though my py skills are almost non-existent. I will have a look into it. Oct 8 '15 at 1:13
• You do not need python at all. You have to understand algorithm (dividing section here) and use it in C. Your code will be a fraction of mine, because it took me 20 lines to replace single var endPoint = GeometryEngine.GeodesicMove(startPoint, 5000, LinearUnits.Meters, 45); Oct 8 '15 at 1:25

Use the .NET interface to GeographicLib. This is callable from C# and the Geodesic::Inverse method allows you to compute the distance and bearing given two points.

Surprisingly one of the other guys in the company needed the same thing on a completely different project in javascript.

His answer is to use Vincenty's Inverse method. There are a few implementations around:

javascript - http://www.movable-type.co.uk/scripts/latlong-vincenty.html scroll right down to the bottom.

There is also an Australian Government validation site in order to check your results http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp