calculate bearing from one MapPoint to another MapPoint

Question

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

Answer 1

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

# your function here
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:

Answer 2

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.

Answer 3

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.

C# - http://www.codeproject.com/Articles/19939/GPS-Receivers-Geodesy-and-Geocaching-Vincenty-s-Fo

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