Comment calculer la boîte englobante pour un emplacement lat / lng donné?

101

J'ai donné un emplacement défini par la latitude et la longitude. Maintenant, je veux calculer une boîte englobante à moins de 10 kilomètres par exemple de ce point.

La boîte englobante doit être définie comme latmin, lngmin et latmax, lngmax.

J'ai besoin de ce matériel pour utiliser l' API panoramio .

Est-ce que quelqu'un connaît la formule pour obtenir ces points?

Edit: Les gars, je cherche une formule / fonction qui prend lat & lng comme entrée et renvoie une boîte englobante comme latmin & lngmin et latmax & latmin. Mysql, php, c #, javascript est très bien mais aussi le pseudocode devrait être correct.

Edit: je ne cherche pas une solution qui me montre la distance de 2 points

Michal
la source
Si vous utilisez une géodatabase quelque part, elles ont sûrement un calcul de boîte englobante intégré. Vous pouvez même aller vérifier la source de PostGIS / GEOS, par exemple.
Vinko Vrsalovic

Réponses:

61

Je suggère d'approcher localement la surface de la Terre comme une sphère de rayon donné par l'ellipsoïde WGS84 à la latitude donnée. Je soupçonne que le calcul exact de latMin et latMax exigerait des fonctions elliptiques et ne donnerait pas une augmentation appréciable de la précision (WGS84 est lui-même une approximation).

Ma mise en œuvre suit (elle est écrite en Python; je ne l'ai pas testée):

# degrees to radians
def deg2rad(degrees):
    return math.pi*degrees/180.0
# radians to degrees
def rad2deg(radians):
    return 180.0*radians/math.pi

# Semi-axes of WGS-84 geoidal reference
WGS84_a = 6378137.0  # Major semiaxis [m]
WGS84_b = 6356752.3  # Minor semiaxis [m]

# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
def WGS84EarthRadius(lat):
    # http://en.wikipedia.org/wiki/Earth_radius
    An = WGS84_a*WGS84_a * math.cos(lat)
    Bn = WGS84_b*WGS84_b * math.sin(lat)
    Ad = WGS84_a * math.cos(lat)
    Bd = WGS84_b * math.sin(lat)
    return math.sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) )

# Bounding box surrounding the point at given coordinates,
# assuming local approximation of Earth surface as a sphere
# of radius given by WGS84
def boundingBox(latitudeInDegrees, longitudeInDegrees, halfSideInKm):
    lat = deg2rad(latitudeInDegrees)
    lon = deg2rad(longitudeInDegrees)
    halfSide = 1000*halfSideInKm

    # Radius of Earth at given latitude
    radius = WGS84EarthRadius(lat)
    # Radius of the parallel at given latitude
    pradius = radius*math.cos(lat)

    latMin = lat - halfSide/radius
    latMax = lat + halfSide/radius
    lonMin = lon - halfSide/pradius
    lonMax = lon + halfSide/pradius

    return (rad2deg(latMin), rad2deg(lonMin), rad2deg(latMax), rad2deg(lonMax))

EDIT: Le code suivant convertit (degrés, nombres premiers, secondes) en degrés + fractions de degré, et vice versa (non testé):

def dps2deg(degrees, primes, seconds):
    return degrees + primes/60.0 + seconds/3600.0

def deg2dps(degrees):
    intdeg = math.floor(degrees)
    primes = (degrees - intdeg)*60.0
    intpri = math.floor(primes)
    seconds = (primes - intpri)*60.0
    intsec = round(seconds)
    return (int(intdeg), int(intpri), int(intsec))
