Cálculo de latitud / longitud X millas desde el punto?

Quiero encontrar un punto de latitud y longitud dado un rumbo, una distancia y una latitud y longitud iniciales.

Esto parece ser lo opuesto a esta pregunta ( Distancia entre puntos lat / long ).

Ya he examinado la fórmula de Haversine y creo que su aproximación al mundo probablemente sea lo suficientemente cercana.

Supongo que necesito resolver la fórmula de Haversine para mi lat / long desconocido, ¿es esto correcto? ¿Hay algún buen sitio web que hable sobre este tipo de cosas? Parece que sería común, pero mi búsqueda en Google solo ha generado preguntas similares a la anterior.

Lo que realmente estoy buscando es solo una fórmula para esto. Me gustaría darle un lat / lng inicial, un rumbo y una distancia (millas o kilómetros) y me gustaría obtener un par lat / lng que represente dónde habría terminado uno si hubieran viajado esa ruta

Me gustaría saber cómo se comparan los resultados de esta fórmula con el pe.dll de Esri .

( cita ).

Un punto {lat, lon} es una distancia d fuera del radio tc desde el punto 1 si:

 lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 IF (cos(lat)=0)
    lon=lon1      // endpoint a pole
 ELSE
    lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
 ENDIF

Este algoritmo se limita a distancias tales que dlon <pi / 2, es decir, aquellas que se extienden alrededor de menos de un cuarto de la circunferencia de la tierra en longitud. Un algoritmo completamente general pero más complicado es necesario si se permiten distancias mayores:

 lat =asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat))
 lon=mod( lon1-dlon +pi,2*pi )-pi

Aquí hay una página html para probar .

Si se va en un plano, entonces el punto de que es r metros de distancia en un rumbo de unos grados al este del norte se desplaza por r * cos (a) en la dirección norte y r * sen (a) en la dirección este. (Estas declaraciones definen más o menos el seno y el coseno).

Aunque no está en un plano, está trabajando en la superficie de un elipsoide curvo que modela la superficie de la Tierra, cualquier distancia inferior a unos pocos cientos de kilómetros cubre una parte tan pequeña de la superficie que, para la mayoría de los propósitos prácticos, puede ser considerado plano La única complicación restante es que un grado de longitud no cubre la misma distancia que un grado de latitud. En un modelo esférico de la Tierra, un grado de longitud es solo cos (latitud) mientras que un grado de latitud. (En un modelo elipsoidal, esta sigue siendo una aproximación excelente, buena para aproximadamente 2.5 cifras significativas).

Finalmente, un grado de latitud es aproximadamente 10 ^ 7/90 = 111,111 metros. Ahora tenemos toda la información necesaria para convertir metros a grados:

El desplazamiento hacia el norte es r * cos (a) / 111111 grados;

El desplazamiento hacia el este es r * sin (a) / cos (latitud) / 111111 grados.

Por ejemplo, a una latitud de -0.31399 grados y una demora de a = 30 grados al este del norte, podemos calcular

cos(a) = cos(30 degrees) = cos(pi/6 radians) = Sqrt(3)/2 = 0.866025.
sin(a) = sin(30 degrees) = sin(pi/6 radians) = 1/2 = 0.5.
cos(latitude) = cos(-0.31399 degrees) = cos(-0.00548016 radian) = 0.999984984.
r = 100 meters.
east displacement = 100 * 0.5 / 0.999984984 / 111111 = 0.000450007 degree.
north displacement = 100 * 0.866025 / 111111 = 0.000779423 degree.

Por lo tanto, a partir de (-78.4437, -0.31399), la nueva ubicación está en (-78.4437 + 0.00045, -0.31399 + 0.0007794) = (-78.4432, -0.313211).

Una respuesta más precisa, en el moderno sistema de referencia ITRF00, es (-78.4433, -0.313207): está a 0.43 metros de la respuesta aproximada, indicando que la aproximación erra en 0.43% en este caso. Para lograr una mayor precisión, debe usar fórmulas de distancia elipsoidal (que son mucho más complicadas) o una proyección conforme de alta fidelidad con divergencia cero (para que el rumbo sea correcto).

