Evaluación del estado de la vegetación a partir de índices espectrales en la zona de compensación de la Hidroeléctrica El Quimbo, Huila mediante imágenes Sentinel 2
1. Introducción
En este cuaderno expondrá el procedimiento y resultado de la evaluación del estado de la vegetación de la zona de compensación de la Central Hidroeléctrica El Quimbo, en el departamento del Huila (Colombia) mediante el uso de índices espectrales con imágenes del sensor remoto Sentinel 2.
El Quimbo es una hidroeléctrica establecida en el centro del departamento del Huila, alrededor de la cual establecieron un área de aproximádamente 11 mil ha de restauración, realizando un proyecto piloto que finalizó en el año 2020. A continuación se muestra el área de restauración que será el área de estudio.
library(leaflet)
library(raster)## Loading required package: spzonaest =shapefile("~/Maestria/Tesis/shp/zona_84.shp")## Warning in rgdal::readOGR(dirname(x), fn, stringsAsFactors = stringsAsFactors, :
## Z-dimension discardedm <- leaflet(data=zonaest) %>% addTiles() %>%
addPolygons(fill = T, stroke = T, color = "#ecc617 ") %>% addLegend("bottomright", colors = "#ecc617 ", labels = "Zona Estudio - El Quimbo")
m2. Datos y métodos
2.1. Descarga y preprocesamiento de imágenes
Las imágenes usadas fueron obtenidas del sitio web de libre acceso earthexplorer.usgs.gov del día 11 de noviembre de 2020, correspondiente al cuadrante 18NVH. El procedimiento de corrección atmosférica se realizó en el software QGIS mediante el Plugin SCP con el método disponible DOS1 para imágenes Sentinel 2, para todas las bandas con excepción de la 1, 9 y 10. Se trabajará principalmente con las librerías de análisis espacial terra y raster. Teniendo los archivos en el directorio de trabajo, se procede a crear los objetos SpatRaster de cada una de las bandas, y se verifica la correcta creación de objeto con la banda 8 de prueba.
library(terra)## terra version 1.0.10#list of bands
b2<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B02.tif")
b3<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B03.tif")
b4<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B04.tif")
b5<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B05.tif")
b6<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B06.tif")
b7<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B07.tif")
b8<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B08.tif")
b8A<-rast("~/Maestria/Tesis/img/Sent/RT_T18NVH_20201109T152639_B8A.tif")
#verifying
b8## class : SpatRaster
## dimensions : 10980, 10980, 1 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : RT_T18NVH_20201109T152639_B08.tif
## name : RT_T18NVH_20201109T152639_B08Luego se realiza una combinación entre bandas para la creación de un archivo raster multibanda, renombrando posteriormente cada una de las bandas a su longitud de onda correspondiente. Se guardó el archivo multibanda en el disco con el nombre “sent_qb.tif”. El archivo tiene las capas correspondientes a las bandas 2,3,4,5,6,7,8 y 8A. Posteriormente se obtuvo un resumen de los datos de cada banda.
#Multiband raster and rename
qb<- c(b2,b3,b4,b5,b6,b7,b8,b8A)
qb## class : SpatRaster
## dimensions : 10980, 10980, 8 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 399960, 509760, 190200, 3e+05 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## sources : RT_T18NVH_20201109T152639_B02.tif
## RT_T18NVH_20201109T152639_B03.tif
## RT_T18NVH_20201109T152639_B04.tif
## ... and 5 more source(s)
## names : RT_T1~9_B02, RT_T1~9_B03, RT_T1~9_B04, RT_T1~9_B05, RT_T1~9_B06, RT_T1~9_B07, ...names(qb)## [1] "RT_T18NVH_20201109T152639_B02" "RT_T18NVH_20201109T152639_B03"
## [3] "RT_T18NVH_20201109T152639_B04" "RT_T18NVH_20201109T152639_B05"
## [5] "RT_T18NVH_20201109T152639_B06" "RT_T18NVH_20201109T152639_B07"
## [7] "RT_T18NVH_20201109T152639_B08" "RT_T18NVH_20201109T152639_B8A"names(qb)<-c("blue","green","red","re1","re2","re3","NIR","NIRn")
#new names
names(qb)## [1] "blue" "green" "red" "re1" "re2" "re3" "NIR" "NIRn"saving results: writeRaster(qb,“~/Maestria/Tesis/img/sent_qb.tif”, overwrite = T)
2.2. Corte de la imagen Sentinel 2 a zona de estudio
Para realizar esta delimitación se realizó el corte de la imagen en el software ArcGIS a la extensión de la zona de estudio para tener un panorama general de la región. De esta manera se obtuvo un archivo raster con las 8 capas de interes.
#creating the spatvector of study zone in corresponding CRS
zona_est<- vect("~/Maestria/Tesis/shp/zona_final.shp")
newqb<-rast("~/Maestria/Tesis/img/new_qb.tif")
#checking CRS and extent
newqb## class : SpatRaster
## dimensions : 3291, 1916, 8 (nrow, ncol, nlyr)
## resolution : 10, 10 (x, y)
## extent : 421310, 440470, 241750, 274660 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : new_qb.tif
## names : blue, green, red, re1, re2, re3, ...
## min values : 0.000000000, 0.010299999, 0.000000000, 0.000000000, 0.004799999, 0.003200000, ...
## max values : 1, 1, 1, 1, 1, 1, ...zona_est## class : SpatVector
## geometry : polygons
## dimensions : 1, 5 (geometries, attributes)
## extent : 421312, 440467.8, 241748.2, 274662.6 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : OBJECTID_1 OBJECTID Shape_Leng Shape_Le_1 Shape_Area
## type : <num> <num> <num> <num> <num>
## values : 1 1 3.111e+05 3.111e+05 1.108e+08A partir de este archivo se procedió a realizar una composición True color RGB432 para la zona .
#Composition True color RGB234 of study area
plotRGB(newqb,r=3,g=2,b=1, stretch="lin",axes=T, main="Zona de compensación El Quimbo - RGB432")
2.3. Estadísticas de la imagen
Para tener un panorama del área de estudio se obtuvieron estadísticas de la zona con la librería raster.
qb_raster<- brick("~/Maestria/Tesis/img/new_qb.tif")
newqb_mean <- cellStats(qb_raster, mean)
newqb_mean## blue green red re1 re2 re3 NIR
## 0.04722916 0.06059206 0.05961452 0.08592595 0.16275797 0.20075244 0.19784697
## NIRn
## 0.23453534Desviación estándar para las bandas
newqb_sd<-cellStats(qb_raster,sd)
newqb_sd## blue green red re1 re2 re3 NIR
## 0.01109552 0.01400979 0.02189057 0.02466336 0.05635303 0.07338374 0.07340452
## NIRn
## 0.08565559Visualización (linear strecht) de todas las bandas y sus correspondientes histogramas
par(mfrow=c(2,4))
tmp<-stretch(qb_raster[[1]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B2 ", axes=F, legend=F )
tmp<-stretch(qb_raster[[2]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B3", axes=F, legend=F )
tmp<-stretch(qb_raster[[3]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B4",axes=F, legend=F )
tmp<-stretch(qb_raster[[4]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B5",axes=F, legend=F )
tmp<-stretch(qb_raster[[5]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B6", axes=F, legend=F )
tmp<-stretch(qb_raster[[6]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B7", axes=F, legend=F )
tmp<-stretch(qb_raster[[7]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B8",axes=F, legend=F )
tmp<-stretch(qb_raster[[8]],minv=0,maxv=255,minq=0,maxq=1)
plot(tmp,col=gray(1:100/100),main="Sent 2 2020 B8A",axes=F, legend=F )
hist(qb_raster, c(1:8), maxpixels=5000000, plot=TRUE, main="Histogramas B2-B8A", xlab= "Surface reflectance", ylab="Frequency", col="#66CC99")
2.4. Obtención índices de vegetación
Este análisis se realizó solo para la zona de restauración. En primer lugar, se obtuvieron los índices correspondientes a NDVI, EVI, GNDVI y SAVI. Los cuales se añadieron como capas al archivo raster de Sentinel-2. A partir de este archivo, se procedió a realizar el corte por máscara en el software QGIS. El resultado final se visualiza posterior a los procedimientos.
#funtion for NDVI
vi <- function(img, nir, red) {
bn <- img[[nir]]
br <- img[[red]]
vi <- (bn - br) / (bn + br)
return(vi)
}
ndvi<-vi(qb_raster, 7,3)
#EVI function
vi2<- function(img,near, red, blue){
bn<-img[[near]]
br<-img[[red]]
bb<--img[[blue]]
vi2 <-2.5*((bn-br)/(bn+6*br-7.5*bb+1))
return(vi2)
}
evi<-vi2(qb_raster,7,3,1)
#GNDVI function
vi3 <- function(img, nir, green) {
bn <- img[[nir]]
bg <- img[[green]]
vi3 <- (bn - bg) / (bn + bg)
return(vi3)
}
gndvi<-vi3(qb_raster,7,2)
#SAVI
vi4<- function(img, nir, red) {
bn <- img[[nir]]
br <- img[[red]]
vi4 <- (bn - br) / (bn + br + 0.428) * (1.428)
return(vi4)
}
savi<-vi4(qb_raster,7,3)Combinando bandas de los índices obtenidos con las bandas Sentinel-2
#Creating a rasterbirck from vegetation index
vi_rast<-brick(ndvi,evi,gndvi,savi)
names(vi_rast) <- c("ndvi","evi","gndvi","savi")
#Adding vegetation index as layers to spatraster file
newqb_brick<-brick(newqb)
qb_vi<-brick(c(newqb_brick,vi_rast))
#checking new rasterbrick
qb_vi## class : RasterBrick
## dimensions : 3291, 1916, 6305556, 12 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 421310, 440470, 241750, 274660 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : C:/Users/juanc/AppData/Local/Temp/RtmpK8Fzic/raster/r_tmp_2021-04-14_212728_18400_56679.grd
## names : blue, green, red, re1, re2, re3, NIR, NIRn, ndvi, evi, gndvi, savi
## min values : 0.000000000, 0.010299999, 0.000000000, 0.000000000, 0.004799999, 0.003200000, 0.007900000, 0.004899999, -0.571045578, -0.076109268, -0.630841076, -0.224628091
## max values : 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 1.0000000, 0.6086919, 0.9079946, 0.7000458##Saving file:
##writeRaster(qb_vi,"~/Maestria/Tesis/img/qbvi.tif")#calling masked raster
masked<- brick("~/Maestria/Tesis/img/maskedqb.tif")
names(masked)<-names(qb_vi)
masked## class : RasterBrick
## dimensions : 3291, 1916, 6305556, 12 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 421311, 440471, 241750, 274660 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : C:/Users/juanc/Documents/Maestria/Tesis/img/maskedqb.tif
## names : blue, green, red, re1, re2, re3, NIR, NIRn, ndvi, evi, gndvi, savi
## min values : 0.02550000, 0.02360000, 0.02160000, 0.02600000, 0.03000000, 0.03210000, 0.02730000, 0.02830000, -0.40947071, -0.05928553, -0.36904758, -0.12195327
## max values : 0.1432000, 0.1760000, 0.2297000, 0.2204000, 0.3266000, 0.4323000, 0.4674000, 0.4903000, 0.8758939, 0.5523337, 0.8155838, 0.6636578#checking a layer
plot(masked[[9]],main="NDVI")
Para el análisis de las regiones con altos índices de vegetación se realizó la reclasificación de estos índices por valores por debajo y sobre la mediana, con valores de 0 y 1 respectivamente.
#Establishing median for each vegetation index
median9<-cellStats(masked[[9]],stat = median)
median10<-cellStats(masked[[10]],stat = median)
median11<-cellStats(masked[[11]],stat = median)
median12<-cellStats(masked[[12]],stat = median)
#subsetting high index pixels (above of median), equaling all these values to 1
hndvi<- reclassify(masked[[9]],c(-Inf,median9,0,median9,Inf,1))
hevi<- reclassify(masked[[10]],c(-Inf,median10,0,median10,Inf,1))
hgndvi<- reclassify(masked[[11]],c(-Inf,median11,0,median11,Inf,1))
hsavi<- reclassify(masked[[12]],c(-Inf,median12,0,median12,Inf,1))Visualización de los nuevos mapas
par(mfrow=c(2,2))
plot(hndvi,legend=F, main="High NDVI")
plot(hevi,legend=F,main="High EVI")
plot(hgndvi,legend=F,main="High GNDVI")
plot(hsavi,legend=F,main="High SAVI")
Posterior a esto, se obtiene un mapa en el que todos los índices se encuentran sobre la mediana mediante álgebra de mapas.
hvi<-hndvi*hevi*hgndvi*hsavi
##Saving results:
##writeRaster(hvi,"~/Maestria/Tesis/img/highvi.tif")De esta manera se obtuvieron las áreas de mejor calidad en términos de vegetación según los 4 índices obtenidos. Para obtener el área de estas zonas se realizó el análisis con la libreria landscapemetrics.Teniendo en cuenta las clases se obtuvieron mediante los índices de vegetación: Clase 0: Vegetación en mal estado Clase 1: Vegetación en buen estado
2.5. Clasificación supervisada
Para realizar la clasificación de las coberturas se tomarán como datos de referencia los obtenidos en el trabajo der Ruiz (sin publicar) para el año 2019, en el que se llegó clases por etapa sucesionales derivadas de cobertura de tierra, como se muestra en la tabla a continuación.
library(imager)## Warning: package 'imager' was built under R version 4.0.5## Loading required package: magrittr##
## Attaching package: 'magrittr'## The following objects are masked from 'package:terra':
##
## extract, inset## The following object is masked from 'package:raster':
##
## extract##
## Attaching package: 'imager'## The following object is masked from 'package:magrittr':
##
## add## The following objects are masked from 'package:terra':
##
## bbox, depth, fill## The following object is masked from 'package:raster':
##
## bbox## The following object is masked from 'package:sp':
##
## bbox## The following objects are masked from 'package:stats':
##
## convolve, spectrum## The following object is masked from 'package:graphics':
##
## frame## The following object is masked from 'package:base':
##
## save.imagetabla_clases<-load.image("~/Maestria/Tesis/CLC/Tabla_convenciones.jpg")
plot(tabla_clases, axes=F)
ref<-brick("~/Maestria/Tesis/CLC/CLC_18N.tif")
#reclassify values to succesion classes
temp<-c(6,7,4,5,2,1,8,10,11,9,3,1,2,3,4,5,6,7,8,9,10,11)
dim(temp)<-c(11,2)
temp## [,1] [,2]
## [1,] 6 1
## [2,] 7 2
## [3,] 4 3
## [4,] 5 4
## [5,] 2 5
## [6,] 1 6
## [7,] 8 7
## [8,] 10 8
## [9,] 11 9
## [10,] 9 10
## [11,] 3 11ref1<-reclassify(ref,temp)
classname<-c("Ta-Av","Ta-In","I-Av","I-In","MS","MD","T-Av","T-In","T-In/De","T-De","OC")
classcolor=c("#FB3636","#FB5A36","#FB7B36","#FBAA36","#FBD736","#DDFB36","#BCFB36","#8DFB36","#4EFB36","#36FB45","#878A87")
classdf <- data.frame(classvalue = c(1,2,3,4,5,6,7,8,9,10,11),
classnames = classname)
#checking reference map
plot(ref1, col=topo.colors(11), legend=F)
legend("right",legend=classname,col=topo.colors(11),pch = 15, bty = "n")
set.seed(100)
ptref<-sampleStratified(ref1,100, na.rm=T,sp=T)
ptref## class : SpatialPointsDataFrame
## features : 1100
## extent : 421721.5, 440251.5, 242856.5, 274521.5 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## variables : 2
## names : cell, CLC_18N
## min values : 110502, 1
## max values : 24414610, 11#extract values for locations sampled
sampvals<-raster::extract(masked,ptref,df=T)
#drop ID col
sampvals<-sampvals[,-1]
sampdata <- data.frame(class = ptref@data$CLC_18N, sampvals)
head(sampdata)Con el cuadro de datos obtenido se entrenará el clasificador
library(rpart)
# Train the model
cart <- rpart(as.factor(class)~., data=sampdata, method = 'class', minsplit = 5)
print(cart)## n=1096 (4 observations deleted due to missingness)
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 1096 996 1 (0.091 0.091 0.091 0.091 0.091 0.091 0.09 0.09 0.091 0.089 0.091)
## 2) ndvi>=0.5427877 520 423 1 (0.19 0.16 0.098 0.15 0.12 0.058 0.054 0.015 0.05 0.037 0.079)
## 4) gndvi>=0.6140694 117 68 1 (0.42 0.19 0.051 0.11 0.12 0.0085 0.0085 0 0 0.026 0.068) *
## 5) gndvi< 0.6140694 403 339 4 (0.12 0.15 0.11 0.16 0.12 0.072 0.067 0.02 0.065 0.04 0.082)
## 10) blue< 0.04775 349 288 4 (0.13 0.16 0.12 0.17 0.12 0.077 0.066 0.014 0.069 0.029 0.037)
## 20) red< 0.04685 121 87 2 (0.083 0.28 0.19 0.22 0.066 0.041 0.025 0.0083 0.041 0.0083 0.033) *
## 21) red>=0.04685 228 192 1 (0.16 0.096 0.088 0.15 0.14 0.096 0.088 0.018 0.083 0.039 0.039)
## 42) re1< 0.08735 83 60 1 (0.28 0.072 0.096 0.096 0.11 0.048 0.096 0.036 0.084 0.024 0.06) *
## 43) re1>=0.08735 145 119 4 (0.09 0.11 0.083 0.18 0.17 0.12 0.083 0.0069 0.083 0.048 0.028) *
## 11) blue>=0.04775 54 34 11 (0.037 0.056 0.037 0.056 0.13 0.037 0.074 0.056 0.037 0.11 0.37) *
## 3) ndvi< 0.5427877 576 485 8 (0.0052 0.033 0.085 0.04 0.066 0.12 0.12 0.16 0.13 0.14 0.1)
## 6) evi>=0.1130017 436 367 10 (0.0069 0.034 0.099 0.053 0.085 0.12 0.11 0.13 0.14 0.16 0.062)
## 12) blue< 0.05774999 392 335 9 (0.0077 0.036 0.1 0.054 0.087 0.13 0.12 0.13 0.15 0.12 0.069)
## 24) blue< 0.05225 284 237 6 (0 0.046 0.11 0.063 0.088 0.17 0.12 0.1 0.13 0.13 0.042) *
## 25) blue>=0.05225 108 86 8 (0.028 0.0093 0.074 0.028 0.083 0.028 0.13 0.2 0.18 0.1 0.14) *
## 13) blue>=0.05774999 44 23 10 (0 0.023 0.068 0.045 0.068 0.068 0.023 0.14 0.091 0.48 0) *
## 7) evi< 0.1130017 140 106 8 (0 0.029 0.043 0 0.0071 0.12 0.16 0.24 0.093 0.071 0.23)
## 14) blue< 0.0567 85 62 8 (0 0.024 0.059 0 0 0.14 0.22 0.27 0.13 0.071 0.082) *
## 15) blue>=0.0567 55 30 11 (0 0.036 0.018 0 0.018 0.091 0.073 0.2 0.036 0.073 0.45) *plot(cart, uniform=TRUE, main="Classification Tree")
text(cart, cex =0.8)
De acuerdo con el modelo obtenido, es posible ver como es posible llegar a 7 de las 11 clases posibles (1:Ta-Av,2:Ta-In,4:I-In,6:MD,8:T-In,10:T-De,11:OC), lo cual se encuentra limitado por los datos de cobertura de tierra que se tienen disponibles para la zona. Sin embargo, se realiza la clasificación de coberturas de la siguiente manera.
#running model to classify masked file
classified <- predict(masked, cart, type="class",na.rm = TRUE)
#checking new file
classified## class : RasterLayer
## dimensions : 3291, 1916, 6305556 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 421311, 440471, 241750, 274660 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 11 (min, max)
## attributes :
## ID value
## from: 1 1
## to : 11 11#checking values of classified
unique(classified)## [1] 1 2 4 6 8 10 11#saving file:
#writeRaster(classified,"~/Maestria/Tesis/img/classified.tif")2.6.Evaluación de exactitud
Se realizó un análisis de exactitud total, de usuario y de productor, con fines prácticos, sabiendo que el modelo con los datos disponibles resultaron siendo de baja calidad. Se realizó de acuerdo a la tecnica de validación cruzada de k-grupos con k=5, uno de los grupos se utiliza para probar el modelo y los demás para entrenar el modelo (de acuerdo a la GFC https://rspatial.org/raster/rs/5-supclassification.html)
#sampling for 5 model testing
library(dismo)## Warning: package 'dismo' was built under R version 4.0.5##
## Attaching package: 'dismo'## The following objects are masked from 'package:imager':
##
## circles, threshold## The following object is masked from 'package:terra':
##
## voronoiset.seed(10)
j <- kfold(sampdata, k = 5, by=sampdata$class)
#training models
x <- list()
for (k in 1:5) {
train <- sampdata[j!= k, ]
test <- sampdata[j == k, ]
cart <- rpart(as.factor(class)~., data=train, method = 'class', minsplit = 5)
pclass <- predict(cart, test, type='class')
# create a data.frame using the reference and prediction
x[[k]] <- cbind(test$class, as.integer(pclass))
}
#combining in a df and computing confusion matrix
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
conmat <- table(y)
conmat## predicted
## observed 1 2 4 5 6 7 8 10 11
## 1 57 11 20 7 0 0 0 3 2
## 2 33 20 20 5 2 1 5 10 4
## 3 15 16 15 7 4 4 14 22 3
## 4 29 17 21 8 4 1 3 13 4
## 5 26 4 22 7 4 2 8 18 9
## 6 7 4 9 10 4 6 24 30 6
## 7 6 4 13 5 4 6 29 28 5
## 8 1 0 4 7 5 9 35 31 8
## 9 5 5 11 6 4 7 25 34 3
## 10 6 2 3 10 5 5 20 39 10
## 11 15 3 2 5 1 4 23 19 28# number of cases
n <- sum(conmat)
n## [1] 1100# number of correctly classified cases per class
#sum of diagonal (not same number of colums and rows)
diag<-c(57,20,0,21,7,4,6,35,0,39,28)
# Overall Accuracy / Exactitud total
OA <- sum(diag) / nTambién se obtiene el estadístico de kappa, la exactitud del usuario y la del productor.
# observed (true) cases per class
rowsums <- apply(conmat, 1, sum)
# predicted cases per class
colsums <- apply(conmat, 2, sum)
colsums## 1 2 4 5 6 7 8 10 11
## 200 86 140 77 37 45 186 247 82#adding 0's for not existing columns
cs<-c(200,86,0,140,77,37,45,186,0,247,82)
p <- rowsums / n
q <- cs / n
#kappa
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
# Producer accuracy
PA <- diag / cs
# User accuracy
UA <- diag / rowsums
outt<- data.frame(classname,producerAccuracy = PA,userAccuracy = UA)3. Resultados
3.1 Índices de vegetación
Visualización de índices de vegetación
#visualization of vegetation index
par(mfrow=c(2,2))
plot(masked[[9]], col=rev(terrain.colors(10)),main="NDVI El Quimbo",axes=F)
plot(masked[[10]], col=rev(terrain.colors(10)),main="EVI El Quimbo",axes=F)
plot(masked[[11]], col=rev(terrain.colors(10)),main="GNDVI El Quimbo",axes=F)
plot(masked[[12]], col=rev(terrain.colors(10)),main="SAVI El Quimbo",axes=F)
#histograms of vegetation index
par(mfrow=c(2,2))
hist(masked[[9]],
main = "Distribución de valores de NDVI",
xlab = "NDVI",
ylab= "Frecuencia",
col = "#66CC99",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2%
## of the raster cells were used. 100000 values used.axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
hist(masked[[10]],
main = "Distribución de valores de EVI",
xlab = "EVI",
ylab= "Frecuencia",
col = "#66CC99",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2%
## of the raster cells were used. 100000 values used.axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
hist(masked[[11]],
main = "Distribución de valores de GNDVI",
xlab = "GNDVI",
ylab= "Frecuencia",
col = "#66CC99",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2%
## of the raster cells were used. 100000 values used.axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
hist(masked[[12]],
main = "Distribución de valores de SAVI",
xlab = "SAVI",
ylab= "Frecuencia",
col = "#66CC99",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')## Warning in .hist1(x, maxpixels = maxpixels, main = main, plot = plot, ...): 2%
## of the raster cells were used. 100000 values used.axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
Matriz de correlación entre índices
#correlation between layers
pairs(masked[[9:12]], hist=TRUE, cor=TRUE, use="pairwise.complete.obs", maxpixels=100000)
Como se observa en la matriz de correlación, los índices más relacionados entre sí son NDVI - GNDVI y EVI-SAVI, con coeficientes de correlación de 0,96 y 0,99 respectivamente. Esto se puede evidenciar también con los mapas y los histogramas de estos.
Zonas con mayores índices de vegetación
plot(hvi, main="High Vegetation Index Areas",legend=F)
Áreas con índices de vegetación altos
library(landscapemetrics)## Warning: package 'landscapemetrics' was built under R version 4.0.5#Class area (ha)
lsm_c_ca(hvi,directions = 8)#porcent of landscape
lsm_c_pland(hvi,directions = 8)#total area (ha)
lsm_l_ta(hvi,directions = 8)De esta manera, es evidente como en el 43% de la zona cuenta con vegetación en buen estado, según los 4 índices de vegetación analizados, correspondiente a 4792 ha.
3.2. Clasificación supervisada
Mapa de clasificación supervisada
#Plotting new classified raster
plot(classified, col=topo.colors(11), legend=F)
legend("right",legend=classname,col=topo.colors(11),pch = 15, bty = "n")
Matriz de confusión, exactitud total, de productor y de usuario.
#confusion matrix
conmat## predicted
## observed 1 2 4 5 6 7 8 10 11
## 1 57 11 20 7 0 0 0 3 2
## 2 33 20 20 5 2 1 5 10 4
## 3 15 16 15 7 4 4 14 22 3
## 4 29 17 21 8 4 1 3 13 4
## 5 26 4 22 7 4 2 8 18 9
## 6 7 4 9 10 4 6 24 30 6
## 7 6 4 13 5 4 6 29 28 5
## 8 1 0 4 7 5 9 35 31 8
## 9 5 5 11 6 4 7 25 34 3
## 10 6 2 3 10 5 5 20 39 10
## 11 15 3 2 5 1 4 23 19 28#kappa
kappa## [1] 0.117#overall accuracy
OA## [1] 0.1972727#Producer and user accuracy
outt#Se ve como el valor de la exactitud total es bajo, teniendo un rango de [0,1], lo cual se esperaba para el modelo como se explicó anteriormente. De igual manera, es posible observar valores bajos para el índice de kappa y las exactitudes de productor y usuario.
4. Discusión y conclusiones
A partir de los resultados obtenidos se observó como una buena parte de la zona de estudio cuenta con vegetación en buen estado, según los índices obtenidos, lo cual es un indicador del estado general de la zona de restauración, con un 43% del total del área. Sin embargo, se es necesario ajustar los datos mediante la verificación en campo, para obtener rangos de índices óptimos para la zona, con el fin de que con esta metodología se pueda realizar el monitoreo del avance de la restauración ecológica de la zona.
La clasificación por encima de la mediana coincidió con clasificaciones como las de Martunuzzi et al. (2008) y Zaitunah et al. (2018), en los cuales las zonas de bosques se encontraban con NDVI>0,5, siendo la obtenida en el presente trabajo de 0.54. Al afinar esta metodología se permitiría comparar el estudio con otros realizados en bosque seco tropical teninendo en cuenta su actual estado de amenaza a nivel mundial (Jenzen, 1988).
En general, se realizó una aproximación adecuada, con imágenes de resolución espacial relativamente buena (20 m), desde la cual podrían realizarse trabajan a mayor escala para tener un panorama del BST en en Valle del Magdalena, por ejemplo, además de servir como indicadores proxy de otras características del bosque como los servicios ecosistémicos que brinda a distintos usuarios, trabajo que tendría que complementarse con datos de campo y otra información espacial.
Respecto a la clasificación supervisada se resalta la importancia de los datos de referencia o entrenamiento, ya que con los que se tenían para la zona contaban con una clasificación que no era la mjeor para realizar este tipo de procesos, por lo cual se recomienda tener mapas de coberturas, no de estados sucesionales como en este caso, ya que se llegan a clasificaciones poco exactas, como fue claro desde la ejecución inicial del modelo de clasificación en donde no se llegó a incluir todas las clases.
5. Referencias
Martinuzzi, S., Gould, W. A., Ramos Gonzalez, O. M., Martinez Robles, A., Calle Maldonado, P., Pérez-Buitrago, N., & Fumero Caban, J. J. (2008). Mapping tropical dry forest habitats integrating Landsat NDVI, Ikonos imagery, and topographic information in the Caribbean Island of Mona. Revista de Biología Tropical, 56(2), 625-639. Janzen, D. H. (1988). Tropical Dry Forests The Most Endangered Major Tropical Ecosystem. En Biodiversity. National Academies Press (US). https://www.ncbi.nlm.nih.gov/books/NBK219281/ Zaitunah, A., Ahmad, A. G., & Safitri, R. A. (2018). Normalized difference vegetation index (ndvi) analysis for land cover types using landsat 8 oli in besitang watershed, Indonesia. In IOP Conference Series: Earth and Environmental Science (Vol. 126, No. 1, p. 012112). IOP Publishing.
6. Reproducibilidad
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252 LC_CTYPE=Spanish_Colombia.1252
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] landscapemetrics_1.5.2 dismo_1.3-3 rpart_4.1-15
## [4] imager_0.42.8 magrittr_2.0.1 terra_1.0-10
## [7] raster_3.4-5 sp_1.4-5 leaflet_2.0.4.1
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 pillar_1.5.0 compiler_4.0.4 highr_0.8
## [5] tools_4.0.4 digest_0.6.27 lifecycle_1.0.0 jsonlite_1.7.2
## [9] evaluate_0.14 tibble_3.0.6 lattice_0.20-41 pkgconfig_2.0.3
## [13] png_0.1-7 rlang_0.4.10 igraph_1.2.6 crosstalk_1.1.1
## [17] yaml_2.2.1 rgdal_1.5-23 xfun_0.21 readbitmap_0.1.5
## [21] bmp_0.3 stringr_1.4.0 knitr_1.31 htmlwidgets_1.5.3
## [25] vctrs_0.3.6 grid_4.0.4 glue_1.4.2 R6_2.5.0
## [29] jpeg_0.1-8.1 fansi_0.4.2 rmarkdown_2.7 purrr_0.3.4
## [33] ellipsis_0.3.1 codetools_0.2-18 htmltools_0.5.1.1 tiff_0.1-8
## [37] utf8_1.1.4 stringi_1.5.3 crayon_1.4.1