Federico A. Ramponi
la source
4
Comme indiqué dans la documentation de la bibliothèque CPAN suggérée, cela n'a de sens que pour halfSide <= 10km.
Federico A. Ramponi
1
Cela fonctionne-t-il près des pôles? Cela ne semble pas le cas, car il semble que cela se termine par latMin <-pi (pour le pôle sud) ou latMax> pi (pour le pôle nord)? Je pense que lorsque vous êtes à moins de la moitié d'un poteau, vous devez renvoyer une boîte englobante qui comprend toutes les longitudes et les latitudes calculées normalement pour le côté éloigné du poteau et au poteau du côté près du poteau.
Doug McClean
1
Voici une implémentation PHP de la spécification trouvée sur JanMatuschek.de: github.com/anthonymartin/GeoLocation.class.php
Anthony Martin
2
J'ai ajouté une implémentation C # de cette réponse ci-dessous.
Ε Г И І И О
2
@ FedericoA.Ramponi quel est le haldSideinKm ici? ne comprends pas ... ce que je dois passer dans ces argyments, le rayon entre deux points de la carte ou quoi?
53

J'ai écrit un article sur la recherche des coordonnées de délimitation:

http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates

L'article explique les formules et fournit également une implémentation Java. (Cela montre également pourquoi la formule de Federico pour la longitude min / max est inexacte.)

Jan Philip Matuschek
la source
4
J'ai créé un port PHP de votre classe GeoLocation. Il peut être trouvé ici: pastie.org/5416584
Anthony Martin
1
Je l'ai téléchargé sur github maintenant: github.com/anthonymartin/GeoLocation.class.php
Anthony Martin
1
Cela répond-il même à la question? Si nous n'avons qu'un seul point de départ, nous ne pouvons pas calculer la distance du grand cercle comme cela est fait dans ce code, qui nécessite deux emplacements latéraux longs.
mdoran3844
il y a un mauvais code dans votre variante C #, par exemple:, public override string ToString()il est très mauvais de remplacer une telle méthode globale uniquement dans un seul but, mieux vaut simplement ajouter une autre méthode, puis remplacer la méthode standard, qui peut être utilisée dans d'autres parties de l'application, pas pour le gis exact ...
Voici un lien mis à jour vers le port PHP de la classe GeoLocaiton de Jan: github.com/anthonymartin/GeoLocation.php
Anthony Martin
34

Ici, j'ai converti la réponse de Federico A. Ramponi en C # pour toute personne intéressée:

public class MapPoint
{
    public double Longitude { get; set; } // In Degrees
    public double Latitude { get; set; } // In Degrees
}

public class BoundingBox
{
    public MapPoint MinPoint { get; set; }
    public MapPoint MaxPoint { get; set; }
}        

// Semi-axes of WGS-84 geoidal reference
private const double WGS84_a = 6378137.0; // Major semiaxis [m]
private const double WGS84_b = 6356752.3; // Minor semiaxis [m]

// 'halfSideInKm' is the half length of the bounding box you want in kilometers.
public static BoundingBox GetBoundingBox(MapPoint point, double halfSideInKm)
{            
    // Bounding box surrounding the point at given coordinates,
    // assuming local approximation of Earth surface as a sphere
    // of radius given by WGS84
    var lat = Deg2rad(point.Latitude);
    var lon = Deg2rad(point.Longitude);
    var halfSide = 1000 * halfSideInKm;

    // Radius of Earth at given latitude
    var radius = WGS84EarthRadius(lat);
    // Radius of the parallel at given latitude
    var pradius = radius * Math.Cos(lat);

    var latMin = lat - halfSide / radius;
    var latMax = lat + halfSide / radius;
    var lonMin = lon - halfSide / pradius;
    var lonMax = lon + halfSide / pradius;

    return new BoundingBox { 
        MinPoint = new MapPoint { Latitude = Rad2deg(latMin), Longitude = Rad2deg(lonMin) },
        MaxPoint = new MapPoint { Latitude = Rad2deg(latMax), Longitude = Rad2deg(lonMax) }
    };            
}

// degrees to radians
private static double Deg2rad(double degrees)
{
    return Math.PI * degrees / 180.0;
}

