Les API de l’IGN sont bien pratiques, plusieurs existent : calcul d’itinéraires, d’isochrones autour d’une position, géocodage d’adresses, ou calcul d’altitude en un point donné. Comme la plupart des API, on les utilise sous contrainte. Pour l’API altimétrie : moins de 5 000 positions « alti-codables » par requête, et moins de 5 requêtes par seconde. Avec plusieurs millions de positions à coder, on peut s’attendre à des heures d’interrogation, et des essais répétés pour tenir compte des requêtes rejetées.
Voyons comment alticoder à coup sûr et en un temps raisonnable de tels volumes avec DuckDB, ce moteur SQL portable et véloce. En débouchant sur quelques lignes de SQL moderne, cela va nous amener à :
• comprendre la doc IGN, qui est encore perfectible 😉;
• savoir mobiliser l’API altimétrie en passant les coordonnées via le protocole POST ;
• discrétiser l’espace avec H3 pour maitriser le volume à alticoder et le ramener à moins de 5 millions de positions ;
• utiliser le SQL récursif dans une variante récente particulièrement élégante ;
• exploiter le nouveau format Raquet (raster en parquet) pour alticoder certaines positions au-delà des frontières, hors du champ de l’IGN.
Cas d'usage : les vols en parapente
J’ai entrepris d’étudier tous les vols en parapente dans les Pyrénées, tels que déclarés et stockés sur le site de la Coupe fédérale de distance de la FFVL. Comme pilote, je m’intéresse surtout aux phases les plus spectaculaires d’un vol, aux fortes ascendances et à leur symétrique, les descendances (ou « dégueulantes »). Je souhaite les examiner de manière statistique, plutôt que de les analyser individuellement.
Il y a plus de 12 000 vols déclarés dans les Pyrénées depuis 2000. Chaque vol est décrit par une trace, la suite des positions GPS et des altitudes mesurées par l’altimètre du pilote.
Voici un exemple d’un joli vol réalisé par un ami pilote ariégeois, avec son tracé en 2D et son profil altimétrique :
Cette trace comprend 13 000 points pour un vol de 3,5 heures, soit une position par seconde.
En compilant dans une base parquet tous les vols pyrénéens, on aboutit à près de 100 millions de triplets coordonnées GPS + altitude du pilote.
J’ai besoin de l’API altimétrie de l’IGN pour déterminer l’altitude du terrain survolé à chaque instant du vol. Par différence, j’en déduis la hauteur sol, soit la distance verticale qui sépare le pilote du relief sous ses pieds. C’est une donnée importante pour apprécier la sécurité d’une phase de vol : prendre une ascendance à proximité du sol est toujours plus exigeant et engageant.
Rappelons que l’API d’altimétrie gère un maximum de 5 000 positions par demande. À l’usage, une telle requête prend 4-5 secondes. Pour alticoder 100 millions de positions, il faudrait près de deux heures. J’aimerais réduire ce temps à quelque chose de plus raisonnable, de l’ordre de la demi-heure.
Comme je connais bien la grille hexagonale H3, je songe à discrétiser les positions de vol selon ce maillage. Je vais chercher la résolution offrant le meilleur compromis entre précision et nombre final d’hexagones. Le niveau 11 de la grille H3 décrit des hexagones d’un rayon de 25 mètres, ce qui correspond à la résolution de la BD Alti de l’IGN.
En remplaçant chaque position de vol à la seconde par le centre de chaque hexagone de niveau 11 la contenant, je réduis mon volume de points à alticoder de 100 millions à 3 millions environ (pour les Pyrénées). Découpés en requête de 5 000 positions, cela me donne en gros 620 appels de l’API.
Dans un monde idéal, il suffirait d’enchainer ces requêtes. Chacune prenant 5 secondes, il est exclu que le quota de 5 requêtes max par seconde soit dépassé. Dans le monde réel, certaines de ces requêtes vont échouer, parfois 2/620, parfois 12/620, ou bien aucune… C’est aléatoire, c’est la vie interne du serveur de l’IGN ; mais c’est comme ça avec la plupart des API.
Il faut donc pouvoir itérer :
- lancer ces 620 requêtes chainées,
- constater celles qui ont échoué,
- relancer ces dernières,
- et recommencer jusqu’à ce qu’il ne reste aucune requête en échec !
Il faudra parfois jusqu’à 5 itérations pour tout alticoder. Là où en programmation classique une boucle « while » traite ce sujet, en SQL récursif « avec clé » on trouvera une alternative élégante et sûre.
Place au code commenté pour les geeks !
Je charge la base des 10 millions de positions et, avec l’extension H3, j’en dérive le centroïde des hexagones de résolution 11 qui les recouvrent.
La table h3_11_pyr_0 n’a plus que 3 millions de lignes (hexagones distincts identifiés par un h3_index).
INSTALL H3 FROM COMMUNITY ;
LOAD H3 ;
CREATE OR REPLACE TABLE h3_11_pyr_0 AS
WITH t1 AS (
FROM 'hf://datasets/ericmauviere/cfd_ffvl/cfd_traces_PYR_2000_2026-03.parquet'
SELECT DISTINCT h3_latlng_to_cell(lat, lon, 11) AS h3_index
)
FROM t1
SELECT h3_index,
h3_cell_to_latlng(h3_index)[1].round(5) AS h3_lat,
h3_cell_to_latlng(h3_index)[2].round(5) AS h3_lon ;
Pour l’API IGN Altimétrie, je dois organiser mes coordonnées en paquets de 5 000 max (par sécurité, je réduis à 4 990). Longitudes et latitudes seront transmises chacune sous forme d’une chaine délimitée par des |.
CREATE OR REPLACE TABLE h3_11_pyr_chunked AS
WITH t0 AS (
FROM h3_11_pyr_0
SELECT *,
c: ((row_number() OVER(ORDER BY h3_index))/4990)::INT
)
FROM t0
SELECT c,
list(h3_index) AS h3_indexes,
lon_delim: array_to_string(list(h3_lon), '|'),
lat_delim: array_to_string(list(h3_lat), '|')
GROUP BY c
ORDER BY c
;
Pour utiliser l’API altimétrie, je respecte le formalisme suivant, décrit dans une macro qui va traiter d’un coup un paquet de coordonnées :
INSTALL http_client FROM community ;
LOAD http_client ;
CREATE OR REPLACE MACRO get_elev_by_chunk(lon_delim, lat_delim) AS
http_post( 'https://data.geopf.fr/altimetrie/1.0/calcul/alti/rest/elevation.json',
headers => MAP {
'accept': 'application/json',
},
params => MAP {
'lon': lon_delim,
'lat': lat_delim,
'resource': 'ign_rge_alti_wld',
'delimiter': '|',
'indent': 'false',
'measures': 'false',
'zonly': 'false'
}
) ;
Testons la macro sur un premier paquet :
FROM h3_11_pyr_chunked
SELECT c,
res: get_elev_by_chunk(lon_delim, lat_delim)::json
WHERE c = 1 ;
Voici un aperçu de la réponse de l’API.
L’array « elevations » me donne bien un z d’altitude (admirez la précision) du terrain.
Status = 200 me confirme que la requête s’est exécutée sans problème.
Mais ce ne sera pas toujours le cas, la requête peut échouer, rarement il est vrai, pour des raisons obscures (Bad Gateway, HTTP POST request failed…).
Plutôt que lancer naïvement l’alticodage de tous les paquets (620) :
FROM h3_11_pyr_chunked
SELECT c,
res: get_elev_by_chunk(lon_delim, lat_delim)::json ;
je vais écrire une requête SQL récursive, avec la nouvelle option (DuckDB) USING KEY :
CREATE OR REPLACE TABLE h3_11_pyr_elev0 AS
WITH RECURSIVE ign_api_elevation(c, h3_indexes, lon_delim, lat_delim, res, status, reason, iter)
USING KEY(c) AS (
FROM h3_11_pyr_chunked
SELECT
c, h3_indexes, lon_delim, lat_delim,
res: get_elev_by_chunk(lon_delim, lat_delim),
status: [res->>'status'],
reason: [res->>'reason'],
iter: 1
UNION ALL
FROM ign_api_elevation
SELECT
c, h3_indexes, lon_delim, lat_delim,
res: get_elev_by_chunk(lon_delim, lat_delim),
status: reason.list_append(res->>'status'),
reason: reason.list_append(res->>'reason'),
iter: iter + 1
WHERE status[-1] <> '200' AND iter < 10
)
FROM ign_api_elevation;
La lecture est ici plus exigeante. Mais pour celles ou ceux qui connaissent déjà le SQL récursif, il s’agit d’une variante plus facile à saisir. Ce n’est plus un process cumulatif où les lignes s’empilent, mais un algorithme récursif de mise à jour d’une table, qui préserve son nombre de lignes en utilisant une clé primaire (USING KEY). UNION ici signifie plutôt UPDATE.
La première étape lance l’alticodage des 620 paquets.
La seconde considère les paquets en échec (status <> ‘200’), mettons 10, et les relance. Le résultat de cette relance met à jour 10 lignes sur cette table de 620 paquets. Puis considère à nouveau les éventuels échecs résiduels (mettons 2). Et poursuit jusqu’à ce qu’il n’y en ait plus.
Par sécurité, le nombre total d’itérations est borné. Dans mes tests, je n’ai constaté qu’un maximum de 3 itérations, dont je peux tracer l’historique jusqu’à résolution :
Je vais ensuite décoder la réponse pour chaque paquet et déplier tout cela en une table de mes 3 millions d’hexagones :
CREATE OR REPLACE TEMP TABLE h3_11_pyr_elev1 AS
FROM h3_11_pyr_elev0
SELECT c,
h3_index: unnest(h3_indexes),
unnest(
(res->>'body'->>'elevations')::STRUCT(lon DOUBLE, lat DOUBLE, z DOUBLE, acc VARCHAR)[],
recursive := TRUE
),
status, reason, iter
ORDER BY c ;
Tant qu’à faire, j’affecte un code département à chaque hexagone par une jointure spatiale avec le fond data.gouv des départements français.
Et je nettoie la donnée d’altitude du terrain : -99999 signifie position non alticodée (typiquement parce que l’hexagone est à l’étranger, hors champ IGN). Je la transforme en NULL et en profite pour arrondir au mètre la donnée d’altitude terrain, renommée pour plus de clarté et pour la distinguer par la suite de l’altitude du pilote.
LOAD SPATIAL ;
CREATE OR REPLACE TABLE h3_11_pyr_elev AS
FROM h3_11_pyr_elev1 AS h3
LEFT JOIN (
FROM st_read('https://object.data.gouv.fr/contours-administratifs/2025/geojson/departements-50m.geojson')
) AS depts ON st_contains(depts.geom, st_point(lon, lat))
SELECT
h3.h3_index, h3.lon, h3.lat,
alt_terrain: NULLIF(round(z)::INT, -99999),
dept_hex: depts.code ;
Il me reste 100 000 (soit 0,3 %) d’altitudes NULL : des pilotes pyrénéens qui ont fait un petit tour en Andorre ou en Espagne, ou qui ont survolé la mer ou l’océan, les chanceux !
Hello Raquet
Pour alticoder à l’étranger, je vais m’appuyer sur une base mondiale de données d’élevation – certes moins précise et homogène que la BD Alti de l’IGN, mais libre et globale.
Ces données sont habituellement stockées dans des images raster (GeoTiff par exemple). En lisant le bon pixel et sa couleur, on peut en déduire une valeur, l’altitude.
Mais, à l’instar de Postgis, il est possible de stocker ces données dans une table, et même, depuis peu, dans un fichier parquet.
Bienvenue ici au projet Raquet (Rasters dans Parquet) !
INSTALL RAQUET FROM community;
LOAD RAQUET;
-- Macro de récupération d'une donnée d'altitude par position - monde entier
CREATE OR REPLACE MACRO getWorldElevation(lon, lat) AS (
SELECT ST_RasterValue(block, band_1, ST_Point(lon,lat), metadata)::int
AS elevation_meters
FROM read_raquet_at(
'https://storage.googleapis.com/raquet_demo_data/world_elevation.parquet', lon, lat)
);
-- Isolement des hexagones non alticodés par l'IGN
CREATE OR REPLACE TABLE missings_elevations AS
FROM h3_11_pyr_elev
SELECT h3_index, getWorldElevation(lon, lat) AS z
WHERE alt_terrain IS NULL ;
-- Mise à jour des données d'élevation manquantes
UPDATE h3_11_pyr_elev h
SET alt_terrain = m.z
FROM missings_elevations m
WHERE h.h3_index = m.h3_index
AND h.alt_terrain IS NULL;
Et voici, en guise d’aboutissement, une carte de ces hexagones H3 coloriés selon l’altitude du terrain. Ces hexagones ont été traversés par des parapentistes. On est ici au sud de l’Ariège, à la frontière avec l’Andorre et l’Espagne. On voit bien que certains territoires sont inlassablement sillonnés, d’autres bien moins ; et que se dessinent certaines « autoroutes » du ciel.
Conclusion
DuckDB est un petit moteur passe-partout, super rapide à mettre en œuvre, mais il brille aussi par son extensibilité : des extensions gérées par l’équipe DuckDB comme SPATIAL, ou des extensions communautaires fiables comme H3, HTTP_CLIENT ou RAQUET.
Un prochain article approfondira l’étude des ascendances et des descendances vécues en parapente dans les Pyrénées.