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

46

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

Jason Whitehorn
fuente
¿Está buscando una herramienta que haga esto (como el pe.dll de Esri) o una fórmula?
Kirk Kuykendall
Lo siento, no fui específico ... Estoy buscando una fórmula. Actualizaré mi pregunta para que sea más específica.
Jason Whitehorn
Hay un montón de ejemplos matemáticos resueltos aquí <a href=" movable-type.co.uk/scripts/latlong.html"> Calcule la distancia, la demora y más entre los puntos de latitud / longitud </a> que incluye la solución para "Destino distancia dada punto y rumbo desde el punto inicial ".
Stephen Quan el
1
Muy relacionado: gis.stackexchange.com/questions/2951/… .
whuber
aquí está la página que enlaza con los cálculos lat / long [Lat / long cálculos] ( movable-type.co.uk/scripts/latlong.html ) también esta página Lat / long cálculos hay un código + calculadora
Abo gaseem

Respuestas:

21

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 .

Kirk Kuykendall
fuente
Gracias por la rápida respuesta. Permíteme digerir parte de esta información y volveré contigo. En la superficie, sin embargo, esto parece perfecto.
Jason Whitehorn
1
Probé el caso directo usando pe.dll (en realidad libpe.so en solaris) después de recuperar la distancia y el acimut directo desde la página html para 32N, 117W a 40N, 82W. Usando 32N, 117W, d = 3255.056515890041, azi = 64.24498012065699, obtuve 40N, 82W (en realidad 82.00000000064).
mkennedy
3
¡Increíble! Muchas gracias por el enlace al artículo del Formulario de Aviación de Ed Williams, no lo había visto antes, pero hasta ahora ha demostrado ser una gran lectura. Solo una nota para cualquiera que vea esto en el futuro, las entradas y salidas de esta fórmula están TODAS en radianes, incluso en la distancia.
Jason Whitehorn
1
¿Cuál es la unidad de distancia en esta fórmula?
Karthik Murugan
1
La introducción de @KarthikMurugan Ed dice que las unidades de distancia están en radianes a lo largo de un gran círculo.
Kirk Kuykendall
18

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

whuber
fuente
2
+1 para comprender adecuadamente el contexto matemático (es decir, su plano local)
Frank Conry
4

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);
Jonatas Walker
fuente
2

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);
}
BigCheesy
fuente
Estoy buscando una solución donde quiero obtener 4 lat, largos desde el punto donde estoy a 50 millas al norte, 50 millas al oeste, 50 millas al este, etc. En la respuesta anterior, ¿puede explicar cómo puedo lograrlo? esta ?
Rahul Vyas
1

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

Dan S.
fuente
Muy buena pregunta, a pesar de la redacción en mi pregunta que puede haber indicado que estaba calculando un destino para un viajero, no lo estoy. Buen punto sin embargo. Esto fue principalmente para poder calcular un área delimitadora (en un pedido pequeño, digamos 50 millas) ... Espero que tenga sentido.
Jason Whitehorn
gis.stackexchange.com/questions/3264/… tenía la misma pregunta (construir áreas de delimitación desde un punto y distancia)
Dan S.
0

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))
Liam Horne
fuente
-2

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&lt;polygon.length; i++) {
    xPts[i-1] = ( polygon[i].lat() - polygon[0].lat() ) * latDiv;
  }
  var yPts=[];
  for (i=1; i&lt;polygon.length; i++) {
    yPts[i-1] = ( polygon[i].lng() - polygon[0].lng() ) * lngDiv;
  }
  var area = 0;
  j = polygon.length-2;
  for (i=0; i&lt;polygon.length-1; i++) {
    area = area +  ( xPts[j] + xPts[i] ) * ( yPts[j] - yPts[i] );
    j = i;
  }
  area = Math.abs(area/2);
  return area;
}
Dante Valenzuela
fuente
1
La fórmula cartesiana para el área de polígono que intenta aplicar aquí no es aplicable a distancias y ángulos calculados en una superficie curva (como un esferoide). Esta fórmula comete un error adicional al usar latitud y longitud como si fueran coordenadas cartesianas. Las únicas circunstancias bajo las cuales podría considerarse su uso serían exactamente aquellas (para polígonos muy pequeños) donde la fórmula de haversina es innecesaria de todos modos. En general, parece que este código funciona demasiado duro sin ninguna ganancia.
whuber