// radians to degrees
private static double Rad2deg(double radians)
{
    return 180.0 * radians / Math.PI;
}

// Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
private static double WGS84EarthRadius(double lat)
{
    // http://en.wikipedia.org/wiki/Earth_radius
    var An = WGS84_a * WGS84_a * Math.Cos(lat);
    var Bn = WGS84_b * WGS84_b * Math.Sin(lat);
    var Ad = WGS84_a * Math.Cos(lat);
    var Bd = WGS84_b * Math.Sin(lat);
    return Math.Sqrt((An*An + Bn*Bn) / (Ad*Ad + Bd*Bd));
}
Ε Г И І И О
la source
1
Merci, ce travail pour moi. J'ai dû tester le code à la main, je ne savais pas comment écrire un test unitaire pour cela, mais cela génère des résultats précis avec le degré de précision dont j'ai besoin
mdoran3844
quel est le haldSideinKm ici? ne comprends pas ... ce que je dois passer dans ces argyments, le rayon entre deux points de la carte ou quoi?
@GeloVolro: C'est la demi-longueur du cadre de sélection que vous voulez.
Ε Г И І И О
1
Vous ne devez pas nécessairement écrire votre propre classe MapPoint. Il existe une classe GeoCoordinate dans System.Device.Location qui prend la latitude et la longitude comme paramètres.
Lawyerson
1
Cela fonctionne à merveille. J'apprécie vraiment le port C #.
Tom Larcher
10

J'ai écrit une fonction JavaScript qui renvoie les quatre coordonnées d'une boîte englobante carrée, étant donné une distance et une paire de coordonnées:

'use strict';

/**
 * @param {number} distance - distance (km) from the point represented by centerPoint
 * @param {array} centerPoint - two-dimensional array containing center coords [latitude, longitude]
 * @description
 *   Computes the bounding coordinates of all points on the surface of a sphere
 *   that has a great circle distance to the point represented by the centerPoint
 *   argument that is less or equal to the distance argument.
 *   Technique from: Jan Matuschek <http://JanMatuschek.de/LatitudeLongitudeBoundingCoordinates>
 * @author Alex Salisbury
*/

getBoundingBox = function (centerPoint, distance) {
  var MIN_LAT, MAX_LAT, MIN_LON, MAX_LON, R, radDist, degLat, degLon, radLat, radLon, minLat, maxLat, minLon, maxLon, deltaLon;
  if (distance < 0) {
    return 'Illegal arguments';
  }
  // helper functions (degrees<–>radians)
  Number.prototype.degToRad = function () {
    return this * (Math.PI / 180);
  };
  Number.prototype.radToDeg = function () {
    return (180 * this) / Math.PI;
  };
  // coordinate limits
  MIN_LAT = (-90).degToRad();
  MAX_LAT = (90).degToRad();
  MIN_LON = (-180).degToRad();
  MAX_LON = (180).degToRad();
  // Earth's radius (km)
  R = 6378.1;
  // angular distance in radians on a great circle
  radDist = distance / R;
  // center point coordinates (deg)
  degLat = centerPoint[0];
  degLon = centerPoint[1];
  // center point coordinates (rad)
  radLat = degLat.degToRad();
  radLon = degLon.degToRad();
  // minimum and maximum latitudes for given distance
  minLat = radLat - radDist;
  maxLat = radLat + radDist;
  // minimum and maximum longitudes for given distance
  minLon = void 0;
  maxLon = void 0;
  // define deltaLon to help determine min and max longitudes
  deltaLon = Math.asin(Math.sin(radDist) / Math.cos(radLat));
  if (minLat > MIN_LAT && maxLat < MAX_LAT) {
    minLon = radLon - deltaLon;
    maxLon = radLon + deltaLon;
    if (minLon < MIN_LON) {
      minLon = minLon + 2 * Math.PI;
    }
    if (maxLon > MAX_LON) {
      maxLon = maxLon - 2 * Math.PI;
    }
  }
  // a pole is within the given distance
  else {
    minLat = Math.max(minLat, MIN_LAT);
    maxLat = Math.min(maxLat, MAX_LAT);
    minLon = MIN_LON;
    maxLon = MAX_LON;
  }
  return [
    minLon.radToDeg(),
    minLat.radToDeg(),
    maxLon.radToDeg(),
    maxLat.radToDeg()
  ];
};
Asalisbury
la source
Ce code ne fonctionne pas du tout. Je veux dire, même après avoir corrigé les erreurs évidentes comme minLon = void 0;et maxLon = MAX_LON;cela ne fonctionne toujours pas.
aroth
1
@aroth, je viens de le tester et je n'ai eu aucun problème. N'oubliez pas que l' centerPointargument est un tableau composé de deux coordonnées. Par exemple, getBoundingBox([42.2, 34.5], 50)- void 0est la sortie CoffeeScript pour "non défini" et n'affectera pas la capacité des codes à s'exécuter.
asalisbury
ce code ne fonctionne pas. degLat.degToRadn'est pas une fonction
user299709
Le code fonctionnait dans Node et Chrome tel quel, jusqu'à ce que je le mette dans un projet sur lequel je travaille et que je commence à avoir des degToRaderreurs " n'est pas une fonction". Je Number.prototype.n'ai jamais découvert pourquoi mais ce n'est pas une bonne idée pour une fonction utilitaire comme celle-ci, alors je les ai converties en fonctions locales normales. Il est également important de noter que la boîte retournée est [LNG, LAT, LNG, LAT] au lieu de [LAT, LNG, LAT, LNG]. J'ai modifié la fonction de retour lorsque j'ai utilisé cela pour éviter toute confusion.
KernelDeimos
9

