El objetivo es procesar los datos registrados en las UM relevadas en el marco del curso: “Teoría del Muestreo e Inventarios Forestales” dictado por la Facultad de Cs. Forestales de la Universidad Nacional de Misiones, y estimar diferentes variables indicadoras a distintos niveles: individual (árbol), parcela (muestra), lote o estrato.

Cargamos los paquetes necesarios:

library(tidyverse)
library(BIOMASS)
library(vegan)
library(BiodiversityR)
library(tidymodels)
library(MASS)
library(fitdistrplus)

Cargamos las bases. Revisar cómo están guardadas y modificar el código según su configuración. “sep” es el separador de valores (columnas) y “dec” el caracter usado para indicar números decimales.

arboles <- read.delim("./bases/arboles.csv", sep = ";", dec = ",")
renovales <- read.delim("./bases/renovales.csv", sep = ";", dec = ",")
sotobosque  <- read.delim("./bases/sotobosque.csv", sep = ";", dec = ",")
especies <- read.delim("./bases/especies_listado_completa.csv", sep=";")
matlenioso <- read.delim("./bases/materiallenioso.csv", sep = ";", dec = ",")

Estrato Arboreo

Miramos un poco la base y el tipo de datos

glimpse(arboles)
## Rows: 495
## Columns: 10
## $ PARQUE      <chr> "RUMG", "RUMG", "RUMG", "RUMG", "RUMG", "RUMG", "RUMG", "R…
## $ PARCELA     <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2…
## $ N_ARBOL     <int> 1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, …
## $ PARCELA_CON <chr> "A", "A", "B", "A", "A", "B", "A", "A", "B", "B", "B", "B"…
## $ ESPECIE     <chr> "PARRIGID", "PARRIGID", "SIDOBTUS", "OCOPUBER", "PARRIGID"…
## $ ESTADO      <int> 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 1, 4, 1, 1, 0, 1, 1, 4, 1, 2…
## $ ALTURA_TOT  <dbl> 16, 16, 7, 17, 21, 16, 15, 18, 10, 11, 12, 7, 22, 14, 4, 1…
## $ N_FUSTES    <int> 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ DAP         <dbl> 24.50986, 38.51550, 19.57606, 56.97747, 29.92113, 17.50704…
## $ COB_LIAN    <int> 1, 1, 2, 1, 1, 3, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…

Descriptiva

summary(arboles)
##     PARQUE             PARCELA         N_ARBOL      PARCELA_CON       
##  Length:495         Min.   : 1.00   Min.   : 1.00   Length:495        
##  Class :character   1st Qu.: 5.00   1st Qu.: 7.00   Class :character  
##  Mode  :character   Median :10.00   Median :13.50   Mode  :character  
##                     Mean   :10.09   Mean   :14.85                     
##                     3rd Qu.:15.00   3rd Qu.:21.00                     
##                     Max.   :20.00   Max.   :45.00                     
##                                     NA's   :19                        
##    ESPECIE              ESTADO        ALTURA_TOT       N_FUSTES     
##  Length:495         Min.   :0.000   Min.   : 1.50   Min.   : 1.000  
##  Class :character   1st Qu.:1.000   1st Qu.:10.00   1st Qu.: 1.000  
##  Mode  :character   Median :1.000   Median :14.00   Median : 1.000  
##                     Mean   :1.101   Mean   :14.14   Mean   : 1.046  
##                     3rd Qu.:1.000   3rd Qu.:18.00   3rd Qu.: 1.000  
##                     Max.   :4.000   Max.   :28.00   Max.   :11.000  
##                                                                     
##       DAP             COB_LIAN    
##  Min.   :  4.138   Min.   :0.000  
##  1st Qu.: 19.417   1st Qu.:0.000  
##  Median : 24.828   Median :1.000  
##  Mean   : 30.000   Mean   :1.578  
##  3rd Qu.: 36.128   3rd Qu.:3.000  
##  Max.   :140.375   Max.   :4.000  
##

