Manual estimación biomasa con drones

library(ForestTools)
library(BIOMASS)
library(ggplot2)
library(readxl)
library(scales)
library(raster)
library(dplyr)
library(tidyr)
library(terra)
library(dplyr)
library(tmap)
library(sf)

Manual para estimación de biomasa con imágenes de dron, utilizando software libre

Presentamos una metodología de 3 componentes: fotogrametría, dendrometría en campo y modelos alométricos, para estimar biomasa en parcelas agrícolas, específicamente del tipo de Milpa Intercalada con Árboles Frutales (MIAF). La metodología consiste en los siguientes pasos:

En campo:

1. Captura de imágenes haciendo un vuelo programado para fotogrametría. Para este caso usamos la aplicación Drone Harmony para la planeación de vuelo y drones Mavic Air 2S y Mavic Mini de la empresa DJI, estos pasos no serán detallados en este manual.
2. Toma de datos de dendrometría en campo: altura y circunferencia a la altura del pecho (para después convertirla en diámetro a la altura del pecho o DAP), tomados con cinta métrica. Este paso no necesita ser detallado en este manual.

En gabinete:

1. Generación de modelos digitales de elevación de terreno y superficie, a partir de las imágenes obtenidas con el dron, en el programa opendronemap OpenDroneMap (2024).
2. Derivación de un Modelo de Altura de Dosel (Canopy Height Model, en inglés) con la librería ForestTools (Plowright 2024) del lenguaje de programación estadística R. En este paso nos interesa el archivo de cúspides (tree tops) generado por el programa.
3. Ajuste de un modelo lineal siguiendo la fórmula DAP ~ altura a los datos tomados en campo.
4. Aplicación del modelo lineal del paso anterior a los datos tomados en dron, para estimar el DAP a partir de la altura de las cúspides; la altura lleva un ajuste previo, que puede ser manual o programático, por la diferencia entre datos medidos en campo y la altura que un modelo puede estimar.
5. Estimación de biomasa aérea (above ground biomass) utilizando el paquete BIOMASS (Rejou-Mechain et al. 2017), el cual usa modelos pantropicales de alometría. A la función BIOMASS::computeAGB le suministra la altura de la cúspide del Modelo de Altura de Dosel, el diámetro estimado con el modelo lineal y la especie de árbol: el paquete BIOMASS tiene a la fecha 16,467 valores de densidad de maderas para distintas especies, incluyendo las monitoreadas en este proyecto (carambolo, mango, chico zapote, durazno y manzano).

OpenDroneMap. 2024. OpenDroneMap Authors ODM – A command line toolkit to generate maps, point clouds, 3D models and DEMs from drone, balloon or kite images. OpenDroneMap. https://github.com/OpenDroneMap/ODM.

Plowright, Andrew. 2024. «ForestTools: Tools for Analyzing Remote Sensing Forest Data». https://CRAN.R-project.org/package=ForestTools.

Rejou-Mechain, Maxime, Ariane Tanguy, Camille Piponiot, Jerome Chave, y Bruno Herault. 2017. «BIOMASS : an R package for estimating above-ground biomass and its uncertainty in tropical forests». Editado por Sarah Goslee. Methods in Ecology and Evolution 8 (9). https://doi.org/10.1111/2041-210X.12753.

Generación de modelos digitales de terreno y de elevación

Para este proceso utilizamos la aplicación para línea de comandos OpenDroneMap, la cual corre en un contenedor Docker. En el siguiente bloque de código se muestra cómo instalar Docker para Ubuntu. (para derivativos de Ubuntu cambiar $VERSION_CODENAME por $UBUNTU_CODENAME)

# Para instalar docker

# Add Docker's official GPG key:
sudo apt-get update
sudo apt-get install ca-certificates curl
sudo install -m 0755 -d /etc/apt/keyrings
sudo curl -fsSL https://download.docker.com/linux/ubuntu/gpg -o /etc/apt/keyrings/docker.asc
sudo chmod a+r /etc/apt/keyrings/docker.asc