Comme j'avais besoin d'une estimation très approximative, donc pour filtrer certains documents inutiles dans une requête elasticsearch, j'ai utilisé la formule ci-dessous:

Min.lat = Given.Lat - (0.009 x N)
Max.lat = Given.Lat + (0.009 x N)
Min.lon = Given.lon - (0.009 x N)
Max.lon = Given.lon + (0.009 x N)

N = kms requis pour l'emplacement donné. Pour votre cas N = 10

Pas précis mais pratique.

Ajay
la source
En effet, pas précis mais toujours utile et très facile à mettre en œuvre.
MV.
6

Vous recherchez une formule ellipsoïde.

Le meilleur endroit que j'ai trouvé pour commencer à coder est basé sur la bibliothèque Geo :: Ellipsoid du CPAN. Il vous donne une base de référence pour créer vos tests et comparer vos résultats avec ses résultats. Je l'ai utilisé comme base pour une bibliothèque similaire pour PHP chez mon ancien employeur.

Géo :: Ellipsoïde

Jetez un œil à la locationméthode. Appelez-le deux fois et vous avez votre bbox.

Vous n'avez pas publié la langue que vous utilisiez. Une bibliothèque de géocodage est peut-être déjà disponible pour vous.

Oh, et si vous ne l'avez pas encore compris, Google Maps utilise l'ellipsoïde WGS84.

jcoby
la source
5

Illustration de @Jan Philip Matuschek excellente explication. (Veuillez voter pour sa réponse, pas celle-ci; j'ajoute ceci car j'ai pris un peu de temps pour comprendre la réponse originale)

La technique de la boîte englobante consistant à optimiser la recherche des voisins les plus proches aurait besoin de dériver les paires de latitude et de longitude minimale et maximale pour un point P à la distance d. Tous les points qui se trouvent en dehors de ceux-ci sont définitivement à une distance supérieure à d du point. Une chose à noter ici est le calcul de la latitude d'intersection comme le souligne l'explication de Jan Philip Matuschek. La latitude d'intersection n'est pas à la latitude du point P mais légèrement décalée de celui-ci. Il s'agit d'une partie souvent manquée mais importante pour déterminer la longitude de délimitation minimale et maximale correcte du point P pour la distance d. Ceci est également utile pour la vérification.