El siguiente bloque es para unificar categorías utilizadas

arboles %>% 
  mutate(COB_LIAN = case_when(COB_LIAN == 0 ~ 1, 
                              T ~ COB_LIAN)) -> arboles
  • Cantidad de árboles vivos:
arboles %>%
  filter(ESTADO != 0) %>% 
  count()
  • Arboles vivos por parcela:
arboles %>%
  filter(ESTADO != 0) %>% 
  group_by(PARCELA) %>% 
  count() %>% 
  arrange(n)
  • Arboles muertos:
arboles %>%
  filter(ESTADO == 0) %>% 
  count()
  • Arboles muertos por parcela:
arboles %>%
  filter(ESTADO == 0) %>% 
  group_by(PARCELA) %>% 
  count() %>% 
  arrange(n)
  • Promedio por parcela:
465/20
## [1] 23.25
30/20
## [1] 1.5

Total arboles = 495, promedio 23.25/parcela, rango 11 a 38 muertos = 30 , promedio 1.5/parcela

Podemos ver la distribución de algunas variables. Para eso nos quedamos solo con árboles DAP > 20 cm

arboles %>% 
  filter(PARCELA_CON == "A") -> arbolesA
  • DAP
ggplot(arboles, aes(DAP)) +
  geom_density() +
  facet_grid(.~PARCELA_CON, scales = "free_x") +
  theme_bw()

ggplot(arbolesA, aes(DAP)) +
  geom_density() +
  theme_bw()

  • Altura
ggplot(arbolesA, aes(ALTURA_TOT)) +
  geom_density() +
  theme_bw()

  • Estado de los árboles
ggplot(arboles, aes(x = ESTADO)) +
  geom_bar(aes(y = (after_stat(count))/sum(after_stat(count)))) + ylab("Porcentaje") +
  theme_bw()

  • Carga de lianas
ggplot(arboles, aes(x = COB_LIAN)) +
  geom_bar(aes(y = 100*(after_stat(count))/sum(after_stat(count)))) + ylab("Porcentaje") +
  facet_grid(.~PARCELA_CON) + 
  theme_bw()

Atributos comunitarios

Vamos a calcular valor de importancia por especie, riqueza y diversidad por parcela usando los paquetes ‘vegan’ y ‘BiodiversityR’.

arboles %>% 
  filter(ESTADO > 0) %>%     # es lo mismo que hacer ESTADO !=0? por qué?
  group_by(PARCELA, PARCELA_CON, ESPECIE) %>% 
  #calculo las variables para cada subparcela, escalo y después las sumo
  summarise(dens = case_when(PARCELA_CON =="A" ~n()*10,
                           T ~ n()*40),
            ab = case_when(PARCELA_CON =="A" ~sum(pi*(DAP/2)^2)/1000,
                           T ~ sum(pi*(DAP/2)^2)/250)) %>% 
  group_by(PARCELA, ESPECIE) %>% 
  summarise(cant = sum(dens),
            abt = sum(ab)) -> datos_abund
## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## `summarise()` has grouped output by 'PARCELA', 'PARCELA_CON', 'ESPECIE'. You
## can override using the `.groups` argument.
## `summarise()` has grouped output by 'PARCELA'. You can override using the
## `.groups` argument.
abund <- as.data.frame(datos_abund)

#vemos como quedó:
head(abund)

El valor de importancia (%) es una medida que ponera la frecuencia, densidad y dominancia (área basal) de cada especie.

ivg <- importancevalue(abund, 
                       site = "PARCELA",
                       species='ESPECIE',
                       count='cant',
                       basal='abt',
                       factor = "")

# algunas operaciones para acomodar los daatos
dfivgd <- as.data.frame(ivg) %>% 
  rownames_to_column("especie") %>% 
  arrange(density.percent) %>%           #ordeno en base a la densidad
  mutate(Especie = as_factor(especie))


