ESTIMACIÓN DE BIOMASA AÉREA FORESTAL
Adriana Tovar
27 de noviembre de 2017
El objetivo de este estudio fue estimar la biomasa aérea de una plantación forestal comercial de Eucalyptus grandis localizada en el departamento del Cauca, Colombia.
Las fuentes utilizadas fueron datos dasométricos de parcelas circulares tomados en campo y datos satelitales multiespectrales Sentinel-2 con reflectancia al tope de la atmósfera (TOA); con estos últimos se generaron los índices espectrales NDVI, EVI, SAVI, OSAVI Y RERNDVI.
La biomasa aérea es uno de los parámetros cruciales para evaluar el almacenamiento y las emisiones de C causadas por la deforestación y la degradación de los bosques. Dado que la biomasa de la vegetación afecta a una gama de procesos ecosistémicos, como el ciclo del carbono y el agua, así como los flujos de energía, se requiere información exacta de BA para el desarrollo de estrategias de manejo forestal sostenible.
Sentinel-2 pueden proporcionar información relevante debido a su resolución espectral, espacial y temporal. El algoritmo Random Forest ha recibido una atención creciente en la última década debido a que esta técnica de aprendizaje automático agrega árboles, aprende la relación entre el conjunto de datos de entrenamiento y respuesta, y calcula la variable de respuesta. Por ende se aplicó esta técnica con el fin de seleccionar el mejor modelo explicativo.
Área de estudio
El estudio se realizó en diversas regiones del departamento del Cauca, en zonas de bosque húmedo tropical (bh-T) y bosque seco tropical (bs-T), en un rango altitudinal de 1000 a 2200 m.s.n.m.
library(knitr)
library(sp)
library(raster)
library(rasterVis)## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewerlibrary(tmap)
setwd("D:/Proyecto_final")
#list.files(".")
bio <- shapefile("bio_cauca_wgs.shp")
pol <- shapefile("pol_cauca_wgs.shp")
# Visualicemos el área de estudio
library(leaflet)
m <- leaflet(pol) %>%
addPolygons() %>%
addTiles()
addMiniMap(m, tiles = providers$Esri.WorldStreetMap,
toggleDisplay = T, width = 120, height=80)De acuerdo con el siguiente histograma, los datos tienden a ser normales.
names(bio)## [1] "OBJECTID_1" "OBJECTID" "ID" "id_invent" "finca"
## [6] "lote" "especie" "dominantes" "volumTCC" "volumTSC"
## [11] "area_basal" "altura_med" "diametro_m" "TonTSC" "ton_Arb_DQ"
## [16] "area_parce" "edad" "Coo_X" "Coo_Y"# Histograma de la biomasa aérea forestal
hist(bio$TonTSC)
# Gráfico de la distribución de los datos
tm_shape(pol) + tm_fill(col="grey") +
tm_shape(bio) +
tm_dots(col="TonTSC", palette="RdYlGn",
auto.palette.mapping =FALSE,
title="Biomasa aérea de Eucalyptus grandis",
size=0.1) +
#
tm_legend(legend.outside=TRUE)
Variables basadas en Sentinel-2
La misión Sentinel-2, de la Agencia Espacial Europea (ESA), cuenta con 13 bandas espectrales de 10, 20 y 60 m de resolución espacial, de las cuales las bandas 5, 6 y 7 corresponden a los rangos espectrales del borde rojo, las cuales pueden proporcionar información clave sobre el estado de la vegetación (Drusch et al., 2012)
Se emplearon algunos índices de vegetación para la estimación de biomasa a partir de las imágenes Sentinel-2 teniendo en cuenta el poder predictivo de cada uno y la capacidad de comparación con otros estudios relacionados.
Índices espectrales. R = Reflectancia en la banda espectral
Descripción de los índices espectrales usados
ndvi = raster("ndvi.tif")
ndvi## class : RasterLayer
## dimensions : 3994, 1775, 7089350 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 312120, 329870, 253320, 293260 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\Proyecto_final\ndvi.tif
## names : ndvirequire(rasterVis)
mapTheme <- rasterTheme(region=brewer.pal(10,"RdYlGn"))
plt <- levelplot(ndvi, margin=F, par.settings=mapTheme, at = do.breaks(c(-1,1),30),main = "NDVI")
plt + layer(sp.lines(pol, col="black", lwd=1))
evi = raster("evi.tif")
evi## class : RasterLayer
## dimensions : 3994, 1775, 7089350 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 312120, 329870, 253320, 293260 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\Proyecto_final\evi.tif
## names : evisummary(evi)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (1.41% of all cells)## evi
## Min. -12475.000000
## 1st Qu. -9.397116
## Median 3.322476
## 3rd Qu. 11.959790
## Max. 13355.000000
## NA's 0.000000plt2 <- levelplot(evi, margin=F, par.settings=mapTheme, at = do.breaks(c(-50,50),30),main = "EVI")
plt2 + layer(sp.lines(pol, col="black", lwd=1))
savi = raster("savi.tif")
savi## class : RasterLayer
## dimensions : 3994, 1775, 7089350 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 312120, 329870, 253320, 293260 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\Proyecto_final\savi.tif
## names : savisummary(savi)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (1.41% of all cells)## savi
## Min. 0.08248472
## 1st Qu. 1.51098725
## Median 1.66519988
## 3rd Qu. 1.76033044
## Max. 2.01862907
## NA's 0.00000000plt3 <- levelplot(savi, margin=F, par.settings=mapTheme, at = do.breaks(c(-0.5,2.5),30),main = "SAVI")
plt3 + layer(sp.lines(pol, col="black", lwd=1))
osavi = raster("osavi.tif")
osavi## class : RasterLayer
## dimensions : 3994, 1775, 7089350 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 312120, 329870, 253320, 293260 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\Proyecto_final\osavi.tif
## names : osavisummary(osavi)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (1.41% of all cells)## osavi
## Min. -0.2746349
## 1st Qu. 0.6059565
## Median 0.6881583
## 3rd Qu. 0.7385939
## Max. 0.8742043
## NA's 0.0000000plt4 <- levelplot(osavi, margin=F, par.settings=mapTheme, at = do.breaks(c(-0.01,1),30),main = "OSAVI")
plt4 + layer(sp.lines(pol, col="black", lwd=1))
rerndvi = raster("rerndvi20.tif")
rerndvi## class : RasterLayer
## dimensions : 1997, 888, 1773336 (nrow, ncol, ncell)
## resolution : 20, 20 (x, y)
## extent : 312120, 329880, 253320, 293260 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : D:\Proyecto_final\rerndvi20.tif
## names : rerndvi20rerndvi<-resample(rerndvi, ndvi, method='bilinear')
summary(rerndvi)## rerndvi20
## Min. -3.543298e-01
## 1st Qu. 0.000000e+00
## Median 0.000000e+00
## 3rd Qu. 0.000000e+00
## Max. 2.946682e-01
## NA's 7.449860e+05plt5 <- levelplot(rerndvi, margin=F, par.settings=mapTheme,
at = do.breaks(c(-0.001,0.001),50),main = "RERNDVI")
plt5 + layer(sp.lines(pol, col="black", lwd=1))
##KRIGGING
Generación de puntos para la predicción con Krigging
setwd("D:/Proyecto_final")
malla<-read.csv2("malla_bio.csv", sep=";",header=TRUE,enc="latin1")
str(malla)## 'data.frame': 254528 obs. of 2 variables:
## $ x: int 1043560 1043610 1043660 1043710 1043760 1043810 1043860 1043910 1043960 1044010 ...
## $ y: int 745526 745526 745526 745526 745526 745526 745526 745526 745526 745526 ...coordinates(malla) <- ~x+y
class(malla)## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"proj4string(bio)=CRS("+init=epsg:3115")## Warning in `proj4string<-`(`*tmp*`, value = <S4 object of class structure("CRS", package = "sp")>): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransformproj4string(malla)=CRS("+init=epsg:3115")
gridded(malla) <- TRUE
library(gstat)
ba = sample(1:nrow(bio), size=0.4*nrow(bio))
ba## [1] 59 55 23 41 47 49 21 69 5 26 38 68 7 72 66 57 22 39 45 43 6 40 13
## [24] 29 67 9 27 11 34test = bio[ba,]
dim(test)## [1] 29 19train = bio[-ba,]
dim(train)#45 29## [1] 45 19class(train)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"proj4string(train)=CRS("+init=epsg:3115")
names(train)## [1] "OBJECTID_1" "OBJECTID" "ID" "id_invent" "finca"
## [6] "lote" "especie" "dominantes" "volumTCC" "volumTSC"
## [11] "area_basal" "altura_med" "diametro_m" "TonTSC" "ton_Arb_DQ"
## [16] "area_parce" "edad" "Coo_X" "Coo_Y"proj4string(test)=CRS("+init=epsg:3115")
names(test)## [1] "OBJECTID_1" "OBJECTID" "ID" "id_invent" "finca"
## [6] "lote" "especie" "dominantes" "volumTCC" "volumTSC"
## [11] "area_basal" "altura_med" "diametro_m" "TonTSC" "ton_Arb_DQ"
## [16] "area_parce" "edad" "Coo_X" "Coo_Y"## Variograma omnidireccional de BA ##
bio.vg1 <- variogram(TonTSC~1, test)
plot(bio.vg1, col='black', main="Variograma omnidireccional para BA")
## Variograma direccional de BA
bio.dvg <- variogram(TonTSC~1, test,alpha=c(0,45,90,135))
plot(bio.dvg, col='black', main="Variograma direccional de BA")
##Variograma usando gstat
bio.vgm <- vgm(psill=900, model="Wav", range=1200, nugget=2300)
plot(bio.vg1, model=bio.vgm, col='black', main="BA: Wave, nugget=2300, sill=900, range=1200")
ba.krige <- krige(TonTSC~1, test, malla, model=bio.vgm, beta=mean(test$TonTSC))## [using simple kriging]summary(ba.krige)## Object of class SpatialPixelsDataFrame
## Coordinates:
## min max
## x 1043560 1059910
## y 745526 784276
## Is projected: TRUE
## proj4string :
## [+init=epsg:3115 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Number of points: 254528
## Grid attributes:
## cellcentre.offset cellsize cells.dim
## x 1043560 50 328
## y 745526 50 776
## Data attributes:
## var1.pred var1.var
## Min. :205.3 Min. :3200
## 1st Qu.:205.3 1st Qu.:3200
## Median :205.3 Median :3200
## Mean :205.3 Mean :3200
## 3rd Qu.:205.3 3rd Qu.:3200
## Max. :205.3 Max. :3200cols <- colorRampPalette(c('red', 'yellow', 'green'))
pts <- list('sp.points', bio, pch=1, cex=0.7, col='black')
spplot(ba.krige, zcol='var1.pred', col.regions=cols(20), main="BA estimada", sp.layout=pts, scales=list(draw = TRUE))
Predicción de biomasa aérea forestal con Random Forest
# Extracción de información de los índices espectrales
library(raster)
drops <- c("OBJECTID_1","OBJECTID","ID","finca","lote","especie","dominantes","volumTCC","area_basal","altura_med" ,"diametro_m","ton_Arb_DQ","area_parce","edad","coords.x1.1","coords.x2.1", "optional" )
#Para eliminar los campos "AREA" y "PERIMETER"
ntrain <- train[,!(names(train) %in% drops)]
names(train)## [1] "OBJECTID_1" "OBJECTID" "ID" "id_invent" "finca"
## [6] "lote" "especie" "dominantes" "volumTCC" "volumTSC"
## [11] "area_basal" "altura_med" "diametro_m" "TonTSC" "ton_Arb_DQ"
## [16] "area_parce" "edad" "Coo_X" "Coo_Y"names(ntrain)## [1] "id_invent" "volumTSC" "TonTSC" "Coo_X" "Coo_Y"summary(ntrain)## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 -76.673174 -76.545347
## coords.x2 2.302318 2.645756
## Is projected: TRUE
## proj4string :
## [+init=epsg:3115 +proj=tmerc +lat_0=4.596200416666666
## +lon_0=-77.07750791666666 +k=1 +x_0=1000000 +y_0=1000000
## +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs]
## Number of points: 45
## Data attributes:
## id_invent volumTSC TonTSC Coo_X
## Length:45 Min. : 82.98 Min. : 70.49 Min. :1044975
## Class :character 1st Qu.:194.48 1st Qu.:163.14 1st Qu.:1047509
## Mode :character Median :224.76 Median :184.74 Median :1055104
## Mean :229.17 Mean :189.98 Mean :1053709
## 3rd Qu.:252.83 3rd Qu.:210.91 3rd Qu.:1058547
## Max. :359.16 Max. :294.31 Max. :1059200
## Coo_Y
## Min. :746351
## 1st Qu.:754919
## Median :760575
## Mean :761371
## 3rd Qu.:763475
## Max. :784328class(ntrain)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"####################library(randomForest)## randomForest 4.6-12## Type rfNews() to see new features/changes/bug fixes.bio2 <- data.frame(ntrain$id_invent, ntrain$TonTSC)