Si necesita una solución de JavaScript, considere estos functionsy este violín :

var gis = {
  /**
  * All coordinates expected EPSG:4326
  * @param {Array} start Expected [lon, lat]
  * @param {Array} end Expected [lon, lat]
  * @return {number} Distance - meter.
  */
  calculateDistance: function(start, end) {
    var lat1 = parseFloat(start[1]),
        lon1 = parseFloat(start[0]),
        lat2 = parseFloat(end[1]),
        lon2 = parseFloat(end[0]);

    return gis.sphericalCosinus(lat1, lon1, lat2, lon2);
  },

  /**
  * All coordinates expected EPSG:4326
  * @param {number} lat1 Start Latitude
  * @param {number} lon1 Start Longitude
  * @param {number} lat2 End Latitude
  * @param {number} lon2 End Longitude
  * @return {number} Distance - meters.
  */
  sphericalCosinus: function(lat1, lon1, lat2, lon2) {
    var radius = 6371e3; // meters
    var dLon = gis.toRad(lon2 - lon1),
        lat1 = gis.toRad(lat1),
        lat2 = gis.toRad(lat2),
        distance = Math.acos(Math.sin(lat1) * Math.sin(lat2) +
            Math.cos(lat1) * Math.cos(lat2) * Math.cos(dLon)) * radius;

    return distance;
  },

  /**
  * @param {Array} coord Expected [lon, lat] EPSG:4326
  * @param {number} bearing Bearing in degrees
  * @param {number} distance Distance in meters
  * @return {Array} Lon-lat coordinate.
  */
  createCoord: function(coord, bearing, distance) {
    /** http://www.movable-type.co.uk/scripts/latlong.html
    * φ is latitude, λ is longitude, 
    * θ is the bearing (clockwise from north), 
    * δ is the angular distance d/R; 
    * d being the distance travelled, R the earth’s radius*
    **/
    var 
      radius = 6371e3, // meters
      δ = Number(distance) / radius, // angular distance in radians
      θ = gis.toRad(Number(bearing));
      φ1 = gis.toRad(coord[1]),
      λ1 = gis.toRad(coord[0]);

    var φ2 = Math.asin(Math.sin(φ1)*Math.cos(δ) + 
      Math.cos(φ1)*Math.sin(δ)*Math.cos(θ));

    var λ2 = λ1 + Math.atan2(Math.sin(θ) * Math.sin(δ)*Math.cos(φ1),
      Math.cos(δ)-Math.sin(φ1)*Math.sin(φ2));
    // normalise to -180..+180°
    λ2 = (λ2 + 3 * Math.PI) % (2 * Math.PI) - Math.PI; 

    return [gis.toDeg(λ2), gis.toDeg(φ2)];
  },
  /**
   * All coordinates expected EPSG:4326
   * @param {Array} start Expected [lon, lat]
   * @param {Array} end Expected [lon, lat]
   * @return {number} Bearing in degrees.
   */
  getBearing: function(start, end){
    var
      startLat = gis.toRad(start[1]),
      startLong = gis.toRad(start[0]),
      endLat = gis.toRad(end[1]),
      endLong = gis.toRad(end[0]),
      dLong = endLong - startLong;

    var dPhi = Math.log(Math.tan(endLat/2.0 + Math.PI/4.0) / 
      Math.tan(startLat/2.0 + Math.PI/4.0));

    if (Math.abs(dLong) > Math.PI) {
      dLong = (dLong > 0.0) ? -(2.0 * Math.PI - dLong) : (2.0 * Math.PI + dLong);
    }

    return (gis.toDeg(Math.atan2(dLong, dPhi)) + 360.0) % 360.0;
  },
  toDeg: function(n) { return n * 180 / Math.PI; },
  toRad: function(n) { return n * Math.PI / 180; }
};