La distance sinusoïdale entre (latitude d'intersection, longitude haute) et (latitude, longitude) de P est égale à la distance d.

Python gist ici https://gist.github.com/alexcpn/f95ae83a7ee0293a5225

entrez la description de l'image ici

Alex Punnen
la source
5

Voici une implémentation simple utilisant javascript qui est basée sur la conversion du degré de latitude en kms où 1 degree latitude ~ 111.2 km.

Je calcule les limites de la carte à partir d'une latitude et d'une longitude données d'une largeur de 10 km.

function getBoundsFromLatLng(lat, lng){
     var lat_change = 10/111.2;
     var lon_change = Math.abs(Math.cos(lat*(Math.PI/180)));
     var bounds = { 
         lat_min : lat - lat_change,
         lon_min : lng - lon_change,
         lat_max : lat + lat_change,
         lon_max : lng + lon_change
     };
     return bounds;
}
Noushad
la source
4

J'ai adapté un script PHP que j'ai trouvé pour faire exactement cela. Vous pouvez l'utiliser pour trouver les coins d'une boîte autour d'un point (par exemple, à 20 km). Mon exemple spécifique concerne l'API Google Maps:

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

Richard
la source
-1 Ce que recherche l'OP est: étant donné un point de référence (lat, lon) et une distance, trouver la plus petite case telle que tous les points qui sont <= "distance" du point de référence ne soient pas en dehors de la case. Votre boîte a ses coins "éloignés" du point de référence et est donc trop petite. Exemple: le point "distance" plein nord est bien en dehors de votre boîte.
John Machin le
Eh bien, par hasard, c'est exactement ce dont j'avais besoin. Alors merci, même si cela ne répond pas tout à fait à cette question :)
Simon Steinberger
Eh bien, je suis content que cela puisse aider quelqu'un!
Richard
1

Je travaillais sur le problème de la boîte englobante en tant que problème secondaire pour trouver tous les points dans le rayon SrcRad d'un point statique LAT, LONG. Il y a eu pas mal de calculs qui utilisent

maxLon = $lon + rad2deg($rad/$R/cos(deg2rad($lat)));
minLon = $lon - rad2deg($rad/$R/cos(deg2rad($lat)));

pour calculer les limites de longitude, mais j'ai trouvé que cela ne donnait pas toutes les réponses nécessaires. Parce que ce que tu veux vraiment faire c'est

(SrcRad/RadEarth)/cos(deg2rad(lat))

Je sais, je sais que la réponse devrait être la même, mais j'ai trouvé que ce n'était pas le cas. Il est apparu qu'en ne m'assurant pas de faire le (SRCrad / RadEarth) d'abord, puis en divisant par la partie Cos, j'omettais certains points de localisation.

Une fois que vous avez obtenu tous les points de votre boîte englobante, si vous avez une fonction qui calcule la distance point à point donnée lat, il est facile d'obtenir uniquement les points qui sont à un certain rayon de distance du point fixe. Voici ce que j'ai fait. Je sais que cela a pris quelques étapes supplémentaires mais cela m'a aidé

-- GLOBAL Constants
gc_pi CONSTANT REAL := 3.14159265359;  -- Pi

-- Conversion Factor Constants
gc_rad_to_degs          CONSTANT NUMBER := 180/gc_pi; -- Conversion for Radians to Degrees 180/pi
gc_deg_to_rads          CONSTANT NUMBER := gc_pi/180; --Conversion of Degrees to Radians

lv_stat_lat    -- The static latitude point that I am searching from 
lv_stat_long   -- The static longitude point that I am searching from 

-- Angular radius ratio in radians
lv_ang_radius := lv_search_radius / lv_earth_radius;
lv_bb_maxlat := lv_stat_lat + (gc_rad_to_deg * lv_ang_radius);
lv_bb_minlat := lv_stat_lat - (gc_rad_to_deg * lv_ang_radius);

