Comment trouver le lat/long qui est x km au nord d'un lat/long donné?

Question

J'ai un code C# qui génère des google maps. Ce code examine tous les points que j'ai besoin de tracer sur la carte, puis calcule les limites d'un rectangle pour inclure ces points. Il transmet ensuite ces limites à l'API Google Maps pour définir correctement le niveau de zoom afin d'afficher tous les points de la carte.

Ce code fonctionne très bien mais j'ai une nouvelle exigence.

Un des points peut avoir une précision associée. Si c'est le cas, je dessine un cercle autour du point avec le rayon défini sur la valeur de précision. Encore une fois, cela fonctionne bien, mais la vérification de mes limites ne fait pas ce que je veux. Je veux que le cadre englobe le cercle complet.

Ceci nécessite un algorithme pour prendre un point x et calculer le point y qui serait z mètres au nord de x et aussi z mètres au sud de x.

Quelqu'un at-il cet algorithme, de préférence en C#. J'ai trouvé un algorithme générique here mais je ne l'ai pas implémenté correctement car les réponses que je reçois sont des milliers de km à la dérive.

Ceci est l'exemple générique

Lat/lon given radial and distance 

A point {lat,lon} is a distance d out on the tc radial from point 1 if: 

    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 

Et voici ma traduction C#.

// Extend a Point North/South by the specified distance 
    public static Point ExtendPoint(Point _pt, int _distance, int _bearing) 
    { 
     Decimal lat = 0.0; 
     Decimal lng = 0.0; 

     lat = Math.Asin(Math.Sin(_pt.Lat) * Math.Cos(_distance) + Math.Cos(_pt.Lat) * 
      Math.Sin(_distance) * Math.Cos(_bearing)); 

     if (Math.Cos(lat) == 0) 
     { 
      lng = _pt.Lng;  // endpoint a pole 
     } 
     else 
     { 
      lng = (
       (_pt.Lng - Math.Asin(Math.Sin(_bearing) * Math.Sin(_distance)/Math.Cos(lat)) 
       + Math.PI) % (2 * Math.PI)) - Math.PI; 
     } 

     ret = new Point(lat,lng); 
     return ret; 
    } 

J'appelle cette fonction avec un palier de 0 à calculer la nouvelle position nord et une valeur de 180 pour calculer la nouvelle position du sud. Est-ce que quelqu'un peut voir ce que j'ai fait de mal ou peut-être fournir un algorithme de travail connu?

Answer 1

Si vous avez une latitude et la longitude, vous pouvez calculer la latitude et la longitude correctes d'un changement x km en latitude comme ceci:

new-lat = ((old-km-north + x-km-change)/40,075) * 360) 
     ^is the ratio of the    ^times the ratio of the circle 
      of the earth the change    by 360 to get the total ratio 
      covers.        covered in degrees. 

La même chose peut appliquer à la longitude. Si vous avez la distance totale plus le changement, vous pouvez calculer le total des degrés d'une manière similaire.

new-long = ((old-km-east + x-km-change)/40,075) * 360) 
     ^is the ratio of the    ^times the ratio of the circle 
      of the earth the change    by 360 to get the total ratio 
      covers.        covered in degrees. 

Encore une fois, ces calculs doivent travailler, mais je suis en cours d'exécution de l'intuition pure, mais la logique ne semble vrai.

Edit: Comme l'a souligné Skizz 40,075 doit être ajusté à la circonférence de la terre à une latitude donnée en utilisant 2.pi.r.cos (LAT) ou 40074.cos (LAT)

Answer 2

Il est plus précis si vous le reprojetez d'abord à UTM puis vérifiez la distance.

J'espère que cela aide

Answer 3

J'ai un morceau de code très similaire. Cela m'a donné des résultats très proches par rapport à une autre implémentation.

Je pense que le problème avec le vôtre est que vous utilisez la "distance" comme distance linéaire en mètres au lieu de la distance angulaire en radians.

/// <summary> 
/// Calculates the end-point from a given source at a given range (meters) and bearing (degrees). 
/// This methods uses simple geometry equations to calculate the end-point. 
/// </summary> 
/// <param name="source">Point of origin</param> 
/// <param name="range">Range in meters</param> 
/// <param name="bearing">Bearing in degrees</param> 
/// <returns>End-point from the source given the desired range and bearing.</returns> 
public static LatLonAlt CalculateDerivedPosition(LatLonAlt source, double range, double bearing) 
{ 
    double latA = source.Latitude * UnitConstants.DegreesToRadians; 
    double lonA = source.Longitude * UnitConstants.DegreesToRadians; 
    double angularDistance = range/GeospatialConstants.EarthRadius; 
    double trueCourse = bearing * UnitConstants.DegreesToRadians; 

    double lat = Math.Asin(
     Math.Sin(latA) * Math.Cos(angularDistance) + 
     Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse)); 

    double dlon = Math.Atan2(
     Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA), 
     Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat)); 

    double lon = ((lonA + dlon + Math.PI) % UnitConstants.TwoPi) - Math.PI; 

    return new LatLonAlt(
     lat * UnitConstants.RadiansToDegrees, 
     lon * UnitConstants.RadiansToDegrees, 
     source.Altitude); 
} 

Où

public const double EarthRadius = 6378137.0; // WGS-84 ellipsoid parameters 

et LatLonAlt est en degrés/mètres (conversion a lieu en interne). Ajuster au besoin.

Je suppose que vous pouvez comprendre ce que la valeur est UnitConstants.DegreesToRadians :)

Answer 4

Je ne sais pas si je manque quelque chose ici, mais je pense que la question pourrait être reformulée comme «J'ai lat/lon point, et je veux trouver le point x mètres nord et x mètres sud de ce point. "

Si c'est la question alors vous n'avez pas besoin de trouver une nouvelle longitude (ce qui rend les choses plus simples), vous avez juste besoin d'une nouvelle latitude. Un degré de latitude est d'environ 60 milles marins partout sur Terre, et un mille nautique est de 1 852 mètres. Ainsi, pour de nouvelles latitudes x mètres au nord et au sud:

north_lat = lat + x/(1852 * 60) 
north_lat = min(north_lat, 90) 

south_lat = lat - x/(1852 * 60) 
south_lat = max(south_lat, -90) 

Ce n'est pas tout à fait exact parce que la Terre est pas une sphère parfaite avec exactement 60 miles nautiques entre chaque degré de latitude. Cependant, les autres réponses supposent que les lignes de latitude sont équidistantes, donc je suppose que vous vous en fichez. Si vous êtes intéressé à combien erreur qui pourrait introduire, il y a une belle table sur Wikipédia qui montre « la distance de la surface par 1 ° changement de latitude » pour différentes latitudes à ce lien:

http://en.wikipedia.org/wiki/Latitude#Degree_length

Answer 5

Il y a problèmes avec les deux équations sur Ed William's rather awesome site ... mais je ne les ai pas analysés pour voir pourquoi.

Une troisième équation qui I found here semble donner des résultats corrects.

Voici le cas de test en php ... la troisième équation est correcte, les deux premiers donnent des valeurs follement incorrectes pour la longitude.

<?php 
      $lon1 = -108.553412; $lat1 = 35.467155; $linDistance = .5; $bearing = 170; 
      $lon1 = deg2rad($lon1); $lat1 = deg2rad($lat1); 
      $distance = $linDistance/6371; // convert dist to angular distance in radians 
      $bearing = deg2rad($bearing); 

      echo "lon1: " . rad2deg($lon1) . " lat1: " . rad2deg($lat1) . "<br>\n"; 

// doesn't work 
      $lat2 = asin(sin($lat1) * cos($distance) + cos($lat1) * sin($distance) * cos($bearing)); 
      $dlon = atan2(sin($bearing) * sin($distance) * cos($lat1), cos($distance) - sin($lat1) * sin($lat2)); 
      $lon2 = (($lon1 - $dlon + M_PI) % (2 * M_PI)) - M_PI; // normalise to -180...+180 

      echo "lon2: " . rad2deg($lon2) . " lat2: " . rad2deg($lat2) . "<br>\n"; 

// same results as above 
      $lat3 = asin((sin($lat1) * cos($distance)) + (cos($lat1) * sin($distance) * cos($bearing))); 
      $lon3 = (($lon1 - (asin(sin($bearing) * sin($distance)/cos($lat3))) + M_PI) % (2 * M_PI)) - M_PI; 

      echo "lon3: " . rad2deg($lon3) . " lat3: " . rad2deg($lat3) . "<br>\n"; 

