Passer au contenu

Tutoriel lidarHD

Jean-Matthieu Monnet

lidarHD.Rmd

Téléchargement des données

Tableau d’assemblage des fichiers LiDARHD IGN

Au 4 décembre 2025, il n’existe pas de moyen simple pour récupérer en masse les urls des fichiers de nuages de points de LiDARHD à jour :

  • les connexions wfs (private ou public) ne semblent pas à jour.
  • la nouvelle interface de téléchargement ne donne les urls que pour des fichiers sélectionnés individuellement ou par polygone dont l’emprise reste restreinte.

Un catalogue des fichiers a cependant été constitué et mis à jour régulièrement sur https://browser.stac.teledetection.fr/collections/lidarhd.

La fonction load_classified_ta() reconstitue le tableau d’assemblage des fichiers LiDARHD classés à partir de ce catalogue. Le téléchargement de l’ensemble du catalogue étant très long il est conseillé :

  • soit de le télécharger en spécifiant une région d’intéret (argument bbox) pour restreindre spatialement la requête,
  • soit de le télécharger en entier puis de le sauvegarder.

Si l’argument save_dir est renseigné, une copie sera stockée dans le dossier indiqué. Elle pourra être chargée avec un prochain appel à load_classified_ta() avec l’option url = NULL.

# creer la region d'interet
roi <- sf::st_as_sfc(
  sf::st_bbox(
    sf::st_as_sf(data.frame(x = c(885900, 886100),
                        y = c(6435600, 6435800)),
             coords = c("x", "y"),
             crs = 2154)))
# ajouter un buffer de 2 km et reprojeter en WGS84
bbox <- sf::st_bbox(sf::st_transform(sf::st_buffer(roi, 4000), 4326))
# charger depuis le catalog stac les emprises des dalles dans la bbox
ta_lidarHD <- load_classified_ta(bbox = bbox)
# charger depuis un fichier present dans un dossier (non exécuté)
ta_lidarHD <- load_classified_ta(url = NULL, save_dir = "./data/sig/")

Le tableau d’assemblage conserve deux champs :

La fonction load_classified_ta() ajoute le nom du bloc dans la table d’attributs, ainsi que le nom du dossier qui contient les fichiers. Celui-ci est le même que le nom de bloc pour les fichiers hébergés sur OVH (première version), et de la forme “NUALHD_1-0__LAZ_LAMB93_OO_2025-03-14” pour les fichiers de la géoplateforme. Dans cet exemple :

  • “OO” est le nom du bloc
  • “_2025-03-14” est une date qui devrait disparaître du nom de dossier lors des prochaines mises à jour.

Les dalles sont des carrés d’un kilomètre de côté. Pour des questions de rapidité de traitement, le tableau d’assemblage est réduit dans le tutoriel aux dalles des blocs OM et PM.

ta_lidarHD <-ta_lidarHD[is.element(ta_lidarHD$bloc, c("OM", "PM")),]

Téléchargement des dalles dans une zone d’intérêt

La fonction download_files() télécharge les dalles qui intersectent une région d’intérêt et qui ne sont pas déjà présentes dans le répertoire de destination. La région d’intérêt peut être un objet spatial au format sf, ou un vecteur comportant les noms des blocs à télécharger.

# creer la region d'interet
roi <- sf::st_as_sfc(
  sf::st_bbox(
    sf::st_as_sf(data.frame(x = c(885900, 886100),
                        y = c(6435600, 6435800)),
             coords = c("x", "y"),
             crs = 2154)))
# dossier pour stocker les données
dossier <- "./data/"
# dossier pour stocker les fichiers LiDARHD
dossier_classe <- paste0(dossier, "/LAZ_classe/")
# dossier pour stocker les fichiers normalises
dossier_norm <- paste0(dossier, "/LAZ_norm/")

L’opération d’intersection est longue du fait du grand nombre de dalles LiDAR.

Le paramère prompt est réglé par défaut pour que l’affichage des dalles et le téléchargement soient confirmés par l’utilisateur (recommandé).

Le paramètre buffer permet de télécharger les dalles adjacentes à celles qui intersectent la région d’intérêt, ce qui peut être utile pour que les traitements ultérieurs sur les dalles d’intérêt ne présentent pas d’effet de bord (notamment calcul de MNT ou de la hauteur des points).