ggplot(dfivgd, aes(Especie, density.percent)) +
  geom_col(fill = "lightblue") +
  theme_bw() +
  scale_x_discrete(limits = rev(levels(dfivgd$Especie))) +
  xlab("Especie") + ylab("Densidad (árboles/ha)") +
  theme(axis.text.x = element_text(angle = 90))

dfivg <- as.data.frame(ivg) %>% 
    rownames_to_column("especie") %>% 
    arrange(importance.value) %>%        #ordeno en base al valor de importancia
    mutate(Especie = as_factor(especie))

dfivg %>% 
  pivot_longer(cols = c(5:7), values_to = "Porcentaje", names_to = "variable")-> dfivg_long
ggplot(dfivg_long, aes(especie, Porcentaje, group = variable, fill = variable)) +
  geom_col() + 
  scale_x_discrete(limits = rev(levels(dfivg$Especie))) +
  theme_bw() +
  scale_fill_discrete(name = "Variable", labels = c("Densidad", "Dominancia", "Frecuencia")) +
  theme(axis.text.x = element_text(angle = 90))

Algunas operaciones para llevar la base a formato de matriz donde cada columna es una especie, las filas corresponden a las parcelas y cada celda indica la densidad de árboles (por ha.).

datos_abund %>% 
  dplyr::select(-abt) %>% 
  pivot_wider(names_from = ESPECIE, values_from = cant, values_fill = 0)  -> mx_sps

Podemos estimar riqueza y diversidad por parcela. Primero creo una base de datos para almacenar las estimaciones

atr_comunitarios <- data.frame(PARCELA = as_factor(c(1:20)))
atr_comunitarios$div  <-  diversity(mx_sps[, 2:68])

atr_comunitarios$riqueza  <-  specnumber(mx_sps[, 2:68])

atr_comunitarios

Estimacion biomasa

Para estimar biomasa necesitamos el DAP, altura y densidad de madera de cada árbol. Para eso usamos el paquete ‘BIOMASS’ que es un paquete creado para estimar biomasa aérea en bosques tropicales en base a la ecuacion alométrica propuesta por Chave et al. (2014) y contiene una variedad de funciones muy útiles:

?`BIOMASS-package`

Nosotros vamos a usar el paquete para obtener la densidad de madera de cada especie de la base global de densidad de madera (Chave et al. 2009, Zanne et al. 2009).

?getWoodDensity

La función requiere que haya una columna con el genero y otra la especie. En nuestro caso tenemos los datos cargados directamente con el ‘CODIGO’ de la especie, así que necesitamos armar un diccionario de codigos de especies y preparar el formato:

especies %>% 
  separate(`NOMBRE.CIENTIFICO`, c("genero", "especie"), sep = " ", remove = F) %>%   
  dplyr::select(-`NOMBRE.VULGAR`) %>% 
  rename(ESPECIE = `CODIGO`) -> dicespecies
## Warning: Expected 2 pieces. Additional pieces discarded in 94 rows [3, 4, 5, 6, 7, 8, 9,
## 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23, ...].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 1 rows [77].
dicespecies

Agregamos las columnas “genero” y “especie” a nuestra base:

arboles %>% 
  left_join(dicespecies, by = "ESPECIE") -> dfarboles

Ahora si tenemos toda la información en la misma base y podemos correr la función para obtener la densidad de madera:

getWoodDensity(dfarboles$genero, dfarboles$especie,
               region = "World", verbose = T) -> dm_biomass
## The reference dataset contains 16467 wood density values
## Your taxonomic table contains 68 taxa
head(dm_biomass)
dm_biomass %>% 
  dplyr::select(2:4) -> dm_biomass_sub

dfarboles %>% 
  left_join(dm_biomass_sub, by = c("genero" = "genus", "especie" = "species")) %>% 
  unique() -> dfarboles
## Warning in left_join(., dm_biomass_sub, by = c(genero = "genus", especie = "species")): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 1 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.

