Diferenças

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

--- lidar:projetos:rcode_invcomlidar [2023/11/27 10:38] – criada lcer
+++ lidar:projetos:rcode_invcomlidar [2024/03/23 20:17] (atual) – edição externa 127.0.0.1
@@ Linha 1: / Linha 1: @@
 [[ lidar:projetos |{{ :lidar:projetos:retornar.png?nolink&25x25 }}]]
 <code>
@@ Linha 8: / Linha 7: @@
 # Autor: Luiz Carlos Estraviz Rodriguez
 #        Departamento de Ciências Florestais
-#        ESALQ/USP - 25/Nov/2023
+#        ESALQ/USP - 02/Dez/2023
 #
 # Inventário florestal da Fazenda Modelo
@@ Linha 14: / Linha 13: @@
 #      - shape files dos talhões florestais
 #      - LiDAR multitemporal (2013 e 2014)
-#   - sugestão de pasta para armazenamento local:
+#   - sugestão de pastas locais
-#        C:/LiDAR/PRJ_Modelo/
+#        C:/LiDAR/PRJ_Modelo/NUVENS/     como pasta para os dados LidAR
+#        C:/GitRepo/PRJ_Modelo/     sendo GitRepo a pasta onde ficam os
+#                                      projetos sincronizados no github
 #
 # Linguagem de programação:
 #       R (v 4.3)
-#       pacote lidR* (v 4.0.3 March 11th, 2023) Jean Romain Roussel
+#       RStudio (v. 2023.09.1 Build 494)
+#       Pacote lidR* (v 4.0.3 March 11th, 2023) Jean Romain Roussel
 #
+# Bibliografia complementar:
+#  https://github.com/r-lidar/lidR/wiki/Area-based-approach-from-A-to-Z
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 rm(list=ls(all=TRUE))                                   # Limpa memória
 gc()
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# 1. Leitura e organização inicial de dados
+# 1. Carrega pacotes, organiza pastas e define funções locais
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-if(!require(tidyverse))  # Para melhor manipulação de dados e funções
+if(!require(tidyverse))    # Para melhor manipulação de dados e funções
   install.packages("tidyverse")
 library(tidyverse)
-# Define diretórios e pastas de trabalho
+if(!require(sf))            # Para manipulação de shapes e outros tipos
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  install.packages("sf")                      # de dados espacializados
-prjNam <- "PRJ_Modelo"
+library(sf)
-wrkDir <- str_c('C:/LiDAR/', prjNam)
-  dir.create(wrkDir, showWarnings = F)
-  setwd(str_c(wrkDir))
-lidDir <- str_c('C:/LiDAR/', prjNam, '/CLOUDS/')
-  dir.create(lidDir, showWarnings = F)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+if(!require(tidyterra))           # Package para processar mapas raster
-# Define os nomes dos tiles LiDAR multitemporais que serão lidos do
+  install.packages("tidyterra")
-# repositório, os respectivos anos e o URL completo do github
+library(tidyterra)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-gitOnde <-"https://github.com/FlorestaR/dados/blob/main/5_LIDARF/Modelo/CLOUDS/"
-gitNome <-c("L0004-C0005.laz", "L0004-C0006.laz",
-            "L0005-C0005.laz", "L0005-C0006.laz",
-            "L0006-C0005.laz", "L0006-C0006.laz")
-# anoData <- c("A13", "A14")
-anoData <- c("A13")                  # Lê apenas os dados LiDAR de 2013
-# Faz o download dos tiles
+if(!require(terra))           # Package para processar mapas raster
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  install.packages("terra")
-options(timeout=1000) # Reset timeout oferecendo mais tempo de download
+library(terra)
-for (ano in anoData) {
-  dirLAZ <- paste0(lidDir, ano)
-  dir.create(dirLAZ, showWarnings = T)
-  for (nome in gitNome){
-    gitFile <- paste0(gitOnde, ano, "/", ano, nome, "?raw=true")
-    localFl <- paste0(dirLAZ, "/", nome)
-    if(!file.exists(localFl))        # garante download binário (wb)
-      download.file(gitFile, mode="wb", destfile = localFl)
-  }
-}
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+if(!require(stars))                  # Package que permite rasterização
-# Download dos shapes da Fazenda Modelo (2 layers: talhoes e parcelas)
+  install.packages("stars")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+library(stars)
-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
+if(!require(tools))           # Package para manipular nomes de aquivos
-zipf <- file.path(tmpd, "shapes.zip")              # arquivo temporário
+  install.packages("tools")
+library(tools)
-options(timeout=1000) # Reset timeout oferecendo mais tempo de download
+if(!require(RColorBrewer))          # Package com mais paletas de cores
-if(!file.exists(zipf))  # garante download de dados binários (wb)
+  install.packages("RColorBrewer")
-  download.file(gitArqv, mode="wb", destfile = zipf)
+library(RColorBrewer)
-shpDir <- str_c(wrkDir, '/SHAPES')           # diretório para os shapes
+if(!require(progress))          # Package com mais paletas de cores
-  dir.create(shpDir, showWarnings = F)
+  install.packages("progress")
-unzip(zipf, exdir = shpDir)                           # shape é unziped
+library(progress)
-unlink(zipf)                                  # deleta o arquivo zipado
-# ----
+if(!require(reshape2))              # Usado na função que gera gráficos
+  install.packages("reshape2")         # mais caprichados de correlação
+library(reshape2)
+if(!require(mapview))             # Package para incluir mapas de fundo
+  install.packages("mapview")
+library(mapview)
-# *********************************************************************
-# 2: A Amostragem Dupla (AD) com variável auxiliar LiDAR começa aqui.
-# **************************************************************** ----
 # devtools::install_github("Jean-Romain/rlas", dependencies=TRUE)
 # devtools::install_github("Jean-Romain/lidR")
 if(!require(lidR))                     # PARA MANPULAÇÃO DE DADOS LiDAR
