Vincent GODARD - V3 - 14/02/2023
Inspiré de :
de
et de
Document Notebook et données (data) à 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 Geopandas ni celle qui va permettre le calcul des distances, Geopy. Il faut les installer si ce n'est fait (fonction Install).
Il faut ensuite les charger (fonction Import). Pour une meilleure mise en relation des bibliothèques et des traitements, les imports apparaissent là où ils sont utiles dans les scripts.
# Pour installer geopy et geopandas, effacez le # (avant le "!" ou le "%").
#!pip install geopy
#%pip install --user geopandas # syntaxe d'installation depuis un Notebook
Nous allons dans un premier temps calculer des distances par la méthode de la distance euclidienne entre deux points. Celles que l'on calcule dans le plan quand on peut faire abstraction de la rotondité de la Terre. Puis, nous comparerons les résultats obtenus par la méthode de la distance euclidienne avec ceux obtenus par la distance géodésique et ceux obtenus par la distance sur le grand cercle (cf. tdpy12sig.ipynb), pour savoir si, sur le territoire de la France métropolitaine, il faut substituer l'une à l'autre en fonction de la latitude (https://pro.arcgis.com/fr/pro-app/latest/tool-reference/spatial-analyst/geodesic-versus-planar-distance.htm).
Nous allons dans un premier temps calculer la distance euclidienne en appliquant la formule "standard", dite "naïve". Puis, ensuite, nous appliquerons quelques formules piochées dans des bibliothèques (de confiance (: !).
Nous allons commencer par la distance euclidienne telle que lue dans les livres !
Nous pouvons aisément calculer la distance entre deux points, même de plus de deux dimensions (x, y, z par exemple) en calculant la différence entre les dimensions des deux points, au carré.
En utilisant le théorème de Pythagore, par exemple. Méthode naïve, qui ne nécessite pas de bibliothèques, juste de se souvenir de la formule :
{d = [(x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 + etc]^(1/2)}.
On se contentera des coordonnées x et y (pas de z ou autres) si la distance est calculée dans un plan, en deux dimensions, donc.
# Calcul de la distance euclidienne entre deux points
# Le couple (lat, long) est un "tuple", une collection ordonnée de plusieurs éléments (on ne permute pas les x, y et z !)
# Cette fois-ci, il est exprimé en coordonnées planes (ici du RGF 1993 - Lambert-93).
# https://datagy.io/python-euclidian-distance/
point_1 = (942174.00, 6858371.20)
point_2 = (956269.10, 6838945.50)
def naive_euclidian_distance(point1, point2):
return sum([(point1[x] - point2[x]) ** 2
for x in range(len(point1))]) ** 0.5 # suppose que tous les points aient la même "longueur" que le point 1
print(naive_euclidian_distance(point_1, point_2))
La bibliothèque Python "math" nous servira de support.
En partie inspiré de How to Calculate Euclidean distance in Python https://pythonsolved.com/how-to-calculate-euclidean-distance-in-python/
# Utilisation de la fonction "math.sqrt"
# import de la bibliothèque
import math
x1 = 942174.00
x2 = 956269.10
y1 = 6858371.20
y2 = 6838945.50
distance = math.sqrt(((x2 - x1) ** 2) + (y2 - y1) ** 2)
# Attention à ne pas mettre des "" autour de 'math.sqrt' mais des '' entres les parenthèses du print !
print("La distance euclidienne, entre les deux points, en utilisant la formule 'math.sqrt', est de : ", str(distance))
# Une écriture plus pratique des coordonnées des deux points
point_1 = (942174.00, 6858371.20) # tuples à deux valeurs
point_2 = (956269.10, 6838945.50)
x1 = point_1[0] # extraction de la position 1 "[0]" du tuple point_1
x2 = point_2[0]
y1 = point_1[1] # extraction de la position 2 "[1]" du tuple point_2
y2 = point_2[1]
distance = math.sqrt(((x2 - x1) ** 2) + (y2 - y1) ** 2)
# Si le nombre de décimales ne convient pas, il est possible de le modifier
print("La distance euclidienne, entre les deux points, utilisant cette formule est de : ", f'{distance:.2f}', 'mètres')
# Si le nombre de décimales ou l'unité ne conviennent pas, il est possible de les modifier
print("La distance euclidienne, entre les deux points, utilisant cette formule est de : ", f'{distance/1000:.5f}', 'kilomètres')
Enfin, il existe plein de fonctions qui calculent directement la distance euclidienne.
La fonction "math.dist" de la bibliothèque "math" semble être la plus rapide (cf. https://datagy.io/python-euclidian-distance/).
Plus d'info sur la bibliothèque "math" ici : cf. https://docs.python.org/3/library/math.html.
# Pour aller plus vite => utilisation de la formule "math.dist" de la bibliothèque "math"
distance = math.dist(point_1, point_2)
# Attention à ne pas mettre des "" autour de 'math.dist' mais des '' entres les parenthèses du print !
print("La distance euclidienne, entre les deux points, utilisant la formule 'math.dist' est de : ", f'{distance:.2f}', 'mètres')
À titre d'exemple, s'il fallait inclure l'altitude dans le calcul de la distance, nous intégrerions les coordonnées z (l'altitude) obtenues sur Google Earth pour les points : point_1 et point_2 qui sont respectivement les coins NW et SE de la coupure (feuille) IGN de St-Nicolas-de-Port, Série bleue 3415E.
Reprenons la formule "naïve" de calcul de la distance euclidienne vue plus haut.
# Ses triplets de coordonnées deviennent :
point_1 = (942174.00, 6858371.20, 252.7014493)
point_2 = (956269.10, 6838945.50, 249.7111228)
def naive_euclidian_distance(point1, point2):
return sum([(point1[x] - point2[x]) ** 2
for x in range(len(point1))]) ** 0.5 # suppose que tous les points aient la même "longueur" que le point 1
print(naive_euclidian_distance(point_1, point_2))
En raison de la faible différence d'altitude, sur le plateau lorrain dans ce secteur, c'est au dixième de millimètres qu'apparaissent les écarts de distances entre : avec et sans l'altitude (z) dans le calcul de la distance euclidienne sur la diagonale NW-SE de la feuille 3415E de la Série bleue de l'IGN.
Sauriez-vous refaire ce calcul, avec le bon "réglage" de décimales, en utilisant une formule comme "math.dist" ?
# Calcul de la distance euclidienne avec la fonction 'math.dist'
Maintenant que l'on sait calculer la distance euclidienne, à quelles échelles peut-on l'utiliser ?
Dans le cas des projections cylindriques directes, il est écrit que les déformations sont plus faibles aux basses latitudes qu'aux hautes latitudes (https://pro.arcgis.com/fr/pro-app/latest/tool-reference/spatial-analyst/geodesic-versus-planar-distance.htm) où l'on a recours aux projections stéréographiques. Dans le cas de la France métropolitaine, avec la projection conique qui caractérise son système légale - celle du RGF 1993 - Lambert-93, est-il préférable d'opter pour une mesure de distance plutôt qu'une autre ?
Nous allons utiliser la Série orange au 1/50 000 de l'IGN (Institut national de l'information géographique et forestière) pour tester le décalage possible entre distances calculées selon les méthodes : euclidiennes ; géodésiques et sur le grand cercle pour répondre à cette question.
Note : pour exécuter les cellules suivantes, stocker les données dans un répertoire nommé "data" placé au même niveau que le Notebook. Sinon, penser à modifier les chemins d'accès aux données avant d'exécuter les scripts.
Deux bibliothèques vont être utilisées pour afficher les secteurs sur lesquels portent les calculs de distances. Il s'agit de Geopandas et de Matplotlib. Ce choix (arbitraire) permet d'afficher directement des couches projetées en RGF 1993 - Lambert-93, sans transformations. Une superposition sur un fond comme Folium nécessiterait une conversion en coordonnées géographiques dans le datum (système géodésique) WGS84.
# Chargement de Geopandas
import geopandas as gpd # Si vous obtenez le message d'erreur : "ModuleNotFoundError: No module named 'geopandas'"
# => redémarrez le noyau (Kernel/Restart).
Le tableau d'assemblage de la Série orange au 1/50 000 utilisé ici est du Type 1950 (https://geoservices.ign.fr/recherche?search=tableau, https://geoservices.ign.fr/sites/default/files/2021-11/GRAPHE-MOSAIQUAGE-SCAN50%C2%AE.zip).
Comme nous le verrons, lorsque nous afficherons les cinq premières lignes du tableau [cf. SC50_cartes_L93.head()], la Série orange est la continuité de la carte topographique de la France au 1:50 000 dite Type 1922 (un type est une édition avec ses caractéristiques, sa légende), voire du type destiné à l'armée (Type M).
L'extension "_L93" précise que la couche est projetée dans le système légal français, le RGF 1993 - Lambert-93, depuis l'an 2000.
# Lecture d'une 1ère couche shape
# Lecture du tableau d'assemblage de la Série orange au 1/50 000
shp = 'data/SC50_cartes_L93.shp'
# Création d'un tableau Geopandas
SC50_cartes_L93 = gpd.read_file(shp) # lire les données de la couche vectorielle définies dans "shp"
# Affichage des 5 premières lignes du tableau "SC50_cartes_L93.dbf"
SC50_cartes_L93.head()
Vous constaterez dans la colonne TYPE_CARTO que la Série orange, selon les coupures, est constituée de Type 1922, de Type M, etc.
# Affichage du nombre de lignes et colonnes du tableau de données
SC50_cartes_L93.shape
# le nombre des lignes indique le nombre de feuilles contenues dans la Série orange
Lecture de la feuille SO_2302 qui se situe au centre-nord du tableau d'assemblage de la Série orange.
La feuille SO_2302, celle de Dunkerque, servira d'étalon pour mesurer la distance de sa base nord.
# Lecture d'une 2ème couche shape
# Lecture de la feuille 3415 de la Série orange au 1/50 000
shp = 'data/SO_2302_L93.shp'
# Création d'un tableau Geopandas
SO_2302_L93 = gpd.read_file(shp) # lire les données de la couche vectorielle définies dans "shp"
# Affichage de la première (et unique ici) ligne (head affiche usuellement les 5 premières lignes)
SO_2302_L93.head()
Lecture de la feuille SO_2350 qui se situe au centre-sud du tableau d'assemblage de la Série orange.
La feuille SO_2350, celle de Prats-de-Mollo-La-Preste, servira d'étalon pour mesurer la distance de sa base sud.
# Lecture d'une 2ème couche shape
# Lecture de la feuille 3415 de la Série orange au 1/50 000
shp = 'data/SO_2350_L93.shp'
# Création d'un tableau Geopandas
SO_2350_L93 = gpd.read_file(shp) # lire les données du geopackage définies dans "shp"
# Affichage de la première (et unique ici) ligne (head affiche usuellement les 5 premières lignes)
SO_2350_L93.head()
Pour visualiser la position des couches (SC50_cartes_L93, SO_2302_L93 et SO_2350_L93), les unes par rapport aux autres, elles vont être superposées avec la fonction ".plot()" de Matplotlib.
# Chargement de Matplotlib
# Chargement des fonctions contenues dans les paquets et sous paquets de "matplotlib"
import matplotlib.pyplot as plt
# Création de la carte avec la fonction "plot" en superposant les couches (overlayed maps)
# Pour en savoir plus sur les paramètres de "fig, ax" : https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html
fig, ax = plt.subplots(figsize=(14, 12)) # figsize=(width, height)
SC50_cartes_L93.plot(color='orange',ax=ax)
SO_2302_L93.plot(markersize=0.2, color='blue',ax=ax) # affiche en rouge la position de la feuille SO_2302_L93
SO_2350_L93.plot(markersize=0.2, color='green',ax=ax) # affiche en rouge la position de la feuille SO_2350_L93
plt.title("Position des feuilles SO_2302_L93 (en bleu) et SO_2350_L93 (en vert)\n dans le tableau d'assemblage de la Série orange")
# Affichage de la carte
plt.show()
# Les coordonnées marginales sont à multiplier par 1 000 000 de mètres
• Rappel : bibliothèques nécessaires pour les calculs dans ce paragraphe
# import de la bibliothèque math pour le calcul de la distance euclidienne
import math
# import de la bibliothèque geopy pour le calcul de la distance géodésique et du grand cercle
import geopy
# from geopy.distance import geodesic
○ Distance euclidienne
# Calcul de la distance euclidienne avec la formule "math.dist" de la bibliothèque "math"
# Coordonnées en RGF 1993 - Lambert-93 de la feuille SO_2302_L93
point_2302_NW_L93 = (7114159.10, 640581.30)
point_2302_NE_L93 = (7114159.10, 666055.46)
distance_eucl = math.dist(point_2302_NW_L93, point_2302_NE_L93)
# Affichage avec deux décimales
print("La distance euclidienne est de :", f'{distance_eucl:.2f}', 'mètres')
○ Distance géodésique
# Calcul de la distance géodésique
from geopy import distance
# Coordonnées angulaires en WGS84 de la feuille SO_3415
point_2302_NW_WGS84 = (51.119981, 2.156407)
point_2302_NE_WGS84 = (51.119981, 2.516486)
distance_geod = distance.distance(point_2302_NW_WGS84, point_2302_NE_WGS84).m
# Affichage avec deux décimales
print("La distance géodésique est de :", f'{distance_geod:.2f}', 'mètres')
○ Distance sur le grand cercle
# Calcul de la distance sur le grand cercle
from geopy.distance import great_circle
# Coordonnées angulaires en WGS84 de la feuille SO_3415
point_2302_NW_WGS84 = (51.119981, 2.156407)
point_2302_NE_WGS84 = (51.119981, 2.516486)
dist_great_circle = great_circle(point_2302_NW_WGS84, point_2302_NE_WGS84).meters
# Affichage avec deux décimales
print("La distance calculée sur le grand cercle est de :", f'{dist_great_circle:.2f}', 'mètres')
○ Évaluation du décalage
• Décalage entre distance euclidienne et distance géodésique
# Évaluation du décalage entre la distance euclidienne et la distance géodésique sur la diagonale NW-SE d'une Série orange
offset = distance_eucl / distance_geod
print("La valeur du facteur d'erreur est de :", f'{offset:.3f}')
• Décalage entre distance euclidienne et distance sur le grand cercle
# Évaluation du décalage entre la distance euclidienne et la distance calculée sur le grand cercle pour la diagonale NW-SE d'une Série orange
offset = distance_eucl / dist_great_circle
print("La valeur du facteur d'erreur est de :", f'{offset:.3f}')
L'ordre de grandeur est d'environ 1 p.100 entre les distances angulaires et la distance euclidienne (calculée dans le plan).
○ Distance euclidienne
# Calcul de la distance euclidienne avec la formule "math.dist" de la bibliothèque "math"
# Coordonnées en RGF 1993 - Lambert-93 de la feuille SO_3415_L93
point_2350_SW_L93 = (6133453.19, 630324.59)
point_2350_SE_L93 = (6133453.19, 660202.02)
distance_eucl = math.dist(point_2350_SW_L93, point_2350_SE_L93)
# Affichage avec deux décimales
print("La distance euclidienne est de :", f'{distance_eucl:.2f}', 'mètres')
Distance à comparer avec celle calculée au point 3.2.4.1. Base supérieure (nord) de la feuille SO_2302_L93.
○ Distance géodésique
# Calcul de la distance géodesique
from geopy import distance
# Coordonnées angulaires en WGS84 de la feuille SV_121
point_2350_SW_WGS84 = (42.300062, 2.156453)
point_2350_SE_WGS84 = (42.300062, 2.516696)
distance_geod = distance.distance(point_2350_SW_WGS84, point_2350_SE_WGS84).m
# Affichage avec deux décimales
print("La distance géodésique est de :", f'{distance_geod:.2f}', 'mètres')
○ Distance sur le grand cercle
# Calcul de la distance sur le grand cercle entre deux points
from geopy.distance import great_circle
# Coordonnées angulaires en WGS84 de la feuille SV_121
point_2350_SW_WGS84 = (42.300062, 2.156453)
point_2350_SE_WGS84 = (42.300062, 2.516696)
dist_great_circle = great_circle(point_2350_SW_WGS84, point_2350_SE_WGS84).meters
# Affichage avec deux décimales
print("La distance calculée sur le grand cercle est de :", f'{dist_great_circle:.2f}', 'mètres')
○ Évaluation du décalage
• Décalage entre distance euclidienne et distance géodésique
# Évaluation du décalage entre la distance euclidienne et la distance géodésique sur la diagonale NW-SE d'une Top100 ex Série bleue
offset = distance_eucl / distance_geod
print("La valeur du facteur d'erreur est de :", f'{offset:.3f}')
• Décalage entre distance euclidienne et distance sur le grand cercle
# Évaluation du décalage entre la distance euclidienne et la distance calculée sur le grand cercle pour la diagonale NW-SE d'une Série bleue
offset = distance_eucl / dist_great_circle
print("La valeur du facteur d'erreur est de :", f'{offset:.3f}')
L'ordre de grandeur est inférieur à 1 p.100 entre les distances angulaires et la distance euclidienne (calculée dans le plan).
En conclusion, en utilisant une projection de type RGF 1993 - Lambert-93 sur le territoire métropolitain de la France, quelque soit la latitude, il semble possible d'utiliser indifféremment les distances angulaires et la distance euclidienne pour mesurer des distances de l'ordre de la trentaine de kilomètres sans introduire de distorsions trop importantes.
En serait-il de même avec une autre projection comme une UTM ?