Diferenças

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

--- lidar:projetos:rcode_invcomlidar [2023/11/28 04:16] – lcer
+++ lidar:projetos:rcode_invcomlidar [2024/03/23 20:17] (atual) – edição externa 127.0.0.1
@@ Linha 7: / Linha 7: @@
 # Autor: Luiz Carlos Estraviz Rodriguez
 #        Departamento de Ciências Florestais
-#        ESALQ/USP - 27/Nov/2023
+#        ESALQ/USP - 02/Dez/2023
 #
 # Inventário florestal da Fazenda Modelo
@@ Linha 23: / Linha 23: @@
 #       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 pastas e local de trabalho (altere de acordo se Linux ou Mac)
+if(!require(sf))            # Para manipulação de shapes e outros tipos
+  install.packages("sf")                      # de dados espacializados
+library(sf)
+if(!require(tidyterra))           # Package para processar mapas raster
+  install.packages("tidyterra")
+library(tidyterra)
+if(!require(terra))           # Package para processar mapas raster
+  install.packages("terra")
+library(terra)
+if(!require(stars))                  # Package que permite rasterização
+  install.packages("stars")
+library(stars)
+if(!require(tools))           # Package para manipular nomes de aquivos
+  install.packages("tools")
+library(tools)
+if(!require(RColorBrewer))          # Package com mais paletas de cores
+  install.packages("RColorBrewer")
+library(RColorBrewer)
+if(!require(progress))          # Package com mais paletas de cores
+  install.packages("progress")
+library(progress)
+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)
+# 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 versão mais recente
+library(lidR)              # e tiver o package devtools instalado, rode
+                                # alternativamente as duas linhas acima
+if(!require(RCSF))    #        Pacote complementar do lidR e necessário
+  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)
+cores  <- as.integer(parallel::detectCores() - 1)
+plan(multisession, workers = cores)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Define pastas e local de trabalho
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 prjNam <- "PRJ_Modelo"                               # Nomeia o projeto
-wrkDir <- str_c('C:/GitRepo/', prjNam)   # Define diretório de trabalho
+dirProjet <- str_c('C:/GitRepo/',prjNam) # Define diretório de trabalho
-if (!dir.exists(wrkDir)) {           # Cria diretórios caso não existam
+if (!dir.exists(dirProjet)) {          # Cria diretório caso não exista
   dir.create('C:/GitRepo/', showWarnings = F)
-  dir.create(wrkDir, showWarnings = F)
+  dir.create(dirProjet, showWarnings = F)
 }
-setwd(wrkDir)                        # Define pasta default de trabalho
+setwd(dirProjet)                     # Define pasta default de trabalho
 datDir <- str_c('C:/LiDAR/', prjNam)    # Define raiz da pasta de dados
@@ Linha 57: / Linha 114: @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Define os nomes dos tiles LiDAR multitemporais que serão lidos do
+# Este bloco cria funções que permitem produzir gráficos esteticamente
-# repositório, os respectivos anos e o URL completo do github
+# caprichados para apresentar os índices de correlação de Pearson
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Parâmetros dos gráficos
+txt.size = 6
+thm.size = (14/5) * txt.size
+Tema <- theme(
+  axis.line    = element_line(size = .5, colour = "black"),
+  axis.text.x  = element_text(angle = 45, size = thm.size, vjust = 1, hjust = 1),
+  axis.title.y = element_blank(),
+  axis.text.y  = element_text(angle = 45, size = thm.size),
+  axis.ticks   = element_blank(),
+  legend.justification = c(1, 0),
+  legend.position  = c(1, 0.5),
+  legend.text      = element_text(size = 12, colour="black"),
+  legend.title = element_blank(),
+  panel.border     = element_blank(),
+  panel.background = element_blank(),
+  panel.grid.major = element_line(colour = "grey"),
+  title = element_text(size = thm.size, colour="black", face="bold"),
+  plot.background = element_rect(colour = "black", fill=NA, size=1)
+)
+# Função para reordenamento da cormat (correlation matrix)
+reorder_cormat <- function(cormat){
+  # Use correlation between variables as distance
+  dd <- as.dist((1-cormat)/2)
+  hc <- hclust(dd)
+  cormat <-cormat[hc$order, hc$order]
+}
+# Função para definição do triângulo superior da cormat
+get_upper_tri <- function(cormat){
+  cormat[lower.tri(cormat)]<- NA
+  return(cormat)
+}
+# Função para embelazamento da cormat
+graphCormat <- function(cormat){
+  # Resets font size for correlation matrix graph
+  txt.size = 6
+  thm.size = (14/5) * txt.size
+  upper_tri <- get_upper_tri(cormat)
+  # Melt the correlation matrix
+  melted_cormat <- melt(upper_tri, na.rm = TRUE)
+  p <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
+    geom_tile(color = "white")+
+    scale_fill_gradient2(low = "blue", high = "red", mid = "white",
+                         midpoint = 0, limit = c(-1,1), space = "Lab",
+                         name="Pearson\nCorrelation") +
+    coord_fixed() +
+    geom_text(aes(Var2, Var1, label = value), color = "black", size = 0.7*txt.size) +
+    theme(
+      axis.line    = element_line(size = .5, colour = "black"),
+      axis.title.x = element_blank(),
+      axis.text.x  = element_text(angle = 45, size = thm.size, vjust = 1, hjust = 1),
+      axis.title.y = element_blank(),
+      axis.text.y  = element_text(angle = 45, size = thm.size),
+      axis.ticks   = element_blank(),
+      legend.direction = "horizontal",
+      legend.justification = c(1, 0),
+      legend.position  = c(0.6, 0.7),
+      legend.text      = element_text(size = 12, colour="black"),
+      panel.border     = element_blank(),
+      panel.background = element_blank(),
+      panel.grid.major = element_line(colour = "grey"),
+      title = element_text(size = thm.size, colour="black", face="bold"),
+      plot.background = element_rect(colour = "black", fill=NA, size=1)
+    ) +
+    guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
+                                 title.position = "top", title.hjust = 0.5))
+  return(p)
+}
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 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
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# 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 <-
@@ Linha 68: / Linha 249: @@
 # anoData <- c("A13", "A14")       # Se quiser baixar os dois conjuntos
-# Faz o download dos tiles
+# Download dos tiles
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 options(timeout=1000) # Reset timeout oferecendo mais tempo de download
 for (ano in anoData) {     # guarda na sub-pasta de dados LiDAR por ano
-  dirLAZ <- str_c(lidDir, ano)
+    dirLAZ <- str_c(lidDir, ano)
   dir.create(dirLAZ, showWarnings = T)
   for (nome in gitNome){         # acrescenta "?raw=true" no fim do URL
@@ Linha 83: / Linha 264: @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Download dos shapes da Fazenda Modelo (2 layers: talhoes e parcelas)
+# 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 <-
@@ Linha 93: / Linha 275: @@
 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)) {
@@ Linha 98: / Linha 282: @@
 }
-shpDir <- str_c(wrkDir, '/SHAPES')    # Define diretório para os shapes
+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)
@@ Linha 107: / Linha 291: @@
 # ----
-# *********************************************************************
-# 2: O processamento da 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
-                           # e tiver o package devtools instalado, rode
-                           # alternativamente as duas linhas acima.
-library(lidR)
-if(!require(future))           # Package que permite ao lidR rodar usar
-  install.packages("future")          # processamento paralelo de dados
-library(future)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Define quantos cores vão ser usados no processamento paralelo
+# 3. Processamento da variável auxiliar LiDAR e estudo da correlação
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+#    das métricas LiDAR com os parâmetros de interesse
-cores  <- as.integer(parallel::detectCores() - 1)
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-if (cores > 4L){cores <- cores - 1L} else {cores <- 3L}
-plan(multisession, workers = cores)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ Linha 137: / Linha 304: @@
 ctg_orig <- readLAScatalog(dirLAZ) # LEITURA DAS NUVENS DE PONTOS LiDAR
 opt_select(ctg_orig) <- "xyzic" # Só atributos xyz intensid. e classif.
-plot(ctg_orig)
-if(!require(mapview))            # Package para incluir mapas de fundo
-  install.packages("mapview")
-library(mapview)
-plot(ctg_orig, mapview = TRUE, map.type = "Esri.WorldImagery")
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Lê shape e atributos dos talhoes
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-if(!require(sf))                           # Para manipulação de shapes
-  install.packages("sf")
-library(sf)
 shpArq  <- str_c(shpDir, "/Modelo_talhoes.shp")      # shape de talhões
 tal_shp <- st_read(shpArq)
 st_is_valid(tal_shp)        # confere se os elementos do shape estão OK
-plot(tal_shp, add = TRUE, col = "white")
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ Linha 169: / Linha 326: @@
 opt_laz_compression(ctg_orig) <- TRUE              # Mantém formato LAZ
 ctg_talh = clip_roi(ctg_orig, tal_shp)       # Usa shape como "tesoura"
-plot(ctg_talh, mapview = TRUE, map.type = "Esri.WorldImagery")
 # Para plotar a nuvem LiDAR de um dos talhões no catálogo ctg_tal,
@@ Linha 182: / Linha 338: @@
 #       ... espere até que o sinal "stop" desapareça na janela de "log"
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-if(!require(RCSF))    # Pacote complementar da função classify_ground()
+dirNoN <- str_c(dirTal, '/NoNORM') # Nuvens c/ ground não normalizadas
-  install.packages("RCSF")
+if (!dir.exists(dirNoN)) {
-library(RCSF)            # A função tempdir() cria uma pasta tempórária
+  dir.create(dirNoN, showWarnings = F)
-opt_output_files(ctg_talh) <- str_c(tempdir(), "/{*}_g")
+}
-ctg_grnd                   <- classify_ground(ctg_talh, csf())
+opt_output_files(ctg_talh) <- str_c(dirNoN, "/{*}_g")
+ctg_gNoN                   <- classify_ground(ctg_talh, csf())
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Normaliza as nuvens de pontos dos talhões          (igualmente lento)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dirNrm <- str_c(dirLAZ, '/NRMLZD')       # Pasta p/ nuvens normalizadas
+dirNrm <- str_c(dirTal, '/SiNORM')       # Pasta p/ nuvens normalizadas
 if (!dir.exists(dirNrm)) {
   dir.create(dirNrm, showWarnings = F)
 }
-dirTaN <- str_c(dirNrm, '/TALHOES') # Sub-pasta p/ talhoes normalizados
+opt_output_files(ctg_gNoN ) <- str_c(dirNrm, "/{*}SiN")
-if (!dir.exists(dirTaN)) {
+ctg_gSiN                    <- normalize_height(ctg_gNoN, tin())
-  dir.create(dirTaN, showWarnings = F)
-}
-opt_output_files(ctg_grnd) <- str_c(dirTaN, "/{*}")
-ctg_taln                   <- normalize_height(ctg_grnd, tin())
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Cria um novo conjunto de tiles normalizados contendo apenas os pontos
+# Retile dos talhões normalizados para geração de um conjunto de núvens
-# LiDAR dentro dos talhões e do buffer de 10m em torno de cada talhão
+# mais apropriado para a correta clipagem das parcelas
-# criado no passo anterior.
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-opt_output_files(ctg_taln) <-     # Onde guardar as nuvens normalizadas
+dirTNrm <- str_c(dirNrm, '/TILES')        # Pasta p/ tiles normalizados
-  str_c(dirNrm, "/FazMod_{XLEFT}_{YBOTTOM}_n")  # renomeadas com coords
+if (!dir.exists(dirTNrm)) {
-opt_chunk_buffer(ctg_taln) <- 0     # sem buffers ao redor de cada tile
+  dir.create(dirTNrm, showWarnings = F)
-opt_chunk_size(ctg_taln) <- 300                     # em tiles de 300 m
+}
-opt_laz_compression(ctg_taln) <- TRUE              # Mantém formato LAZ
+# ctg_gSiN <- readLAScatalog(dirTNrm)
-ctg_tile <- catalog_retile(ctg_taln)                        # e executa
+opt_output_files(ctg_gSiN) <-     # Onde guardar as nuvens normalizadas
-plot(ctg_tile, chunk = TRUE)
+  str_c(dirTNrm, "/FzMod_{XLEFT}_{YBOTTOM}")    # renomeadas com coords
-plot(ctg_tile, mapview = TRUE, map.type = "Esri.WorldImagery")
+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
+# Lê shape e atributos das parcelas p/ acessar parâmetros de interesse
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Lê os centróides das parcelas de inventário convencional
+shpPar  <- str_c(shpDir, "/Modelo_parcentr.shp")    # shape de parcelas
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+par_shp <- st_read(shpPar)
-shpArq  <- str_c(shpDir,'/Modelo_parcentr.shp')   # shape c/ centroídes
-par_shp <- st_read(shpArq)
 st_is_valid(par_shp)        # confere se os elementos do shape estão OK