--Here's the tricky part, accounting for the Longitude getting smaller as we move up the latitiude scale
-- I seperated the parts of the equation to make it easier to debug and understand
-- I may not be a smart man but I know what the right answer is... :-)

lv_int_calc := gc_deg_to_rads * lv_stat_lat;
lv_int_calc := COS(lv_int_calc);
lv_int_calc := lv_ang_radius/lv_int_calc;
lv_int_calc := gc_rad_to_degs*lv_int_calc;

lv_bb_maxlong := lv_stat_long + lv_int_calc;
lv_bb_minlong := lv_stat_long - lv_int_calc;

-- Now select the values from your location datatable 
SELECT *  FROM (
SELECT cityaliasname, city, state, zipcode, latitude, longitude, 
-- The actual distance in miles
spherecos_pnttopntdist(lv_stat_lat, lv_stat_long, latitude, longitude, 'M') as miles_dist    
FROM Location_Table 
WHERE latitude between lv_bb_minlat AND lv_bb_maxlat
AND   longitude between lv_bb_minlong and lv_bb_maxlong)
WHERE miles_dist <= lv_limit_distance_miles
order by miles_dist
;
Greg Beck
la source
0

C'est très simple, il suffit d'aller sur le site Web de panoramio, puis d'ouvrir la carte du monde à partir du site Web de panoramio.Ensuite, allez à l'emplacement spécifié dont la latitude et la longitude sont requises.

Ensuite, vous avez trouvé la latitude et la longitude dans la barre d'adresse par exemple dans cette adresse.

http://www.panoramio.com/map#lt=32.739485&ln=70.491211&z=9&k=1&a=1&tab=1&pl=all

lt = 32,739485 => latitude ln = 70,491211 => longitude

ce widget API JavaScript Panoramio crée un cadre de délimitation autour d'une paire lat / longue, puis renvoie toutes les photos avec dans ces limites.

Un autre type de widget API JavaScript Panoramio dans lequel vous pouvez également changer la couleur d'arrière-plan avec l' exemple et le code est ici .

Cela n'apparaît pas dans l'ambiance de composition, mais après la publication.

<div dir="ltr" style="text-align: center;" trbidi="on">
<script src="https://ssl.panoramio.com/wapi/wapi.js?v=1&amp;hl=en"></script>
<div id="wapiblock" style="float: right; margin: 10px 15px"></div>
<script type="text/javascript">
var myRequest = {
  'tag': 'kahna',
  'rect': {'sw': {'lat': -30, 'lng': 10.5}, 'ne': {'lat': 50.5, 'lng': 30}}
};
  var myOptions = {
  'width': 300,
  'height': 200
};
var wapiblock = document.getElementById('wapiblock');
var photo_widget = new panoramio.PhotoWidget('wapiblock', myRequest, myOptions);
photo_widget.setPosition(0);
</script>
</div>
kahna9
la source
0

Ici, j'ai converti la réponse de Federico A. Ramponi en PHP si quelqu'un est intéressé:

<?php
# deg2rad and rad2deg are already within PHP

# Semi-axes of WGS-84 geoidal reference
$WGS84_a = 6378137.0;  # Major semiaxis [m]
$WGS84_b = 6356752.3;  # Minor semiaxis [m]

# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
function WGS84EarthRadius($lat)
{
    global $WGS84_a, $WGS84_b;

    $an = $WGS84_a * $WGS84_a * cos($lat);
    $bn = $WGS84_b * $WGS84_b * sin($lat);
    $ad = $WGS84_a * cos($lat);
    $bd = $WGS84_b * sin($lat);

    return sqrt(($an*$an + $bn*$bn)/($ad*$ad + $bd*$bd));
}