echo \
  "deb [arch=$(dpkg --print-architecture) signed-by=/etc/apt/keyrings/docker.asc] https://download.docker.com/linux/ubuntu \
  $(. /etc/os-release && echo "$VERSION_CODENAME") stable" | \
  sudo tee /etc/apt/sources.list.d/docker.list > /dev/null
sudo apt-get update

# instalar los paquetes de docker:
sudo apt-get install docker-ce docker-ce-cli containerd.io docker-buildx-plugin docker-compose-plugin

# probamos la instalación
sudo docker run hello-world

Teniendo la instalación de Docker ahora podemos correr OpenDroneMap. Lo más sencillo para usarlo es crear una carpeta con el nombre de nuestro sitio y adentro de este una carpeta llamada images- que tenga las imágenes del vuelo. Una vez lista la carpeta con las imágenes nos posicionamos en la carpeta con el nombre de nuestro sitio podemos ejecutae el comando odm en su forma más elemental:

usuario@computadora:~/sitio/images$ docker run -ti --rm -v /my/project:$(pwd) opendronemap/odm --project-path /datasets

Para el caso de modelos de elevación, frecuentemente debemos de ajustar el comando con los siguientes parámetros:

1. --use-fixed-camera-params Para evitar efectos de “bowling” en el terreno. Este parámetro no lleva más argumentos.
2. --pc-rectify Hacer rectificación del suelo en la nube de puntos. Esto busca que los puntos de suelo clasificados erróneamente serán re-clasificados, rellenándose así los vacíos. Este parámetro no lleva más argumentos.
3. --smrf-scalar Valor de escalamiento: incrementar este parámetro para sitios con mucha variación de alturas. Default: 1.25
4. --smrf-slope Parámetro de tolerancia: una medida de la “tolerancia a pendientes”. Incrementar este parámetro para terrenos con mucha variación de alturas. Entre 0.1 y 1.2. Default: 0.15
5. --smrf-threshold Umbral de elevación. Ajustar este parámetro a la altura mínima, en metros, que esperamos para los objetos para diferenciarlos del terreno. Default: 0.5
6. --smrf-window Parámetro de radio de ventana, en metros, que corresponde al tamaño del objeto más grange (edificio, árboles, etc.) que removeremos de ser considerado como suelo. Debe ser un valor mayor a 10. Default: 18.

Generación del modelo de altura de dosel

Para generar el MAD, el primer paso es restar los modelos de superficie y de terreno (Figura 1). En el código siguiente además lo recortamos a un área específica, trazada sobre el MIAF, eliminando el resto de la imagen que no nos interesa, y usamos la función terra::resample para reducir la resolución y hacer más rápidos los cálculos subsecuentes.

fincas = read_sf("ejemplo/polys_miaf.gpkg") |> 
  st_zm() |> vect() |> project("epsg:32614")

i = 1
orto  = rast("ejemplo/odm_orthophoto.tif")
dtm   = rast("ejemplo/dtm.tif")
dsm   = rast("ejemplo/dsm.tif")

dsm = terra::crop(dsm, fincas[i,]) |> terra::mask(fincas[i,])
dtm = terra::crop(dtm, fincas[i,]) |> terra::mask(fincas[i,])
c_height = dsm - dtm 

ea_5 = rast(ext(c_height), resolution = 0.14, crs = crs(c_height))
ea_5 = resample(c_height, ea_5)
orto = resample(orto, ea_5)
ea_5 = raster(ea_5) # foresttools no usa terra, sino raster


plot(ea_5)

Figura 1: Diferencia DSM-DTM, metros, resampleada a 0.14 cm de resolución

Para hacer un modelo de altura de dosel efectivo tenemos que ajustar los parámetros de la función de la ventana móvil (winFun) y la altura mínima (minHeight) de la función filtro de ventana variable vwf. De acuerdo a la altura y tamaños de copa de nuestros frutales. Por lo general, a mayor altura usaremos ventanas más amplias y alturas mínimas (minHeight) más grandes. Este es un ejemplo Tabla 1 de distintas parcelas.

Tabla 1: Parámetros utilizados para distintas parcelas

Zona	Frutal	Función
Huejotzingo	Manzana y durazno	lin <- function(x){x * 0.05 + .4}
		ttops <- vwf(CHM = ea_5, winFun = lin, minHeight = 0.35)
