Diferenças

Aqui você vê as diferenças entre duas revisões dessa página.

--- lidar:projetos:rcode_alsmodelo [2023/11/26 12:34] – lcer
+++ lidar:projetos:rcode_alsmodelo [2023/11/27 10:41] (atual) – removida lcer
@@ Linha 1: / Linha 1: @@
-[[ lidar:projetos |{{ :lidar:projetos:retornar.png?nolink&25x25 }}]]
-<code>
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Inventário com dados LiDAR na Fazenda Modelo ~~~~~~~~~~~~~~~~~~~~~~~~
-#
-# Autor: Luiz Carlos Estraviz Rodriguez
-#        Departamento de Ciências Florestais
-#        ESALQ/USP - 25/Nov/2023
-#
-# Inventário florestal da Fazenda Modelo
-#   - download dos dados mantidos em um repositório github público
-#      - shape files dos talhões florestais
-#      - LiDAR multitemporal (2013 e 2014)
-#   - sugestão de pasta para armazenamento local:
-#        C:/LiDAR/PRJ_Modelo/
-#
-# Linguagem de programação:
-#       R (v 4.3)
-#       pacote lidR* (v 4.0.3 March 11th, 2023) Jean Romain Roussel
-#
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-rm(list=ls(all=TRUE))                                   # Limpa memória
-gc()
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# 1. Leitura e organização de dados
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-if(!require(tidyverse))  # Para melhor manipulação de dados e funções
-  install.packages("tidyverse")
-library(tidyverse)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Download do shape da Fazenda Modelo (2 layers: talhoes e parcelas)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-gitOnde <- "https://github.com/FlorestaR/dados/blob/main/5_LIDARF/Modelo/SHAPES"
-gitNome <- "fazmodelo.zip"
-gitArqv <- file.path(gitOnde, gitNome) %>% paste0("?raw=true")
-tmpd <- tempdir(check = TRUE)                    # diretório temporário
-zipf <- file.path(tmpd, "shapes.zip")              # arquivo temporário
-options(timeout=1000) # Reset timeout oferecendo mais tempo de download
-if(!file.exists(zipf))  # garante download de dados binários (wb)
-  download.file(gitArqv, mode="wb", destfile = zipf)
-unzip(zipf, exdir = tmpd)     # shape é unziped no diretório temporário
-unlink(zipf)                                  # deleta o arquivo zipado
-if(!require(sf))                           # Para manipulação de shapes
-  install.packages("sf")
-library(sf)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Lê atributos dos talhoes
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-shpArq <- paste0(tmpd, "/Modelo_talhoes.shp")       # shape com talhões
-talhoesComGeo <- read_sf(shpArq)                    # completo com geom
-talhoesSemGeo <- tibble(sf::st_drop_geometry(talhoesComGeo))  # s/ geom
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Lê atributos das parcelas
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-shpArq <- paste0(tmpd, "/Modelo_parcelas.shp")     # shape com parcelas
-parcelasComGeo <- read_sf(shpArq)               # completo com geom
-parcelasSemGeo <- tibble(sf::st_drop_geometry(parcelasComGeo)) # s/ geom
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Cria coluna IDINV na tabela "talhões", extraída da tabela "parcelas"
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-talhoes <-
-  parcelasSemGeo %>%
-  group_by(SUBTALHAO) %>%
-  summarise(IDINV = unique(IDINV)) %>%
-  left_join(talhoesSemGeo) %>%
-  select("SUBTALHAO", "IDINV", "AREA") %>%
-  arrange(SUBTALHAO) %>% as.data.frame
-talhoesComGeo <- inner_join(talhoesComGeo, talhoes, by="SUBTALHAO") %>%
-  select("OBJECTID", "SUBTALHAO", "IDINV", "geometry")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Reorganiza colunas da tabela "parcelas"
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-parcelas <-
-  parcelasSemGeo %>%
-  select(SUBTALHAO, CHAVE2, DATAREALIZ, IDINV, AREAPARCEL, MHDOM, VTCC, AB) %>%
-  arrange(SUBTALHAO) %>% as.data.frame
-# Cria diretórios e pastas para onde os tiles LiDAR serão copiados
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-prjNome <- 'PRJ_Modelo'
-dirNome <- paste0('C:/LiDAR/', prjNome)
-dir.create(dirNome, showWarnings = F)
-# Salva e imprime dados reorganizados
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dirNome <- paste0(dirNome, '/DADOS/')
-dir.create(dirNome, showWarnings = F)
-# Leitura da mais atual versão rio para facilmente ler e salvar Excel
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-if (!require("remotes")){install.packages("remotes")}
-if (!require(rio))      {remotes::install_github("gesistsa/rio")}
-library(rio)
-# Salva dados em planilha Excel
-arqNome <- paste0(dirNome, prjNome, '.xlsx')
-export(list(talhoes = talhoes,
-            parcelas  = parcelas),
-       file = arqNome)
-# Mostra tabelas no Viewer
-# Use o "Export" do Viewer para salvar em diferentes formatos
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-if(!require(kableExtra))
-  install.packages("kableExtra")
-library(kableExtra)
-# Mostra talhões
-AreaTotal <- talhoes$AREA %>% sum
-NotaDeRodape <- paste0(": ", AreaTotal)
-talhoes %>%
-  kbl(caption = "Talhões da Fazenda Modelo", align = "r") %>%
-  kable_classic(full_width = F) %>%
-  footnote(general = NotaDeRodape,
-           general_title = "Área total",
-           footnote_as_chunk = T)
-# Mostra parcelas
-parcelas %>%
-  kbl(caption = "Parcelas da Fazenda Modelo", align = "r") %>%
-  kable_classic(full_width = F)
-# Cria diretório para salvar mapas, figuras e gráficos
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dirNome <- paste0('C:/LiDAR/', prjNome)
-  dir.create(dirNome, showWarnings = F)
-dirNome <- paste0(dirNome, '/GRAPH/')
-  dir.create(dirNome, showWarnings = F)
-# Cria mapa dos talhões com localização das parcelas (EPSG: 31983)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-arqNome <- paste0(dirNome, prjNome, '.png')
-png(arqNome, 30, 20, 'cm', res = 200)                # abre "impressão"
-ggplot() +        #         plot dos talhões e parcelas (cor por idade)
-  geom_sf(data = talhoesComGeo, colour = "black", fill="white") +
-  geom_sf(data = parcelasComGeo, aes(fill = factor(IDINV))) +
-  scale_fill_discrete(name = "Idade") +
-  guides(fill = guide_legend(reverse=F)) +
-  coord_sf(datum=st_crs(31983)) +        # Especifica sistema de coord.
-  scale_y_continuous(breaks = seq(from=7356500,to=7359000, by=200)) +
-  scale_x_continuous(breaks = seq(from=206200, to=207600,  by=200))
-dev.off()                                # fecha "impressão" do aquivo
-# ----
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# 2. Amostragem Casual Simples (ACS)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-# Função para cálculo do Intervalo de Confiança
-calcCI = function(err, n, alpha=.05){
-  return(
-    qt(1 - alpha/2, n-1) * err #/ sqrt(n)
-  )
-}
-# Função para cálculo dos parâmetros de estimação e inferência
-calcPars = function(df, N, alpha=.05){
-  means = apply(df, 2, mean)
-  vars   = apply(df, 2, function(y){
-    (var(y) / length(y)) * ((N - length(y)) / N)
-  })
-  err = sqrt(vars)
-  cis = calcCI(err, nrow(df), alpha)
-  err_pc = 100*cis/means
-  out = rbind(means, vars, err, cis, err_pc, nrow(df)) %>% as.data.frame
-  row.names(out) = c('media', 'var', 'dp', 'ic', 'erro', 'n')
-  names(out) = names(df)
-  return(out)
-}
-# Cria lista com variáveis de interesse
-variaveisDeInteresse <- c("MHDOM", "VTCC", "AB")
-# Cria dataframe com valores observados das variáveis de interesse
-df <- parcelas %>% select(all_of(variaveisDeInteresse))
-# Calcula o número máximo de unidades amostrais possíveis (N)
-parcelaAreaMedia <- mean(parcelas$AREAPARCEL) / 10000          # por ha
-N_ACS <- round(AreaTotal / parcelaAreaMedia , 0)    # N
-# Resultados do inventário (estimação + inferência)
-#             pelo método Amostragem Casual Simples (ACS)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-paramEstatisticosACS  <- df %>% calcPars(N_ACS)
-# Função para cálculo da intensidade amostral recomendada
-# (ha por parcela) da ACS que garante ~10% de erro (altere se desejar)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-erro <- 10
-tamanhoIdealACS = function(y, N, errDesired=erro/100, alpha=erro/100){
-  B  = errDesired * mean(y)
-  qt = qt(1 - alpha/2, length(y)-1)
-  n  = N*var(y)*qt^2 / (N * B^2 + qt^2 * var(y))
-  return(n)
-}
-# Intensidade amostral usada
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-IAS <- paste0(" usada: 1 parc/",
-              round(AreaTotal/paramEstatisticosACS$VTCC[6],0), " ha.")
-# Intensidade amostral que seria necessária para ACS
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-IA <- round(AreaTotal/tamanhoIdealACS(df$VTCC,N_ACS),0)
-erroTexto <- as.character(erro)
-IAI <- paste0("\n Necessárias p/ erro de ", erroTexto,
-              "%: 1 parc/", IA, " ha.")
-# Mostra resultados da ACS no Viewer
-# Use o "Export" do Viewer para salvar em diferentes formatos
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-NotaDeRodape <- paste0(IAS, IAI)
-paramEstatisticosACS %>%
-  kbl(caption = "Amostragem Casual Simples (ACS)", align = "r") %>%
-  kable_classic(full_width = F) %>%
-  footnote(general = NotaDeRodape,
-           general_title = "Intensidade amostral",
-           footnote_as_chunk = T)
-# ----
-# *********************************************************************
-# 3. Amostragem Casual Estratificada (ACE)
-# **************************************************************** ----
-# Escolhe variável de estratificação
-grupos = c("IDINV")
-# Box-plot da variável de interesse por nível de estratificação
-par(mfrow=c( length(grupos), length(variaveisDeInteresse) ))
-attach(parcelas)
-for(i in grupos){
-  for(j in variaveisDeInteresse){
-    plot( get(i) %>% factor , get(j), main = j, sub=i)
-  }
-}
-detach(parcelas)
-# Estimativa da população para amostragem estratificada
-popFromStrata = function(factorStrataList){
-  popEstimates = foreach(i = factorStrataList, .combine = 'c') %do% {
-    gpMeans = lapply(i, function(x) x[1,,drop=F]) %>% do.call(what = rbind)
-    gpVars  = lapply(i, function(x) x[2,,drop=F]) %>% do.call(what = rbind)
-    cols = 1:(ncol(gpMeans)-4)
-    popMean   = apply(gpMeans[,cols], 2,
-                      function(x) sum(x*gpMeans$N) ) /N_ACS
-    popVar    = apply(gpVars[,cols], 2,
-                      function(x) sum( x * (gpVars$N/N_ACS)^2 ) )
-    popStdErr = sqrt(popVar)
-    popCI     = calcCI(popStdErr, sum(gpMeans$n))
-    popPars = data.frame(
-      media = popMean,
-      var   = popVar,
-      dp    = popStdErr,
-      ic    = popCI,
-      erro  = 100 * popCI / popMean,
-      n     = sum( sapply(i, function(x) mean(x$n)) )
-    )
-    return(list(popPars))
-  }
-  names(popEstimates) = names(factorStrataList)
-  return(popEstimates)
-}
-if(!require(foreach))                         # Para loops inteligentes
-  install.packages("foreach")
-library(foreach)
-# Cria vetor de (True ou False) para marcar dados que estejam dentro
-#      do intervalo de idades de interesse (entre 2 e 6)
-rightAges = talhoes$IDINV > 2 & talhoes$IDINV < 6
-# Resultados do inventário (estimação + inferência)
-#             pelo método Amostragem Casual Estratificada (ACE)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-paramEstatisticosACE = foreach(
-  grp = grupos, .combine = 'rbind') %:%
-  foreach(
-    niv = parcelas[,grp] %>% unique %>%
-      as.character %>% sort, .combine = 'rbind'
-  ) %do% {
-    inGroup = talhoes[,grp] == niv
-    popSize = round(
-      sum(talhoes[rightAges & !is.na(inGroup) & inGroup,]$AREA) /
-        parcelaAreaMedia)
-    inGroup = parcelas[,grp] == niv
-    tempPars = parcelas[inGroup, variaveisDeInteresse, drop=F]
-    inventory = calcPars(tempPars, popSize)
-    inventory$grupo = grp
-    inventory$nivel = niv
-    inventory$n     = nrow(tempPars)
-    inventory$N     = popSize
-    return(inventory)
-  }
-paramEstatisticosACE %<>% base::split(f=paramEstatisticosACE$grupo) %>%
-  lapply(function(x) split(x, x$nivel))
-globalparamEstatisticosACE = popFromStrata(paramEstatisticosACE)
-# Função para cálculo da intensidade amostral recomendada
-# (ha por parcela) da ACE que garante ~10% de erro (altere se desejar)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-erro <- 10
-tamanhoIdealACE = function(y, g, Nh, errDesired=erro/100){
-  vars = by(y, g, stats::var)
-  Wh   = by(y, g, length) / length(y)
-  Nh   = Nh[ names(Nh) %in% g ]
-  B    = errDesired * mean(y)
-  n    = sum( Nh^2 * vars / Wh ) /
-         ((B^2 * (sum(Nh)^2))/4 + sum(Nh * vars))
-  return(n)
-}
-# Calcula o número possível de unidades amostrais (N) por estrato
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-parcPorEstrato <- round(                                           # Nh
-  by(st_area(talhoesComGeo), talhoes$IDINV, sum) /
-    (mean(parcelas$AREAPARCEL)), 0)
-IA <- round(AreaTotal/tamanhoIdealACE(y = parcelas$VTCC,
-                                      g = parcelas$IDINV,
-                                      Nh = parcPorEstrato), 0)
-# Intensidade amostral usada por estrato
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-IAE <- round(by(st_area(talhoesComGeo)/10000, talhoes$IDINV, sum) /
-               (parcelas %>% count(IDINV))[2])[1]
-IAS <- paste0(" usada: 1 parc / ", as.character(IAE) , " ha.")
-# Intensidade amostral que seria necessária para ACE
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-erroTexto <- as.character(erro)
-IAI <- paste0("\n Necessárias p/ erro de ", erroTexto,
-              "%: 1 parc / ", IA, " ha.")
-NotaDeRodape <- paste0(IAS, IAI)
-globalparamEstatisticosACE[[grupos]] %>% t %>%  # transpõe o data frame
-  kbl(caption = paste0("Amostragem Casual Estratificada (ACE)"),
-      align = "r") %>%
-  kable_classic(full_width = F) %>%
-  footnote(general = NotaDeRodape,
-           general_title = "Intensidade amostral",
-           footnote_as_chunk = T)
-# Gráfico de número de amostras necessárias por erro amostral aceitável
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-prjNome <- 'PRJ_Modelo'
-dirNome <- paste0('C:/LiDAR/', prjNome)
-dir.create(dirNome, showWarnings = F)
-dirNome <- paste0(dirNome, '/GRAPH/')
-dir.create(dirNome, showWarnings = F)
-arqNome <- paste0(dirNome, prjNome, '_IntAmostral_ACSvsACE.png')
-erro = 1:20
-scsPlotn <- tamanhoIdealACS(parcelas$VTCC, N_ACS,
-                            errDesired = erro/100)
-cssPlotn <- tamanhoIdealACE(parcelas$VTCC,  parcelas$IDINV,
-                            parcPorEstrato, errDesired = erro/100)
-# dsPlotn = dubleSamplePlotNumber(plotMetrics$PC_VOL_TOT,
-#                                 plotMetrics$avgCrownHeight,
-#                                 cropMetrics$avgCrownHeight, 300, erro/100)
-png( arqNome, 30, 20, 'cm', res = 200)           # abre "impressão" PNG
-# jpeg(arqNome, 30, 20, 'cm', res = 200)           # abre "impressão" JPG
-plot( scsPlotn ~ erro, type='l', col='red', lwd=2,
-      ylim=c(0,200), ylab='Número de parcelas', xlab='Erro amostral (%)',
-      axes=F)
-lines(erro, cssPlotn, col='green', lwd=2)
-# lines(erro, dsPlotn, col='blue', lwd=2)
-box()
-axis(1)
-axis(2, at = c(25, 50, 75, 100, 125, 150, 175, 200))
-abline(v=10, col='black', lwd=2, lty=2)
-lines(c(0,15), rep(100, 2), col='black', lty=2, lwd=2)
-# lines(c(0,5), rep(dsPlotn[5], 2), col='blue', lty=2, lwd=2)
-# lines(c(0,5), rep(cssPlotn[5], 2), col='green', lty=2, lwd=2)
-legend('topright', col=rainbow(3), lwd=2,
-       legend = c('ACS', 'ACE'))
-# legend = c('ACS', 'ACE', 'AD'))
-dev.off()                                           # fecha "impressão"
-# Salva resultados em planilha Excel
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-prjNome <- 'PRJ_Modelo'
-dirNome <- paste0('C:/LiDAR/', prjNome)
-dir.create(dirNome, showWarnings = F)
-dirNome <- paste0(dirNome, '/RSLTS/')
-dir.create(dirNome, showWarnings = F)
-arqNome <- paste0(dirNome, prjNome, '_Results.xlsx')
-export(list(ACS = paramEstatisticosACS,
-            ACE37 = paramEstatisticosACE$IDINV$'3.7',
-            ACE52 = paramEstatisticosACE$IDINV$'5.2'),
-       file = arqNome)
-# save.image(file='Modelo_inventario.RData')       # faz backup dos dados
-# load('Modelo_inventario.RData')                       # recupera backup
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Este script foi atualizado em Novembro/2023 depois de instalar,
-# na seguinte sequência, as mais recentes versões do pacote lidR*:
-#
-#   devtools::install_github("Jean-Romain/rlas", dependencies=TRUE)
-#   devtools::install_github("Jean-Romain/lidR")
-#
-#   * lidR latest development version
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#  ----
-</code>
-----
-Referências sobre o pacote //lidR// para processamento de dados LiDAR:
-  * [[https://r-lidar.github.io/lidRbook/ | The lidR package (eBook)]]
-  * [[https://hal.inrae.fr/hal-02972284/file/roussel_2020.pdf | Roussel et al.(2020) lidR: An R package for analysis of ALS data]]
-  * [[https://www.rdocumentation.org/packages/lidR/versions/4.0.1 | RDocumentation (lidR)]]
-  * [[https://cran.rstudio.com/web/packages/lidR/vignettes/ | Vignettes]]