-# ----
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirParc <- str_c(dirNrm, '/PARCELAS')      # Pasta p/ parcelas clipadas
-# Este bloco cria funções que permitem produzir gráficos esteticamente
+if (!dir.exists(dirParc)) {
-# caprichados para apresentar os índices de correlação de Pearson
+  dir.create(dirParc, showWarnings = F)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-# Parâmetros dos gráficos
-txt.size = 6
-thm.size = (14/5) * txt.size
-Tema <- theme(
-  axis.line    = element_line(size = .5, colour = "black"),
-  axis.text.x  = element_text(angle = 45, size = thm.size, vjust = 1, hjust = 1),
-  axis.title.y = element_blank(),
-  axis.text.y  = element_text(angle = 45, size = thm.size),
-  axis.ticks   = element_blank(),
-  legend.justification = c(1, 0),
-  legend.position  = c(1, 0.5),
-  legend.text      = element_text(size = 12, colour="black"),
-  legend.title = element_blank(),
-  panel.border     = element_blank(),
-  panel.background = element_blank(),
-  panel.grid.major = element_line(colour = "grey"),
-  title = element_text(size = thm.size, colour="black", face="bold"),
-  plot.background = element_rect(colour = "black", fill=NA, size=1)
-)
-# Função para reordenamento da cormat (correlation matrix)
-reorder_cormat <- function(cormat){
-  # Use correlation between variables as distance
-  dd <- as.dist((1-cormat)/2)
-  hc <- hclust(dd)
-  cormat <-cormat[hc$order, hc$order]
 }
-# Função para definição do triângulo superior da cormat
+# Exibe localização dos tiles, talhoes e parcelas
-get_upper_tri <- function(cormat){
-  cormat[lower.tri(cormat)]<- NA
-  return(cormat)
-}
-# Função para embelazamento da cormat
-if(!require(reshape2))
-  install.packages("reshape2")
-library(reshape2)
-graphCormat <- function(cormat){
-  # Resets font size for correlation matrix graph
-  txt.size = 6
-  thm.size = (14/5) * txt.size
-  upper_tri <- get_upper_tri(cormat)
-  # Melt the correlation matrix
-  melted_cormat <- melt(upper_tri, na.rm = TRUE)
-  p <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
-    geom_tile(color = "white")+
-    scale_fill_gradient2(low = "blue", high = "red", mid = "white",
-                         midpoint = 0, limit = c(-1,1), space = "Lab",
-                         name="Pearson\nCorrelation") +
-    coord_fixed() +
-    geom_text(aes(Var2, Var1, label = value), color = "black", size = 0.7*txt.size) +
-    theme(
-      axis.line    = element_line(size = .5, colour = "black"),
-      axis.title.x = element_blank(),
-      axis.text.x  = element_text(angle = 45, size = thm.size, vjust = 1, hjust = 1),
-      axis.title.y = element_blank(),
-      axis.text.y  = element_text(angle = 45, size = thm.size),
-      axis.ticks   = element_blank(),
-      legend.direction = "horizontal",
-      legend.justification = c(1, 0),
-      legend.position  = c(0.6, 0.7),
-      legend.text      = element_text(size = 12, colour="black"),
-      panel.border     = element_blank(),
-      panel.background = element_blank(),
-      panel.grid.major = element_line(colour = "grey"),
-      title = element_text(size = thm.size, colour="black", face="bold"),
-      plot.background = element_rect(colour = "black", fill=NA, size=1)
-    ) +
-    guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
-                                 title.position = "top", title.hjust = 0.5))
-  return(p)
-}
-# -----
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Estudo da correlação das métricas LiDAR com parâmetros de interesse
+plot(ctg_tile, mapview = TRUE, map.type = "Esri.WorldImagery")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
+plot(ctg_tile)
-opt_filter(ctg_tile) <-    "-drop_z_below 0"  # Ignora pontos com z < 0
+plot(tal_shp, add = TRUE, col="black")
+plot(par_shp, add = TRUE, col="white")
-# Cálculo "padrão" de métricas LiDAR (veja manuais do lidR)
+# Clipa e salva as núvens das parcelas, e calcula
-# para cada "voxel" (um por parcela)
+# 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)
@@ Linha 318: / Linha 403: @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 X <- tibble(D) %>% select(VTCC, MHDOM, IDINV, zmean,
-                          zq90, zq95, pzabovezmean, pzabove2)
+                          zq45, zq75, zq95, zpcum2, zpcum4, zpcum6,
+                          pzabovezmean, pzabove2)
 # Prepara o gráfico para análise das correlações