Xigüipilincan	Mango	lin <- function(x){x * 0.05 + 1}
		ttops <- vwf(CHM = ea_5, winFun = lin, minHeight = 0.5)
Axochío	Carambolo	lin <- function(x){x * 0.05 + .8}
		ttops <- vwf(CHM = ea_5, winFun = lin, minHeight = 0.5)
Axochío	Chico zapote	lin <- function(x){x * 0.05 + .9}
		ttops <- vwf(CHM = ea_5, winFun = lin, minHeight = 0.5)

Ahora correremos las funciones de ForestTools::vwf (filtro de ancho de ventana variable) para detectar las cúspides de los árboles y el ForestTools:mcws para detectar el área de copas el cual es un algoritmo parecido a los de delimitación de cuencas, pero invertido. El resultado de las cúspides y demarcaciónes de copas está en Figura 2. Podemos ver que el modelo funcionó correctamente porque están en la división de las copas se puede ver el corte en tatura de la parcela MIAF.

lin <- function(x){x * 0.05 + .4}
ttops <- vwf(CHM = ea_5, winFun = lin, minHeight = 0.35)
crownsPoly <- mcws(treetops = ttops, CHM = ea_5, format = "polygons",
                   minHeight = 0.35) |> vect() #, verbose = FALSE)

ttops$layer = paste("Cúspides, total:", nrow(ttops))
crownsPoly$layer = paste("Copas, área total:", comma( sum( round( expanse( crownsPoly )  ,2))), " m²")

mapa_rgb =  
  tm_layout(main.title = sprintf("Modelo de altura de dosel, Finca %s", fincas$sitio[i]),
            main.title.size = 1.2, 
            title = paste0("No. cúspides: ", nrow(crownsPoly), "\n",
                           "Cobertura arbórea: ", round( sum(expanse(crownsPoly) / expanse(fincas[i,]), 2 )),
                            "%\n",
                            "Área predio: ", round(expanse(fincas[i,])/10000,1), " ha" ) 
            ) +
  tm_shape(orto, unit = "m") +
  tm_rgb() +
  tm_shape(ttops) + 
  tm_symbols(col = "height", size = .1, 
              palette = viridis::viridis(8), title.col = "Altura, m", 
              border.col = "black", border.lwd = .8, legend.hist = F,
              legend.col.reverse = T, legend.format = list(text.separator = "a")) +
   tm_shape(ttops) + 
   tm_symbols(col = "height", size = .1, 
              palette = viridis::viridis(8), title.col = "Altura, m", legend.col.show =  F,
              border.col = "black", border.lwd = .8, legend.hist = T, 
              legend.hist.title = "Frec. alturas, m", 
              legend.col.reverse = F, legend.format = list(text.separator = "a")) +  
  tm_shape(st_as_sf(crownsPoly)) + 
  tm_polygons(col = "layer", alpha = 0, border.col = "black", legend.show = T, title = "") + 
  tm_layout(legend.outside = T, legend.hist.width = 4 ) +
  tm_scale_bar(text.color = "black", text.size = .8, breaks = c(0,10,50), position = c(0.04, 0.01),
               bg.color = "white") +
  tm_graticules(projection = 32614, labels.cardinal = F, n.x = 4, n.y = 4, 
                labels.format = list(fun = function(x) x), labels.rot = c(0,90), col = "gray70",
                lwd = .7) 

mapa_rgb

Figura 2: Modelo de altura de dosel

Estimación de biomasa

Ahora cargamos un archivo con las lecturas de circunferencia (la cual transformaremos a diámetro) y altura, con esto obtenemos un modelo lineal, ya sea generalizado o estándar (funciones glmy lm de R) para cada frutal (Figura 3), el cual ajustaremos al objeto de cúspides (ttops) para estimar un diámetro para cada cúspide.

genero = c("Averrhoa", "Manilkara", "Prunus", "Mangifera", "Malus")
especie = c("carambola", "zapota", "persica", "indica", "domestica")
names(genero) = c("Carambola", "Chico zapote", "Durazno", "Mango", "Manzana")
names(especie) = names(genero)

