calcular el rumbo de un MapPoint a otro MapPoint

Question

Estoy utilizando el paquete Esri.ArcGISRuntime dentro de un programa C# que estamos escribiendo.

Necesito crear formas específicas a partir de un punto de referencia y proyectarlas en diferentes direcciones. Para ello es fácil crear un nuevo punto a 5.000 metros de distancia en un rumbo de 45 grados. es decir

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

Más adelante, cuando necesite hacer ingeniería inversa de los puntos, puedo recalcular fácilmente la distancia entre estos puntos utilizando

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

Sin embargo, no puedo encontrar ninguna forma de determinar cuál era el ángulo o la orientación original entre estos 2 puntos.

He intentado usar Arctan2, pero las coordenadas de MapPoint están sesgadas por la curvatura de la tierra (creo) y por lo tanto el cálculo directo es incorrecto.

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

Lo anterior devuelve una cifra de aproximadamente 39 grados (PS He intentado cambiar X e Y alrededor también, no es mejor).

Answer 1

Tome 39 como primera estimación, 29 como límite inferior, 49 como límite superior. Intente minimizar la distancia entre el punto final conocido y el nuevo punto final en el rumbo dado usando la primera ecuación, vea cómo funciona la sección dorada https://en.wikipedia.org/wiki/Golden_section_search

Actualización de la respuesta original No tengo C# etc, así que esta es la solución usando 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"

Resultado:

Answer 2

Utilice el Interfaz .NET a GeographicLib . Esto se puede llamar desde C# y el Geodésica::Inversa permite calcular la distancia y el rumbo dados dos puntos.

Answer 3

Sorprendentemente uno de los otros chicos de la empresa necesitaba lo mismo en un proyecto completamente diferente en javascript.

Su respuesta es utilizar el método inverso de Vincenty. Hay algunas implementaciones por ahí:

javascript - http://www.movable-type.co.uk/scripts/latlong-vincenty.html desplácese hasta el final.

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

También hay un sitio de validación del Gobierno australiano para comprobar sus resultados http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp