lidar:projetos:rcode_invcomlidar
Essa é uma revisão anterior do documento!
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Inventário com dados LiDAR na Fazenda Modelo ~~~~~~~~~~~~~~~~~~~~~~~~
#
# Autor: Luiz Carlos Estraviz Rodriguez
# Departamento de Ciências Florestais
# ESALQ/USP - 27/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 pastas locais
# 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)
# RStudio (v. 2023.09.1 Build 494)
# 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 pastas e local de trabalho (altere de acordo se Linux ou Mac)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
prjNam <- "PRJ_Modelo" # Nomeia o projeto
wrkDir <- str_c('C:/GitRepo/', prjNam) # Define diretório de trabalho
if (!dir.exists(wrkDir)) { # Cria diretórios caso não existam
dir.create('C:/GitRepo/', showWarnings = F)
dir.create(wrkDir, showWarnings = F)
}
setwd(wrkDir) # Define pasta default de trabalho
datDir <- str_c('C:/LiDAR/', prjNam) # Define raiz da pasta de dados
if (!dir.exists(datDir)) { # Cria diretórios caso não existam
dir.create('C:/LiDAR/', showWarnings = F)
dir.create(datDir, showWarnings = F)
}
lidDir <- str_c(datDir, '/NUVENS/') # Cria sub-pasta p/ dados LiDAR
if (!dir.exists(lidDir)) {
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", # 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
# Faz o 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)
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)
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 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" # shapes da fazenda modelo
gitArqv <- file.path(gitOnde, gitNome) %>% str_c("?raw=true")
tmpd <- tempdir(check = TRUE) # diretório temporário
zipf <- file.path(tmpd, "shapes.zip") # arquivo temporário
options(timeout=1000)
if(!file.exists(zipf)) {
download.file(gitArqv, mode="wb", destfile = zipf)
}
shpDir <- str_c(wrkDir, '/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
# ----
# *********************************************************************
# 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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
cores <- as.integer(parallel::detectCores() - 1)
if (cores > 4L){cores <- cores - 1L} else {cores <- 3L}
plan(multisession, workers = cores)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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.
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")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Os tiles de dados originais abrangem uma área maior do que a
# necessária. Para reduzir a informação LiDAR que será usada, uma
# alternativa é usar os limites dos talhões mais um buffer de 10m
# como "tesoura" para clipar as nuvens.
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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"
plot(ctg_talh, mapview = TRUE, map.type = "Esri.WorldImagery")
# 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:
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# las <- readLAS(ctg_tal$filename[3])
# plot(las)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Classifica pontos de solo (este passo pode ser bem lento)
# ... espere até que o sinal "stop" desapareça na janela de "log"
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if(!require(RCSF)) # Pacote complementar da função classify_ground()
install.packages("RCSF")
library(RCSF) # A função tempdir() cria uma pasta tempórária
opt_output_files(ctg_talh) <- str_c(tempdir(), "/{*}_g")
ctg_grnd <- classify_ground(ctg_talh, csf())
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Normaliza as nuvens de pontos dos talhões (igualmente lento)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
dirNrm <- str_c(dirLAZ, '/NRMLZD') # Pasta p/ nuvens normalizadas
if (!dir.exists(dirNrm)) {
dir.create(dirNrm, showWarnings = F)
}
dirTaN <- str_c(dirNrm, '/TALHOES') # Sub-pasta p/ talhoes normalizados
if (!dir.exists(dirTaN)) {
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
# LiDAR dentro dos talhões e do buffer de 10m em torno de cada talhão
# criado no passo anterior.
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
opt_output_files(ctg_taln) <- # Onde guardar as nuvens normalizadas
str_c(dirNrm, "/FazMod_{XLEFT}_{YBOTTOM}_n") # renomeadas com coords
opt_chunk_buffer(ctg_taln) <- 0 # sem buffers ao redor de cada tile
opt_chunk_size(ctg_taln) <- 300 # em tiles de 300 m
opt_laz_compression(ctg_taln) <- TRUE # Mantém formato LAZ
ctg_tile <- catalog_retile(ctg_taln) # e executa
plot(ctg_tile, chunk = TRUE)
plot(ctg_tile, mapview = TRUE, map.type = "Esri.WorldImagery")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Lê os centróides das parcelas de inventário convencional
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
# ----
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# 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_tile) <- "-drop_z_below 0" # Ignora pontos com z < 0
# Cálculo "padrão" de métricas LiDAR (veja manuais do lidR)
# para cada "voxel" (um por parcela)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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, '/RESULTADOS/') # Define pasta para 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
nomeGrf <- "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)
plot(p + labs(title=nomeGrf))
dev.off() # Fecha a "impressora"
# 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) # 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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# É possível gerar uma função que calcula qualquer métrica desejada.
# A função abaixo mostra um exemplo dessa estratégia.
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
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)
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Produz mapas "raster" dos valores calculados pela função de métricas
# para cada talhão
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Mapa raster com pixels hexagonais: Talhão 1 do catálogo normalizado
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
nuvem <- ctg_taln$filename[1]
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 com pixels hexagonais: Talhão 2 do catálogo normalizado
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
nuvem <- ctg_taln$filename[2]
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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
nuvem <- ctg_taln$filename[3]
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)
# Mapa raster de pixels 10x10m: Talhão 4 do catálogo normalizado
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
nuvem <- ctg_taln$filename[4]
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)
# ----
# +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# 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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ----
Referências sobre o pacote lidR para processamento de dados LiDAR:
e sobre os pacotes terra e sf para processamento de dados espaciais:
lidar/projetos/rcode_invcomlidar.1701574361.txt.gz · Última modificação: 2024/03/23 20:17 (edição externa)