Le script affiche le nombre de dalles selon les croisements de catégories suivantes :

  • intersecte ou pas la région d’intérêt (dedans / dehors),
  • déjà téléchargé ou pas (fait / a_faire).

Seules les dalles a_faire seront téléchargées.

# telecharger les fichiers qui intersectent la region d'interet 
download_files(ta_lidarHD, roi = roi, folder = dossier_classe, prompt = TRUE)

Lors du téléchargement, le script place automatiquement les fichiers dans des sous-dossiers correspondant aux “blocs”. Le nom d’origine est conservé. Au cas où des sous-dossiers sont déjà présents avec des fichiers, la fonction propose uniquement le téléchargement des fichiers non présents localement.

Si des fichiers ont été mis à jour (dalles hébergées sur la géoplateforme alors que des fichiers existaient auparavant sur le serveur OVH), ils seront re-téléchargés dans un nouveau dossier avec la nomenclature “géoplateforme”. Il faudra manuellement supprimer ou déplacer les anciens dossiers en doublon. Pour identifier les blocs à traiter, on peut s’aider de la fonction check_folder() et il est important de mettre à jour le catalogue global avec catalog_folder() pour que les anciens fichiers ne soient plus recensés.

# sous-dossiers
list.dirs(dossier_classe, full.names = FALSE, recursive = FALSE)
## [1] "NUALHD_1-0__LAZ_LAMB93_OM_2025-03-25"
## [2] "NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25"

Dans chaque sous-dossier de bloc, un catalogue des dalles LiDAR présentes est créé sous la format d’un objet au format LAS-catalog enregistré une archive R nommée “cata.rda”. Un catalogue de l’ensemble des dalles est créé dans une archive de même nom mais enregistré à la racine du répertoire de destination. Ces catalogues sont mis à jour lorsque des dalles supplémentaires sont ajoutées.

# fichiers dans le sous-dossier du bloc "PM"
dir(paste0(dossier_classe,
           "NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25"))
## [1] "cata.rda"                                   
## [2] "LHD_FXX_0886_6435_PTS_LAMB93_IGN69.copc.laz"
## [3] "LHD_FXX_0886_6436_PTS_LAMB93_IGN69.copc.laz"
## [4] "LHD_FXX_0886_6437_PTS_LAMB93_IGN69.copc.laz"
## [5] "LHD_FXX_0887_6435_PTS_LAMB93_IGN69.copc.laz"
## [6] "LHD_FXX_0887_6436_PTS_LAMB93_IGN69.copc.laz"
## [7] "LHD_FXX_0887_6437_PTS_LAMB93_IGN69.copc.laz"

Le code ci-dessous charge le catalogue, l’affiche avec une couleur par bloc, et ajoute la région d’intérêt avec un contour rouge. Les dalles adjacentes à celles qui intersectent la région d’intérêt ont été également téléchargées.

# mettre a jour les catalogues si besoin
# for (i in list.dirs(dossier_classe, recursive = FALSE, full.names = FALSE))
#   cata <- lidarHD::catalog_bloc(dossier_classe, i)
# cata <- lidarHD::catalog_folder(dossier_classe)
# charger le catalog
load(file = paste0(dossier_classe, "cata.rda"))
lidR::plot(cata,
           col = ifelse(cata$bloc == "NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25",
                        "blue", "green"))
plot(roi, add = TRUE, border = "red")

La fonction check_folder() vérifie que les fichiers et dossiers sont organisés conformément au tableau d’assemblage. Il est utile de faire cette vérification régulièrement au cas où des fichiers locaux seraient rendus obsolètes par des mises à jour sur le serveur IGN.

lidarHD::check_folder(dossier_classe, ta = ta_lidarHD)

Téléchargement des produits dérivés

L’IGN met également à disposition des produits dérivés : modèles numériques de terrain (MNT), de hauteur (MNH) et de surface (MNS).

La fonction load_classified_ta() peut être utilisée pour télécharger le tableau d’assemblage des produits dérivés à partir du serveur WFS de la géoplateforme IGN en changeant l’url par défaut (couche “IGNF_LIDAR-HD_TA:mnt-dalle”, “IGNF_LIDAR-HD_TA:mnh-dalle” ou “IGNF_LIDAR-HD_TA:mns-dalle” pour les MNT, MNH ou MNS).

