Comment rechercher les POI autour d’un point GPS avec une requête SQL

Jean-Baptiste Guerraz

16 février 2009 | Catégorie: MySQL |

En travaillant sur une base essentiellement composée de POI (Point Of Interest) je me suis demandé comment j’allais pouvoir rechercher dans celle-ci tous les points se trouvant à moins de X kms d’une coordonnée GPS donnée. L’idée par exemple, est de rechercher tous les restaurants se trouvant autour d’un musée ou à coté de votre hôtel.

Ma table en question est assez simple et contient les coordonnées GPS sortis de Google Maps de chaque point. Pour simplifier, je donnerai la définition suivante de la structure de ma table :

poi_id : int(10), PK, identifiant unique
latitude : float(12,8), GMAP latitude ou phi
longitude : float(12,8), GMAP longitude ou lambda

Les coordonnées GMAP stockés sont des angles exprimés en degrés, or il me faut calculer la distance entre tous les points renseignés en ma coordonnée de référence.

Rien de fort compliqué, un peu de trigo et le tour est joué :

R = Rayon de la terre ( = 6,371km) Δlat = lat2− lat1 Δlong = long2− long1 a = sin²(Δlat/2) + cos(lat1)*cos(lat2)*sin²(Δlong/2) c = 2*atan2(√a, √(1−a)) d = R*c

La distance entre deux points dont nous connaissons la latitude et la longitude est donc donné par cette formule. Maintenant reste plus qu’à la transformer en requête SQL (attention, par rapport à la formule précédente, j’ai rajouté le passage en radians) :

SET @lat = $lat;
SET @long = $long;
SET @radius = $radius;
 
SELECT 
  poi_id 
FROM poi AS p 
WHERE 6371 * 2 * ATAN2 ( SQRT ( ( SIN( ( RADIANS(p.latitude - @lat) / 2 ) * SIN( RADIANS(p.latitude - @lat) / 2 ) + COS ( RADIANS (@lat )) * COS ( RADIANS ( p.latitude ) ) * SIN ( RADIANS(p.longitude - @long) / 2 ) * SIN ( RADIANS(p.longitude - @long) / 2 ) ) ) ,  SQRT ( 1 - (SIN( RADIANS(p.latitude - @lat) / 2 ) * SIN( RADIANS(p.latitude - @lat) / 2 ) + COS ( RADIANS (@lat) ) * COS ( RADIANS (p.latitude) ) * SIN ( RADIANS(p.longitude - @long) / 2 ) * SIN ( RADIANS(p.longitude - @long) / 2 ) ) ) ) < @radius;

Il faut remplacer $lat, $lon par les coordonnées de notre point de référence et $radius par le rayon dans lequel on cherche.

Pour finir, si vous avez beaucoup de points (ce qui est souvent le cas), je vous conseille de partitionner la table par RANGE en découpant par tranches de longitudes par exemple (ou même en carrés). MySQL branchera le “prunning” car la clause WHERE, même si elle semble compliquée sera bien évaluée par MySQL et il saura donc dans quelle partition(s) chercher les points.

Très souvent, les coordonnées disponibles en France sont au format Lambert II étendu. Or Google Maps utilise le WGS84, voici donc un petit script PHP permettant le passage de Lambert vers GMAP (WGS84) :

<?php
function ALG0001($phi,$e)
{
    $temp = ( 1 - ( $e * sin( $phi ) ) ) / ( 1 + ( $e * sin( $phi ) ) );
 
    $L = log ( tan ( (pi()/4) + ($phi/2) ) * pow ($temp, ($e/2) ));
 
    return $L;
}
 
function ALG0002($L,$e,$epsilon)
{
    $phi[0] = 2 * atan(exp($L)) - (pi()/2);
 
    $i=0;
    do
    {
        $i++;
        $temp = ( 1 + ( $e * sin( $phi[$i-1] ) ) ) / ( 1 - ( $e * sin( $phi[$i-1] ) ) );
        $phi[$i] = 2 * atan ( pow ($temp, ($e/2)) * exp ($L) ) - pi()/2;
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);
 
    return $phi[$i];
}
 
function ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon)
{
    $R = sqrt( pow(($X - $Xs),2) + pow(($Y - $Ys),2) );
    $gamma = atan(($X - $Xs)/($Ys - $Y));
 
    $lambda = $lambdac + ($gamma / $n);
 
    $L = (-1 / $n) * log(abs($R/$c));
 
    $phi = ALG0002($L,$e,$epsilon);
 
    $coords['lambda']=$lambda;
    $coords['phi']=$phi;
 
    return $coords;
}
 
function ALG0009($lambda,$phi,$he,$a,$e)
{
    $N = ALG0021($phi,$a,$e);
 
    $X = ($N + $he) * cos($phi) * cos($lambda);
 
    $Y = ($N + $he) * cos($phi) * sin($lambda);
 
    $Z = ($N * (1 - $e*$e)  + $he) * sin ($phi);
 
    $coords['X']=$X;
    $coords['Y']=$Y;
    $coords['Z']=$Z;
 
    return $coords;
}
 