df = read_excel("ejemplo/mediciones_arbol_fruto.xlsx", sheet = "Mediciones") |>
  dplyr::select(1:6) |> dplyr::filter(!if_all(.cols = everything(), .fns = is.na)) |>
  fill(c(parcela, fruta, sitio), .direction = "down") 

df$diametro_cm = df$diametro_cm/pi
  
df |> ggplot(aes(diametro_cm, altura_m)) +
  geom_point(pch = 17) + 
  geom_smooth(method = "lm", se = FALSE, color = 1, lty = "12", lwd = 0.3) +
  facet_wrap(~fruta, scales = "free") + theme_minimal() + theme(panel.background = element_rect(fill = "cornsilk"))

Figura 3: Dendrometría y modelos lineales por frutal analizado

cuspides = ttops |> st_transform(4326)

temp = df[df$parcela == sprintf("Parcela %s", fincas$unidad[i]) & df$fruta == fincas$frutal[i],]
modelo = glm(diametro_cm ~ altura_m, data = temp)

AIC(modelo) # the smaller the AIC or BIC, the better the fit {help(AIC)}

[1] 48.42084

# para comparar: 
# lm1 <- lm(Fertility ~ . , data = swiss); AIC(lm1)
# nos da AIC: 326

# necesitamos un ajuste por la variación de lo que mide el personal en campo y lo que registra el dron

cuspides_parcela_fruta = ttops #st_filter(cuspides, parcelas[i,])
ajuste = if( max(temp$altura_m) - max(cuspides_parcela_fruta$height)< 0){ 
  0} else { max(temp$altura_m) - max(cuspides_parcela_fruta$height)} 
diametros = predict.glm(modelo, data.frame(altura_m = cuspides_parcela_fruta$height + ajuste), type = "response")

biomasa = BIOMASS::computeAGB(D = diametros, H = cuspides_parcela_fruta$height + ajuste, 
                    WD = rep(getWoodDensity(genero[[fincas$frutal[i]]], especie[[fincas$frutal[i]]],
                                            region = "NorthAmerica")$meanWD, 
                    length(diametros) ))
# Toneladas

par(mar=c(5, 4, 4, 8), bg = "cornsilk", xpd = FALSE) # xpd es para que imprima fuera del margen, se activa en legend

with(temp,
     {plot(diametro_cm,  altura_m, 
           main = sprintf("Parcela %s, %s,  %s", fincas$unidad[i], fincas$sitio[i], fincas$frutal[i]),
           sub = paste("Ajuste del modelo lineal; AIC:", 
                       round( AIC(modelo),1) )) }
  )

lines(diametros, cuspides_parcela_fruta$height + ajuste)
legend( "bottomleft",  cex = 0.8,y.intersp = 1.3, xpd = NA, 
        c("Mediciones\nen campo.", "Modelo lineal\nde dron."), pch = c(1, NA), lty = c(NA, 1))
mtext(paste("Biomasa:\n", round(sum(biomasa),1), "tons\n", 
            "Superficie:\n", round(expanse(fincas[i,]),1), "m²"), cex = 0.8, side = 4, line = 2, las = 1)

Figura 4

Nuesto modelo lineal tiene un AIC de 48.4, lo cual es un buen ajust. La estimación para esta parcela nos da 4.9 toneladas de biomasa de los árboles de manzana, la cual tiene una densidad de madera promedio de 0.6017 g/cm³, según el paquete BIOMASS para la región de Norteamérica.

Conclusiones

El método mixto de dendrometría manual, vuelo de dron y modelos alométricos resultó ser útil para la estimación de biomasa en 5 cultivos distintos en dos regiones de México (en este ejercicio sólo presentamos un frutal de una región). Hemos presentado un método plenamente replicable que hace uso exclusivamente de software libre, que permite la estimación para extensiones de varias hectáreas con una sola visita al sitio. Por otra parte, la metodología requiere un mínimo capacidades técnicas del personal: realizar vuelos programados, generar ortomosaicos, ajustar parámetros para el modelo de altura de dosel y leer y reemplazar valores en un archivo de R; estas capacidades son alcanzables y este manual ayuda en cada parte del proceso, a excepción de hacer vuelos programados.