AGB forest sampling calculations
This post shows an example of how to calculate some common plot-level variables from a forest sampling.
library(BIOMASS)
library(stringr)
library(tidyverse)Read data
# Field data with individual tree measures
df <- read.csv("D:/Drive/Jonathan_trabaggio/Doctorado/R/Ayuquila_Degradation/CleanData/df_all.csv",
na.strings = "NA")
# Coordinates of each site.
coords <- read.csv("D:/Drive/Jonathan_trabaggio/Doctorado/R/Ayuquila_Degradation/Data/gpscoords.csv")How the headers of the data look like the following for the df object.
Nombre Especie Observaciones DAP1 DAP2 DAP3 DAP4 DAP5 DAP6 DAP7 DAP8 DAP9 DAP10 DAP11
1 <NA> <NA> <NA> 3.44 NA NA NA NA NA NA NA NA NA NA
2 <NA> <NA> <NA> 6.81 NA NA NA NA NA NA NA NA NA NA
3 lysiloma Lysiloma divaricatum <NA> 3.12 NA NA NA NA NA NA NA NA NA NA
4 leocarpus Heliocarpus sp. <NA> 8.72 NA NA NA NA NA NA NA NA NA NA
5 muerto <NA> <NA> 4.04 NA NA NA NA NA NA NA NA NA NA
6 leocarpus Heliocarpus sp. <NA> 3.34 NA NA NA NA NA NA NA NA NA NA
DAP12 DAP13 DAP14 DAP15 DAP16 DAP17 DAP18 DAP19 DAP20 DAP21 DAP22 DAP23 DAP24 Altura parcela x
1 NA NA NA NA NA NA NA NA NA NA NA NA NA 1.66 Amacuau1 19°53.904
2 NA NA NA NA NA NA NA NA NA NA NA NA NA 4.56 Amacuau1 19°53.904
3 NA NA NA NA NA NA NA NA NA NA NA NA NA 3.65 Amacuau1 19°53.904
4 NA NA NA NA NA NA NA NA NA NA NA NA NA 5.29 Amacuau1 19°53.904
5 NA NA NA NA NA NA NA NA NA NA NA NA NA 2.69 Amacuau1 19°53.904
6 NA NA NA NA NA NA NA NA NA NA NA NA NA 3.02 Amacuau1 19°53.904
y cobertura register Observaciones.sitio id DAP1.BA DAP2.BA DAP3.BA DAP4.BA DAP5.BA DAP6.BA
1 104°06.428 70-80 Yan <NA> 1 9.294088 NA NA NA NA NA
2 104°06.428 70-80 Yan <NA> 2 36.423704 NA NA NA NA NA
3 104°06.428 70-80 Yan <NA> 3 7.645380 NA NA NA NA NA
4 104°06.428 70-80 Yan <NA> 4 59.720420 NA NA NA NA NA
5 104°06.428 70-80 Yan <NA> 5 12.818955 NA NA NA NA NA
6 104°06.428 70-80 Yan <NA> 6 8.761588 NA NA NA NA NA
DAP7.BA DAP8.BA DAP9.BA DAP10.BA DAP11.BA DAP12.BA DAP13.BA DAP14.BA DAP15.BA DAP16.BA DAP17.BA DAP18.BA
1 NA NA NA NA NA NA NA NA NA NA NA NA
2 NA NA NA NA NA NA NA NA NA NA NA NA
3 NA NA NA NA NA NA NA NA NA NA NA NA
4 NA NA NA NA NA NA NA NA NA NA NA NA
5 NA NA NA NA NA NA NA NA NA NA NA NA
6 NA NA NA NA NA NA NA NA NA NA NA NA
DAP19.BA DAP20.BA DAP21.BA DAP22.BA DAP23.BA DAP24.BA
1 NA NA NA NA NA NA
2 NA NA NA NA NA NA
3 NA NA NA NA NA NA
4 NA NA NA NA NA NA
5 NA NA NA NA NA NA
6 NA NA NA NA NA NAand like the following for the coords file.
name elevation date x y xutm ytum
1 LIMON1 1118.1133 2022/05/12 -104.1652 19.83815 587412.5 2193788
2 AMACUAU1 998.5761 2022/05/13 -104.1071 19.89840 593460.8 2200486Fix taxonomy and get wood density for AGB calculation
The first step is to get the corrected names of the species registered in the field. This information will be used to get the wood density by species, genus or family, depending on what info is available in the global woodensity database. Then this wood density is going to be used to calculate AGB as
# Separate genus and species
df <- df |>
mutate(genus = str_extract(Especie, "[A-z]+(?= )"),
species = str_extract(Especie, "(?<= )[A-z]+(?=)")) |>
mutate(across(species, ~gsub("sp|sp.","NA",.x)))
# Correct taxo names
taxo <- correctTaxo(genus = df$genus,
species = df$species,
useCache = F,
verbose = FALSE)
# Add as new columns
df <- df |>
mutate(genuscorr = taxo$genusCorrected,
speciescorr = taxo$speciesCorrected)
# Get family according to APG 3
APG <- getTaxonomy(df$genuscorr, findOrder = TRUE)
# Add to original df
df <- df |>
mutate(family = APG$family)
# Get wood density
dataWD <- getWoodDensity(
genus = df$genuscorr,
species = df$speciescorr,
family = df$family,
region = "World",
stand = df$parcela
)
# add as columns to df
df <- df |>
mutate(meanwood = dataWD$meanWD,
sdwood = dataWD$sdWD,
levelwood = dataWD$levelWD,
nInd = dataWD$nInd)
# Compute AGB for all DAP
df_agb <- df |>
mutate(across(matches("DAP[0-9]+$"), ~ computeAGB(D = .x,
WD = meanwood,
H = Altura))) |>
select(id, matches("DAP[0-9]+$"))
colnames(df_agb) <- c("id", paste0(("AGB"), seq(1,24)))
# Join AGB by id
df <- df |>
left_join(df_agb, "id")After these steps we now got the individual AGB for each stem and individual, as new columns of the original df. So the next step is going to calculate these same metrics for each individual (e.g., individual AGB sum of all of its stems). Watch that you might need to change some parameters of the subplot_size definition according to your own sampling design.
# Max DAP, plot extent
subplot_size <- df|>
select(id, matches("DAP[0-9]+$")) |>
pivot_longer(cols = -id,
names_to = c("AGB")) |>
drop_na(value) |>
group_by(id) |>
summarise(DAP_max = max(value)) |>
mutate(subplot_size = case_when(DAP_max >= 5 ~ 500,
DAP_max >= 2.5 & DAP_max < 5 ~ 29,
DAP_max < 2.5 ~ 0))
# AGB sum
AGB <- df|>
select(id, starts_with("AGB")) |>
pivot_longer(cols = -id,
names_to = c("AGB")) |>
drop_na(value) |>
group_by(id) |>
summarise(AGB_sum = sum(value))
# number of stems
Stems <- df|>
select(id, starts_with("AGB")) |>
pivot_longer(cols = -id,
names_to = c("AGB")) |>
drop_na(value) |>
group_by(id) |>
summarise(Stem_sum = n())
# Basal area
BA <- df|>
select(id, ends_with("BA")) |>
pivot_longer(cols = -id,
names_to = c("BA")) |>
drop_na(value) |>
group_by(id) |>
summarise(BA_sum = sum(value))
# Join previous calculation to original df
df <- df |>
left_join(AGB, "id") |>
left_join(BA, "id") |>
left_join(Stems, "id") |>
left_join(subplot_size, "id")
# Select variables of interest
df <- df |>
select(parcela, x, y, cobertura, Observaciones.sitio,
id, subplot_size, AGB_sum, BA_sum, Stem_sum, Altura,
genuscorr, speciescorr, meanwood, levelwood,
register) Now that we have the individual measures we need to calculate the per plot variables: AGBplot, BAplot, Dplot, Stemplot, Hmplot, H10plot. Since some of these measures are sums, extrapolated to 1 ha sums, means, etc, each one is summarised using the most common function to calculate it (e.g., AGB sum, height mean). Finally, we assume that the best registered coordinates are located in the coords file; thus, these coordinates are pasted on to the final result.
# -------------------Per plot variables------------------------------
# Summarise variables that need to be extrapolated
vars1 <- df |>
group_by(parcela, subplot_size) |>
# Sums by subplot_size
summarise(AGBsubplot = sum(AGB_sum),
BAsubplot = sum(BA_sum) / 10000,
Dsubplot = n(),
Stemsubplot = sum(Stem_sum)) |>
# Extrapolate to 1 ha according to the subplot size
mutate(across(c(AGBsubplot, BAsubplot, Dsubplot, Stemsubplot), ~.x * 10000 / subplot_size)) |>
ungroup() |>
group_by(parcela) |>
# Sum both subplot estimates
summarise(AGBplot = sum(AGBsubplot),
BAplot = sum(BAsubplot),
Dplot = sum(Dsubplot),
Stemplot = sum(Stemsubplot))
# Calculate top 10 mean height
vars2 <- df |>
group_by(parcela) |>
slice_max(Altura, n = 10) |>
summarise(H10plot = mean(Altura))
# Cover
vars3 <- df |>
select(parcela, cobertura) |>
separate(col = cobertura,
sep = "-",
into = c("cob1", "cob2")) |>
mutate(across(starts_with("cob"), ~as.numeric(.x))) |>
group_by(parcela) |>
summarise(cob = mean(c(cob1, cob2), na.rm = T))
# Mean height
vars4 <- df |>
group_by(parcela) |>
summarise(Hmplot = mean(Altura))
# Join all calculations
resul <- vars1 |>
left_join(vars2, "parcela") |>
left_join(vars3, "parcela") |>
left_join(df |>
select(c(x, y, parcela)) |>
distinct(parcela, .keep_all = T),
"parcela") |>
left_join(vars4, "parcela") |>
select(parcela, x, y, cob, AGBplot, BAplot, Dplot, Stemplot, Hmplot, H10plot)
# Rename columns
colnames(resul) <- c("Plot", "Lat", "Long",
"Cob(%)", "AGB(Mgha-1)", "BA(m2ha-1)",
"Dplot(indha-1)", "Stemplot(stemha-1)",
"Hmean(m)", "H10mean(m)")
resul <- resul |>
select(-c(Lat, Long))
coords <- coords |>
rename("yutm" = "ytum") |>
select(name, elevation, date, xutm, yutm) |>
mutate(across(name, ~str_to_title(.x))) |>
rename("Plot" = "name")
resul <- resul |>
left_join(coords, "Plot") |>
select(Plot, xutm, yutm, everything()) |>
mutate(year = 2022)Results look like the following.
# A tibble: 2 × 13
Plot xutm yutm `Cob(%)` `AGB(Mgha-1)` `BA(m2ha-1)` `Dplot(indha-1)` `Stemplot(stemha-1)` `Hmean(m)`
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Amacu… 5.93e5 2.20e6 75 37.6 18.2 2919. 5159. 4.46
2 Limon1 5.87e5 2.19e6 75 52.8 29.4 2130. 4410. 3.77
# ℹ 4 more variables: `H10mean(m)` <dbl>, elevation <dbl>, date <chr>, year <dbl>