# Bounding box surrounding the point at given coordinates,
# assuming local approximation of Earth surface as a sphere
# of radius given by WGS84
function boundingBox($latitudeInDegrees, $longitudeInDegrees, $halfSideInKm)
{
    $lat = deg2rad($latitudeInDegrees);
    $lon = deg2rad($longitudeInDegrees);
    $halfSide = 1000 * $halfSideInKm;

    # Radius of Earth at given latitude
    $radius = WGS84EarthRadius($lat);
    # Radius of the parallel at given latitude
    $pradius = $radius*cos($lat);

    $latMin = $lat - $halfSide / $radius;
    $latMax = $lat + $halfSide / $radius;
    $lonMin = $lon - $halfSide / $pradius;
    $lonMax = $lon + $halfSide / $pradius;

    return array(rad2deg($latMin), rad2deg($lonMin), rad2deg($latMax), rad2deg($lonMax));
}
?>
Joe Black
la source
0

Merci @Fedrico A. pour l'implémentation de Phyton, je l'ai porté dans une classe de catégorie Objective C. Voici:

#import "LocationService+Bounds.h"

//Semi-axes of WGS-84 geoidal reference
const double WGS84_a = 6378137.0; //Major semiaxis [m]
const double WGS84_b = 6356752.3; //Minor semiaxis [m]

@implementation LocationService (Bounds)

struct BoundsLocation {
    double maxLatitude;
    double minLatitude;
    double maxLongitude;
    double minLongitude;
};

+ (struct BoundsLocation)locationBoundsWithLatitude:(double)aLatitude longitude:(double)aLongitude maxDistanceKm:(NSInteger)aMaxKmDistance {
    return [self boundingBoxWithLatitude:aLatitude longitude:aLongitude halfDistanceKm:aMaxKmDistance/2];
}

#pragma mark - Algorithm 

+ (struct BoundsLocation)boundingBoxWithLatitude:(double)aLatitude longitude:(double)aLongitude halfDistanceKm:(double)aDistanceKm {
    double radianLatitude = [self degreesToRadians:aLatitude];
    double radianLongitude = [self degreesToRadians:aLongitude];
    double halfDistanceMeters = aDistanceKm*1000;


    double earthRadius = [self earthRadiusAtLatitude:radianLatitude];
    double parallelRadius = earthRadius*cosl(radianLatitude);

    double radianMinLatitude = radianLatitude - halfDistanceMeters/earthRadius;
    double radianMaxLatitude = radianLatitude + halfDistanceMeters/earthRadius;
    double radianMinLongitude = radianLongitude - halfDistanceMeters/parallelRadius;
    double radianMaxLongitude = radianLongitude + halfDistanceMeters/parallelRadius;

    struct BoundsLocation bounds;
    bounds.minLatitude = [self radiansToDegrees:radianMinLatitude];
    bounds.maxLatitude = [self radiansToDegrees:radianMaxLatitude];
    bounds.minLongitude = [self radiansToDegrees:radianMinLongitude];
    bounds.maxLongitude = [self radiansToDegrees:radianMaxLongitude];

    return bounds;
}

+ (double)earthRadiusAtLatitude:(double)aRadianLatitude {
    double An = WGS84_a * WGS84_a * cosl(aRadianLatitude);
    double Bn = WGS84_b * WGS84_b * sinl(aRadianLatitude);
    double Ad = WGS84_a * cosl(aRadianLatitude);
    double Bd = WGS84_b * sinl(aRadianLatitude);
    return sqrtl( ((An * An) + (Bn * Bn))/((Ad * Ad) + (Bd * Bd)) );
}

+ (double)degreesToRadians:(double)aDegrees {
    return M_PI*aDegrees/180.0;
}

+ (double)radiansToDegrees:(double)aRadians {
    return 180.0*aRadians/M_PI;
}



@end

Je l'ai testé et semble fonctionner correctement. Struct BoundsLocation devrait être remplacé par une classe, je l'ai utilisé juste pour le partager ici.

Jésuslg123
la source
0

