Comment obtenir l’altitude de surface pour une latitude et une longitude spécifiques avec Python

Introduction

Une manière simple d’obtenir l’altitude de surface pour un point défini par une latitude et une longitude consiste à échantillonner un modèle numérique d’élévation. Ici, nous utilisons le Copernicus DEM GLO-30, un modèle numérique de surface global à 30 m, disponible sous forme de fichiers Cloud Optimized GeoTIFF sur AWS. Il représente la surface terrestre, incluant le terrain, les bâtiments et la végétation. Les valeurs d’altitude sont exprimées en mètres.

Installer les packages

1
conda install -c conda-forge rasterio numpy

Exemple avec un point

Dans cet exemple, nous récupérons l’altitude autour du Mont Blanc, situé dans les Alpes françaises, près de la frontière entre la France et l’Italie.

1
2
3
4
5
6
7
8
import os
import rasterio
import numpy as np

os.environ["AWS_NO_SIGN_REQUEST"] = "YES"

lon = 6.8652
lat = 45.8326

Ouvrir la bonne tuile Copernicus DEM

Les tuiles Copernicus DEM sont organisées sur AWS dans des dossiers de 1° × 1°.

Pour la latitude 45.8326 et la longitude 6.8652, la tuile correspondante est :

1
2
3
4
5
6
7
tile = "Copernicus_DSM_COG_10_N45_00_E006_00_DEM"

url = (
    f"/vsis3/copernicus-dem-30m/"
    f"{tile}/"
    f"{tile}.tif"
)

Ouvrir la tuile :

1
2
with rasterio.open(url) as src:
    print(src.profile)

Extraire l’altitude en un point

1
2
3
4
with rasterio.open(url) as src:
    elevation_m = list(src.sample([(lon, lat)]))[0][0]

print(elevation_m)

Exemple de sortie :

1
4797.0034

Cela signifie :

1
2
4797.0034 mètres au-dessus du niveau de la mer
 4.78 km

Rasterio peut ouvrir des fichiers locaux, des URLs ou des chemins virtuels compatibles avec GDAL, comme /vsis3/. La méthode sample() permet d’extraire les valeurs d’un raster à partir de coordonnées géographiques.

Visualiser le point sur une carte satellite avec Folium

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
import folium

# Coordonnées du Mont Blanc
lon = 6.8652
lat = 45.8326

# Créer la carte
m = folium.Map(
    location=[lat, lon],
    zoom_start=12,
    tiles=None
)

# Ajouter l’imagerie satellite Esri
folium.TileLayer(
    tiles=(
        "https://server.arcgisonline.com/ArcGIS/rest/services/"
        "World_Imagery/MapServer/tile/{z}/{y}/{x}"
    ),
    attr="Esri World Imagery",
    name="Esri Satellite",
    overlay=False,
    control=True,
).add_to(m)

# Ajouter un marqueur
folium.CircleMarker(
    location=[lat, lon],
    radius=8,
    color="red",
    fill=True,
    fill_color="red",
    fill_opacity=0.9,
    popup=(
        f"<b>Mont Blanc</b><br>"
        f"Latitude: {lat}<br>"
        f"Longitude: {lon}"
    ),
).add_to(m)

# Contrôle des couches, optionnel
folium.LayerControl().add_to(m)

# Afficher la carte
m

Comment obtenir l’altitude de surface pour une latitude et une longitude spécifiques avec Python
Comment obtenir l’altitude de surface pour une latitude et une longitude spécifiques avec Python

Échantillonner plusieurs points

1
2
3
4
5
6
7
8
9
lons = np.array([-104.98, -104.50])
lats = np.array([40.20, 40.80])

coords = list(zip(lons, lats))

with rasterio.open(url) as src:
    elevations_m = np.array([v[0] for v in src.sample(coords)])

print(elevations_m)

Exemple de sortie :

1
[1515.6829 1583.732 ]

Ajouter l’altitude à un DataFrame

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
import pandas as pd

df = pd.DataFrame({
    "longitude": [-104.98, -104.50],
    "latitude": [40.20, 40.80],
})

coords = list(zip(df["longitude"], df["latitude"]))

with rasterio.open(url) as src:
    df["elevation_m"] = [v[0] for v in src.sample(coords)]

df["elevation_km"] = df["elevation_m"] / 1000

print(df)

Fonction pour une seule tuile

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
def sample_elevation_from_tile(lons, lats, tile):
    url = (
        f"/vsis3/copernicus-dem-30m/"
        f"{tile}/"
        f"{tile}.tif"
    )

    coords = list(zip(lons, lats))

    with rasterio.open(url) as src:
        elevations = np.array([v[0] for v in src.sample(coords)])

    return elevations

Utilisation :

1
2
3
4
5
6
7
tile = "Copernicus_DSM_COG_10_N40_00_W105_00_DEM"

df["elevation_m"] = sample_elevation_from_tile(
    df["longitude"].values,
    df["latitude"].values,
    tile
)

Note importante :

Cet exemple fonctionne uniquement si tous les points se trouvent dans la même tuile DEM. Pour des jeux de données plus grands, il faut d’abord regrouper les points par nom de tuile, puis échantillonner chaque tuile séparément.

Visualiser plusieurs points sur une carte satellite avec Folium

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import folium

# Centre de la carte
center_lat = lats.mean()
center_lon = lons.mean()