function ALG0012($X,$Y,$Z,$a,$e,$epsilon)
{
    $lambda = atan ($Y/$X);
 
    $P = sqrt($X*$X + $Y*$Y);
    $phi[0] = atan ($Z/ ($P * (1 - ( ($a*$e*$e)/sqrt($X*$X + $Y*$Y + $Z*$Z) ) ) ) );
 
    $i = 0;
    do
    {
        $i++;
        $temp =  pow((1 - ( $a * $e*$e * cos($phi[$i - 1] )/( $P * sqrt(1 - $e*$e * sin($phi[$i - 1])*sin($phi[$i - 1]))))),-1);
        $phi[$i] = atan( $temp * $Z / $P );
    }
    while (abs($phi[$i] - $phi[$i - 1]) >= $epsilon);
 
    $phix = $phi[$i];
 
    $he = ($P/cos($phix)) - ($a/sqrt(1 - $e*$e * sin($phix)*sin($phix)));
 
    $coords['lambda']=$lambda;
    $coords['phi']=$phix;
    $coords['he']=$he;
 
    return $coords;
}
 
function ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$U)
{
    $V['X'] = $Tx + $U['X'] * (1 + $D) + $U['Z'] * $Ry - $U['Y'] * $Rz;
    $V['Y'] = $Ty + $U['Y'] * (1 + $D) + $U['X'] * $Rz - $U['Z'] * $Rx;
    $V['Z'] = $Tz + $U['Z'] * (1 + $D) + $U['Y'] * $Rx - $U['X'] * $Ry;
 
    return $V;
}
 
function ALG0019($lambda0,$phi0,$k0,$X0,$Y0,$a,$e)
{
    $lambdac = $lambda0;
    $n = sin($phi0);
    $C = $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) * exp ( $n * ALG0001($phi0,$e) );
    $Xs = $X0;
    $Ys = $Y0 + $k0 * ALG0021($phi0,$a,$e) * tan (pi()/2 - $phi0) ;
 
    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;
 
    return $tab;
 
}
 
function ALG0021($phi,$a,$e)
{
    $N = $a/sqrt( 1 - $e * $e * sin($phi) * sin($phi) );
    return $N;
 
}
 
function ALG0054($lambda0,$phi0,$X0,$Y0,$phi1,$phi2,$a,$e)
{
    $lambdac = $lambda0;
    $n = ( (log( (ALG0021($phi2,$a,$e)*cos($phi2))/(ALG0021($phi1,$a,$e)*cos($phi1)) ))/(ALG0001($phi1,$e) - ALG0001($phi2,$e) ));
    $C = ((ALG0021($phi1,$a,$e)* cos($phi1))/$n) * exp($n * ALG0001($phi1,$e));
 
    if ($phi0 == (pi()/2))
    {
        $Xs = $X0;
        $Ys = $Y0;
    }
    else
    {
            echo ('coucou');
        $Xs = $X0;
        $Ys = $Y0 + $C * exp(-1 * $n * ALG0001($phi0,$e));
    }
 
    $tab ['e'] = $e;
    $tab ['n'] = $n;
    $tab ['C'] = $C;
    $tab ['lambdac'] = $lambdac;
    $tab ['Xs'] = $Xs;
    $tab ['Ys'] = $Ys;
 
    return $tab;
 
}
 
function Lambert2WGS84($orig,$X,$Y)
{
    $epsilon = 0.00000000001;
 
    switch ($orig)
    {
        case 'LII' :
            $n = 0.7289686274;
            $c = 11745793.39;
            $Xs = 600000;
            $Ys = 6199695.768;
            $lambdac = 0.04079234433; // pour greenwich
 
            $e = 0.08248325676; //(première excentricité de l’ellipsoïde Clarke 1880 français)
 
            $he = 100;
            $a = 6378249.2;
 
            $Tx = -168;
            $Ty = -60;
            $Tz = +320;
            $D  = 0;
            $Rx = $Ry = $Rz = 0;
        break;
        case 'LIIe' :
            $n = 0.7289686274;
            $c = 11745793.39;
            $Xs = 600000;
            $Ys = 8199695.768;
            $lambdac = 0.04079234433; // for greenwich
 
            $e = 0.08248325676; 
 
            $he = 100;
            $a = 6378249.2;
 
            $Tx = -168;
            $Ty = -60;
            $Tz = +320;
            $D  = 0;
            $Rx = $Ry = $Rz = 0;
        break;
        case 'L93' :
            $n = 0.7256077650;
            $c = 11745255.426;
            $Xs = 700000;
            $Ys = 12655612.050;
            $lambdac = 0.04079234433; // for greenwich
 
            $e = 0.08248325676; 
 
            $he = 100;
            $a = 6378249.2;
 
            $Tx = -168;
            $Ty = -60;
            $Tz = +320;
            $D  = 0;
            $Rx = $Ry = $Rz = 0;
        break;
    }
 
    $coords = ALG0004($X,$Y,$n,$c,$Xs,$Ys,$lambdac,$e,$epsilon);
 
    $coords = ALG0009($coords['lambda'],$coords['phi'],$he,$a,$e);
 
    $coords = ALG0013($Tx,$Ty,$Tz,$D,$Rx,$Ry,$Rz,$coords);
 
    $a = 6378137.0;    // ellipsoïdes WGS84
    $f = 1/298.257223563;
    $b = $a*(1-$f);
    $e = sqrt(($a*$a - $b*$b)/($a*$a));
 
    $X = $coords['X'];
    $Y = $coords['Y'];
    $Z = $coords['Z'];
    $coords = ALG0012($X,$Y,$Z,$a,$e,$epsilon);
 
    $coords['lambda']=rad2deg($coords['lambda']);
    $coords['phi']   =rad2deg($coords['phi']);
    return $coords;
}
 
