Passer au contenu

Téléchargement des données

Tableau d’assemblage des fichiers LiDARHD IGN

La fonction load_classified_ta() télécharge le tableau d’assemblage des fichiers LiDARHD classés à partir d’un serveur WFS de la géoplateforme IGN (couche “IGNF_LIDAR-HD_TA:nuage-dalle”). Le téchargement est assez long car il existe au 20 décembre 2024 plus de 337000 dalles, qui sont chargées par plusieurs requêtes successives car le nombre maximal d’entités retournées est limité à 5000. 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.

# charger depuis le serveur WFS (non exécuté)
ta_lidarHD <- load_classified_ta()
# charger depuis un fichier presente dans un dossier
ta_lidarHD <- load_classified_ta(url = NULL, save_dir = "./data/sig/")
## [1] "Reading file from 2023-12-18"
## Reading layer `TA_diff_pkk_lidarhd_classe_OM_PM' from data source 
##   `/home/jean-matthieu_monnet/R/lidar/lidarHD/vignettes/data/sig/ta_231218.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 5000 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 836000 ymin: 6433000 xmax: 936000 ymax: 6483000
## Projected CRS: RGF93 v1 / Lambert-93

Le tableau d’assemblage contient deux champs :

  • nom_pkk avec le nom du fichier sur le modèle “LHD_FXX_0377_6228_PTS_C_LAMB93_IGN69.copc.laz”,
    • LHD: campagne LiDAR HD,
    • FXX: ?
    • 0377 et 6228: les deux nombres à quatre chiffres sont les coordonnées du coin haut et gauche de la dalle, en milliers de kilomètres,
    • PTS_C: nuage de points classés,
    • LAMB93_IGN69 : système de coordonnées Lambert 93, altitude IGN69,
    • copc : les fichiers sont au format copc (spatialement indexés, ils sont légèrement moins compressés que des fichiers laz classiques, mais il peuvent être requêtés plus efficacement),
  • url_telech est le lien de téléchargement du fichier, l’arborescence contient en dernier le nom du fichier, et en avant dernier le nom du bloc. La fonction load_classified_ta() ajoute le nom du bloc dans la table d’attributs.

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.

# sous-dossiers
list.dirs(dossier_classe, full.names = FALSE, recursive = FALSE)
## [1] "OM" "PM"

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, "PM"))
## [1] "cata.rda"                                     
## [2] "LHD_FXX_0886_6435_PTS_O_LAMB93_IGN69.copc.laz"
## [3] "LHD_FXX_0886_6436_PTS_O_LAMB93_IGN69.copc.laz"
## [4] "LHD_FXX_0886_6437_PTS_O_LAMB93_IGN69.copc.laz"
## [5] "LHD_FXX_0887_6435_PTS_O_LAMB93_IGN69.copc.laz"
## [6] "LHD_FXX_0887_6436_PTS_O_LAMB93_IGN69.copc.laz"
## [7] "LHD_FXX_0887_6437_PTS_O_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 == "PM", "blue", "green"))
plot(roi, add = TRUE, border = "red")

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.

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 à télécharger. 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 téléchargées.

# telecharger les fichiers qui intersectent la region d'interet 
normalize_files(dossier_classe, dossier_norm, roi = roi, prompt = TRUE, ncores = 2L)

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] "OM" "PM"

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, "PM"))
## [1] "cata.rda"                                          
## [2] "LHD_FXX_0886_6436_PTS_O_LAMB93_IGN69.copc.norm.lax"
## [3] "LHD_FXX_0886_6436_PTS_O_LAMB93_IGN69.copc.norm.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 == "PM", "blue", "green"))
plot(roi, add = TRUE, border = "red")

# Calcul de modèles numériques de terrain —-

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à téléchargé 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)

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/PM"))
## [1] "dtm1m_LHD_FXX_0886_6436_PTS_O_LAMB93_IGN69.copc.tif"
## [2] "rasterize_terrain.vrt"

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èle 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é pour le prendre en compte. L’option save permet de sauvegarder le tableau modifié dans un fichier dont l’argument ta_lidarHD 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/OM/dtm1m_LHD_FXX_0885_6436_PTS_C_LAMB93_IGN69.copc.tif"
## [2] "/DTM_1m/PM/dtm1m_LHD_FXX_0886_6436_PTS_O_LAMB93_IGN69.copc.tif"