Diferenças

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

--- lidar:projetos:rcode_alsmodelo [2023/11/27 10:33] – criada 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 inicial de dados
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-if(!require(tidyverse))  # Para melhor manipulação de dados e funções
-  install.packages("tidyverse")
-library(tidyverse)
-# Define diretórios e pastas de trabalho
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-prjNam <- "PRJ_Modelo"
-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)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Define os nomes dos tiles LiDAR multitemporais que serão lidos do
-# repositório, os respectivos anos e o URL completo do github
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-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
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-options(timeout=1000) # Reset timeout oferecendo mais tempo de download
-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)
-  }
-}
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Download dos shapes 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)
-shpDir <- str_c(wrkDir, '/SHAPES')           # diretório para os shapes
-  dir.create(shpDir, showWarnings = F)
-unzip(zipf, exdir = shpDir)                           # shape é unziped
-unlink(zipf)                                  # deleta o arquivo zipado
-# ----
-# *********************************************************************
-# 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
-                           # e tiver o package devtools instalado, rode
-                           # alternativamente as duas linhas acima.
-library(lidR)
-if(!require(future))        # Package que permite ao lidR rodar algumas
-  install.packages("future")           # com processmamento em paralelo
-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)
-# Exibe conteúdo do catálogo e confere se os mapas de tiles são iguais
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-dirLAZ <- paste0(lidDir, anoData[1])                             # 2013
-ctg    <- readLAScatalog(dirLAZ)
-opt_select(ctg) <- "xyzic"
-plot(ctg)
-if(!require(mapview))          # Permite criar plots com mapas de fundo
-  install.packages("mapview")
-library(mapview)
-plot(ctg, mapview = TRUE, map.type = "Esri.WorldImagery")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Recorta nuvens usando os limites dos talhões com buffer de 10m
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-opt_chunk_buffer(ctg) <- 10                               # 10 m buffer
-opt_output_files(ctg) <- paste0(tempfile(), "_{SUBTALHAO}")
-ctg_tal = clip_roi(ctg, tal_shp)
-plot(ctg_tal, chunk=TRUE)
-plot(ctg_tal, mapview = TRUE, map.type = "Esri.WorldImagery")
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# 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
-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
-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
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ----
-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)
-# Escolhe um subgrupo de métricas e dados para estudo da correlação
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-X <- tibble(D) %>% select(VTCC, MHDOM, IDINV, zmean,
-                          zq90, zq95, 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(wrkDir, '/RSLTS/') dir.create(grfDir, showWarnings = F)
-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)
-summary(m)
-plot(X$VTCC, predict(m))
-abline(0,1)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Gera métricas para Área Total
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# É possível gerar uma função que calcula qualquer métrica desejada
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-metriLiDAR <-function(z){             # Criação de métricas específicas
-  if (length(z)<= 0) return(NULL)
-  metris <- list(
-    n        <- length(z),
-    zmean    <- mean(z),
-    zq90     <- quantile(z, 0.90),
-    zq95     <- quantile(z, 0.95),
-    zrelief  <- mean(z)/max(z),
-    zabvmean <- (length(z[z>mean(z)])/length(z))*100
-  )
-  names(metris) <- c("n","zmean","zq90","zq95",
-                     "zrelief","zabvmean")
-  return(metris)
-}
-# Aplica a função de métricas específicas em cada talhão separadamente
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-prjNam <- "PRJ_Modelo"
-wrkDir <- str_c('C:/LiDAR/', prjNam)
-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')
-arq2a <- paste0(lidDir, anoData[1], '/TALHAO/302a.laz')
-arq2c <- paste0(lidDir, anoData[1], '/TALHAO/302c.laz')
-las <- readLAS(arq1a)
-d1a <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
-plot(d1a, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-las <- readLAS(arq1d)
-d1d <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
-plot(d1d, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-las <- readLAS(arq2a)
-d2a <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
-plot(d2a, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-las <- readLAS(arq2c)
-d2c <- hexagon_metrics(las, metriLiDAR(Z), area = 400)
-plot(d2c, pal = heat.colors, axes = TRUE, key.pos = NULL, reset = FALSE)
-# ----
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-#    AGORA É SUA VEZ
-#    ESTUDO DIRIGIDO ...
-#    PRODUZA O MAPA DE ESTIMATIVAS DE VOLUME
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-#    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){
-  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)
-}
-# 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,
-                                 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)
-}
-spatialPlotMetrics  <- spatialPlotMetrics %>%
-  mutate(VTCC = ifelse(VTCC <50, VTCC*10,VTCC))
-lims = spatialPlotMetrics %>% names %in% c('zmean', 'DNT') %>% which
-lidarOnly = spatialPlotMetrics[,lims[1]:lims[2]]
-interestVars = c('MHDOM','VTCC', 'AB')
-corrMat = cor(spatialPlotMetrics[,interestVars], lidarOnly)
-corrMat = corrMat[,apply(corrMat, 2, function(x) !any(is.na(x)))]
-handPick = c('zq99', 'zq99', 'zq99')
-vars     = c('MHDOM','VTCC', 'AB')
-handp    = as.data.frame(cbind(vars, handPick))
-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()
-doubleSamplePars %>%
-  kbl(caption = "Amostragem Dupla", align = "r") %>%
-  kable_classic(full_width = F)
-# Gráfico para análise da precisão dos métodos amostrais
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-erro = 1:20
-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)
-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)
-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', 'AD'))
-dev.off()
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# 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>