Y el último paso, estimar biomasa aérea utilizando la ecuación correspondiente. Para árboles, estamos utilizando la ecuación propuesta para bosques tropicales semi-húmedos propuesta por Chave et al. (2005).

dfarboles %>% 
  mutate(biomasa = case_when(ESPECIE %in% c("SYAROMAN", "EUTEDULI") ~ 6.6666+ (12.826*(DAP^0.5))*log(DAP),
                             TRUE~ exp(-2.977+log(meanWD*ALTURA_TOT*DAP^2))),
         AREA_BASAL_i = pi*((DAP/2)/100)^2) -> dfarboles

dfarboles

Atributos estructurales

Hasta acá venimos trabajando a escala o nivel de cada árbol. Lo siguiente es estimar las variables de interés a nivel de parcela y finalmente, hacer las inferencias a nivel de lote o estrato.

Por un lado vamos a separar los registros que corresponden a árboles vivos y árboles muertos en pie. Además, necesitamos llevar todo a la misma unidad (hectárea^-1).

dfarboles %>%
  filter(ESTADO != 0) %>%                     
  group_by(PARQUE, PARCELA, PARCELA_CON) %>%
  summarise(H_mean = mean(ALTURA_TOT), 
            AB = sum(AREA_BASAL_i),
            dens = n(),
            biomasat = sum(biomasa)/1000) %>% 
  mutate(across(c(AB, dens, biomasat), ~(case_when(PARCELA_CON == "A" ~ .*10,
                                               PARCELA_CON == "B" ~ .*40)))) -> arb_vivos 
## `summarise()` has grouped output by 'PARQUE', 'PARCELA'. You can override using
## the `.groups` argument.
arb_vivos
dfarboles %>% 
  filter(ESTADO == 0) %>% 
  group_by(PARQUE, PARCELA, PARCELA_CON) %>%
  summarise(ABm = sum(AREA_BASAL_i),
            dens_m = n(), 
            biomasat_m = sum(biomasa)/1000) %>% 
  mutate(across(c(ABm, dens_m, biomasat_m), ~(case_when(PARCELA_CON == "A" ~ .*10,
                                               PARCELA_CON == "B" ~ .*40)))) -> arb_muertos
## `summarise()` has grouped output by 'PARQUE', 'PARCELA'. You can override using
## the `.groups` argument.
arb_muertos
# ponemos todo en una misma base

arb_vivos %>% 
  left_join(arb_muertos, by = c("PARQUE", "PARCELA", "PARCELA_CON")) %>% 
  replace(is.na(.), 0) -> estructura_arborea

estructura_arborea
estructura_arborea %>% 
  pivot_longer(4:10, names_to = "variable", values_to = "valor") %>% 
  pivot_wider(names_from = PARCELA_CON, values_from = valor, values_fill = 0) %>% 
  mutate(TOTAL = case_when(variable == "H_mean" ~ A,
                            TRUE ~ A+B)) %>% 
  dplyr::select(-c(A, B)) %>% 
  pivot_wider(names_from = variable, values_from = TOTAL) %>% 
  mutate(PARCELA = as_factor(PARCELA))-> estructura_arborea_tots

estructura_arborea_tots

Podemos visualizar las variables estimadas:

ggplot(estructura_arborea_tots, aes(as_factor(PARCELA), biomasat))+
  geom_col() + xlab("Parcela") +
  theme_bw()

ggplot(estructura_arborea_tots, aes(as_factor(PARCELA), dens))+
  geom_col() + xlab("Parcela") +
  theme_bw()

Por último, le agregamos los atributos comunitarios estimados antes a nuestra base:

estructura_arborea_tots %>% 
  left_join(atr_comunitarios, by = "PARCELA") -> estructura_arborea_tots

Coberturas estratos

La cobertura de estratos se mide a lo largo de una linea transecta que recorre la subparcela A de N a S. Los datos registrados en m lineales hay que llevarlos a %.

