Um inventário florestal é uma prática que busca a quantificação, e por conseguinte, a qualificação das florestas em diferentes aspectos, como a produção de madeira entre outros possíveis produtos, através do uso de técnicas estatísticas de amostragem. Sendo assim, o inventário visa obter o máximo de informações de uma floresta com a máxima precisão e o mínimo custo. Esses três fatores desempenham papel importante no planejamento de um inventário florestal (JÚNIOR et al., 2020).
Realizar o inventário florestal de um povoamento com um cultivar de Eucalyptus grandis W. Hill ex Maiden x Eucalyptus urophylla S. T. Blake na região noroeste do Estado de Minas Gerais, para a estimar a produção madeireira dos seguintes produtos de forma independente:
PRODUTO 1: diâmetro mínimo de 5, 5cm sem casca e comprimento de tora de 5, 6m.
PRODUTO 2: diâmetro mínimo de 7, 5cm com casca e comprimento de tora de 3, 3m.
VOLUME TOTAL COM CASCA.
arv <- read.csv2( "cubagem.csv" )
arv$dap <- arv$cap / pi
arv$disc <- arv$dicc - ( 2 * arv$espcasca )
par( mfrow = c( 1, 2 ),
pch = 20,
bg = "#171717",
col.axis = "white",
col.lab = "white",
col.main = "white",
fg = "white" )
with( arv, plot( hi/ht, dicc, pch = 20 ) )
with( arv, plot( hi, dicc, pch = 20 ) )
Figura 1. Relação altura com o diâmetro com casca
Quatro modelos foram ajustados considerando:
\(\beta_1,\beta_2 ...,\beta_n\) = parâmetros do modelo;
\(h_i\) = altura até seção i do fuste;
\(d_i\) = diâmetro na seção i do fuste;
\(d\) = diâmetro à altura do peito;
\(h\) = altura da árvore.
\(p\) = primeiro ponto de inflexão calculado no modelo segmentado de Max & Burkhart (1976).
A partir disso foram feitas as estimativas dos diâmetros sem casca e com casca para cada modelo.
Função de Schöepfer (1966)
#' #### Ajuste com casca
#' ##### Ajuste linear
arv$hih <- arv$hi / arv$ht
arv$did <- arv$dicc / arv$dap
ajlin <- lm( did ~ hih + I( hih^2 ) + I( hih^3 ) + I( hih^4 ) + I( hih^5 ), arv )
#' ##### Erro padrão residual
syx <- summary( ajlin )$sigma
syxperc <- syx/mean( arv$dicc ) * 100 * mean( arv$dap )
#' ##### RMSE
arv$diccest5grau <- predict( ajlin ) * arv$dap
RMSE_policc <- arv %>%
summarise( RMSE = sqrt( sum ( ( diccest5grau - dicc )^2 ) / mean( diccest5grau ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
#' #### Ajuste sem casca
#' ##### Ajuste linear
arv$didsc <- arv$disc / arv$dap
ajlinsc <- lm( didsc ~ hih + I( hih^2 ) + I( hih^3 ) + I( hih^4 ) + I( hih^5 ), arv )
#' ##### Erro padrão residual
syxsc <- summary( ajlinsc )$sigma
syxpercsc <- syxsc/mean( arv$disc ) * 100 * mean( arv$dap )
#' ##### RMSE
arv$discest5grau <- predict( ajlinsc ) * arv$dap
RMSE_polisc <- arv %>%
summarise( RMSE = sqrt( sum ( ( discest5grau - disc )^2 ) / mean( discest5grau ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
Função trigonométrica de forma variável de Bi (2000).
#' #### Ajuste com casca
#' ##### Ajuste não linear
arv$hih <- arv$hi / arv$ht
bi <- nlsLM( dicc ~ taper_bi( dap, ht, hih, b0, b1, b2, b3, b4, b5, b6 ),
data = arv,
start = list( b0 = 1.8,
b1 = -0.2,
b2 = -0.04,
b3 = -0.9,
b4 = -0.0006,
b5 = 0.07,
b6 = -0.14 ) )
#' ##### Erro padrão residual
bisyx <- summary( bi )$sigma
bisyxperc <- bisyx/mean( arv$dicc ) * 100
#' ##### RMSE
arv$diccestbi <- predict( bi )
RMSE_bi <- arv %>%
summarise( RMSE = sqrt( sum ( ( diccestbi - dicc )^2 ) / mean( diccestbi ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
#' #### Ajuste sem casca
#' ##### Ajuste não linear
bisc <- nlsLM( disc ~ taper_bi( dap, ht, hih, b0, b1, b2, b3, b4, b5, b6 ),
data = arv,
start = list( b0 = 1.8,
b1 = -0.2,
b2 = -0.04,
b3 = -0.9,
b4 = -0.0006,
b5 = 0.07,
b6 = -0.14 ) )
#' ##### Erro padrão residual
bisyxsc <- summary( bisc )$sigma
bisyxpercsc <- bisyxsc/mean( arv$disc ) * 100
#' ##### RMSE
arv$discestbi <- predict( bisc )
RMSE_bisc <- arv %>%
summarise( RMSE = sqrt( sum ( ( discestbi - disc )^2 ) / mean( discestbi ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
Função de Kozak (2004)
#' #### Ajuste com casca
#' ##### Ajuste não linear
kozak <- nlsLM( dicc ~ taper_kozak( dap, ht, hih, b0, b1, b2, b3, b4, b5, b6, b7, b8, p ),
start = list( b0 = 1.00, b1 = 0.97, b2 = 0.03, b3 = 0.49, b4 = -0.87, b5 = 0.50, b6 = 3.88, b7 = 0.03, b8 = -0.19, p = 0.1 ),
data = arv,
control = nls.lm.control( maxiter = 1000, maxfev = 2000 )
)
#' ##### Erro padrão residual
kosyx <- summary( kozak )$sigma
kosyxperc <- kosyx/mean( arv$dicc ) * 100
#' ##### RMSE
arv$diccestkoz <- predict( kozak )
RMSE_kozcc <- arv %>%
summarise( RMSE = sqrt( sum ( ( diccestkoz - dicc )^2 ) / mean( diccestkoz ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
#' #### Ajuste sem casca
#' ##### Ajuste não linear
kozaksc <- nlsLM( disc ~ taper_kozak( dap, ht, hih, b0, b1, b2, b3, b4, b5, b6, b7, b8, p ),
start = list( b0 = 1.00, b1 = 0.97, b2 = 0.03, b3 = 0.49, b4 = -0.87, b5 = 0.50, b6 = 3.88, b7 = 0.03, b8 = -0.19, p = 0.1 ),
data = arv,
control = nls.lm.control( maxiter = 1000, maxfev = 2000 )
)
#' ##### Erro padrão residual
kosyxsc <- summary( kozaksc )$sigma
kosyxpercsc <- kosyxsc/mean( arv$disc ) * 100
#' ##### RMSE
arv$discestkoz <- predict( kozaksc )
RMSE_kozsc <- arv %>%
summarise( RMSE = sqrt( sum ( ( discestkoz - disc )^2 ) / mean( discestkoz ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
Função de Ormerod (1973)
#' #### Ajuste com casca
#' ##### Ajuste não linear
ormerod <- nls("dicc~dap*((ht-hi)/(ht-1.3))^b1",
arv,
start = list( b1 = 0.1 ) )
#' ##### Erro padrão residual
orsyx <- summary( ormerod )$sigma
orsyxperc <- orsyx/mean( arv$dicc ) * 100
#' ##### RMSE
arv$diccestorm <- predict( ormerod )
RMSE_ormcc <- arv %>%
summarise( RMSE = sqrt( sum ( ( diccestorm - dicc )^2 ) / mean( diccestorm ) ) ) %>%
pull( RMSE ) %>%
round( 2 )
#' #### Ajuste sem casca
ormerodsc <- nls("disc~dap*((ht-hi)/(ht-1.3))^b1",
arv,
start = list( b1 = 0.1 ) )
#' ##### Erro padrão residual
orsyxsc <- summary( ormerodsc )$sigma
orsyxpercsc <- orsyxsc/mean( arv$disc ) * 100
#' ##### RMSE
arv$discestorm <- predict( ormerodsc )
RMSE_ormsc <- arv %>%
summarise( RMSE = sqrt( sum ( ( discestorm - disc )^2 ) / mean( discestorm ) ) ) %>%
pull( RMSE ) %>%
round( 2 )| Modelos | Erro Padrão Residual (% com casca) | Erro Padrão Residual (% sem casca) |
|---|---|---|
| Schöepfer (1966) | 5.17 | 5.22 |
| Bi (2000) | 4.24 | 6.56 |
| Kozak (2004) | 4.24 | 4.39 |
| Ormerod (1973) | 6.66 | 7.41 |
Tabela 1.
Quanto menor o erro padrão residual mais distante do ideal o modelo se encontra, ou seja, os resíduos estão mais distantes do valor real. Desse modo, comparando os modelos de forma isolada a partir dessa métrica estatística os modelos de Kozak e Bi obtiveram os melhores resultados para as estimativas de diâmetro com casca. Já para as estimativas de diâmetro sem casca, o modelo de Kozak foi o melhor. No entanto, vale ressaltar que uma métrica não é suficiente para escolhermos o melhor modelo de acordo com a nossa base, é necessário analisarmos outros fatores, que serão discutidos nesse relatório.
| Modelos | RMSE (com casca) | RMSE (sem casca) |
|---|---|---|
| Schöepfer (1966) | 8.51 | 8.18 |
| Bi (2000) | 7.49 | 11.05 |
| Kozak (2004) | 7.47 | 7.47 |
| Ormerod (1973) | 11.84 | 12.47 |
Tabela 2.
A raiz quadrada do erro médio, RMSE (Root Mean Squared Error) é a raiz quadrada do erro médio, que nos auxilia a verificar a acurácia dos modelos, penalizando erros maiores, ou seja, é uma métrica mais sensível aos erros devido ao processo de quadratura que a produziu. Sendo assim, analisando os resultados dessa métrica quanto maior o RMSE pior será o modelo, logo, os modelos de Kozak e Bi novamente foram os melhores para as estimativas dos diâmetros com casca. E para as estimativas dos diâmetros sem casca, o modelo de Kozak também foi o melhor, se destacando entre os modelos testados para duas métricas estatísticas.
## Resíduos
#' ##### 5° grau
arv$res5grau <- residuals( ajlin ) * arv$dap
arv$res5grausc <- residuals( ajlinsc ) * arv$dap
#' ##### methol
arv$resbi <- residuals( bi )
arv$resbisc <- residuals( bisc )
#' ##### kozak
arv$reskoz <- residuals( kozak )
arv$reskozsc <- residuals( kozaksc )
#' ##### ormerod
arv$resorm <- residuals( ormerod )
arv$resormsc <- residuals( ormerodsc )
par( mfrow = c( 2, 2 ),
pch = 20,
bg = "#171717",
col.axis = "white",
col.lab = "white",
col.main = "white",
fg = "white" )
fBasics::qqnormPlot( ( arv$res5grausc ), title = FALSE, main = "Schöepfer" )
fBasics::qqnormPlot( ( arv$resbisc ), title = FALSE, main = "Bi" )
fBasics::qqnormPlot( ( arv$reskozsc ), title = FALSE, main = "Kozak" )
fBasics::qqnormPlot( ( arv$resormsc ), title = FALSE, main = "Ormerod" )
Figura 2. Gráfico de quantis de todos modelos testados - Diâmetros com casca.
par( mfrow = c( 2, 2 ),
pch = 20,
bg = "#171717",
col.axis = "white",
col.lab = "white",
col.main = "white",
fg = "white" )
graficos_res <- function( hi, res, modelo, xlab = "hi(m)" ) {
plot( hi, res,
pch = 20,
ylim = c( -5, 5 ),
ylab = "Resíduo (m³)",
main = modelo )
abline( h = 0, col = "red" )
}
graficos_res( arv$hi, arv$res5grau, "Schöepfer" )
graficos_res( arv$hi, arv$resbi, "Bi" )
graficos_res( arv$hi, arv$reskoz, "Kozak" )
graficos_res( arv$hi, arv$resorm, "Ormerod" )
Figura 3. Gráfico de resíduos de todos modelos testados - Diâmetros com casca.
Quando analisados os gráficos de resíduos (Figura 3) e os gráficos de normalidade dos resíduos QQ-plot (Figura 2), foi possível notar que os melhores modelos que apresentaram distribuições mais agrupadas em torno do zero, foram os modelos de Kozak e Bi. Sendo o modelo de Schöepfer aparentemente mais próximo do ideal nas ponteiras das árvores, no entanto esse modelo apresentou erros maiores (mais distantes de zero), principalmente na base das árvores, em relação aos modelos de Kozak e Bi.
par( mfrow = c( 2, 2 ),
pch = 20,
bg = "#171717",
col.axis = "white",
col.lab = "white",
col.main = "white",
fg = "white" )
fBasics::qqnormPlot((arv$res5grausc),title = FALSE, main = "Schöepfer")
fBasics::qqnormPlot((arv$resbisc),title = FALSE, main = "Bi")
fBasics::qqnormPlot((arv$reskozsc),title = FALSE, main = "Kozak")
fBasics::qqnormPlot((arv$resormsc),title = FALSE, main = "Ormerod")
Figura 4. Gráfico de quantis de todos modelos testados - Diâmetros sem casca.
par( mfrow = c( 2, 2 ),
pch = 20,
bg = "#171717",
col.axis = "white",
col.lab = "white",
col.main = "white",
fg = "white" )
graficos_res( arv$hi, arv$res5grausc, "Schöepfer" )
graficos_res( arv$hi, arv$resbisc, "Bi" )
graficos_res( arv$hi, arv$reskozsc, "Kozak" )
graficos_res( arv$hi, arv$resormsc, "Ormerod" )
Figura 5. Gráfico de resíduos de todos modelos testados - Diâmetros sem casca.
A partir da Figura 4, a qual representa o gráfico de quantis é possível observar que os resíduos afirmam a normalidade. Além disso, com a análise dos gráficos de resíduos dispostos na Figura 5, foi possível notar que os melhores modelos, os quais apresentaram as melhores dispersões de resíduos devido às suas distribuições serem mais agrupadas em torno do zero, foram os modelos de Kozak e Schöepfer. Em contraponto, os demais modelos apresentaram erros maiores (mais distantes de zero), principalmente na base das árvores, apresentando alta heterogeneidade.
Desse modo, com o que foi observado nas análises gráficas e estatísticas o melhor modelo testado, e que será usado para os cálculos dos produtos para produção madeireira é o modelo de Kozak (2004).
# coeficientes para obtenção dos volumes
coef_kozak <- coef( kozak )
coef_kozaksc <- coef( kozaksc )
###### volume produto com casca --------------------------------------------------------------------
# A tabela de sortimentos deve conter cinco colunas, que podem ter qualquer nome, mas devem estar na seguinte ordem:
#diâmetro da tora na ponta fina, em centímetros;
#comprimento mínimo da tora (líquido), em metros; comprimento máximo da tora (líquido), em metros;
#e a perda resultante do traçamento de cada tora, em centímetros. Um ponto de atenção é que as linhas da tabela representam os produtos da madeira, e estes devem estar ordenados do mais valioso para o menos valioso, de modo que o algoritmo de extração de toras priorize aos produtos de maior valor comercial.
# afilamento usando timbeR::poly5_logs
fustes <- read.csv2( "fustes.csv" )
sortimentoscc <- data.frame(
NAME = c( "cc" ),
SED = c( 7.5 ),
min = c( 2 ),
max = c( 3.3 ),
LOSS = c( 10 )
)
# testando para produtos com casca
vticc <- mapply( kozak_logs,
dbh = fustes$cap/pi,
h = fustes$htre,
MoreArgs = list(
coef = coef_kozak[ -10 ],
p = coef_kozak[ 10 ],
assortments = sortimentoscc,
total_volume = TRUE,
only_vol = TRUE ) )
vticc <- t( vticc )
vticc <- as.data.frame( vticc )
vticc$cc <- as.numeric( vticc$cc )
vticc$Total <- as.numeric( vticc$Total )
###### volume produto sem casca --------------------------------------------------------------------
sortimentos <- data.frame(
NAME = c( "sc" ),
SED = c( 5.5 ),
min = c( 2 ),
max = c( 5.6 ),
LOSS = c( 10 )
)
# testando para produtos sem casca
vtisc <- mapply( kozak_logs,
dbh = fustes$cap/pi,
h = fustes$htre,
MoreArgs = list(
coef = coef_kozaksc[ -10 ],
p = coef_kozaksc[ 10 ],
assortments = sortimentos,
total_volume = TRUE,
only_vol = TRUE ) )
vtisc <- t( vtisc )
vtisc <- as.data.frame( vtisc )
vtisc$sc <- as.numeric( vtisc$sc )
###### Juntando resultados para fazer inventário de cada sortimento.
fustes <- read.csv2( "fustes.csv" )
fustes$vprodisc <- vtisc$sc
fustes$vprodicc <- vticc$cc
fustes$vti <- vticc$Total##### produto 1 --------------------------------------------------------------------
parcela <- aggregate(list( vprod1 = fustes$vprodisc ),
list( talhao = fustes$talhao, parcela = fustes$parcela ), sum )
#[m³/parcela]
media <- mean(parcela$vprod1)
# area da parcela
areaparc <- mean(((fustes$raio^2)*pi))
# m³/parc para m³/ha.
parcela$vprod1_ha<-parcela$vprod1*10000/areaparc
# para calcular o volume medio m³/ha
vmed <- mean(parcela$vprod1_ha)
# variancia do volume m6/ha
vvar <- var(parcela$vprod1_ha)
# desvio padrao do volume m³/ha
vdesv <- sd(parcela$vprod1_ha)
n <- length(parcela$vprod1_ha)
talhao <- subset(fustes, duplicated( talhao ) == F,
c( "talhao", "area") )
#area da fazenda em ha
area_fazenda <- sum(talhao$area)
#numero de parcelas que cabem na fazenda
N <- area_fazenda * 10000/areaparc
erro_padrao_media <- sqrt( vvar/n ) * ( 1 - n / N )
#erro do inv
erro_inv_unid <- qt( 0.975, n - 1 ) * erro_padrao_media
#erro percentual
erro_inv_perc_prod1 <- erro_inv_unid/vmed * 100
#Intervalo de confiança em m³/ha
# cat(paste(round(vmed-erro_inv_unid,2)," a ", round( vmed+erro_inv_unid, 2 ),"m³/ha com 95% de confiança (produto 1)", sep = "" ) )
prod1_inf <- round(( vmed - erro_inv_unid) * area_fazenda, 2 )
prod1_sup <- round(( vmed + erro_inv_unid) * area_fazenda, 2 )
prod1_media <- round( vmed * area_fazenda, 0 )De acordo com os dados obtidos no processamento do inventário, pode-se esperar uma produção média do produto 1 de 98104m³, com um intervalo de confiança de 95115.29m³ a 101093.6m³. Sendo o erro percentual amostral de 3.05%.
##### produto 2 --------------------------------------------------------------------
parcela <- aggregate(list( vprod2 = fustes$vprodicc ),
list( talhao = fustes$talhao, parcela = fustes$parcela), sum )
#[m³/parcela]
media <- mean( parcela$vprod2 )
#area da parcela
areaparc <- mean( ( ( fustes$raio^2 ) * pi ) )
# m³/parc para m³/ha.
parcela$vprod2_ha <- parcela$vprod2 * 10000/areaparc
#para calcular o volume medio m³/ha
vmed <- mean( parcela$vprod2_ha )
#variancia do volume m6/ha
vvar <- var( parcela$vprod2_ha )
#desvio padrao do volume m?/ha
vdesv <- sd( parcela$vprod2_ha )
n <- length( parcela$vprod2_ha )
talhao <- subset( fustes, duplicated( talhao ) == F,
c("talhao", "area" ) )
#area da fazenda em ha
area_fazenda <- sum(talhao$area)
#numero de parcelas que cabem na fazenda
N <- area_fazenda * 10000 / areaparc
erro_padrao_media <- sqrt( vvar/ n ) * ( 1 - n/ N )
#erro do inv
erro_inv_unid <- qt( 0.975, n-1) * erro_padrao_media
#erro percentual
erro_inv_perc_prod2 <- erro_inv_unid/vmed * 100
#Intervalo de confiança em m³/ha
# cat(paste(round( vmed - erro_inv_unid, 2 ), " a ", round( vmed + erro_inv_unid, 2 ),"m³/ha com 95% de confiança (produto2)", sep = "" ) )
prod2_inf <- round(( vmed - erro_inv_unid) * area_fazenda, 2 )
prod2_sup <- round(( vmed + erro_inv_unid) * area_fazenda, 2 )
prod2_media <- round( vmed * area_fazenda, 0 )De acordo com os dados obtidos no processamento do inventário, pode-se esperar uma produção média do produto 2 de 108568m³, com um intervalo de confiança de 105160.1m³ a 111975m³. Sendo o erro percentual amostral de 3.14%.
##### vtcc --------------------------------------------------------------------
parcela <- aggregate(list( vtcc = fustes$vti ),
list( talhao = fustes$talhao, parcela = fustes$parcela ), sum )
#[m³/parcela]
media <- mean( parcela$vtcc )
#area da parcela
areaparc <- mean( ( ( fustes$raio^2 ) * pi ) )
# m³/parc para m³/ha.
parcela$vtcc_ha <- parcela$vtcc * 10000/areaparc
#para calcular o volume medio m³/ha
vmed <- mean( parcela$vtcc_ha )
#variancia do volume m6/ha
vvar <- var( parcela$vtcc_ha )
#desvio padrao do volume m³/ha
vdesv <- sd( parcela$vtcc_ha )
n <- length( parcela$vtcc_ha )
talhao <- subset( fustes, duplicated( talhao ) == F,
c("talhao", "area" ) )
#area da fazenda em ha
area_fazenda <- sum( talhao$area )
#numero de parcelas que cabem na fazenda
N <- area_fazenda * 10000 / areaparc
erro_padrao_media <- sqrt( vvar /n ) * ( 1 - n / N )
#erro do inv
erro_inv_unid <- qt( 0.975, n - 1 ) * erro_padrao_media
#erro percentual
erro_inv_perc_prod3 <- erro_inv_unid / vmed * 100
#Intervalo de confiança em m³/ha
# cat(paste(round( vmed - erro_inv_unid, 2 )," a ",round( vmed + erro_inv_unid, 2 ),
# "m³/ha com 95% de confiança para volume total com casca", sep = "" ) )
prod3_inf <- round(( vmed - erro_inv_unid) * area_fazenda, 2 )
prod3_sup <- round(( vmed + erro_inv_unid) * area_fazenda, 2 )
prod3_media <- round( vmed * area_fazenda, 0 )De acordo com os dados obtidos no processamento do inventário, pode-se esperar uma produção média do produto 3 de 116896m³, com um intervalo de confiança de 113424.6m³ a 120367.8m³. Sendo o erro percentual amostral de 2.97%.
O volume total com casca apresentou maior potencial de produção segundo as estimativas presentes neste relatório.
JÚNIOR, C. F. D. et al. Measure: sistema de coleta de informações para inventário florestal. Tecno-Lógica, 2020.