print_r(Lambert2WGS84('LIIe',591647.56,2426659.65));
?>

5 commentaires sur “Comment rechercher les POI autour d’un point GPS avec une requête SQL”

Plop dit :

29 mars 2009 à 11 h 19 min

La conversion lambert II étendu / wgs84 ne semble pas fonctionner.

J’ai pour ma part des résultats très différents en utilisant ce script php et Circé (ou tout autre converstisseur).

Et pas moyen de trouver d’où provient le problème. Pb de référentiel peut etre auquel s’ajoute un second il me semble…

Répondre
Chandon dit :

5 mai 2009 à 15 h 05 min

Bonjour,
très intéressé par ce script, pour trouver les membres géolocalisés de mon site autour de soi, je me heurte lors de son application à une erreur de syntaxe sql que je ne trouve pas.
Voilà ma requete :

SELECT id_geo
FROM mabase.matable AS p
WHERE 6371 * 2 * ATAN2 ( SQRT ( ( SIN( ( RADIANS(p.lat - 43.27292469899955000) / 2 ) * SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) + COS ( RADIANS (43.27292469899955000 )) * COS ( RADIANS ( p.lat ) ) * SIN ( RADIANS(p.lon - 5.36235809326171900) / 2 ) * SIN ( RADIANS(p.lon - 5.36235809326171900) / 2 ) ) ) , SQRT ( 1 - (SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) * SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) + COS ( RADIANS (43.27292469899955000) ) * COS ( RADIANS (p.lat) ) * SIN ( RADIANS(p.lon - 5.36235809326171900) / 2 ) * SIN ( RADIANS(p.lon - 5.36235809326171900) / 2 ) ) ) ) < 1

J’ai bien vérifié les parenthèses, tout est bon (lat et lon sont mes champs de la table ‘matable’
cela pourrait il venir de ma version de php (5.2.4) ou mysql (5.0.24) ?
Merci de votre aide

Répondre
Maxime Topolov dit :

6 mai 2009 à 4 h 11 min

Bonjour Chandon !

En fait votre requête présente plusieurs erreurs….

1. “…WHERE 6371 * 2 * ATAN2 ( SQRT ( ( SIN( (…” >> Après le SIN il y a une parenthèse ouvrante en trop…

2. J’ai remarqué par ailleurs que MySQL n’accepte pas les espaces après les fonctions… mea culpa, notre exemple on a inséré des espaces après le copier-coller pour plus de clarté.

Voici la requête qui devrait marcher

SELECT id_geo
FROM mabase.matable AS p
WHERE 6371 * 2 *
ATAN2(
SQRT((SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) * SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) + COS( RADIANS(43.27292469899955000 )) * COS( RADIANS( p.lat ) ) * SIN( RADIANS(p.lon - 5.36235809326171900) / 2 ) * SIN( RADIANS(p.lon - 5.36235809326171900) / 2 ) ) ) ,

SQRT( 1 - (SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) * SIN( RADIANS(p.lat - 43.27292469899955000) / 2 ) + COS( RADIANS(43.27292469899955000) ) * COS( RADIANS(p.lat) ) * SIN(RADIANS(p.lon - 5.36235809326171900) / 2 ) * SIN( RADIANS(p.lon - 5.36235809326171900) / 2 ) ) ) ) < 1

Répondre
Chandon dit :

6 mai 2009 à 7 h 30 min

Ahhhh merci !!!!!
Entre temps j’avais trouvé la parenthèse de trop, mais je me heurait justement à ‘ATAN2 fonction inconnue’.
ça marche, merci infiment, vous allez faire le bonheur de mon intranet pro et de mon site web perso

Petite précision, la valeur $radius est en km, c’est bien ça ?

Répondre
Kapweb dit :

16 décembre 2009 à 16 h 59 min

Bonjour,

merci pour ce script fort utile.
Pour répondre à Chandon, oui la valeur de Radius est en Km.

Bonne journée.

Répondre

Laisser une réponse

« Lenteur de MySQL, Drupal, Views & User Relationship [English] Drupal Tutorial: Organizing content on pages »

Méthodes agiles, Drupal, Magento, J2ME, Ergonomie et Utilisabilité

Mots-clefs

Pages

Catégories

Suivez nous sur Twitter