sotobosque %>%
  arrange(PARCELA) %>%
  dplyr::select(1:5) %>% 
  pivot_longer(3:5, names_to = "Estrato", values_to = "metrosl") %>%
  mutate(Porcentaje = metrosl/35.5*100,
         PARCELA = as_factor(PARCELA)) %>% 
  pivot_wider(id_cols = PARCELA, names_from = Estrato, values_from = Porcentaje) -> sotobosque_porc
colnames(sotobosque_porc) <- c("PARCELA", "COB_ARBOREA_7", "COB_ARBOREA_3A7", "COB_BAMBU")

sotobosque_porc

Agregamos los datos a la base

estructura_arborea_tots %>% 
  left_join(sotobosque_porc, by = "PARCELA") -> estructura_arborea_tots

Podemos visualizar la cobertura por estrato.

sotobosque_porc %>% 
  pivot_longer(2:4, names_to = "Estrato", values_to = "Porcentaje") %>% 
  mutate(Estrato = as_factor(Estrato),
         Parcela = as_factor(PARCELA)) -> sotol

ggplot(sotol, aes(Parcela, Porcentaje)) +
  geom_col() + 
  facet_wrap(vars(Estrato), ncol = 1) +
  theme_bw()

Regeneración

Al igual que con los individuos adultos, con los renovales podemos trabajar a nivel de especie o en general. Se pueden calcular los mismos atributos comunitarios que calculamos, aunque no el valor de importancia porque no tenemos medidas de dominancia.

renovales %>% 
  mutate(PARCELA = as.numeric(str_sub(PARCELA, start = 5, end = -1L))) %>%
  group_by(PARCELA, CLASE, ESPECIE) %>%
  summarise(Renovales = sum(CANTIDAD)*200) %>% 
  pivot_wider(names_from = CLASE, values_from = Renovales, values_fill = 0) %>% 
  rename(R1 = 3, R2 = 4) %>% 
  mutate(TOTAL = R1 + R2,
         PARCELA = as_factor(PARCELA)) -> renovs_tot_especies
## `summarise()` has grouped output by 'PARCELA', 'CLASE'. You can override using
## the `.groups` argument.
renovs_tot_especies
renovs_tot_especies %>%
  group_by(PARCELA) %>%
  summarise(TOTAL = sum(TOTAL),
            R1 = sum(R1),
            R2 = sum(R2)) %>% 
  arrange(as.numeric(PARCELA)) -> renovs_tot

renovs_tot
estructura_arborea_tots %>%
  mutate(PARCELA = as_factor(PARCELA)) %>% 
  left_join(renovs_tot, by = "PARCELA") %>% 
  replace(is.na(where(is.numeric)), 0) -> estructura_arborea_tots
## Warning in is.na(where(is.numeric)): is.na() aplicado a un objeto que no es
## (lista o vector) de tipo 'closure
head(estructura_arborea_tots)
ggplot(estructura_arborea_tots, aes(as_factor(PARCELA), TOTAL))+
  geom_col() + xlab("Parcela") +
  theme_bw()

Podemos separar los árboles por clase diamétrica y combinarla con los registros de renovales para evaluar la distribución de tamaños.

dfarboles %>% 
  mutate(CLASE_D = as_factor(cut_interval(DAP, length = 10, right = F))) %>% 
  group_by(PARCELA, CLASE_D) %>% 
  summarise(Cant = n()) %>% 
  mutate(Densidad = case_when(CLASE_D == "[10,20)" ~ Cant*40,
                             T ~ Cant * 10),
         PARCELA = as_factor(PARCELA)) %>% 
  dplyr::select(PARCELA, CLASE_D, Densidad) -> dfclasesd