Entonces, si desea calcular una nueva coordenada, puede ser así:

var start = [15, 38.70250];
var end = [21.54967, 38.70250];
var total_distance = gis.calculateDistance(start, end); // meters
var percent = 10;
// this can be also meters
var distance = (percent / 100) * total_distance;
var bearing = gis.getBearing(start, end);
var new_coord = gis.createCoord(icon_coord, bearing, distance);

Obtuve esto trabajando en ObjectiveC. La clave aquí es saber que necesita puntos lat / lng en radianes y debe convertirlos nuevamente a lat / lng después de aplicar la ecuación. Además, sepa que la distancia y tc están en radianes.

Aquí está la ecuación original:

 lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
 IF (cos(lat)=0)
    lon=lon1      // endpoint a pole
 ELSE
    lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
 ENDIF

Aquí se implementa en ObjC donde radian es un radian medido en sentido antihorario desde N (por ejemplo, PI / 2 es W, PI es S, 2 PI / 3 es E) y la distancia es en kilómetros.

+ (CLLocationCoordinate2D)displaceLatLng:(CLLocationCoordinate2D)coordinate2D withRadian:(double)radian
                            withDistance:(CGFloat)distance {
  double lat1Radians = coordinate2D.latitude * (M_PI / 180);
  double lon1Radians = coordinate2D.longitude * (M_PI / 180);
  double distanceRadians = distance / 6371;
  double lat = asin(sin(lat1Radians)*cos(distanceRadians)+cos(lat1Radians)*sin(distanceRadians)
      *cos(radian));
  double lon;
  if (cos(lat) == 0) {
    lon = lon1Radians;
  } else {
    lon = fmodf((float) (lon1Radians - asin(sin(radian) * sin(distanceRadians) / cos(lat)) + M_PI),
        (float) (2 * M_PI)) - M_PI;
  }
  double newLat = lat * (180 / M_PI);
  double newLon = lon * (180 / M_PI);
  return CLLocationCoordinate2DMake(newLat, newLon);
}

Si le interesa una mayor precisión, está Vincenty . (El enlace es a la forma 'directa', que es exactamente lo que buscas).

Hay bastantes implementaciones existentes, si buscas código.

Además, una pregunta: no estás asumiendo que el viajero mantiene la misma orientación durante todo el viaje, ¿verdad? Si es así, estos métodos no responden la pregunta correcta. (Sería mejor proyectar en mercator, dibujar una línea recta y luego des-proyectar el resultado).

Aquí hay una solución de Python:

    def displace(self, theta, distance):
    """
    Displace a LatLng theta degrees counterclockwise and some
    meters in that direction.
    Notes:
        http://www.movable-type.co.uk/scripts/latlong.html
        0 DEGREES IS THE VERTICAL Y AXIS! IMPORTANT!
    Args:
        theta:    A number in degrees.
        distance: A number in meters.
    Returns:
        A new LatLng.
    """
    theta = np.float32(theta)

    delta = np.divide(np.float32(distance), np.float32(E_RADIUS))

    def to_radians(theta):
        return np.divide(np.dot(theta, np.pi), np.float32(180.0))

    def to_degrees(theta):
        return np.divide(np.dot(theta, np.float32(180.0)), np.pi)

    theta = to_radians(theta)
    lat1 = to_radians(self.lat)
    lng1 = to_radians(self.lng)

    lat2 = np.arcsin( np.sin(lat1) * np.cos(delta) +
                      np.cos(lat1) * np.sin(delta) * np.cos(theta) )

    lng2 = lng1 + np.arctan2( np.sin(theta) * np.sin(delta) * np.cos(lat1),
                              np.cos(delta) - np.sin(lat1) * np.sin(lat2))

    lng2 = (lng2 + 3 * np.pi) % (2 * np.pi) - np.pi

    return LatLng(to_degrees(lat2), to_degrees(lng2))

