Vincent GODARD - V2 - 07/12/2022
Inspiré de :
de
et de
Document Notebook à télécharger => ici.
Dans le Jupyter Notebook ou le JupyterLab, la bibliothèque Pandas est déjà installée mais peut-être pas la bibliothèque folium ni celle qui va permettre le géocodage sur les données OSM, geopy. Il faut les installer.
Il faut ensuite les charger (fonction import).
# Pour installer folium et geopy, effacez le # (avant le "!").
#!pip install folium
#!pip install geopy
import geopy
from geopy.distance import geodesic
import pandas as pd
import folium
Nous pourrions programmer un calcul de distance, mais cette fonctionalité existe déjà dans de nombreuses bibliothèques.
Dans la bibliothèqye GeoPy, nous utiliserons la fonction "geodesic" pour calculer la distance entre deux points à la surface de la Terre.
La consultation de l'aide en ligne de la bibliothèque Geopy (https://pypi.org/project/geopy/), nous apprend qu'elle permet deux calculs : l'un fournit la distance géodésique, encore appelée distance géographique (https://en.wikipedia.org/wiki/Geographical_distance), le deuxième, la distance sphérique ou orthodromique, celle des grands cercles (https://en.wikipedia.org/wiki/Great-circle_distance).
La distance géodésique est la plus petite distance à la surface d'un ellipsoïde terrestre. C'est celle de l'algorithme par défaut.
Nous allons commencer par la distance géodésique. Celle qui est proposée par défaut dans Geopy.
# Calcul de la distance entre deux points
# Le couple (lat, long) est un "tuple", une collection ordonnée de plusieurs éléments (on ne permutent pas les x, y et z !)
point1 = (49.916081731637746, -6.317862543498883)
point2 = (49.95567132388271, -6.339694772850285)
distance = geodesic(point1, point2)
La distance est calculée, mais rien ne s'affiche. Il faut demander l'impression du résultat !
# Affichage de la distance en mètres.
print(distance.meters, ' meters')
Il y a peut-être trop de décimales ?
# Affichage en kilomètres avec un choix de décimales égal à 2 plus un texte de d'accompagnement
print(f"Résultats : La distance géodésique entre le point 1 et le point 2 est de {distance.kilometers:.2f}", 'kilomètres')
D'autres unités sont possibles :
# En miles impériaux (ceux des anglais)
print(f"Résultats : La distance géodésique entre le point 1 et le point 2 est de {distance.miles:.2f}", 'en imperial miles')
# En milles marins (mille nautique ou nautique)
print(f"Résultats : La distance géodésique entre le point 1 et le point 2 est de {distance.nautical:.2f}", 'nautiques')
Mais où sommes nous ?
Une carte serait peut-être utiles pour faire notre choix d'unité ?
# Indiquer les coordonnées lat, long en DD du point 1.
lat1 = 49.916081731637746
long1 = -6.317862543498883
lat2 = 49.95567132388271
long2 = -6.339694772850285
Map = folium.Map(location = [lat1, long1], zoom_start = 12)
Map
# Indiquer les coordonnées lat, long en DD des points 1 et 2 avec un marqueur.
folium.Marker([lat1, long1], popup = "Point mystère n°1").add_to(Map)
folium.Marker([lat2, long2], popup = "Point mystère n°2").add_to(Map)
Map
Where is it ?
Utilise un modèle sphérique de la Terre, où le rayon moyen de la Terre est d'environ 6371.009 km (pour du WGS84). Se référer à https://geopy.readthedocs.io/en/stable/#module-geopy.distance pour les sources algorithmiques.
Correspond à la distance orthodromique des navigateurs.
# Calcul de la distance sur le grand cercle entre deux points
from geopy.distance import great_circle
point1 = (49.916081731637746, -6.317862543498883)
point2 = (49.95567132388271, -6.339694772850285)
dist_great_circle = great_circle(point1, point2).kilometers
print(dist_great_circle)
Est-ce différent de la distance géodésique ?
# Comparaison de deux distances calculées entre deux points
# Distance sur le grand cercle
print('Distance sur le grand cercle =', dist_great_circle, 'km')
# Distance géodésique
print('Distance géodésique =', distance.kilometers, 'km')
Pour améliorer la lisibilité au niveau décimale :
# Comparaison de deux distances calculées entre deux points
# Distance sur le grand cercle
print(dist_great_circle, 'km (Distance sur le grand cercle)')
# Distance géodésique
print(distance.kilometers, 'km (Distance géodésique)')
Cette différence est-elle significative au regard de la façon de la calculer (comment avons-nous positionné les pointeurs par exemple ?) ?
C'est sans doute sur les grandes distances que cette différence gagne en importance !
# Calcul de la distance entre Paris (Notre-Dame) et New York (Statue de la Liberté)
from geopy.distance import great_circle
Paris = (48.85395617692726, 2.3497107002361703)
New_York = (40.68958997693833, -74.04456826872635)
# Distance sur le grand cercle
dist_great_circle = great_circle(Paris, New_York).kilometers
print(dist_great_circle, 'km (Distance sur le grand cercle)')
# Distance géodésique
distance = geodesic(Paris, New_York)
print(distance.kilometers, 'km (Distance géodésique)')
La distance sur le grand cercle est nettement inférieure à la distance géodésique. C'est plus net que dans l'exemple précédent sur les Sorlingues (nom français des Îles Scilly) !
Pour adapter les calculs de distances à différents systèmes géodésiques (par exemple aux systèmes légaux en vigueur), il est possible d'utiliser différents ellipsoïdes disponibles dans Geopy ou d'en paramétrer à la demande (cf. https://geopy.readthedocs.io/en/stable/#module-geopy.distance).
En premier lieu, il faut passer en revue le dictionnaire d'ellipsoïdes de Geopy.
Un dictionnaire dans Python est une liste d'éléments auxquels on accède par une clef et non plus par un indice comme dans un dans les "listes" ou les "tuples" vus précédemment. Dans Geopy, ils sont fixes, mais on peut quand même en utiliser d'autres.
On les identifie par : {clé: valeur}.
Pour compléter cette notion : https://courspython.com/dictionnaire.html# Dictionnaire fixe indiquant quelques ellipsoïdes {'nom de la clef': (demi grand axe, demi petit axe, aplatissement)}
ELLIPSOIDS = {'WGS-84': (6378.137, 6356.7523142, 1 / 298.257223563),
'GRS-80': (6378.137, 6356.7523141, 1 / 298.257222101),
'Airy (1830)': (6377.563396, 6356.256909, 1 / 299.3249646),
'Intl 1924': (6378.388, 6356.911946, 1 / 297.0),
'Clarke (1880)': (6378.249145, 6356.51486955, 1 / 293.465),
'GRS-67': (6378.1600, 6356.774719, 1 / 298.25),
}
Pour information ou rappel :
• Le système géodésique NTF (ancien système légal français jusqu'au tournant des années 2000) est associé à l'ellipsoïde Clarke 1880. C'est un système local et pas mondial issu de géodésie spatiale.
• Le système géodésique ED50 (ancien système des alliés pendant la Guerre froide) est associé à l'ellipsoïde Hayford 1909 (ici Intl 1924).
• Le système géodésique RGF93 (système légal français actuel) est associé à l'ellipsoïde IAG-GRS80.
• Le système géodésique WGS-84 (système utilisé par les GPS) est très proche de l'ellipsoïde IAG-GRS80, seul le demi petit axe varie.
Les deux derniers sont des systèmes spatiaux (à vocation mondiale) alors que les premiers sont qualifiés de locaux (https://fr.wikipedia.org/wiki/Syst%C3%A8me_g%C3%A9od%C3%A9sique).
○ Pour mémoire : distance sur le grand cercle vs distance géodésique
# Calcul de la distance sur le grand cercle entre Brest (Le château) et Strasbourg (La cathédrale)
from geopy.distance import great_circle
Brest = (48.381236144334025, -4.494679502022013)
Strasbourg = (48.58197325516147, 7.751035828319029)
# Distance sur le grand cercle
dist_great_circle = great_circle(Brest, Strasbourg).kilometers
print(dist_great_circle, 'km (Distance sur le grand cercle)')
# Distance géodésique
distance = geodesic(Brest, Strasbourg)
print(distance.kilometers, 'km (Distance géodésique)')
Par défaut "geopy.distance.distance" utilise la distance géodésique.
Celle-ci se calcule sur des tuples de (lat, long).
# distance.distance pour geodesic
from geopy import distance
Brest = (48.381236144334025, -4.494679502022013)
Strasbourg = (48.58197325516147, 7.751035828319029)
print(distance.distance(Brest, Strasbourg).km)
L'ellipsoïde par défaut est le WGS-84.
# Ellipsoïde WGS-84
print(distance.geodesic(Brest, Strasbourg, ellipsoid='WGS-84').km)
# Ellipsoïde GRS-80 du système géodésique RGF93 (en fait IAG-GRS80 du système européen ETRS89)
print(distance.geodesic(Brest, Strasbourg, ellipsoid='GRS-80').km)
À quelle décimale apparaissent les différences ?
# Reprise des deux résultats précédents
print(f"Résultats : La distance géodésique entre Brest et Strasbourg est de {distance.geodesic(Brest, Strasbourg, ellipsoid='WGS-84').km:.8f}", 'kilomètres selon le WGS-84')
print(f"Résultats : La distance géodésique entre Brest et Strasbourg est de {distance.geodesic(Brest, Strasbourg, ellipsoid='GRS-80').km:.8f}", 'kilomètres selon le GRS-80')
Soit, combien de micromètres ?
Et avec l'ancien système géodésique français, la nouvelle triangulation française (NTF) ?
# Avec la NTF
print(distance.geodesic(Brest, Strasbourg, ellipsoid='Clarke (1880)').km)
On est à une cinquantaine de mètres d'écart avec le RGF93.
○ En s'inspirant de la page wikipédia Earth ellipsoid : https://en.wikipedia.org/wiki/Earth_ellipsoid
On peut calculer la même distance en paramétrant des ellipsoïdes non contenus dans le disctionnaire.
Comme celui de Maupertuis (1738).# Maupertuis (1738)
distance.geodesic(Brest, Strasbourg, ellipsoid=(6397.300, 6363.806283, 1 / 191.)).km
On en conclut que les distances étaient plus courtes au XVIIIème siècle (: !
Saurez-vous faire les calculs en reprenant les scripts de la section 3.2 ?
À venir !