m = folium.Map(
    location=[center_lat, center_lon],
    zoom_start=12,
    tiles=None
)

# Imagerie satellite
folium.TileLayer(
    tiles=(
        "https://server.arcgisonline.com/ArcGIS/rest/services/"
        "World_Imagery/MapServer/tile/{z}/{y}/{x}"
    ),
    attr="Esri World Imagery",
    name="Esri Satellite",
    overlay=False,
    control=True,
).add_to(m)

# Ajouter les points
for i, (lat, lon, elev) in enumerate(
    zip(lats, lons, elevations_m)
):

    popup_text = (
        f"<b>Point {i}</b><br>"
        f"Latitude: {lat:.4f}<br>"
        f"Longitude: {lon:.4f}<br>"
        f"Altitude: {elev:.1f} m"
    )

    folium.CircleMarker(
        location=[lat, lon],
        radius=7,
        color="red",
        fill=True,
        fill_color="red",
        fill_opacity=0.9,
        popup=popup_text,
    ).add_to(m)

folium.LayerControl().add_to(m)

m

Enregistrer la carte au format HTML :

1
m.save("mont_blanc_elevation_map.html")

Trouver la bonne tuile Copernicus DEM pour une latitude et une longitude

Le Copernicus DEM est organisé en tuiles géographiques de 1° × 1°.

Chaque nom de tuile suit cette structure :

1
Copernicus_DSM_COG_10_<LAT>_00_<LON>_00_DEM

Exemples :

1
2
Copernicus_DSM_COG_10_N45_00_E006_00_DEM
Copernicus_DSM_COG_10_N40_00_W105_00_DEM

Étape 1 — Déterminer les indices de tuile

Pour une coordonnée :

1
2
Latitude  = 45.8326
Longitude = 6.8652

on prend la partie entière inférieure :

1
2
3
4
5
6
7
8
9
import numpy as np

lat = 45.8326
lon = 6.8652

lat_floor = int(np.floor(lat))
lon_floor = int(np.floor(lon))

print(lat_floor, lon_floor)

Sortie :

1
45 6

Étape 2 — Convertir au format de nommage Copernicus

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
if lat_floor >= 0:
    lat_str = f"N{lat_floor:02d}"
else:
    lat_str = f"S{abs(lat_floor):02d}"

if lon_floor >= 0:
    lon_str = f"E{lon_floor:03d}"
else:
    lon_str = f"W{abs(lon_floor):03d}"

print(lat_str, lon_str)

Sortie :

1
N45 E006

Étape 3 — Construire le nom de la tuile

1
2
3
4
5
6
7
tile = (
    f"Copernicus_DSM_COG_10_"
    f"{lat_str}_00_"
    f"{lon_str}_00_DEM"
)

print(tile)

Sortie :

1
Copernicus_DSM_COG_10_N45_00_E006_00_DEM

Fonction complète

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
import numpy as np

def get_copernicus_tile(lat, lon):

    lat_floor = int(np.floor(lat))
    lon_floor = int(np.floor(lon))

    # Chaîne pour la latitude
    if lat_floor >= 0:
        lat_str = f"N{lat_floor:02d}"
    else:
        lat_str = f"S{abs(lat_floor):02d}"

    # Chaîne pour la longitude
    if lon_floor >= 0:
        lon_str = f"E{lon_floor:03d}"
    else:
        lon_str = f"W{abs(lon_floor):03d}"

    tile = (
        f"Copernicus_DSM_COG_10_"
        f"{lat_str}_00_"
        f"{lon_str}_00_DEM"
    )

    return tile

Exemple : Mont Blanc

1
2
3
4
5
6
lat = 45.8326
lon = 6.8652

tile = get_copernicus_tile(lat, lon)

print(tile)

Sortie :

1
Copernicus_DSM_COG_10_N45_00_E006_00_DEM

Exemple : Colorado

1
2
3
4
5
6
lat = 40.20
lon = -104.98

tile = get_copernicus_tile(lat, lon)

print(tile)

Sortie :

1
Copernicus_DSM_COG_10_N40_00_W105_00_DEM

Construire automatiquement l’URL AWS

1
2
3
4
5
6
7
8
9
tile = get_copernicus_tile(lat, lon)

url = (
    f"/vsis3/copernicus-dem-30m/"
    f"{tile}/"
    f"{tile}.tif"
)

print(url)

Sortie :

1
2
3
/vsis3/copernicus-dem-30m/
Copernicus_DSM_COG_10_N45_00_E006_00_DEM/
Copernicus_DSM_COG_10_N45_00_E006_00_DEM.tif

Chaque tuile DEM :

  • couvre 1 degré de latitude ;
  • couvre 1 degré de longitude ;
  • contient une grille d’altitude à 30 m de résolution.

Note importante concernant les coordonnées négatives :

Coordonnée Préfixe
Latitude positive N
Latitude négative S
Longitude positive E
Longitude négative W

Exemples :

Localisation Tuile
France N45 E006
Colorado N40 W105
Chili S33 W071

Références

Ressource Pourquoi c’est utile
Copernicus DEM sur AWS Accès au jeu de données et structure S3
Manuel produit Copernicus DEM Référence verticale et détails du produit
Documentation Rasterio Lecture de GeoTIFFs et échantillonnage raster
Image

of