Calculando latitude / longitude X milhas a partir do ponto?

46

Estou querendo encontrar um ponto de latitude e longitude, dado um rumo, uma distância e uma latitude e longitude inicial.

Parece ser o oposto desta pergunta ( distância entre pontos lat / long ).

Eu já examinei a fórmula do haversine e acho que a aproximação do mundo provavelmente está próxima o suficiente.

Estou assumindo que preciso resolver a fórmula de Haversine para meu lat / long desconhecido, isso está correto? Existem bons sites que falam sobre esse tipo de coisa? Parece que isso seria comum, mas meu Google apenas apresentou perguntas semelhantes à acima.

O que realmente estou procurando é apenas uma fórmula para isso. Gostaria de fornecer um lat / lng inicial, um rumo e uma distância (milhas ou quilômetros) e gostaria de extrair dele um par de lat / lng que representa onde um teria terminado se eles viajassem ao longo essa rota.

Jason Whitehorn
fonte
Você está procurando uma ferramenta que faça isso (como o pe.dll da Esri) ou uma fórmula?
Kirk Kuykendall
Desculpe, não fui específico ... Estou procurando uma fórmula. Vou atualizar minha pergunta para ser mais específico.
Jason Whitehorn
Muitas amostras de matemática elaboradas estão aqui <a href=" movable-type.co.uk/scripts/latlong.html"> Calcule a distância, o rumo e muito mais entre os pontos de latitude / longitude </a>, que inclui a solução para "Destino ponto dado a distância e rumo do ponto inicial ".
Stephen Quan
1
Intimamente relacionado: gis.stackexchange.com/questions/2951/… .
whuber
aqui está a página que apontam para latitude / longitude cálculos [latitude / longitude cálculos] ( movable-type.co.uk/scripts/latlong.html ) Também esta página Lat / cálculos longos , há um código + calculadora
Abo gaseem

Respostas:

21

Gostaria de saber como os resultados desta fórmula se comparam ao pe.dll da Esri .

( citação ).

Um ponto {lat, lon} é uma distância d no radial tc do ponto 1 se:

 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

Esse algoritmo é limitado a distâncias tais que dlon <pi / 2, ou seja, aquelas que se estendem por menos de um quarto da circunferência da Terra em longitude. Um algoritmo completamente geral, mas mais complicado, é necessário se maiores distâncias forem permitidas:

 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

Aqui está uma página html para teste .

Kirk Kuykendall
fonte
Obrigado pela resposta rápida. Deixe-me digerir algumas dessas informações e eu voltarei com você. Na superfície, porém, isso parece perfeito.
Jason Whitehorn
1
Tentei o caso direto usando pe.dll (na verdade libpe.so no solaris) depois de recuperar a distância e o azimute direto da página html para 32N, 117W a 40N, 82W. Usando 32N, 117W, d = 3255.056515890041, azi = 64.24498012065699, obtive 40N, 82W (na verdade 82.00000000064).
Mcknedy #
3
Impressionante! Muito obrigado pelo link para o artigo do Formulário de Aviação de Ed Williams, eu não tinha visto isso antes, mas até agora provou ser uma ótima leitura. Apenas uma observação para quem olhar para isso no futuro, as entradas e saídas desta fórmula são TODAS em radianos, até a distância.
Jason Whitehorn
1
Qual é a unidade de distância nesta fórmula?
Karthik Murugan
1
A introdução de @KarthikMurugan Ed diz que as unidades de distância estão em radianos ao longo de um grande círculo.
Kirk Kuykendall
18

Se você estivesse em um avião, então o ponto que é r metros de distância em um rolamento de uma graus a leste de North é deslocado por r * cos (a) na direção norte e r * sin (a) na direção leste. (Essas declarações definem mais ou menos o seno e o cosseno.)

Embora você não esteja em um avião - você está trabalhando na superfície de um elipsóide curvado que modela a superfície da Terra - qualquer distância inferior a algumas centenas de quilômetros cobre uma parte tão pequena da superfície que, para fins mais práticos, pode ser considerado plano. A única complicação restante é que um grau de longitude não cobre a mesma distância que um grau de latitude. Em um modelo esférico da Terra, um grau de longitude é apenas cos (latitude) contanto que um grau de latitude. (Em um modelo elipsóide, essa ainda é uma excelente aproximação, boa para cerca de 2,5 algarismos significativos.)