-  install.packages("lidR") # Se preferir instalar a mais recente versão
+  install.packages("lidR") # Se preferir instalar a versão mais recente
-                           # e tiver o package devtools instalado, rode
+library(lidR)              # e tiver o package devtools instalado, rode
-                           # alternativamente as duas linhas acima.
+                                # alternativamente as duas linhas acima
-library(lidR)
-if(!require(future))        # Package que permite ao lidR rodar algumas
+if(!require(RCSF))    #        Pacote complementar do lidR e necessário
-  install.packages("future")           # com processmamento em paralelo
+  install.packages("RCSF")           # rodar a função classify_ground()
+library(RCSF)
+if(!require(future))         # Package que permite ao lidR rodar usando
+  install.packages("future")          # processamento paralelo de dados
 library(future)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Define quantos cores vão ser usados no processamento paralelo
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 cores  <- as.integer(parallel::detectCores() - 1)
-if (cores > 4L){cores <- cores - 1L} else {cores <- 3L}
 plan(multisession, workers = cores)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Leitura dos catálogos de dados LiDAR (2013 e 2014)
+# Define pastas e local de trabalho
-# Exibe conteúdo do catálogo e confere se os mapas de tiles são iguais
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dirLAZ <- paste0(lidDir, anoData[1])                             # 2013
+prjNam <- "PRJ_Modelo"                               # Nomeia o projeto
-ctg    <- readLAScatalog(dirLAZ)
-opt_select(ctg) <- "xyzic"
-plot(ctg)
-if(!require(mapview))          # Permite criar plots com mapas de fundo
+dirProjet <- str_c('C:/GitRepo/',prjNam) # Define diretório de trabalho
-  install.packages("mapview")
+if (!dir.exists(dirProjet)) {          # Cria diretório caso não exista
-library(mapview)
+  dir.create('C:/GitRepo/', showWarnings = F)
-plot(ctg, mapview = TRUE, map.type = "Esri.WorldImagery")
+  dir.create(dirProjet, showWarnings = F)
+}
+setwd(dirProjet)                     # Define pasta default de trabalho
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+datDir <- str_c('C:/LiDAR/', prjNam)    # Define raiz da pasta de dados
-# Recorta nuvens usando os limites dos talhões com buffer de 10m
+if (!dir.exists(datDir)) {           # Cria diretórios caso não existam
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  dir.create('C:/LiDAR/', showWarnings = F)
-opt_chunk_buffer(ctg) <- 10                               # 10 m buffer
+  dir.create(datDir, showWarnings = F)
-opt_output_files(ctg) <- paste0(tempfile(), "_{SUBTALHAO}")
+}
-ctg_tal = clip_roi(ctg, tal_shp)
-plot(ctg_tal, chunk=TRUE)
+lidDir <- str_c(datDir, '/NUVENS/')     # Cria sub-pasta p/ dados LiDAR
-plot(ctg_tal, mapview = TRUE, map.type = "Esri.WorldImagery")
+if (!dir.exists(lidDir)) {
+  dir.create(lidDir, showWarnings = F)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+}
-# Classifica pontos de solo
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-if(!require(RCSF))    # Pacote complementar da função classify_ground()
-  install.packages("RCSF")
-library(RCSF)
-opt_output_files(ctg_tal) <- paste0(tempdir(), "/{*}_clas")
-ctg                       <- classify_ground(ctg_tal, csf())
-rm(ctg_tal)                                             # limpa memória
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Normaliza as nuvens de pontos dos talhões
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-opt_output_files(ctg) <- paste0(tempdir(), "/{*}_norm")
-ctg                   <- normalize_height(ctg, tin())
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Converte e salva o novo conjunto de tiles normalizado
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dirLAZ <- paste0(dirLAZ, '/NRMLZD')
-opt_output_files(ctg) <-
-  paste0(dirLAZ, "/Modelo_{XLEFT}_{YBOTTOM}")    # renomeia com coords
-opt_chunk_buffer(ctg) <- 20        # com buffers ao redor de cada tile
-opt_chunk_size(ctg) <- 300                # retile para tiles de 300 m
-opt_laz_compression(ctg) <- TRUE                  # mantém formato LAZ
-ctg <- catalog_retile(ctg)                                   # executa
-plot(ctg, mapview = TRUE, map.type = "Esri.WorldImagery")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Lê o shape e atributos dos talhoes
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-if(!require(sf))                           # Para manipulação de shapes
-  install.packages("sf")
-library(sf)
-shpArq  <- paste0(shpDir, "/Modelo_talhoes.shp")     # shape de talhões
-tal_shp <- st_read(shpArq)
-st_is_valid(tal_shp)
-plot(tal_shp, add = TRUE, col = "green")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Lê os centróides das parcelas de inventário convencional
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-shpArq  <- paste0(shpDir,'/Modelo_parcentr.shp')  # centro das parcelas
-par_shp <- st_read(shpArq)
-st_is_valid(par_shp)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Este bloco cria funções que permitem produzir gráficos esteticamente
 # caprichados para apresentar os índices de correlação de Pearson
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Parâmetros dos gráficos
 txt.size = 6
@@ Linha 218: / Linha 153: @@
 # Função para embelazamento da cormat
-if(!require(reshape2))
-  install.packages("reshape2")
-library(reshape2)
 graphCormat <- function(cormat){
   # Resets font size for correlation matrix graph
@@ Linha 256: / Linha 188: @@
   return(p)
 }