# charger depuis le serveur WFS (non exécuté)
ta_mnt <- load_classified_ta(url = "https://data.geopf.fr/private/wfs/wfs?apikey=interface_catalogue&SERVICE=WFS&REQUEST=GetFeature&VERSION=2.0.0&TYPENAMES=IGNF_LIDAR-HD_TA:mnt-dalle")
# sauvegarder en gpkg
sf::st_write(ta_mnt, dsn = "./data/sig/ta_mnt.gpkg", layer = "ta_mnt")
# charger depuis un fichier present dans un dossier
ta_mnt <- sf::st_read(dsn = "./data/sig/ta_mnt.gpkg")
## Reading layer `ta_mnt' from data source 
##   `/home/jean-matthieu_monnet/R/lidar/lidarHD/vignettes/data/sig/ta_mnt.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 482740 features and 3 fields (with 4383 geometries empty)
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 98000 ymin: 6046000 xmax: 1243000 ymax: 7111000
## Projected CRS: RGF93 v1 / Lambert-93

La fonction download_mnx() télécharge les dalles qui intersectent une région d’intérêt (objet spatial au format sf) et qui ne sont pas déjà présentes dans le répertoire de destination.

L’opération d’intersection est longue du fait du grand nombre de dalles LiDAR.

L’option prompt est réglée par défaut pour que l’affichage des dalles soit confirmé par l’utilisateur (recommandé). L’option compress permet de comprimer les fichiers.

Le paramètre buffer permet de télécharger les dalles adjacentes à celles qui intersectent la région d’intérêt, ce qui peut être utile pour que les traitements ultérieurs sur les dalles d’intérêt ne présentent pas d’effet de bord.

Le script affiche le nombre de dalles selon les croisements de catégories suivantes :

  • déjà téléchargé ou pas (fait / a_faire).

Seules les dalles a_faire seront téléchargées.

dossier_mnt <- paste0(dossier, "/DTM_0.5m/")
# telecharger les fichiers qui intersectent la region d'interet 
ta_mnt <- download_mnx(ta_mnt, roi = roi, folder = dossier_mnt, prompt = TRUE,
             buffer = FALSE, compress = TRUE, workers = 2)

Contenu des fichiers LAZ

Un tutoriel pour la prise en main des nuages de points LiDAR est disponible sur le site du package lidaRtRee. Il détaille notamment les aspects suivants :

  • contenu d’un fichier LAZ,
  • vérification des paramètres d’acquisition,
  • calculs sur les nuages de points :
    • modèles numériques d’élévation,
    • normalisation du nuage de points.

Les commandes R pour réaliser des opérations similaires sur les nuages de points LiDAR HD seront ajoutées prochainement dans ce tutoriel.

Calcul de produits dérivés

Les paragraphes suivants présentent comment calculer, avec les outils du package lidR, les produits dérivés suivants :

  • nuages de points normalisés avec normalize_files(),
  • modèle numérique de terrain avec compute_dtm(), au cas où ceux-ci ne seraient pas déjà disponible sur le site de l’IGN.

Une fonction alternative basée sur le package lasR permet de réaliser l’une, l’autre ou ces deux opérations de manière plus efficace que ce soit en temps de calcul et en utilisation de la mémoire vive. Les méthodes et paramètres internes sont similaires à ceux des fonctions basées sur lidR.

# produire les fichiers normalisés et les MNT à résolution 1m, en calculant sur 
# deux fichiers à la fois, 5 coeurs max par fichier
process_files(dossier_classe, dossier_norm, res = 1, roi = roi, prompt = TRUE,
              lasR_options = list(ncores = lasR::nested(2, 5)))

Normalisation des fichiers

On appelle couramment “normalisation” d’un nuage de points LiDAR le remplacement de la valeur d’altitude des points par leur hauteur par rapport au sol. La fonction normalize_files() normalise les dalles contenues dans le dossier original qui intersectent une région d’intérêt roi et qui ne sont pas déjà présentes dans le répertoire de destination normalized. La région d’intérêt peut être un objet spatial au format sf, ou un vecteur comportant les noms des blocs à traiter. Si elle n’est pas spécifiée, tous les fichiers sont normalisés (non recommandé). Il faut s’assurer que les dalles adjacentes à celles qui vont être normalisées sont également disponibles dans le dossier original, afin d’éviter les effets de bords.

