workflow.normalisation2021.R

# workflow associated with package lidaRtRee
# Copyright Irstea
# Author(s): Jean-Matthieu Monnet
# Licence: CeCill 2.1
###################################
# NORMALISATION D'UN ENSEMBLE DE FICHIER LAS / LAZ AVEC lidR, méthode sans discrétisation (triangulation de Delaunay des points sol)
#
library(foreach)
# number of cores to use
doParallel::registerDoParallel(cores=4)
# 
# faire le catalog des fichiers las : fonctionne avec des dalles rectangulaires alignées sur les coordonnées
# cata <- lidR::catalog("~/test_normalisation/")
cata <- lidR::catalog("~/Bureau/NCaledonie/ncaledonie/LAS_retiled/")

sp::plot(cata)
#
# au cas où les dalles ne soient pas bien alignees: refaire le dallage
reshape_tiles <- FALSE
tile.size <- 100
reshape_folder <- "~/test_normalisation/reshape/"
reshape_prefix <- "CM2015_"
if (reshape_tiles)
{
  # reshape tiles
  cataR <- catalog_reshape(cata, tile.size, reshape_folder, prefix=reshape_prefix, ext="laz")
  # replace catalog
  cata <- cataR
}
#
# taille du tampon (m) pour eviter les effets de bord lors de la triangulation
bsize <-5
#
catad <- cata@data
#################
# calcul parallelise du nuage normalise
# ecriture dans un fichier suffixe par "norm"
# renvoi un data.frame avec nom de la dalles, nb points avant, nb points après (perte possible en bordure d'acquisition si les points ne sont pas dans le convex hull des points sol)
stats <- foreach::foreach (i=1:nrow(catad), .combine=rbind) %dopar%
{
  # load tile extent plus buffer
  a <- try(lidR::clip_rectangle(cata, catad$Min.X[i]-bsize, catad$Min.Y[i]-bsize, catad$Max.X[i]+bsize, catad$Max.Y[i]+bsize)) # using only one core to load files because the for loop is already parallelised
  if (class(a)=="try-error")
  {
    # print("no data")
    temp <- data.frame(id=catad$filename[i], nb.init=NA,nb.norm=NA)
  } else {
    if (!is.null(a))
    {
      # normalize
      a <- lidR::normalize_height(a, lidR::tin())
      # remove buffer points
      a <- lidR::filter_poi(a, X >= catad$Min.X[i] & Y >= catad$Min.Y[i] & X < catad$Max.X[i] & Y < catad$Max.Y[i])
      if (nrow(a@data)>0)
      {
        # write normalized tile
        # convert to las if relevant
        # new name
        filename <- strsplit(catad$filename[i], "/")
        filename <- paste0(paste(filename[[1]][-length(filename[[1]])], collapse="/"), "/norm_",gsub(".las",".laz", filename[[1]][length(filename[[1]])]))
        lidR::writeLAS(a, file=filename)
        temp <- data.frame(id=catad$filename[i], nb.init=catad$Number.of.point.records[i], nb.norm=nrow(a@data))
      } else {
        temp <- data.frame(id=catad$filename[i], nb.init=0,nb.norm=0)
      }
    } else {
      temp <- data.frame(id=catad$filename[i], nb.init=NA,nb.norm=NA)
    }
  }
  return(temp)