Índice de uniformidade (PV50)
O PV50 é hoje o índice mais utilizado quando queremos expressar a uniformidade de um plantio florestal. Hakamada (2012) apresentou um estudo detalhado sobre diversos índices e concluiu que o PV50 é o índice mais indicado para explicar a relação entre uniformidade, qualidade silvicultural e produtividade em plantios homogêneos de Eucalyptus.
O objetivo deste post é mostrar, passo a passo, como calcular este índice no R e fazer uma breve análise de seus resultados.
O PV50 é a porcentagem de volume acumulado das 50% menores árvores do seu conjunto de dados, considerando as falhas de plantio e árvores mortas (Hakamada et al. 2015). A expressão do índice é dada da seguinte forma:
Onde: PV50 = porcentagem acumulada do volume das 50% menores árvores plantadas; V = volume da árvore i; n = número de árvores plantadas ordenadas (da menor para a maior).
Primeiro vamos entender os cálculos do índice, considerando apenas 10 árvores hipotéticas com 0,1 metros cúbicos de volume.
# carrega os pacotes necessários
if (!require("pacman")) install.packages("pacman")
pacman::p_load(readr, dplyr, ggplot2, forcats)# exemplo com número par
arv10 <- rep(0.1, 10)
str(arv10)## num [1:10] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1Este é o referencial teórico de uniformidade, todas as árvores do mesmo tamanho. Sem precisar fazer conta, sabemos que o volume das 50% menores árvores é igual a 50% do volume total, o que equivale a um PV50 = 50.
# identifica a metade do numero de árvores
metade <- length(arv10)/2
metade## [1] 5# soma todas as árvores
soma_todas <- sum(arv10, na.rm = TRUE)
soma_todas## [1] 1# soma o valor de metade das árvores em ordem crescente
soma_metade <- sum(sort(arv10)[1:metade], na.rm = TRUE)
soma_metade## [1] 0.5# calcula o PV50
PV50 <- soma_metade / soma_todas * 100
PV50## [1] 50Agora vamos simular 11 árvores com o mesmo volume, veja o que acontece.
# exemplo com número impar
arv11 <- rep(0.1, 11)
str(arv11)## num [1:11] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 ...# metade do numero de árvores
metade <- length(arv11)/2
metade## [1] 5.5# soma todas as árvores
soma_todas <- sum(arv11, na.rm = TRUE)
soma_todas## [1] 1.1# soma o valor de metade das árvores em ordem crescente
soma_metade <- sum(sort(arv11)[1:metade], na.rm = TRUE)
soma_metade## [1] 0.5# calcula o PV50
PV50 <- soma_metade / soma_todas * 100
PV50## [1] 45.45455O resultado deveria ser 50, mas como o número de árvores é ímpar, o R arredonda a posição 5,5 para 5 e pega até a quinta árvore no momento em que queremos somar as 50% menores. Para contornar isso, vamos calcular a soma das 50% menores árvores de uma forma diferente. Primeiro calculamos a soma acumulada e depois extraímos a média (semelhante ao modo de se calcular uma mediana).
# vetor de soma acumulada
soma_acumulada <- cumsum(sort(arv11))
soma_acumulada## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1# soma a metade do vetor de soma acumulada
soma_metade <- mean(soma_acumulada[metade + 0L:1L], na.rm = TRUE)
soma_metade## [1] 0.55# calcula o PV50
PV50 <- soma_metade / soma_todas * 100
PV50## [1] 50Agora que a questão do número de árvores foi superada, podemos incluir árvores mortas, o que equivale a elementos do tipo NA no R. Veja que o resultado não está consistente pois a soma_acumulada ignorou as árvores mortas.
# exemplo com valores perdidos
arv11_na <- rep(0.1, 11)
arv11_na[c(3,4)] <- NA
str(arv11_na)## num [1:11] 0.1 0.1 NA NA 0.1 0.1 0.1 0.1 0.1 0.1 ...# metade do numero de árvores
metade <- length(arv11_na)/2
metade## [1] 5.5# soma todas as árvores
soma_todas <- sum(arv11_na, na.rm = TRUE)
soma_todas## [1] 0.9# vetor de soma acumulada
soma_acumulada <- cumsum(sort(arv11_na))
soma_acumulada## [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9# soma a metade do vetor de soma acumulada
soma_metade <- mean(soma_acumulada[metade + 0L:1L], na.rm = TRUE)
soma_metade## [1] 0.55# calcula o PV50
PV50 <- soma_metade / soma_todas * 100
PV50## [1] 61.11111Para corrigir este o erro, temos de incluir manualmente as árvores mortas na sequência. Veja que agora o resultado está de acordo com o esperado.
# vetor de valores perdidos
mortas <- arv11_na[is.na(arv11_na)]
mortas## [1] NA NA# metade do numero de árvores
metade <- length(arv11_na)/2
metade## [1] 5.5# soma todas as árvores
soma_todas <- sum(arv11_na, na.rm = TRUE)
soma_todas## [1] 0.9# vetor de soma acumulada com valores perdidos
soma_acumulada <- c(mortas, cumsum(sort(arv11_na)))
soma_acumulada## [1] NA NA 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9# soma a metade do vetor de soma acumulada
soma_metade <- mean(soma_acumulada[metade + 0L:1L], na.rm = TRUE)
soma_metade## [1] 0.35# calcula o PV50
PV50 <- soma_metade / soma_todas * 100
PV50## [1] 38.88889Agora eu vou dar aquele passo mágico dos livros de matemática e física, em que o autor diz “é fácil notar que o resultado leva a …” e apresentar uma função que lida com as questões que mostramos acima e retorna o PV50 do nosso conjunto de dados de forma correta.
pv50 <- function(x) {
mortas <- x[is.na(x)]
metade <- length(x)/2
soma_todas <- sum(x, na.rm = TRUE)
soma_acumulada <- c(mortas, cumsum(sort(x)))
if (metade%%2L == 1L)
soma_metade <- mean(soma_acumulada[metade], na.rm = TRUE)
else
soma_metade <- mean(soma_acumulada[metade + 0L:1L], na.rm = TRUE)
z <- soma_metade / soma_todas * 100
return(z)
}Podemos rapidamente verificar se os resultados estão consistentes fazendo alguns testes.
a <- rep(10, 10)
str(a)## num [1:10] 10 10 10 10 10 10 10 10 10 10pv50(a) # Ok!## [1] 50a1 <- rep(10 ,11)
str(a1)## num [1:11] 10 10 10 10 10 10 10 10 10 10 ...pv50(a1) # Ok!## [1] 50b <- a
b[c(3, 7)] <- NA
str(b)## num [1:10] 10 10 NA 10 10 10 NA 10 10 10pv50(b) # Ok!## [1] 37.5b1 <- a1
b1[c(3, 7)] <- NA
str(b1)## num [1:11] 10 10 NA 10 10 10 NA 10 10 10 ...pv50(b1) # Ok!## [1] 38.88889Boa, já temos uma função para calcular o PV50 e podemos aplicá-la em um conjunto de dados para podermos interpretar. Utilizaremos mais uma vez os dados do Projeto TUME, referente ao TUME 134 plantado em Piracicaba-SP. O volume individual foi calculado arbitrariamente utilizando o fator de forma 0,5.
# importa o arquivo tume_55.csv
dados <- read_csv2(
"https://github.com/italocegatta/italocegatta.github.io_source/raw/master/content/dados/tume_55.csv"
)
glimpse(dados)## Observations: 1,222
## Variables: 7
## $ Esp <chr> "E_camaldulensis", "E_camaldulensis", "E_camaldulensis...
## $ I_meses <dbl> 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34, 34...
## $ Parc_m2 <dbl> 288, 288, 288, 288, 288, 288, 288, 288, 288, 288, 288,...
## $ N_arv <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,...
## $ DAP_cm <dbl> 5.411268, 12.254931, 3.978874, 6.429860, 9.676621, 5.6...
## $ H_m <dbl> 7.651490, 11.424046, 5.909205, 8.572873, 10.498957, 7....
## $ Vol <dbl> 0.008798406, 0.067375427, 0.003673747, 0.013918399, 0....Iremos calcular o PV50 e o volume por hectare para cada fator Esp e I_meses e em seguida ordenar as espécies pelo PV50.
# agrupa os dados em função de espécie e idade para
# calcular o pv50 e o volume
dados_pv50 <- dados %>%
group_by(Esp, I_meses) %>%
summarise(
Parc_m2 = mean( Parc_m2),
PV50 = pv50(Vol),
Vol_ha = sum(Vol, na.rm = TRUE) * (10000/Parc_m2)
) %>%
ungroup() %>%
# ordena o fator de espécies de forma decrescente em função do pv50
mutate(Esp = fct_reorder(Esp, -PV50))
dados_pv50## # A tibble: 20 x 5
## Esp I_meses Parc_m2 PV50 Vol_ha
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 E_camaldulensis 34 288 21.4 47.5
## 2 E_camaldulensis 46 288 14.3 79.1
## 3 E_camaldulensis 60 600 13.2 110.
## 4 E_camaldulensis 85 288 12.7 204.
## 5 E_citriodora 34 288 17.8 46.6
## 6 E_citriodora 46 288 16.4 85.0
## 7 E_citriodora 60 600 12.5 97.0
## 8 E_citriodora 85 288 12.7 205.
## 9 E_dunnii 34 288 28.6 103.
## 10 E_dunnii 46 288 29.2 161.
## 11 E_dunnii 60 600 27.5 198.
## 12 E_dunnii 85 288 26.4 350.
## 13 E_paniculata 34 288 27.4 46.4
## 14 E_paniculata 46 288 24.1 84.4
## 15 E_paniculata 60 600 19.6 115.
## 16 E_paniculata 85 288 18.2 195.
## 17 E_urophylla_grandis 34 288 26.3 85.7
## 18 E_urophylla_grandis 46 288 24.2 157.
## 19 E_urophylla_grandis 60 600 20.3 217.
## 20 E_urophylla_grandis 85 288 17.9 277.Para entendermos os dados, vamos primeiro ver o crescimento em volume de cada espécies em função do tempo (Figura 1). Note que E. dunnii e E. urophylla x E. grandis tinham crescimento muito parecido até os 60 meses de idade.
ggplot(dados_pv50, aes(I_meses, Vol_ha, color = Esp)) +
geom_point() +
geom_line() +
labs(
color = "Espécies",
x = "Idade (meses)",
y = Volume~m^3~ha^-1
) +
scale_color_brewer(palette = "Set1") +
theme_bw(16) +
theme(legend.justification = "top")
Figura 1: Crescimento em volume por hectare em função da idade.
Agora podemos construir um gráfico que relaciona o PV50 e a idade (Figura 2). A interpretação do índice é simples, o PV50 representa a porcentagem em volume que as 50% menores árvores contribuem para o volume total. Em nossos dados, E. dunnii, aos 85 meses de idade, tem um PV50 de aproximadamente 26. Isso quer dizer que aos 7 anos, as 50% menores árvores da parcela de E. dunnii representam apenas 26% do volume total. Ou seja, 50% das árvores contribuem muito pouco para o volume total da parcela e isso tem um impacto direto na produtividade.
ggplot(dados_pv50, aes(I_meses, PV50, color = Esp)) +
geom_point() +
geom_line() +
labs(color = "Espécies", x = "Idade (meses)", y = "PV50") +
scale_color_brewer(palette = "Set1") +
scale_y_continuous(breaks = seq(10, 30, 2)) +
theme_bw(16) +
theme(legend.justification = "top")
Figura 2: Variação do PV50 por espécies em função da idade.
A Figura 3 mostra claramente a relação direta que há entre produção de madeira e a uniformidade ao longo do crescimento da floresta. Note também que na medida em que a idade avança, a uniformidade diminui, uma vez que a dominância das árvores maiores sobre as menores fica cada vez mais forte.
ggplot(dados_pv50, aes(Vol_ha, PV50)) +
geom_point(aes(color = factor(I_meses))) +
geom_smooth(method = "lm", formula = y ~x, se = FALSE) +
facet_wrap(~Esp, dir = "v") +
labs(color = "Idade (meses)", x = Volume~m^3~ha^-1, y = "PV50") +
scale_color_brewer(palette = "Dark2") +
theme_bw(16) +
theme(legend.justification = "top")
Figura 3: Relação entre o PV50 e volume por hectare em função da idade.
Por fim, para colocar tudo em um só gráfico, podemos adicionar ao gráfico de crescimento em volume a informação do PV50 para evidenciar que as espécies mais produtivas tem PV50 elevado e que este índice consegue explicar muito bem essa relação (Figura 4).
Um comentário interessante é que dentre as espécies que estamos estudando, todas são de origem seminal, com exceção do E. dunnii, que é um clone. Este fator explica sua produtividade e alta homogeneidade, principalmente frente ao hibrido de E. urophylla x E. grandis, que é seu concorrente direto. Quando estivermos analisando dados de plantios clonais, o PV50 vai expressar a qualidade silvicultural do plantio, uma vez que a base genética é a mesma em todas as plantas.
ggplot(dados_pv50, aes(I_meses, Vol_ha, color = Esp)) +
geom_point(aes(size = PV50), alpha = 0.4) +
geom_line() +
labs(
color = "Espécies",
x = "Idade (meses)",
y = Volume~m^3~ha^-1
) +
scale_color_brewer(palette = "Set1") +
theme_bw(16) +
theme(legend.justification = "top")
Figura 4: Crescimento do volume em função da idade, com informação do PV50 no tamanho do ponto.
Caso tenha alguma dúvida ou sugestão sobre o post, fique à vontade para fazer um comentário ou me contatar por E-mail.
sessioninfo::session_info(c("readr", "dplyr", "ggplot2", "forcats"))## - Session info ----------------------------------------------------------
## setting value
## version R version 3.5.3 (2019-03-11)
## os Windows 10 x64
## system x86_64, mingw32
## ui RTerm
## language (EN)
## collate Portuguese_Brazil.1252
## ctype Portuguese_Brazil.1252
## tz America/Sao_Paulo
## date 2019-08-25
##
## - Packages --------------------------------------------------------------
## package * version date lib source
## assertthat 0.2.1 2019-03-21 [1] CRAN (R 3.5.3)
## BH 1.69.0-1 2019-01-07 [1] CRAN (R 3.5.2)
## cli 1.1.0 2019-03-19 [1] CRAN (R 3.5.3)
## clipr 0.7.0 2019-07-23 [1] CRAN (R 3.5.3)
## crayon 1.3.4 2017-09-16 [1] CRAN (R 3.5.1)
## digest 0.6.20 2019-07-04 [1] CRAN (R 3.5.3)
## dplyr * 0.8.3 2019-07-04 [1] CRAN (R 3.5.3)
## ellipsis 0.2.0.1 2019-07-02 [1] CRAN (R 3.5.3)
## fansi 0.4.0 2018-10-05 [1] CRAN (R 3.5.1)
## forcats * 0.4.0 2019-02-17 [1] CRAN (R 3.5.2)
## ggplot2 * 3.2.1 2019-08-10 [1] CRAN (R 3.5.3)
## glue 1.3.1 2019-03-12 [1] CRAN (R 3.5.3)
## gtable 0.3.0 2019-03-25 [1] CRAN (R 3.5.3)
## hms 0.5.0 2019-07-09 [1] CRAN (R 3.5.3)
## labeling 0.3 2014-08-23 [1] CRAN (R 3.5.0)
## lazyeval 0.2.2 2019-03-15 [1] CRAN (R 3.5.3)
## magrittr 1.5 2014-11-22 [1] CRAN (R 3.5.1)
## MASS 7.3-51.1 2018-11-01 [2] CRAN (R 3.5.3)
## Matrix 1.2-17 2019-03-22 [1] CRAN (R 3.5.3)
## mgcv 1.8-28 2019-03-21 [1] CRAN (R 3.5.3)
## munsell 0.5.0 2018-06-12 [1] CRAN (R 3.5.1)
## nlme 3.1-137 2018-04-07 [2] CRAN (R 3.5.3)
## pillar 1.4.2 2019-06-29 [1] CRAN (R 3.5.3)
## pkgconfig 2.0.2 2018-08-16 [1] CRAN (R 3.5.1)
## plogr 0.2.0 2018-03-25 [1] CRAN (R 3.5.1)
## plyr 1.8.4 2016-06-08 [1] CRAN (R 3.5.1)
## purrr 0.3.2 2019-03-15 [1] CRAN (R 3.5.3)
## R6 2.4.0 2019-02-14 [1] CRAN (R 3.5.2)
## RColorBrewer 1.1-2 2014-12-07 [1] CRAN (R 3.5.0)
## Rcpp 1.0.2 2019-07-25 [1] CRAN (R 3.5.3)
## readr * 1.3.1 2018-12-21 [1] CRAN (R 3.5.2)
## reshape2 1.4.3 2017-12-11 [1] CRAN (R 3.5.1)
## rlang 0.4.0 2019-06-25 [1] CRAN (R 3.5.3)
## scales 1.0.0 2018-08-09 [1] CRAN (R 3.5.1)
## stringr 1.4.0 2019-02-10 [1] CRAN (R 3.5.2)
## tibble 2.1.3 2019-06-06 [1] CRAN (R 3.5.3)
## tidyselect 0.2.5 2018-10-11 [1] CRAN (R 3.5.1)
## vctrs 0.2.0 2019-07-05 [1] CRAN (R 3.5.3)
## viridisLite 0.3.0 2018-02-01 [1] CRAN (R 3.5.1)
## withr 2.1.2 2018-03-15 [1] CRAN (R 3.5.1)
##
## [1] C:/Users/Italo/Documents/R/win-library/3.5
## [2] C:/Program Files/R/R-3.5.3/libraryReferências
Hakamada, Rodrigo Eiji. 2012. “Uso do inventário florestal como ferramenta de monitoramento da qualidade silvicultura em povoamentos clonais de Eucalyptus.” PhD thesis, Piracicaba: Universidade de São Paulo; Biblioteca Digital de Teses e Dissertações da Universidade de São Paulo. https://doi.org/10.11606/D.11.2012.tde-05072012-100431.
Hakamada, Rodrigo Eiji, José Luiz Stape, Cristiane Camargo Zani de Lemos, Adriano Emanuel Amaral Almeida, and Luis Fernando Silva. 2015. “Uniformidade entre árvores durante uma rotação e sua relação com a produtividade em Eucalyptus clonais.” CERNE 21 (3): 465–72. https://doi.org/10.1590/01047760201521031716.