-# -----
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Estudo da correlação das métricas LiDAR com parâmetros de interesse
+# Estimativas. e estatísticas de qualidade da inferência, resultantes
+# da amostragem dupla (double sampling)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dsStats = function(y, x, xLarge, alpha=.05){
+  n = length(y)
+  beta = ( sum(y*x, na.rm=T) - ( sum(x, na.rm=T)*sum(y, na.rm=T) / n ))  / ( sum(x^2, na.rm=T) - (sum(x, na.rm=T)^2 / n) )
+  rho = cor(y,x)
+  N = length(xLarge)
+  ydsr = mean(y, na.rm=T) + beta * ( mean(xLarge, na.rm=T) - mean(x, na.rm=T) )
+  vardsr = (var(y, na.rm=T)/n)*(1 - (rho^2)*(N-n)/N)
+  stderr = sqrt(vardsr)
+  ci     = calcCI(stderr, n, alpha)
+  out = c(media = ydsr, var = vardsr, dp = stderr, ic = ci, erro = 100*ci/ydsr, n=n, rho=rho)
+  return(out)
+}
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Cálculo do número ideal de unidades amostrais em double sampling, ou
+# seja, da intensidade amostral necessária para gerar o erro admissível
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dsNumberOfPlots = function(y, x, xLarge, Cpg = 300,
+                           errDesired = .05, alpha = .05){
+  rho = cor(x,y)
+  a = var(y) * (1 - rho^2)
+  b = var(y) * rho^2
+  B = errDesired * mean(y)
+  qt = qt(1 - alpha/2, length(y)-1)
+  nG = (sqrt( a*b*Cpg ) + b) / (B^2 / qt^2)
+  nP = (sqrt( a*b/Cpg ) + a) / (B^2 / qt^2)
+  return(nP)
+}
+# ----
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 2. Leitura de dados
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-opt_filter(ctg) <- "-drop_z_below 0"          # Ignora pontos com z < 0
-# Cálculo "padrão" de métricas LiDAR (veja manuais do lidR)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-D <- plot_metrics(ctg, .stdmetrics_z, par_shp, radius = 11.28)
+# Define o nome do arquivo com os dados espacializados disponíveis para
+# a área de estudo, faz o download e salva esses dados na devida pasta
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+gitOnde <-
+ "https://github.com/FlorestaR/dados/blob/main/5_LIDARF/Modelo/CLOUDS/"
+gitNome <-c("L0004-C0005.laz", "L0004-C0006.laz",
+            "L0005-C0005.laz", "L0005-C0006.laz",      # arquivos LiDAR
+            "L0006-C0005.laz", "L0006-C0006.laz")
+anoData <- c("A13")                  # Lê apenas os dados LiDAR de 2013
+# anoData <- c("A13", "A14")       # Se quiser baixar os dois conjuntos
-# Escolhe um subgrupo de métricas e dados para estudo da correlação
+# Download dos tiles
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-X <- tibble(D) %>% select(VTCC, MHDOM, IDINV, zmean,
+options(timeout=1000) # Reset timeout oferecendo mais tempo de download
-                          zq90, zq95, pzabovezmean, pzabove2)
+for (ano in anoData) {     # guarda na sub-pasta de dados LiDAR por ano
+    dirLAZ <- str_c(lidDir, ano)
+  dir.create(dirLAZ, showWarnings = T)
+  for (nome in gitNome){         # acrescenta "?raw=true" no fim do URL
+    gitFile <- str_c(gitOnde, ano, "/", ano, nome, "?raw=true")
+    localFl <- str_c(dirLAZ, "/", nome)
+    if(!file.exists(localFl))           # garante download binário (wb)
+      download.file(gitFile, mode="wb", destfile = localFl)
+  }
+}
-# Prepara o gráfico para análise das correlações
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-cormat <- X %>% cor()
+# Define o nome do arquivo com os dados espacializados disponíveis para
-cormat <- round(cor(cormat), 2) %>% reorder_cormat()
+# a área de estudo, faz o download e salva esses dados na devida pasta
-p      <- graphCormat(cormat)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+gitOnde <-
+  "https://github.com/FlorestaR/dados/blob/main/5_LIDARF/Modelo/SHAPES"
+gitNome <- "fazmodelo.zip"                   # shapes da fazenda modelo
+gitArqv <- file.path(gitOnde, gitNome) %>% str_c("?raw=true")
-grfDir <- str_c(wrkDir, '/RSLTS/') dir.create(grfDir, showWarnings = F)
+tmpd <- tempdir(check = TRUE)                    # diretório temporário
+zipf <- file.path(tmpd, "shapes.zip")              # arquivo temporário
+# Download e unzip do coonteúdo do arquivo zipado
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+options(timeout=1000)
+if(!file.exists(zipf)) {
+  download.file(gitArqv, mode="wb", destfile = zipf)
+}
+shpDir <- str_c(dirProjet, '/SHAPES') # Define diretório para os shapes
+if (!dir.exists(shpDir)) {             # Cria diretório caso não exista
+  dir.create(shpDir, showWarnings = F)
+}
+unzip(zipf, exdir = shpDir)   # unzipa shps e guarda na pasta de shapes
+unlink(zipf)                                  # deleta o arquivo zipado
+# ----
-fileGrf <- paste0(grfDir, "CorPearson.jpg")
-nomeGrf <- "Matriz de Correlações de Pearson"
-jpeg(fileGrf, width = 1000, height = 1000, units = "px", quality = 100)
-plot(p + labs(title=nomeGrf))
-dev.off()
-# Análise de regressão por quadrados mínimos ordinários (OLS) para
-# determinação do modelo de predição.  Testando: VTCC ~f(zq90,IDINV)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-m <- lm(VTCC ~ zq90 + IDINV, data = X)
+# 3. Processamento da variável auxiliar LiDAR e estudo da correlação
-summary(m)
+#    das métricas LiDAR com os parâmetros de interesse
-plot(X$VTCC, predict(m))
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-abline(0,1)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Gera métricas para Área Total
+# Leitura do catálogo de dados LiDAR                     (2013 ou 2014)
+# Exibe conteúdo do catálogo e confere se os mapas de tiles são iguais
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# É possível gerar uma função que calcula qualquer métrica desejada
+dirLAZ  <- str_c(lidDir, anoData[1])   #  Define local das nuvens LiDAR
+ctg_orig <- readLAScatalog(dirLAZ) # LEITURA DAS NUVENS DE PONTOS LiDAR
+opt_select(ctg_orig) <- "xyzic" # Só atributos xyz intensid. e classif.
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-metriLiDAR <-function(z){             # Criação de métricas específicas
+# Lê shape e atributos dos talhoes
-  if (length(z)<= 0) return(NULL)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  metris <- list(
+shpArq  <- str_c(shpDir, "/Modelo_talhoes.shp")      # shape de talhões
-    n        <- length(z),
+tal_shp <- st_read(shpArq)
-    zmean    <- mean(z),
+st_is_valid(tal_shp)        # confere se os elementos do shape estão OK
-    zq90     <- quantile(z, 0.90),
-    zq95     <- quantile(z, 0.95),
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    zrelief  <- mean(z)/max(z),
+# Os tiles de dados originais abrangem uma área maior do que a
-    zabvmean <- (length(z[z>mean(z)])/length(z))*100
+# necessária. Para reduzir a informação LiDAR que será usada, uma
-  )
+# alternativa é usar os limites dos talhões mais um buffer de 10m
-  names(metris) <- c("n","zmean","zq90","zq95",
+# como "tesoura" para clipar as nuvens.
-                     "zrelief","zabvmean")
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  return(metris)
+dirTal <- str_c(dirLAZ, '/TALHOES')        # Pasta p/ nuvens por talhao
+if (!dir.exists(dirTal)) {
+  dir.create(dirTal, showWarnings = F)
 }
+opt_chunk_buffer(ctg_orig) <- 10                          # 10 m buffer
+opt_output_files(ctg_orig) <- str_c(dirTal, "/{SUBTALHAO}")
+opt_laz_compression(ctg_orig) <- TRUE              # Mantém formato LAZ
+ctg_talh = clip_roi(ctg_orig, tal_shp)       # Usa shape como "tesoura"
-# Aplica a função de métricas específicas em cada talhão separadamente
+# Para plotar a nuvem LiDAR de um dos talhões no catálogo ctg_tal,
+# basta ler o respectivo arquivo salvo pela função clip_roi(),
+# como no exemplo abaixo:
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-prjNam <- "PRJ_Modelo"
+# las <- readLAS(ctg_tal$filename[3])
-wrkDir <- str_c('C:/LiDAR/', prjNam)
+# plot(las)
-lidDir <- str_c(wrkDir, '/CLOUDS/')
-ctgDir <- str_c(lidDir, anoData[1], '/TALHAO')
-arq1a <- paste0(lidDir, anoData[1], '/TALHAO/301a.laz')
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-arq1d <- paste0(lidDir, anoData[1], '/TALHAO/301d.laz')
+# Classifica pontos de solo             (este passo pode ser bem lento)
-arq2a <- paste0(lidDir, anoData[1], '/TALHAO/302a.laz')
+#       ... espere até que o sinal "stop" desapareça na janela de "log"
-arq2c <- paste0(lidDir, anoData[1], '/TALHAO/302c.laz')
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirNoN <- str_c(dirTal, '/NoNORM') # Nuvens c/ ground não normalizadas
+if (!dir.exists(dirNoN)) {
+  dir.create(dirNoN, showWarnings = F)
+}
+opt_output_files(ctg_talh) <- str_c(dirNoN, "/{*}_g")
+ctg_gNoN                   <- classify_ground(ctg_talh, csf())
-las <- readLAS(arq1a)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-d1a <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
+# Normaliza as nuvens de pontos dos talhões          (igualmente lento)
-plot(d1a, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirNrm <- str_c(dirTal, '/SiNORM')       # Pasta p/ nuvens normalizadas
+if (!dir.exists(dirNrm)) {
+  dir.create(dirNrm, showWarnings = F)
+}
+opt_output_files(ctg_gNoN ) <- str_c(dirNrm, "/{*}SiN")
+ctg_gSiN                    <- normalize_height(ctg_gNoN, tin())
-las <- readLAS(arq1d)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-d1d <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
+# Retile dos talhões normalizados para geração de um conjunto de núvens
-plot(d1d, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
+# mais apropriado para a correta clipagem das parcelas
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirTNrm <- str_c(dirNrm, '/TILES')        # Pasta p/ tiles normalizados
+if (!dir.exists(dirTNrm)) {
+  dir.create(dirTNrm, showWarnings = F)
+}
+# ctg_gSiN <- readLAScatalog(dirTNrm)
+opt_output_files(ctg_gSiN) <-     # Onde guardar as nuvens normalizadas
+  str_c(dirTNrm, "/FzMod_{XLEFT}_{YBOTTOM}")    # renomeadas com coords
+opt_chunk_buffer(ctg_gSiN) <- 0     # sem buffers ao redor de cada tile
+opt_chunk_size(ctg_gSiN) <- 300                     # em tiles de 300 m
+opt_laz_compression(ctg_gSiN) <- TRUE              # Mantém formato LAZ
+ctg_tile <- catalog_retile(ctg_gSiN)                        # e executa
+rm(ctg_gNoN, ctg_gSiN, ctg_orig, ctg_talh) # limpa da memória catálogos
-las <- readLAS(arq2a)
+# Lê shape e atributos das parcelas p/ acessar parâmetros de interesse
-d2a <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-plot(d2a, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
+shpPar  <- str_c(shpDir, "/Modelo_parcentr.shp")    # shape de parcelas
+par_shp <- st_read(shpPar)
+st_is_valid(par_shp)        # confere se os elementos do shape estão OK
-las <- readLAS(arq2c)
+dirParc <- str_c(dirNrm, '/PARCELAS')      # Pasta p/ parcelas clipadas
-d2c <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
+if (!dir.exists(dirParc)) {
-plot(d2c, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
+  dir.create(dirParc, showWarnings = F)
-# ----
+}
+# Exibe localização dos tiles, talhoes e parcelas
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+plot(ctg_tile, mapview = TRUE, map.type = "Esri.WorldImagery")
+plot(ctg_tile)
+plot(tal_shp, add = TRUE, col="black")
+plot(par_shp, add = TRUE, col="white")
+# Clipa e salva as núvens das parcelas, e calcula
+# métricas para cada "voxel" clipado com o shape de parcelas no tiles
+# normalizados usando a função .stdmetrics_z de métricas genéricas
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+opt_output_files(ctg_tile) <-     # Onde guardar as nuvens das parcelas
+  str_c(dirParc, "/P_{NUMPARCELA}") # renomeadas pelo respectivo número
+opt_select(ctg_tile) <- "xyz"   # Carrega na memória apenas coordenadas
+opt_filter(ctg_tile) <- "-drop_z_below 0"     # Ignora pontos com z < 0
+D <- plot_metrics(ctg_tile, .stdmetrics_z, par_shp, radius = 11.28)
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# Escolhe um subgrupo de métricas e dados para estudo da correlação
-#    AGORA É SUA VEZ
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#    ESTUDO DIRIGIDO ...
+X <- tibble(D) %>% select(VTCC, MHDOM, IDINV, zmean,
-#    PRODUZA O MAPA DE ESTIMATIVAS DE VOLUME
+                          zq45, zq75, zq95, zpcum2, zpcum4, zpcum6,
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+                          pzabovezmean, pzabove2)
+# Prepara o gráfico para análise das correlações
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+cormat <- X %>% cor()
+cormat <- round(cor(cormat), 2) %>% reorder_cormat()
+p      <- graphCormat(cormat)
+grfDir <- str_c(dirProjet, '/RESULTADOS/') # Define pasta p/ resultados
+if (!dir.exists(grfDir)) {                # Cria pasta, caso não exista
+  dir.create(grfDir, showWarnings = F)
+}
+fileGrf <- str_c(grfDir, "MatrizDeCorrelacoes.jpg")   # Arq para matriz
+tituGrf <- "Matriz de Correlações de Pearson"        # Define um título
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# Abre "impressora" para gerar e salvar o gráfico no formato JPEG
-#    O BLOCO A SEGUIR SERÁ UTIL SE QUISER "ESPECULAR" MAIS
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#    CALCULA AS INFORMAÇÕES NECESSÁRIAS PARA GERAR O QUADRO
+jpeg(fileGrf, width = 1000, height = 1000, units = "px", quality = 200)
-#    DE COMPARAÇÃO ENTRE ACS, ACE E AD !!
+plot(p + labs(title=tituGrf))
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -----
+dev.off()                                        # Fecha a "impressora"
-# Funções para a
+# Análise de regressão por quadrados mínimos ordinários (OLS) para
-# Amostragem dupla
+# determinação do modelo de predição.  Testando: VTCC ~f(zq90,IDINV)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-doubleSampleRegPars = function(y, x, xLarge, alpha=.05){
+m <- lm(VTCC ~ zq95 + IDINV, data = X)    # Análise de Regressão Linear
-  n = length(y)
+summary(m)                          # Mostra os resultados da regressão
-  beta = ( sum(y*x, na.rm=T) - ( sum(x, na.rm=T)*sum(y, na.rm=T) / n ))
+plot(X$VTCC, predict(m))              # Gráfico de observado vs predito
-       / ( sum(x^2, na.rm=T) - (sum(x, na.rm=T)^2 / n) )
+abline(0,1)
+# ----
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 4. Criação de mapas raster com as métricas mais promissoras
+#    É possível processar a partir deste bloco, desde que já tenham
+#    sido processadas as linhas 1 a 229 (Bloco 1.)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
+prjNam    <- "PRJ_Modelo"                             # Nome do projeto
+dirRaizNuvens <- str_c('C:/LiDAR/',prjNam, '/NUVENS') # Raiz das nuvens
+dirDadosLiDAR <- str_c(dirRaizNuvens, '/A13/TALHOES/SiNORM')
+# Lê catálogo de com as núvens por talhão normalizadas
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ctg_taln <- readLAScatalog(dirDadosLiDAR)
+nNuvens  <- length(ctg_taln$filename)    # Número de núvens no catálogo
+# Criação de mapas raster formato tif para as métricas de interesse
+#   Atribua à variável "interesse" uma das métricas calculadas pela
+#   função padrão ".stdmetrics_z", dentre as várias: zmax, zmean,
+#   zsd, zskew, zkurt, zentropy, pzabovezmean, pzabove2, zq5, zq10,
+#   zq15, zq20, zq25, zq30, zq35, zq40, zq45, zq50, zq55, zq60, zq65,
+#   zq70, zq75, zq80, zq85, zq90, zq95, zpcum1, zpcum2, zpcum3, zpcum4,
+#   zpcum5, zpcum6, zpcum7, zpcum8, zpcum9
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+interesse <- "zq95"
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Criação de mapas raster para a métrica de interesse  (PIXEL QUADRADO)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirDadosRaster <- str_c(dirDadosLiDAR, '/RSTR_qua/')        # Quadrados
+if (!dir.exists(dirDadosRaster)) {
+  dir.create(dirDadosRaster, showWarnings = F)
+}
+pb <- progress_bar$new(total = 2*nNuvens) # Reset da barra de progresso
+for (i in 1:nNuvens) {
+  nuvem   <- ctg_taln$filename[i]
+  talhao  <- ctg_taln$filename[i] %>% basename() %>% file_path_sans_ext()
-  rho = cor(y,x)
+  las <- readLAS(nuvem, select = "xyz",
-  N = length(xLarge)
+                 filter = "-drop_z_below 2")    # lê apenas pontos > 2m
-  ydsr = mean(y, na.rm=T) + beta * ( mean(xLarge, na.rm=T) - mean(x, na.rm=T) )
+  pb$tick()                                 # Avança barra de progresso
-  vardsr = (var(y, na.rm=T)/n)*(1 - (rho^2)*(N-n)/N)
+  inteRaster <- (las %>%   # Cria objeto raster da métrica de interesse
-  stderr = sqrt(vardsr)
+                   pixel_metrics(.stdmetrics_z, res = 20))[[interesse]]
-  ci     = calcCI(stderr, n, alpha)
-  out = c(media = ydsr, var = vardsr, dp = stderr, ic = ci, erro = 100*ci/ydsr, n=n, rho=rho)
+  rstNome <- str_c(dirDadosRaster, talhao,'_sqr_', interesse, '.tif')
+  writeRaster(inteRaster, rstNome, overwrite=TRUE)  # Salva arquivo tif
-  return(out)
+  pb$tick()                                 # Avança barra de progresso
+  rstNome <- str_c(dirDadosRaster, talhao,'_sqr_', interesse, '.png')
+  # Cria parâmetros para melhorar a legenda do gráfico raster
+  minL = round(minmax(inteRaster)[1])-1         # Menor valor no raster
+  maxL = round(minmax(inteRaster)[2])+1         # Maior valor no raster
+  brkL = round((maxL - minL) / 5)            # Legenda com cinco breaks
+  intL = c(minL + brkL, minL + 2*brkL, minL + 3*brkL, minL + 4*brkL)
+  limL = c(minL, maxL)
+  titulo  <- str_c('Talhão: ', talhao, ' | Métrica: ', interesse,
+                   ' | Mín: ', minL, ' | Máx: ', maxL)
+  p <- ggplot() + geom_spatraster(data = inteRaster, aes()) +
+    guides(fill = guide_legend(reverse=F)) +
+    coord_sf(datum=st_crs(31983)) +       # Gráfico numa certa projeção
+    scale_fill_whitebox_c(palette = "atlas", direction = -1,
+                          breaks = intL, limits = limL)+ggtitle(titulo)
+  png(rstNome, 30, 20, 'cm', res = 200)              # Abre "impressão"
+  print(p)
+  dev.off()                               # Fecha "impressão" do aquivo
 }
-# Função para determinação da quantidade ideal de parcelas amostrais,
-# isto é, da intensidade amostral que gera erro admissível
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dubleSamplePlotNumber = function(y, x, xLarge, Cpg = 300,
+# !!! Método alternativo:
-                                 errDesired = .05, alpha = .05){
+# Criação de mapas raster para a métrica de interesse (PIXEL HEXAGONAL)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# dirDadosRaster <- str_c(dirDadosLiDAR, '/RSTR_hex/')       # Hexagonais
+# if (!dir.exists(dirDadosRaster)) {
+#   dir.create(dirDadosRaster, showWarnings = F)
+# }
+#
+# pb <- progress_bar$new(total = 2*nNuvens) # Reset da barra de progresso
+# for (i in 1:nNuvens) {
+#   nuvem   <- ctg_taln$filename[i]
+#   talhao  <- ctg_taln$filename[i] %>% basename() %>% file_path_sans_ext()
+#
+#   las <- readLAS(nuvem, select = "xyz",
+#                  filter = "-drop_z_below 2")    # lê apenas pontos > 2m
+#
+#   pb$tick()                                 # Avança barra de progresso
+#   inteRaster <- las %>%
+#     hexagon_metrics(.stdmetrics_z, area = 400) %>%
+#     select(contains(c(interesse, "geometry"))) %>%
+#     st_rasterize()
+#
+#   rstNome <- str_c(dirDadosRaster, talhao,'_hex_', interesse, '.tif')
+#   write_stars(inteRaster, rstNome)
+#
+#   pb$tick()                                 # Avança barra de progresso
+#   rstNome <- str_c(dirDadosRaster, talhao,'_hex_', interesse, '.png')
+#   titulo  <- str_c('Talhão: ', talhao, '  /  Métrica: ', interesse)
+#   inteRaster <- as(inteRaster, "SpatRaster")
+#   # Cria parâmetros para melhorar a legenda do gráfico raster
+#   minL = round(minmax(inteRaster)[1])-1         # Menor valor no raster
+#   maxL = round(minmax(inteRaster)[2])+1         # Maior valor no raster
+#   brkL = round((maxL - minL) / 5)            # Legenda com cinco breaks
+#   intL = c(minL + brkL, minL + 2*brkL, minL + 3*brkL, minL + 4*brkL)
+#   limL = c(minL, maxL)
+#   titulo  <- str_c('Talhão: ', talhao, ' | Métrica: ', interesse,
+#                    ' | Mín: ', minL, ' | Máx: ', maxL)
+#   p <- ggplot() + geom_spatraster(data = inteRaster, aes()) +
+#     guides(fill = guide_legend(reverse=F)) +
+#     coord_sf(datum=st_crs(31983)) +       # Gráfico numa certa projeção
+#     scale_fill_whitebox_c(palette = "atlas", direction = -1,
+#                           breaks = intL, limits = limL)+ggtitle(titulo)
+#   png(rstNome, 30, 20, 'cm', res = 200)              # Abre "impressão"
+#   print(p)
+#   dev.off()                               # Fecha "impressão" do aquivo
+# }
+# Cálculo de mapas raster como função dos mapas rasters de métricas
+# criados nos passos anteriores.  Por exemplo, MAPA DE VOLUME
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ctg_taln <- readLAScatalog(dirDadosLiDAR)
+nNuvens  <- length(ctg_taln$filename)    # Número de núvens no catálogo
+# Criação de mapas raster formato tif para as métricas de interesse
+#   Atribua à variável "interesse" uma das métricas calculadas pela
+#   função .stdmetrics_z: zmax, zmean, zsd, zskew, zkurt, zentropy,
+#   pzabovezmean, pzabove2, zq5, zq10, zq15, zq20, zq25, zq30, zq35,
+#   zq40, zq45, zq50, zq55, zq60, zq65, zq70, zq75, zq80, zq85, zq90,
+#   zq95, zpcum1, zpcum2, zpcum3, zpcum4, zpcum5, zpcum6, zpcum7,
+#   zpcum8, zpcum9
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+interesse <- "zq95"
+# Criação de lista de mapas raster criados nos passos anteriores
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+listaRasters <- list()
+for (i in 1:nNuvens) {
+  nuvem   <- ctg_taln$filename[i]
+  talhao  <- ctg_taln$filename[i] %>% basename() %>% file_path_sans_ext()
+  rstNome <- str_c(dirDadosRaster, talhao,'_sqr_', interesse, '.tif')
-  rho = cor(x,y)
+  listaRasters[[i]] <- rast(rstNome)
-  a = var(y) * (1 - rho^2)
-  b = var(y) * rho^2
-  B = errDesired * mean(y)
-  qt = qt(1 - alpha/2, length(y)-1)
-  nG = (sqrt( a*b*Cpg ) + b) / (B^2 / qt^2)
-  nP = (sqrt( a*b/Cpg ) + a) / (B^2 / qt^2)
-  return(nP)
 }
+# Junta os rasters para gerar um mapa único com a métrica de interesse
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rastersJuntos <- do.call(merge, listaRasters)
-spatialPlotMetrics  <- spatialPlotMetrics %>%
+# Cria gráfico caprichado da métrica de interesse com todos os rasters
-  mutate(VTCC = ifelse(VTCC <50, VTCC*10,VTCC))
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+minL = round(minmax(rastersJuntos)[1])-1        # Menor valor no raster
+maxL = round(minmax(rastersJuntos)[2])+1        # Maior valor no raster
+brkL = round((maxL - minL) / 5)              # Legenda com cinco breaks
+intL = c(minL + brkL, minL + 2*brkL, minL + 3*brkL, minL + 4*brkL)
+limL = c(minL, maxL)
+titulo  <- str_c('Métrica: ', names(rastersJuntos),
+                 ' | Mín: ', minL, ' | Máx: ', maxL)
+p <- ggplot() + geom_spatraster(data = rastersJuntos, aes()) +
+  guides(fill = guide_legend(reverse=F)) +
+  coord_sf(datum=st_crs(31983)) +       # Gráfico numa certa projeção
+  scale_fill_whitebox_c(palette = "atlas", direction = -1,
+                        breaks = intL, limits = limL)+ggtitle(titulo)
+rstNome <- str_c(dirDadosRaster, 'talhoes_', names(rastersJuntos), '.png')
+png(rstNome, 30, 20, 'cm', res = 200)              # Abre "impressão"
+print(p)
+dev.off()                               # Fecha "impressão" do aquivo
-lims = spatialPlotMetrics %>% names %in% c('zmean', 'DNT') %>% which
+# Lê o arquivo shape com os limites e atributos dos talhoes
-lidarOnly = spatialPlotMetrics[,lims[1]:lims[2]]
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirShape <- str_c('C:/GitRepo/',prjNam, '/SHAPES')   # Pasta com shapes
+arqShape <- str_c(dirShape, '/Modelo_talhoes.shp')      # Nome do shape
-interestVars = c('MHDOM','VTCC', 'AB')
+# Cria estrutura de dados (data frame) com os atributos dos talhoes
-corrMat = cor(spatialPlotMetrics[,interestVars], lidarOnly)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-corrMat = corrMat[,apply(corrMat, 2, function(x) !any(is.na(x)))]
+poligonos <- vect(arqShape)
-handPick = c('zq99', 'zq99', 'zq99')
+# Plota mapa de rasters juntados sobreposto com os limites dos talhões
-vars     = c('MHDOM','VTCC', 'AB')
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-handp    = as.data.frame(cbind(vars, handPick))
+plot(rastersJuntos)
-doubleSamplePars = foreach(i = interestVars, .combine = 'rbind') %do% {
+plot(poligonos, add = TRUE, col = "white")
-  aux = corrMat[i,]
+# ----
-  #pick = which( aux == max(aux) ) %>% names
-  pick = handp$handPick[vars == i]
-  XL = metrics[,pick]
-  est = doubleSampleRegPars(spatialPlotMetrics[,i],
-                            spatialPlotMetrics[,pick], XL) %>% data.frame()
-  #est = doubleSampleRatioPars(spatialPlotMetrics[,i], spatialPlotMetrics[,pick], XL, populationSize) %>% data.frame()
-  names(est) = i
-  est %<>% t %>% as.data.frame
-  est$aux = pick
-  return(est)
-}
-doubleSamplePars <- doubleSamplePars %>% select(-rho) %>% t()
-doubleSamplePars %>%
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-  kbl(caption = "Amostragem Dupla", align = "r") %>%
+# 5 . Cria um novo raster com as estimativas de VOLUME para produção
-  kable_classic(full_width = F)
+#     do mapa de VTCC. A esrimativa será obtida a partir do modelo
+#     ajustado na fase anterior:
+#                           VTCC = -444.142 + 24.63 zq95 + 20.133 IDINV
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-# Gráfico para análise da precisão dos métodos amostrais
+# A lista poligonos$SUBTALHAO contém os nomes dos talhões da Fazenda
+# Modelo:  "302c" "301a" "301d" "302a". Criação de um vetor com as
+# correspondentes idades na mesma ordem, para inclusão no objeto
+# com informações sobres os poligonos
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-erro = 1:20
+Idades <- as.numeric(c(5.2, 3.7, 3.7, 5.2))
-scsPlotn <- tamanhoIdealACS(parcelas$VTCC, tamanhoPopulacao,
+poligonos$IDINV <- Idades
-                            errDesired = erro/100)
-tamanhoIdealACE(y = parcelas$VTCC,
-                g = parcelas$IDINV,
-                Nh = parcPorEstrato,
-                erro/100)
-cssPlotn <- tamanhoIdealACE(parcelas$VTCC, parcelas$IDINV, parcPorEstrato,
+valoresAtrib <- poligonos$IDINV
-                            errDesired = erro/100)
-dsPlotn = dubleSamplePlotNumber(spatialPlotMetrics$VTCC,
+# Cria um novo raster com a idade do respectivo polígono atribuída ao
-                                spatialPlotMetrics$zq99,
+# respectivo pixel sobreposto
-                                metrics$zq99, 300, erro/100)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rastersIdade <- rasterize(poligonos, rastersJuntos,
+                          field=valoresAtrib, ext=ext(rastersJuntos),
+                          res=res(rastersJuntos)) %>% rename(IDINV = layer)
-png('nParcelas.png', 15, 10, 'cm', res = 200)
+# Junta rasters, mantendo os respecivos valores em layers separados,
-plot( scsPlotn ~ erro, type='l', col='red', lwd=2,
+# renomeia as colunas e assim, agora rastersJuntos tem as varíaveis
-      ylim=c(0,200), ylab='Número de parcelas', xlab='Erro amostral (%)',
+# necessárias para estimar o parâmetro de interesse, uma em cada layer.
-      axes=F)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-lines(erro, cssPlotn, col='green', lwd=2)
+rastersFinal <- c(rastersJuntos, rastersIdade)
-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', 'AD'))
-dev.off()
+# Usa a equação ajustada (final do bloco 3 deste script) para estimar
+# o parâmetro de interesse
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Este script foi atualizado em Novembro/2023 depois de instalar,
+rastersVolume <- -444.142 + 24.630 * rastersFinal["zq95"] +
+.133 * rastersFinal["IDINV"]
+rastersVolume <- rastersVolume %>% rename(VTCC = zq95)
+# Junta layer VTCC como novo layer no raster final
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rastersFinal = c(rastersFinal, rastersVolume)
+# Cria novo raster de volumes mudando para zero valores negativos
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+vtcc <- ifel(rastersFinal[["VTCC"]] < 0, 0, rastersFinal)[["VTCC"]]
+minmax(vtcc)[1]
+minmax(vtcc)[2]
+# Apaga rasters intermediários
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rm(inteRaster, rastersIdade, rastersJuntos, rastersVolume)
+# /////////////////////////////////////////////////////////////////////
+#       MAPA DE ESTIMATIVAS DE VOLUME TOTAL COM CASCA (VTCC)
+# /////////////////////////////////////////////////////////////////////
+minL = round(minmax(vtcc)[1])                   # Menor valor no raster
+maxL = round(minmax(vtcc)[2])                   # Maior valor no raster
+brkL = round((maxL - minL) / 6)              # Legenda com cinco breaks
+intL = c(minL+brkL, minL+2*brkL, minL+3*brkL, minL+4*brkL, minL+5*brkL)
+limL = c(minL, maxL)
+titulo  <- str_c('Mapa de estimativas de VTCC | Mín: ', minL,
+                 ' | Máx: ', maxL)
+p <- ggplot() + geom_spatraster(data = vtcc, aes()) +
+  guides(fill = guide_legend(reverse=F)) +
+  coord_sf(datum=st_crs(31983)) +       # Gráfico numa certa projeção
+  scale_fill_whitebox_c(palette = "atlas", direction = -1,
+                        breaks = intL, limits = limL)+ggtitle(titulo)
+rstNome <- str_c(dirDadosRaster, 'talhoes_VTCC.png')
+png(rstNome, 30, 20, 'cm', res = 200)              # Abre "impressão"
+print(p)
+dev.off()                               # Fecha "impressão" do aquivo
+# ----
+# /////////////////////////////////////////////////////////////////////
+#    TENTE AGORA CRIAR A TABELA COM A ESTIMATIVA DE VOLUME CALCULADA
+#    PELA AMOSTRAGEM DUPLA (SEMELHANTE À QUE FIZEMOS PARA ACS E ACE)
+# /////////////////////////////////////////////////////////////////////
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Este script foi atualizado em Dezembro/2023 depois de instalar,
 # na seguinte sequência, as mais recentes versões do pacote lidR*:
 #
@@ Linha 481: / Linha 708: @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #  ----
 </code>
@@ Linha 488: / Linha 716: @@
   * [[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://github.com/r-lidar/lidR/wiki/Area-based-approach-from-A-to-Z | Inventário passo a passo com lidR]]
   * [[https://www.rdocumentation.org/packages/lidR/versions/4.0.1 | RDocumentation (lidR)]]
   * [[https://cran.rstudio.com/web/packages/lidR/vignettes/ | Vignettes]]
+e sobre os pacotes //terra// e //sf// para processamento de dados espaciais:
+  * [[https://rspatial.org/terra/pkg/index.html | The terra package]]
+  * [[https://r-spatial.github.io/sf/articles/sf1.html | Simple Features for R ]]