1

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

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:

enter image description here

  • @simo.379209 have a look at update in my answer – FelixIP 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. – simo.37920 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); – FelixIP Oct 8 '15 at 1:25
1

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.

0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.