Utilizo el enfoque descrito a continuación para determinar la siguiente coordenada dada la demora y la distancia desde la coordenada anterior. Tengo un problema de precisión con otro enfoque que leí en Internet.

Lo uso para determinar el área de tierra, que es un polígono, y trazo ese polígono en Google Earth. Un título de tierra tiene orientaciones y distancias escritas de esta manera: "Norte o Sur x grados y minutos Este u Oeste, z metros hasta el punto n;".

Entonces, a partir de las coordenadas dadas del punto de referencia, primero calculo la distancia por un grado de latitud y un grado de longitud usando la fórmula de Haversine porque esto varía según la ubicación. Luego determino la siguiente coordenada a partir de la fórmula de trigonometría seno y coseno.

A continuación se muestra el javascript:

var mapCenter = new google.maps.LatLng(referencePointLatitude, referencePointLongitude); //the ref point lat and lon must be given, usually a land mark (BLLM)
var latDiv = latDiv(mapCenter); //distance per one degree latitude in this location
var lngDiv = lngDiv(mapCenter); //distance per one degree longitude in this location
var LatLng2 = NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z); //next coordinate given the bearing and distance from previous coordinate
var Lat2 = LatLng2.lat(); //next coord latitude in degrees
var Lng2 = LatLng2.lng(); //next coord longitude in degrees
var polygon=[p1,p2,...,pn-1,pn,p1]; //p1,p2,etc. are the coordinates of points of the polygon, i.e. the land Title. Be sure to close the polygon to the point of beginning p1
var area = Area(polygon); //area of the polygon in sq.m.
function NextCoord(PrevCoord,NorthOrSouth,x,y,EastOrWest,z) {
  var angle = ( x + ( y / 60 ) ) * Math.PI / 180;
  var a = 1;
  var b = 1;
  if (NorthOrSouth == 'South') { a = -1; }
  if (EastOrWest == 'West') { b = -1; }
  var nextLat = PrevCoord.lat() +  ( a * z * Math.cos(angle) / latDiv );
  var nextLng = PrevCoord.lng() +  ( b * z * Math.sin(angle) / lngDiv );
  var nextCoord = new google.maps.LatLng(nextLat, nextLng);
  return nextCoord;
}
function latDiv(mapCenter) {
  var p1 = new google.maps.LatLng(mapCenter.lat()-0.5, mapCenter.lng());
  var p2 = new google.maps.LatLng(mapCenter.lat()+0.5, mapCenter.lng());
  return dist(p1,p2);
}
function lngDiv(mapCenter) {
  var p3 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()-0.5);
  var p4 = new google.maps.LatLng(mapCenter.lat(), mapCenter.lng()+0.5);
  return dist(p3,p4);
}
function dist(pt1, pt2) {
    var dLat  = ( pt2.lat() - pt1.lat() ) * Math.PI / 180;
    var dLng = ( pt2.lng() - pt1.lng() ) * Math.PI / 180;
    var a = Math.sin(dLat/2) * Math.sin(dLat/2) +                 
            Math.cos(rad(pt1.lat())) * Math.cos(rad(pt2.lat())) *
            Math.sin(dLng/2) * Math.sin(dLng/2);
    var R = 6372800; //earth's radius
    var distance = R * 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    return distance;
}
function Area(polygon) {
  var xPts=[];
  for (i=1; i<polygon.length; i++) {
    xPts[i-1] = ( polygon[i].lat() - polygon[0].lat() ) * latDiv;
  }
  var yPts=[];
  for (i=1; i<polygon.length; i++) {
    yPts[i-1] = ( polygon[i].lng() - polygon[0].lng() ) * lngDiv;
  }
  var area = 0;
  j = polygon.length-2;
  for (i=0; i<polygon.length-1; i++) {
    area = area +  ( xPts[j] + xPts[i] ) * ( yPts[j] - yPts[i] );
    j = i;
  }
  area = Math.abs(area/2);
  return area;
}