lidar:projetos:rcode_alsmodelo
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 - 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
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ----
lidar/projetos/rcode_alsmodelo.1701081237.txt.gz · Última modificação: 2024/03/23 20:17 (edição externa)