@@ Linha 326: / Linha 412: @@
 p      <- graphCormat(cormat)
-grfDir <- str_c(wrkDir, '/RESULTADOS/')  # Define pasta para resultados
+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, "CorPearson.jpg")    # Nomeia arq para gráfico
+fileGrf <- str_c(grfDir, "MatrizDeCorrelacoes.jpg")   # Arq para matriz
-nomeGrf <- "Matriz de Correlações de Pearson"        # Define um título
+tituGrf <- "Matriz de Correlações de Pearson"        # Define um título
 # Abre "impressora" para gerar e salvar o gráfico no formato JPEG
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-jpeg(fileGrf, width = 1000, height = 1000, units = "px", quality = 100)
+jpeg(fileGrf, width = 1000, height = 1000, units = "px", quality = 200)
-plot(p + labs(title=nomeGrf))
+plot(p + labs(title=tituGrf))
 dev.off()                                        # Fecha a "impressora"
@@ Linha 343: / Linha 429: @@
 # determinação do modelo de predição.  Testando: VTCC ~f(zq90,IDINV)
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-m <- lm(VTCC ~ zq90 + IDINV, data = X)    # Análise de Regressão Linear
+m <- lm(VTCC ~ zq95 + IDINV, data = X)    # Análise de Regressão Linear
 summary(m)                          # Mostra os resultados da regressão
 plot(X$VTCC, predict(m))              # Gráfico de observado vs predito
 abline(0,1)