// gives correct answer... go figure 
      $lat4 = asin(sin($lat1) * cos($linDistance/6371) + cos($lat1) * sin($linDistance/6371) * cos($bearing)); 
      $lon4 = $lon1 + atan2((sin($bearing) * sin($linDistance/6371) * cos($lat1)), (cos($linDistance/6371) - sin($lat1) * sin($lat2))); 

      echo "lon4: " . rad2deg($lon4) . " lat4: " . rad2deg($lat4) . "<br>\n"; 
?> 

Remarque J'ai reçu par e-mail de l'auteur (Ed Williams) des deux premières équations:

De mes "notes de mise en œuvre":

Note sur la fonction mod. Cela semble être mis en œuvre différemment dans différentes langues, avec des conventions différentes sur si le signe du résultat suit le signe du diviseur ou le dividende. (Nous voulons que le signe suive le diviseur ou soit euclidien.) Fmod C et Java% ne fonctionnent pas.) Dans ce document, Mod (y, x) est le reste en divisant y par x et toujours se trouve dans la gamme 0 < = mod < x. Par exemple: mod (. 2.3,2) = 0,3 et mod (. -2.3,2) = 1,7

Si vous avez une fonction de plancher (int dans Excel), qui retourne plancher (x) = « le plus grand nombre entier inférieur ou égal à x "par exemple plancher (-2,3) = - 3 et plancher (2.3) = 2

mod(y,x) = y - x*floor(y/x) 

Ce qui suit devrait fonctionner en l'absence d'un plancher Fonction-, peu importe si "int" ou tronque tours vers le bas:

mod=y - x * int(y/x) 
if (mod < 0) mod = mod + x 

php est comme fmod en C et le fait "mal" pour mes fins.

Answer 6

Pour ce que ça vaut, j'ai un exemple en PHP qui peut faire ce que l'OP demande. Dans mon exemple, il dessine une boîte autour d'une coordonnée lat/long de départ, mais le code peut facilement être utilisé pour obtenir un seul point, X nombre de km ou miles.

http://www.richardpeacock.com/blog/2011/11/draw-box-around-coordinate-google-maps-based-miles-or-kilometers

Answer 7

Pour les paresseux, (comme moi;)) une solution de copier-coller, la version de Erich Mirabal avec des changements mineurs:

using System.Device.Location; // add reference to System.Device.dll 
public static class GeoUtils 
{ 
    /// <summary> 
    /// Calculates the end-point from a given source at a given range (meters) and bearing (degrees). 
    /// This methods uses simple geometry equations to calculate the end-point. 
    /// </summary> 
    /// <param name="source">Point of origin</param> 
    /// <param name="range">Range in meters</param> 
    /// <param name="bearing">Bearing in degrees</param> 
    /// <returns>End-point from the source given the desired range and bearing.</returns> 
    public static GeoCoordinate CalculateDerivedPosition(this GeoCoordinate source, double range, double bearing) 
    { 
     var latA = source.Latitude * DegreesToRadians; 
     var lonA = source.Longitude * DegreesToRadians; 
     var angularDistance = range/EarthRadius; 
     var trueCourse = bearing * DegreesToRadians; 

     var lat = Math.Asin(
      Math.Sin(latA) * Math.Cos(angularDistance) + 
      Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse)); 

     var dlon = Math.Atan2(
      Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA), 
      Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat)); 

     var lon = ((lonA + dlon + Math.PI) % (Math.PI*2)) - Math.PI; 

     return new GeoCoordinate(
      lat * RadiansToDegrees, 
      lon * RadiansToDegrees, 
      source.Altitude); 
    } 

    private const double DegreesToRadians = Math.PI/180.0; 
    private const double RadiansToDegrees = 180.0/ Math.PI; 
    private const double EarthRadius = 6378137.0; 
} 

Utilisation:

[TestClass] 
public class CalculateDerivedPositionUnitTest 
{ 
    [TestMethod] 
    public void OneDegreeSquareAtEquator() 
    { 
     var center = new GeoCoordinate(0, 0); 
     var radius = 111320; 
     var southBound = center.CalculateDerivedPosition(radius, -180); 
     var westBound = center.CalculateDerivedPosition(radius, -90); 
     var eastBound = center.CalculateDerivedPosition(radius, 90); 
     var northBound = center.CalculateDerivedPosition(radius, 0); 

     Console.Write($"leftBottom: {southBound.Latitude} , {westBound.Longitude} rightTop: {northBound.Latitude} , {eastBound.Longitude}"); 
    } 
}