Finalmente, um grau de latitude é de aproximadamente 10 ^ 7/90 = 111.111 metros. Agora temos todas as informações necessárias para converter medidores em graus:

O deslocamento para o norte é r * cos (a) / 111111 graus;

O deslocamento para o leste é r * sin (a) / cos (latitude) / 111111 graus.

Por exemplo, a uma latitude de -0,31399 graus e um rolamento de a = 30 graus a leste do 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.

De onde, começando em (-78.4437, -0.31399), o novo local está em (-78.4437 + 0.00045, -0.31399 + 0.0007794) = (-78.4432, -0.313211).

Uma resposta mais precisa, no moderno sistema de referência ITRF00, é (-78.4433, -0.313207): está a 0,43 metros da resposta aproximada, indicando que a aproximação erra 0,43% neste caso. Para obter maior precisão, você deve usar fórmulas de distância elipsoidal (que são muito mais complicadas) ou uma projeção conforme de alta fidelidade com divergência zero (para que o rolamento esteja correto).

whuber
fonte
2
+1 para compreender corretamente o contexto matemático (ou seja, seu plano localmente)
Frank Conry
4

Se você precisar de uma solução JavaScript, considere estes functionse este violino :

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; }
};

Portanto, se você deseja calcular uma nova coordenada, pode ser assim:

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
fonte
2

Eu consegui isso trabalhando no ObjectiveC. A chave aqui é saber que você precisa de pontos lat / lng em radianos e precisa convertê-los novamente em lat / lng após aplicar a equação. Além disso, saiba que a distância e tc estão em radianos.

Aqui está a equação 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

Aqui, ele é implementado na ObjC, onde radiano é um radiano medido no sentido anti-horário de N (por exemplo, PI / 2 é W, PI é S, 2 PI / 3 é E) e a distância é em quilô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
fonte
Estou procurando uma solução em que eu queira obter 4 lat, longos desde o ponto em que estou a 50 milhas ao norte, 50 milhas a oeste, 50 milhas a leste e assim por diante ... Na resposta acima, você pode explicar como posso alcançar isto ?
Rahul Vyas
1

Se você estiver interessado em uma melhor precisão, existe o Vincenty . (O link é para o formulário "direto", que é exatamente o que você procura).

Existem algumas implementações existentes, se você estiver atrás do código.

Além disso, uma pergunta: você não está assumindo que o viajante mantém o mesmo rumo durante toda a viagem, não é? Nesse caso, esses métodos não estão respondendo à pergunta certa. (É melhor você projetar para o mercator, desenhar uma linha reta e depois não projetar o resultado.)

Dan S.
fonte
Muito boa pergunta, apesar da redação da minha pergunta que possa ter indicado que eu estava calculando um destino para um viajante, não sou. Bom ponto embora. Isso foi principalmente para que eu pudesse calcular uma área delimitadora (em uma ordem pequena, digamos 80 quilômetros) ... Espero que faça sentido.
Jason Whitehorn
gis.stackexchange.com/questions/3264/... tinha a mesma pergunta (construção de áreas delimitadas a partir de um ponto de distância &)
Dan S.
0

Aqui está uma solução 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
fonte
-2

Utilizo a abordagem descrita abaixo para determinar a próxima coordenada, considerando o rumo e a distância da coordenada anterior. Tenho problemas de precisão com outra abordagem que li na internet.

Eu uso isso para determinar a área da terra, que é um polígono, e traçar esse polígono no google earth. Um título de terra possui rolamentos e distâncias escritos da seguinte maneira: "NorthOrSouth x graus e minutos EastOrWest, z metros ao ponto n;".

Portanto, partindo das coordenadas fornecidas do ponto de referência, primeiro calculo a distância por um grau de latitude e um grau de longitude usando a fórmula haversine, pois isso varia dependendo da localização. Em seguida, determino a próxima coordenada a partir da fórmula trigonométrica seno e cosseno.

Abaixo está o 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
fonte
1
A fórmula cartesiana para a área de polígono que você tenta aplicar aqui não é aplicável a distâncias e ângulos calculados em uma superfície curva (como um esferóide). Essa fórmula comete um erro adicional usando latitude e longitude como se fossem coordenadas cartesianas. As únicas circunstâncias sob as quais seu uso poderia ser considerado seriam exatamente aquelas (para polígonos muito pequenos) em que a fórmula haversine é desnecessária de qualquer maneira. No geral, parece que esse código funciona muito demais para nenhum ganho.
whuber