Tutoriel lidarHD
Jean-Matthieu Monnet
lidarHD.Rmd
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
et6228
: 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 fichierslaz
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 fonctionload_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.
## [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.
## [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_folder
sur 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”.
## [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"