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

43

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

Esto parece ser lo contrario de esta pregunta ( Distancia entre lat / puntos ).

Ya he examinado la fórmula de haversine y creo que es probable que la aproximación del mundo sea lo suficientemente cercana.

Supongo que debo resolver la fórmula de haversine para mi latitud / longitud desconocida, ¿es correcto? ¿Hay buenos sitios web que hablan de este tipo de cosas? Parece que sería común, pero mi Google solo ha encontrado preguntas similares a la anterior.

Lo que realmente estoy buscando es solo una fórmula para esto. Me gustaría darle una latitud inicial, una marcación y una distancia (millas o kilómetros) y me gustaría sacar de ella una pareja lat / nng que represente donde uno hubiera terminado si hubiera viajado esa ruta.

    
pregunta Jason Whitehorn 04.02.2011 - 20:03

7 respuestas

20

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

( cita ).

  

Un punto {lat, lon} es una distancia d fuera   en el tc radial 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 está limitado a distancias   tal que dlon < pi / 2, es decir, aquellos que   extenderse alrededor de menos de un cuarto de   la circunferencia de la tierra en   longitud. Un completamente general, pero   algoritmo más complicado es   necesario si son mayores distancias   permitido:

 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 realizar pruebas .

    
respondido por el Kirk Kuykendall 04.02.2011 - 20:35
14

Si estaba en un avión, entonces el punto que está r a una distancia de a grados al este del norte es desplazado por r * cos (a) en la dirección norte y r * sin (a) en la dirección este. (Estas declaraciones más o menos definen 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 casos, es práctica. Para fines puede ser considerado plano. La única complicación que queda 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 excelente aproximación, de 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 en 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 marcación 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): esto está a 0.43 metros de la respuesta aproximada, lo que indica que la aproximación es de 0.43% en este caso. Para lograr una mayor precisión, debe utilizar fórmulas de distancia elipsoidales (que son mucho más complicadas) o una proyección conforme de alta fidelidad con divergencia cero (para que el rumbo sea correcto).

    
respondido por el whuber 21.07.2012 - 19:06
3

Si necesita una solución de JavaScript, considere estos functions y 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 quieres 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);
    
respondido por el Jonatas Walker 17.06.2015 - 16:26
2

Conseguí esto trabajando en ObjectiveC. La clave aquí es saber que necesitas puntos de latencia / latencia en radianes y debes convertirlos de nuevo a latencia / lng después de aplicar la ecuación. Además, debes saber 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 el radián es un radián medido en sentido contrario a las agujas del reloj 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);
}
    
respondido por el BigCheesy 15.07.2015 - 18:32
1

Si está interesado en una mayor precisión, hay Vincenty . (El enlace está en la forma "directa", que es exactamente lo que está buscando).

Hay bastantes implementaciones existentes, si te interesa el código.

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

    
respondido por el Dan S. 05.02.2011 - 23:02
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))
    
respondido por el Liam Horne 08.07.2015 - 17:27
-2

Utilizo el enfoque que se describe a continuación para determinar la siguiente coordenada dado el rumbo y la distancia desde la coordenada anterior. Tengo problemas de precisión con otro enfoque que leo de Internet.

Utilizo esto para determinar el área de tierra, que es un polígono, y graficar ese polígono en google earth. El Título de un terreno tiene marcadas las orientaciones y las distancias de esta manera: "NorthOrSouth x grados y minutos EastOrWest, z metros al punto n;".

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

Debajo está 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;
}
    
respondido por el Dante Valenzuela 28.06.2014 - 14:34

Lea otras preguntas en las etiquetas