Toutes les réponses ci-dessus ne sont que partiellement correctes . Surtout dans une région comme l'Australie, ils incluent toujours le pôle et calculent un très grand rectangle même pour 10 km.

En particulier, l'algorithme de Jan Philip Matuschek à http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates#UsingIndex incluait un très grand rectangle de (-37, -90, -180, 180) pour presque tous les points en Australie. Cela touche un grand nombre d'utilisateurs dans la base de données et la distance doit être calculée pour tous les utilisateurs dans près de la moitié du pays.

J'ai trouvé que l' algorithme Drupal API Earth de Rochester Institute of Technology fonctionne mieux autour du pôle ainsi qu'ailleurs et est beaucoup plus facile à mettre en œuvre.

https://www.rit.edu/drupal/api/drupal/sites%21all%21modules%21location%21earth.inc/7.54

Utilisez earth_latitude_rangeet à earth_longitude_rangepartir de l'algorithme ci-dessus pour calculer le rectangle englobant

Et utilisez la formule de calcul de distance documentée par google maps pour calculer la distance

https://developers.google.com/maps/solutions/store-locator/clothing-store-locator#outputting-data-as-xml-using-php

Pour effectuer une recherche par kilomètres au lieu de miles, remplacez 3959 par 6371. Pour (Lat, Lng) = (37, -122) et une table de marqueurs avec les colonnes lat et lng , la formule est:

SELECT id, ( 3959 * acos( cos( radians(37) ) * cos( radians( lat ) ) * cos( radians( lng ) - radians(-122) ) + sin( radians(37) ) * sin( radians( lat ) ) ) ) AS distance FROM markers HAVING distance < 25 ORDER BY distance LIMIT 0 , 20;

Lisez ma réponse détaillée sur https://stackoverflow.com/a/45950426/5076414

Sacky San
la source
0

Voici la réponse de Federico Ramponi dans Go. Remarque: pas de vérification des erreurs :(

import (
    "math"
)

// Semi-axes of WGS-84 geoidal reference
const (
    // Major semiaxis (meters)
    WGS84A = 6378137.0
    // Minor semiaxis (meters)
    WGS84B = 6356752.3
)

// BoundingBox represents the geo-polygon that encompasses the given point and radius
type BoundingBox struct {
    LatMin float64
    LatMax float64
    LonMin float64
    LonMax float64
}

// Convert a degree value to radians
func deg2Rad(deg float64) float64 {
    return math.Pi * deg / 180.0
}

// Convert a radian value to degrees
func rad2Deg(rad float64) float64 {
    return 180.0 * rad / math.Pi
}

// Get the Earth's radius in meters at a given latitude based on the WGS84 ellipsoid
func getWgs84EarthRadius(lat float64) float64 {
    an := WGS84A * WGS84A * math.Cos(lat)
    bn := WGS84B * WGS84B * math.Sin(lat)

    ad := WGS84A * math.Cos(lat)
    bd := WGS84B * math.Sin(lat)

    return math.Sqrt((an*an + bn*bn) / (ad*ad + bd*bd))
}

// GetBoundingBox returns a BoundingBox encompassing the given lat/long point and radius
func GetBoundingBox(latDeg float64, longDeg float64, radiusKm float64) BoundingBox {
    lat := deg2Rad(latDeg)
    lon := deg2Rad(longDeg)
    halfSide := 1000 * radiusKm

    // Radius of Earth at given latitude
    radius := getWgs84EarthRadius(lat)

    pradius := radius * math.Cos(lat)

    latMin := lat - halfSide/radius
    latMax := lat + halfSide/radius
    lonMin := lon - halfSide/pradius
    lonMax := lon + halfSide/pradius

    return BoundingBox{
        LatMin: rad2Deg(latMin),
        LatMax: rad2Deg(latMax),
        LonMin: rad2Deg(lonMin),
        LonMax: rad2Deg(lonMax),
    }
}
sma
la source