Le paramètre prompt est réglé par défaut pour que l’affichage des dalles à normaliser et la normalisation soient confirmés par l’utilisateur (recommandé).

Le calcul est parallélisé, le paramètre ncores permet de spécifier combien de process sont lancés simultanément, le facteur limitant sera plutôt la mémoire vive de la machine.

Le script affiche le nombre de dalles selon les catégories suivantes :

  • hors région d’intérêt (hors_zone),
  • dans la région d’intérêt, déjà téléchargé ou pas (fait / a_faire).

Seules les dalles a_faire seront normalisées.

# normaliser les fichiers qui intersectent la region d'interet 
normalize_files(dossier_classe, dossier_norm, roi = roi, prompt = TRUE, ncores = 2L)
# alternative avec les fonctions du package lasR, 2 fichiers à la fois
process_files(dossier_classe, dossier_norm, roi = roi, prompt = TRUE,
              lasR_options = list(ncores = lasR::concurrent_files(2)))

Lors de la normalisation, le script place automatiquement les fichiers dans des sous-dossiers correspondant aux “blocs”. Le suffixe “norm” est ajouté au nom d’origine. Au cas où des sous-dossiers sont déjà présents avec des fichiers, la fonction propose uniquement la normalisation des fichiers non présents localement. Les fichiers normalisés ne sont pas au format copc, un fichier d’indexation au format lax est créé automatiquement pour rendre les requêtes spatiales ultérieures plus rapides.

# sous-dossiers
list.dirs(dossier_norm, full.names = FALSE, recursive = FALSE)
## [1] "NUALHD_1-0__LAZ_LAMB93_OM_2025-03-25"
## [2] "NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25"

Dans chaque sous-dossier de bloc, un catalogue des dalles LiDAR normalisées présentes est créé sous la format d’un objet au format LAS-catalog enregistré dans une archive R nommée “cata.rda”. Un catalogue de l’ensemble des dalles est créé dans une archive de même nom mais enregistré à la racine du répertoire de destination (dossier_norm). Ces catalogues sont mis à jour lorsque des dalles supplémentaires sont calculées.

# fichiers dans le sous-dossier du bloc "PM"
dir(paste0(dossier_norm,
           "NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25"))
## [1] "cata.rda"                                   
## [2] "LHD_FXX_0886_6436_PTS_LAMB93_IGN69.copc.laz"

Le code ci-dessous charge le catalogue, l’affiche avec une couleur par bloc, et ajoute la région d’intérêt avec un contour rouge. Les dalles adjacentes à celles qui intersectent la région d’intérêt n’ont pas été calculées.

# mettre a jour les catalogues si besoin
# for (i in list.dirs(dossier_norm, recursive = FALSE, full.names = FALSE))
#   cata <- lidarHD::catalog_bloc(dossier_norm, i)
# cata <- lidarHD::catalog_folder(dossier_norm)
# charger le catalog
load(file = paste0(dossier_norm, "cata.rda"))
lidR::plot(cata,
           col = ifelse(cata$bloc == "NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25",
                        "blue", "green"))
plot(roi, add = TRUE, border = "red")

La fonction check_folder() vérifie que les fichiers et dossiers sont organisés conformément au tableau d’assemblage. Il est utile de faire cette vérification régulièrement au cas où des fichiers locaux seraient rendus obsolètes par des mises à jour sur le serveur IGN.

lidarHD::check_folder(dossier_norm, ta = ta_lidarHD)

Calcul de modèles numériques

L’IGN fournit déjà des modèles numériques d’élévation au pas de 50 cm. L’inconvénient est qu’ils sont alignés sur un grille différente des orthophotos et des dalles LiDARHD. Il est possible de ré-échantillonner ces modèles, ou de les calculer à partir des nuages de points. L’exemple est donné ci-dessous pour les MNT.

Calcul de MNT

