Calculer la latitude / longitude à X milles du point?

46

Je veux trouver un point de latitude et de longitude étant donné un relèvement, une distance et une latitude et une longitude de départ.

Cela semble être le contraire de cette question ( distance entre les points lat / long ).

J'ai déjà examiné la formule d'Hversine et pense que son approximation du monde est probablement assez proche.

Je suppose que je dois résoudre la formule de haversine pour mon lat / long inconnu, est-ce correct? Existe-t-il de bons sites Web qui parlent de ce genre de chose? Cela semble être courant, mais ma recherche sur Google n'a posé que des questions semblables à celle ci-dessus.

Ce que je recherche vraiment n’est qu’une formule pour cela. J'aimerais lui donner un point de départ, un relèvement et une distance (milles ou kilomètres) et j'aimerais en sortir une paire de points qui représente le lieu où l'on se serait retrouvé cette route.

Jason Whitehorn
la source
Vous recherchez un outil qui fait cela (comme peri d'Esri) ou une formule?
Kirk Kuykendall
Désolé je n'ai pas été spécifique ... Je cherche une formule. Je vais mettre à jour ma question pour être plus précis.
Jason Whitehorn
Voici quelques- uns des exemples mathématiques élaborés: <a href=" movable-type.co.uk/scripts/latlong.html"> Calculer la distance, le relèvement et bien plus entre les points de latitude / longitude </a>, ce qui inclut la solution "Destination point donné distance et relèvement du point de départ ".
Stephen Quan le
1
Étroitement liée: gis.stackexchange.com/questions/2951/… .
whuber
voici la page qui renvoie aux calculs lat / long [calculs lat / long] ( movable-type.co.uk/scripts/latlong.html ) également cette page calculs lat / long il y a un code + calculatrice
Abo gaseem

Réponses:

21

Je serais curieux de savoir comment les résultats de cette formule se comparent à pe.dll d' Esri .

( citation ).

Un point {lat, lon} est une distance d sur la radiale tc à partir du point 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

Cet algorithme est limité aux distances telles que dlon <pi / 2, c'est-à-dire celles qui s'étendent sur moins d'un quart de la circonférence de la Terre en longitude. Un algorithme tout à fait général, mais plus compliqué, est nécessaire si de plus grandes distances sont autorisées:

 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

Voici une page HTML pour les tests .

Kirk Kuykendall
la source
Merci pour la réponse rapide. Laissez-moi digérer certaines de ces informations et je reviendrai avec vous. En surface, cependant, cela semble parfait.
Jason Whitehorn
1
J'ai essayé le cas direct en utilisant pe.dll (en fait libpe.so sur solaris) après avoir récupéré la distance et l'azimut avant à partir de la page html pour 32N, 117W à 40N, 82W. En utilisant 32N, 117W, d = 3255.056515890041, azi = 64.24498012065699, je suis 40N, 82W (en réalité 82.00000000064).
Mkennedy
3
Impressionnant! Merci beaucoup pour le lien vers l'article de Ed Williams sur le formulaire Aviation, je ne l'avais jamais vu auparavant, mais il s'est jusqu'à présent avéré être une excellente lecture. Remarque pour ceux qui regarderont cela à l'avenir, les entrées et les sorties de cette formule sont toutes exprimées en radians, même la distance.
Jason Whitehorn
1
Quelle est l'unité de distance dans cette formule?
Karthik Murugan
1
Dans l' intro de @KarthikMurugan Ed, les unités de distance sont exprimées en radians le long d'un grand cercle.
Kirk Kuykendall
18

Si vous étiez dans un plan, le point est r mètres à un palier d' un degrés à l' est du nord est déplacé par r * cos (a) dans la direction du nord et r * sin (a) dans la direction est. (Ces déclarations définissent plus ou moins le sinus et le cosinus.)

Bien que vous ne soyez pas dans un avion (vous travaillez à la surface d'un ellipsoïde incurvé qui modélise la surface de la Terre), toute distance inférieure à quelques centaines de kilomètres couvre une si petite partie de la surface que, dans la plupart des cas, elle peut être considéré comme plat. La seule complication qui reste est qu'un degré de longitude ne couvre pas la même distance qu'un degré de latitude. Dans un modèle terrestre sphérique, un degré de longitude est seulement cos (latitude) aussi longtemps qu'un degré de latitude. (Dans un modèle ellipsoïdal, il s'agit toujours d'une excellente approximation, bonne à environ 2,5 chiffres significatifs.)

Enfin, un degré de latitude correspond à environ 10 ^ 7/90 = 111 111 mètres. Nous avons maintenant toutes les informations nécessaires pour convertir les mètres en degrés:

Le déplacement vers le nord est r * cos (a) / 111111 degrés;

Le déplacement vers l' est est r * sin (a) / cos (latitude) / 111111 degrés.

Par exemple, à une latitude de -0,31399 degrés et un cap a = 30 degrés nord-est, nous pouvons calculer

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 là, à partir de (-78,4437, -0,31399), le nouvel emplacement est situé à (-78,4437 + 0,00045, -0,31399 + 0,0007794) = (-78,4432, -0,313211).

Une réponse plus précise, dans le système de référence ITRF00 moderne, est (-78,4433, -0,313207): elle se situe à 0,43 mètre de la réponse approximative, ce qui indique que l’approximation est erronée de 0,43% dans ce cas. Pour obtenir une précision supérieure, vous devez utiliser une formule de distance ellipsoïdale (beaucoup plus compliquée) ou une projection conforme haute fidélité avec divergence nulle (pour que le relèvement soit correct).

whuber
la source
2
+1 pour bien comprendre le contexte mathématique (c'est-à-dire son plan local)
Frank Conry
4

Si vous avez besoin d'une solution JavaScript, tenez compte de ces éléments functionset de ce violon :

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

Donc, si vous voulez calculer une nouvelle coordonnée, cela peut ressembler à:

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
la source
2

Je travaille dans ObjectiveC. La clé ici est de savoir que vous avez besoin de points de lat / lng en radians et que vous devez les reconvertir en lat / lng après l’application de l’équation. Sachez également que la distance et tc sont exprimés en radians.

Voici l'équation originale:

 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

Ici, il est implémenté dans ObjC où radian est un radian mesuré dans le sens anti-horaire à partir de N (par exemple, PI / 2 est W, PI est S, 2 PI / 3 est E) et la distance est en kilomètres.

+ (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
la source
Je suis à la recherche d'une solution pour laquelle je souhaite obtenir 4 lat, à partir de l'endroit où je me trouve à 50 milles au nord, 50 milles à l'ouest, 50 milles à l'est, etc. ... Dans la réponse ci-dessus, pouvez-vous expliquer comment je peux atteindre cette ?
Rahul Vyas
1

Si vous êtes intéressé par une meilleure précision, il y a Vincenty . (Lien est à la forme «directe», qui est exactement ce que vous recherchez).

Il y a pas mal d'implémentations existantes, si vous voulez du code.

En outre, une question: vous ne supposez pas que le voyageur conserve le même cap pendant tout le voyage, n'est-ce pas? Si tel est le cas, ces méthodes ne répondent pas à la bonne question. (Vous feriez mieux de projeter sur mercator, de tracer une ligne droite, puis de projeter le résultat.)

Dan S.
la source
Très bonne question, malgré le libellé de ma question qui aurait pu indiquer que je calculais une destination pour un voyageur, je ne le suis pas. Bon point cependant. C'était principalement pour que je puisse calculer une zone de délimitation (sur une petite commande, disons 50 miles) ... J'espère que cela a du sens.
Jason Whitehorn
gis.stackexchange.com/questions/3264/ avait la même question (construction de zones de délimitation à partir d'un point et d'une distance)
Dan S.
0

Voici une solution 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
la source
-2

J'utilise l'approche décrite ci-dessous pour déterminer la coordonnée suivante en fonction du relèvement et de la distance par rapport à la coordonnée précédente. J'ai un problème de précision avec une autre approche que je lis sur Internet.

J'utilise ceci pour déterminer la superficie de la terre, qui est un polygone, et trace ce polygone dans Google Earth. Un titre foncier a des relèvements et des distances écrits de cette manière: "NordOuSud x degrés y minutes EstOuOuest, z mètres au point n;".

Donc, en partant des coordonnées données du point de référence, je calcule d’abord la distance par degré de latitude et de longitude en utilisant la formule de haversine car elle varie en fonction de l’emplacement. Ensuite, je détermine la coordonnée suivante à partir des formules de trigonométrie sinus et cosinus.

Voici le 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
la source
1
La formule cartésienne pour la surface de polygone que vous tentez d'appliquer ici ne s'applique pas aux distances et aux angles calculés sur une surface courbe (telle qu'un sphéroïde). Cette formule commet une erreur supplémentaire en utilisant la latitude et la longitude comme s'il s'agissait de coordonnées cartésiennes. Les seules circonstances dans lesquelles son utilisation pourrait être envisagée seraient exactement celles (pour les très petits polygones) où la formule de haversine est de toute façon inutile. Globalement, il semble que ce code fonctionne beaucoup trop durement sans gain.
whuber