+# ----
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Gera métricas para Área Total
+# 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
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# É possível gerar uma função que calcula qualquer métrica desejada.
+ctg_taln <- readLAScatalog(dirDadosLiDAR)
-# A função abaixo mostra um exemplo dessa estratégia.
+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
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-metriLiDAR <-function(z){             # Criação de métricas específicas
+interesse <- "zq95"
-  if (length(z)<= 0) return(NULL)
-  metris <- list(
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    n        <- length(z),
+# Criação de mapas raster para a métrica de interesse  (PIXEL QUADRADO)
-    zmean    <- mean(z),
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    zq90     <- quantile(z, 0.90),
+dirDadosRaster <- str_c(dirDadosLiDAR, '/RSTR_qua/')        # Quadrados
-    zq95     <- quantile(z, 0.95),
+if (!dir.exists(dirDadosRaster)) {
-    zrelief  <- mean(z)/max(z),
+  dir.create(dirDadosRaster, showWarnings = F)
-    zabvmean <- (length(z[z>mean(z)])/length(z))*100
+}
-  )
-  names(metris) <- c("n","zmean","zq90","zq95",
+pb <- progress_bar$new(total = 2*nNuvens) # Reset da barra de progresso
-                     "zrelief","zabvmean")
+for (i in 1:nNuvens) {
-  return(metris)
+  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 %>%   # Cria objeto raster da métrica de interesse
+                   pixel_metrics(.stdmetrics_z, res = 20))[[interesse]]
+  rstNome <- str_c(dirDadosRaster, talhao,'_sqr_', interesse, '.tif')
+  writeRaster(inteRaster, rstNome, overwrite=TRUE)  # Salva arquivo tif
+  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
 }
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Produz mapas "raster" dos valores calculados pela função de métricas
+# !!! Método alternativo:
-# para cada talhão
+# 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
+# }
-# Mapa raster com pixels hexagonais: Talhão 1 do catálogo normalizado
+# Cálculo de mapas raster como função dos mapas rasters de métricas
+# criados nos passos anteriores.  Por exemplo, MAPA DE VOLUME
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-nuvem <- ctg_taln$filename[1]
+ctg_taln <- readLAScatalog(dirDadosLiDAR)
-las <- readLAS(nuvem, filter = "-drop_z_below 2")      # lê pontos > 2m
+nNuvens  <- length(ctg_taln$filename)    # Número de núvens no catálogo
-map <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
-plot(map, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-# Mapa raster com pixels hexagonais: Talhão 2 do catálogo normalizado
+# 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
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-nuvem <- ctg_taln$filename[2]
+interesse <- "zq95"
-las <- readLAS(nuvem, filter = "-drop_z_below 2")      # lê pontos > 2m
-map <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
-plot(map, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-# Mapa raster de pixels 10x10m: Talhão 3 do catálogo normalizado
+# Criação de lista de mapas raster criados nos passos anteriores
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-nuvem <- ctg_taln$filename[3]
+listaRasters <- list()
-las <- readLAS(nuvem, filter = "-drop_z_below 2")      # lê pontos > 2m
+for (i in 1:nNuvens) {
-map <- pixel_metrics(las, metriLiDAR(Z), res = 10)
+  nuvem   <- ctg_taln$filename[i]
-plot(map, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
+  talhao  <- ctg_taln$filename[i] %>% basename() %>% file_path_sans_ext()
+  rstNome <- str_c(dirDadosRaster, talhao,'_sqr_', interesse, '.tif')
+  listaRasters[[i]] <- rast(rstNome)
+}
-# Mapa raster de pixels 10x10m: Talhão 4 do catálogo normalizado
+# Junta os rasters para gerar um mapa único com a métrica de interesse
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-nuvem <- ctg_taln$filename[4]
+rastersJuntos <- do.call(merge, listaRasters)
-las <- readLAS(nuvem, filter = "-drop_z_below 2")      # lê pontos > 2m
-map <- pixel_metrics(las, metriLiDAR(Z), res = 10)
-plot(map, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-# ----
+# Cria gráfico caprichado da métrica de interesse com todos os rasters
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+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
+# Lê o arquivo shape com os limites e atributos dos talhoes
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+dirShape <- str_c('C:/GitRepo/',prjNam, '/SHAPES')   # Pasta com shapes
+arqShape <- str_c(dirShape, '/Modelo_talhoes.shp')      # Nome do shape
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# Cria estrutura de dados (data frame) com os atributos dos talhoes
-#    AGORA É SUA VEZ
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#    ESTUDO DIRIGIDO ...
+poligonos <- vect(arqShape)
-#    PRODUZA O MAPA DE ESTIMATIVAS DE VOLUME
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# Plota mapa de rasters juntados sobreposto com os limites dos talhões
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+plot(rastersJuntos)
+plot(poligonos, add = TRUE, col = "white")
+# ----
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-#    O BLOCO A SEGUIR SERÁ UTIL SE QUISER "ESPECULAR" MAIS
-#    CALCULA AS INFORMAÇÕES NECESSÁRIAS PARA GERAR O QUADRO
-#    DE COMPARAÇÃO ENTRE ACS, ACE E AD !!
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -----
-# Funções para a
-# Amostragem dupla
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-doubleSampleRegPars = function(y, x, xLarge, alpha=.05){
+# 5 . Cria um novo raster com as estimativas de VOLUME para produção
-  n = length(y)
+#     do mapa de VTCC. A esrimativa será obtida a partir do modelo
-  beta = ( sum(y*x, na.rm=T) - ( sum(x, na.rm=T)*sum(y, na.rm=T) / n ))
+#     ajustado na fase anterior:
-       / ( sum(x^2, na.rm=T) - (sum(x, na.rm=T)^2 / n) )
+#                           VTCC = -444.142 + 24.63 zq95 + 20.133 IDINV
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-  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)
-}
-# Função para determinação da quantidade ideal de parcelas amostrais,
+# A lista poligonos$SUBTALHAO contém os nomes dos talhões da Fazenda
-# isto é, da intensidade amostral que gera erro admissível
+# 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
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dubleSamplePlotNumber = function(y, x, xLarge, Cpg = 300,
+Idades <- as.numeric(c(5.2, 3.7, 3.7, 5.2))
-                                 errDesired = .05, alpha = .05){
+poligonos$IDINV <- Idades
-  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)
-}
+valoresAtrib <- poligonos$IDINV
-spatialPlotMetrics  <- spatialPlotMetrics %>%
+# Cria um novo raster com a idade do respectivo polígono atribuída ao
-  mutate(VTCC = ifelse(VTCC <50, VTCC*10,VTCC))
+# respectivo pixel sobreposto
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rastersIdade <- rasterize(poligonos, rastersJuntos,
+                          field=valoresAtrib, ext=ext(rastersJuntos),
+                          res=res(rastersJuntos)) %>% rename(IDINV = layer)
-lims = spatialPlotMetrics %>% names %in% c('zmean', 'DNT') %>% which
+# Junta rasters, mantendo os respecivos valores em layers separados,
-lidarOnly = spatialPlotMetrics[,lims[1]:lims[2]]
+# renomeia as colunas e assim, agora rastersJuntos tem as varíaveis
+# necessárias para estimar o parâmetro de interesse, uma em cada layer.
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rastersFinal <- c(rastersJuntos, rastersIdade)
-interestVars = c('MHDOM','VTCC', 'AB')
+# Usa a equação ajustada (final do bloco 3 deste script) para estimar
-corrMat = cor(spatialPlotMetrics[,interestVars], lidarOnly)
+# o parâmetro de interesse
-corrMat = corrMat[,apply(corrMat, 2, function(x) !any(is.na(x)))]
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+rastersVolume <- -444.142 + 24.630 * rastersFinal["zq95"] +
+.133 * rastersFinal["IDINV"]
+rastersVolume <- rastersVolume %>% rename(VTCC = zq95)
-handPick = c('zq99', 'zq99', 'zq99')
+# Junta layer VTCC como novo layer no raster final
-vars     = c('MHDOM','VTCC', 'AB')
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-handp    = as.data.frame(cbind(vars, handPick))
+rastersFinal = c(rastersFinal, rastersVolume)
-doubleSamplePars = foreach(i = interestVars, .combine = 'rbind') %do% {
-  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()
+# Cria novo raster de volumes mudando para zero valores negativos
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-doubleSamplePars %>%
+vtcc <- ifel(rastersFinal[["VTCC"]] < 0, 0, rastersFinal)[["VTCC"]]
-  kbl(caption = "Amostragem Dupla", align = "r") %>%
+minmax(vtcc)[1]
-  kable_classic(full_width = F)
+minmax(vtcc)[2]
-# Gráfico para análise da precisão dos métodos amostrais
+# Apaga rasters intermediários
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-erro = 1:20
+rm(inteRaster, rastersIdade, rastersJuntos, rastersVolume)
-scsPlotn <- tamanhoIdealACS(parcelas$VTCC, tamanhoPopulacao,
-                            errDesired = erro/100)
-tamanhoIdealACE(y = parcelas$VTCC,
-                g = parcelas$IDINV,
-                Nh = parcPorEstrato,
-                erro/100)
-cssPlotn <- tamanhoIdealACE(parcelas$VTCC, parcelas$IDINV, parcPorEstrato,
+# /////////////////////////////////////////////////////////////////////
-                            errDesired = erro/100)
+#       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
-dsPlotn = dubleSamplePlotNumber(spatialPlotMetrics$VTCC,
-                                spatialPlotMetrics$zq99,
-                                metrics$zq99, 300, erro/100)
-png('nParcelas.png', 15, 10, 'cm', res = 200)
+# ----
-plot( scsPlotn ~ erro, type='l', col='red', lwd=2,
-      ylim=c(0,200), ylab='Número de parcelas', xlab='Erro amostral (%)',
+# /////////////////////////////////////////////////////////////////////
-      axes=F)
+#    TENTE AGORA CRIAR A TABELA COM A ESTIMATIVA DE VOLUME CALCULADA
-lines(erro, cssPlotn, col='green', lwd=2)
+#    PELA AMOSTRAGEM DUPLA (SEMELHANTE À QUE FIZEMOS PARA ACS E ACE)
-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()
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Este script foi atualizado em Novembro/2023 depois de instalar,
+# Este script foi atualizado em Dezembro/2023 depois de instalar,
 # na seguinte sequência, as mais recentes versões do pacote lidR*:
 #
@@ Linha 540: / Linha 708: @@
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 #  ----
 </code>
@@ Linha 547: / 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 ]]