La fonction compute_dtm() calcul le MNT pour les dalles contenues dans le dossier laz_folder qui intersectent une région d’intérêt roi et qui ne sont pas déjà présentes dans le répertoire de destination dtm_folder. Si il n’est pas spécifié, le répertoire dtm_folder est créé dans le même dossier que laz_foldersur le format “DTM_resm” avec res la résolution cible res. La région d’intérêt peut être un objet spatial au format sf, ou un vecteur comportant les noms des blocs à traiter. Il faut s’assurer que les dalles adjacentes à celles qui vont être caclulées sont également disponibles dans le dossier laz_folder, afin d’éviter les effets de bords. C’est le cas si elles ont été téléchargées avec la fonction download_files().

Le MNT est calculé avec la fonction lidR::rasterize_terrain() et l’algorithme lidR::tin(). Les classes de points utilisées sont par défaut 2 et 9 (sol et eau).

Les MNT sont enregistrés dans des fichiers selon le format “dtmresm_nom.tif”, selon la même arborescence de blocs que les fichiers LAZ téléchargés depuis le site IGN.

Le paramètre prompt est réglé par défaut pour que l’affichage des dalles à calculer et le calcul soient confirmés par l’utilisateur (recommandé).

Le calcul est parallélisé, le paramètre ncores permet de spécifier combien de process sont lancés simultanément, le facteur limitant sera plutôt la mémoire vive de la machine.

Le script affiche le nombre de dalles selon les croisements de catégories suivantes :

  • intersecte ou pas la région d’intérêt (dedans / hors_zone),
  • déjà calculé ou pas (fait / a_faire).

Seules les dalles a faire seront calculées.

# calculer les MNT qui intersectent la region d'interet 
compute_dtm(dossier_classe, res = 1, roi = roi, prompt = TRUE, ncores = 2L)
# alternative basée sur les fonctions lasR, 1 fichier à la fois avec 2 coeurs
process_files(dossier_classe, res = 1, roi = roi, prompt = TRUE,
              lasR_options = list(ncores = lasR::concurrent_points(2)))

Lors du calcul, le script place automatiquement les fichiers dans des sous-dossiers correspondant aux “blocs”. Au cas où des sous-dossiers sont déjà présents avec des fichiers, la fonction propose uniquement le calcul des fichiers non présents localement. Un raster virtuel consitué des rasters qui viennent d’être calculés est créé pour chaque bloc dans le fichier “rasterize_terrain.vrt”.

# fichiers dans le sous-dossier du bloc "PM"
dir(paste0(dossier, "/DTM_1m/NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25"))
## [1] "dtm1m_LHD_FXX_0886_6436_PTS_LAMB93_IGN69.copc.tif"

Ajout des produits dérivés dans le tableau d’assemblage

Afin de recenser dans un fichier spatial unique l’emprise des dalles LiDAR HD disponibles, le chemin local vers les fichiers téléchargés ainsi que celui vers les produits dérivés (nuages normalisés et modèles de terrain), la fonction add_attributes_ta() peut être utilisée. À partir du tableau d’assemblage des dalles LiDARHD, elle ajoute des attributs qui indiquent le chemin vers les fichiers dérivés :

  • dossier “LAZ_classe” : fichiers LiDAR HD téléchargés,
  • dossier “LAZ_norm” : fichiers normalisés,
  • dossier “DTM_1m” : modèles numériques de terrain à 1 m de résolution,
  • dossier “DTM_5m” : modèles numériques de terrain à 5 m de résolution,
  • dossier “CHM_1m” : modèles numériques de canopée à 1 m de résolution.

Au cas où d’autres produits sont calculés, ou avec des noms de répertoire différents, la fonction pourra être modifiée pour le prendre en compte. L’option save permet de sauvegarder le tableau modifié dans un fichier dont l’argument dossier serait le chemin d’accès.

ta_lidarHD <- lidarHD::add_attributes_ta(ta_lidarHD, dossier, save = FALSE)
# produced dtm at 1 m
ta_lidarHD$DTM_1m[!is.na(ta_lidarHD$DTM_1m)]
## [1] "/DTM_1m/NUALHD_1-0__LAZ_LAMB93_OM_2025-03-25/dtm1m_LHD_FXX_0885_6436_PTS_LAMB93_IGN69.copc.tif"
## [2] "/DTM_1m/NUALHD_1-0__LAZ_LAMB93_PM_2025-03-25/dtm1m_LHD_FXX_0886_6436_PTS_LAMB93_IGN69.copc.tif"