## `summarise()` has grouped output by 'PARCELA'. You can override using the
## `.groups` argument.
dfclasesd
renovs_tot %>% 
  dplyr::select(1:2)  %>% 
  rename(Densidad = TOTAL) %>% 
  mutate(CLASE_D = "[0,10)") %>% 
  dplyr::select(PARCELA, CLASE_D, Densidad) %>% 
  bind_rows(dfclasesd) %>% 
  arrange(PARCELA) %>% 
  mutate(CLASE_D = as_factor(CLASE_D))-> dfclasesd

dfclasesd
ggplot(dfclasesd, aes(CLASE_D, Densidad)) +
  geom_col() +
  facet_wrap(vars(PARCELA)) +
  theme_bw()

Material leñoso

El muestreo de material leñoso caído se utiliza para estimar la biomasa muerta en descomposición y como indicador de las emisiones de C que se producirán.

Primero calculamos la densidad de madera promedio para utilizar de referencia:

mean(dfarboles$meanWD) #0.605
## [1] 0.604998

Necesitamos estimar cómo decae el C almacenado en la madera a medida que se descompone. A campo, el material leñoso es clasificado en 4 categorías según el grado de descomposición. Vamos a suponer:

Categoria MO remanente (%)
1 100
2 90
3 75
4 50

Con esto estimamos la necromasa caída:

matlenioso %>% 
  mutate(PARCELA = as_factor(str_sub(Parcela, start = 5, end = -1L)),
         DM = case_when(Estado == 1 ~ 0.605,
                        Estado == 2 ~ 0.605*0.9,
                        Estado == 3 ~ 0.605*0.75,
                        Estado == 4 ~ 0.605*0.5),
         masa_i = pi*(Diametro..cm./2)^2*DM) %>% 
  group_by(PARCELA) %>% 
  summarise(Necromasa = (sum(masa_i)/35*100)*(10^-2)*0.1) -> necromasa_caida #en tn/ha

y esto lo usamos en combinación con los árboles muertos en pie para estimar la necromasa total

estructura_arborea_tots %>% 
  left_join(necromasa_caida, by="PARCELA") %>% 
  mutate(necromasat = replace_na(biomasat_m + Necromasa, 0)) -> estructura_arborea_tots

y podemos visualizar:

ggplot(estructura_arborea_tots, aes(as_factor(PARCELA), necromasat))+
  geom_col() + xlab("Parcela") +
  theme_bw()

Estimacion a escala de predio/lote/estrato

El último paso es hacer inferencias a partir de las estimaciones muestrales. Vamos a calcular la media (+/- ES) de biomasa, necromasa, área basal y densidad de árboles de la Reserva.

estructura_arborea_tots %>% 
  dplyr::select(AB, dens, biomasat, necromasat) %>%  
  pivot_longer(AB:necromasat, names_to = "variable", values_to = "valor") %>% 
  group_by(variable) %>% 
  summarise(media = mean(valor),
            es = sd(valor)/sqrt(19))
## Adding missing grouping variables: `PARQUE`, `PARCELA`

También podemos estudiar la regeneración a escala de la reserva. Primero, promediamos la densidad de árboles por clase diamétrica:

orden <- c("[0,10)", "[10,20)", "[20,30)","[30,40)", "[40,50)", "[50,60)", "[60,70)","[70,80)", "[80,90)", "[90,100)", "[100,110)", "[110,120)",">=120")
dfclasesd %>% 
  group_by(CLASE_D) %>% 
  summarise(Densm = mean(Densidad)) -> dfclasesm

dfclasesm$CLASE_D = factor(dfclasesm$CLASE_D, levels = orden)

Vemos cómo quedó:

ggplot(dfclasesm, aes(CLASE_D, Densm)) +
  scale_x_discrete(limits = orden)+
  geom_col() +
  theme_bw()
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_col()`).

Podemos estimar los parámetros característicos de la distribución de Weibull usando los paquetes ‘MASS’ y ‘fitdistrplus’

fitdist(dfclasesm$Densm, distr = "weibull")
## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters:
##         estimate  Std. Error
## shape  0.4945837  0.09242295
## scale 73.5073431 44.00866443