Clasificación de coberturas en la Reserva Natural Thomas van der Hammen
Luis H. Romero
14/10/2019
Introducción
En el marco de la Agenda 2030 sobre el Desarrollo Sostenible, se plantearon 17 objetivos de Desarrollo Sostenible con los que se busca mejorar la calidad de vida de la sociedad humana, sin importar sus condiciones socioeconómicas, ubicación, creencias, etc. Los objetivos van desde la erradicación de la pobreza hasta el combate contra el cambio climático, la educación, la igualdad de la mujer, la defensa del medio ambiente y el diseño de nuestras ciudades. - ONU
En lo que respecta a la calidad de vida de las poblaciones urbanas, pueden orientarse varios de los objetivos del desarrollo sostenible para que se conviertan en directrices de crecimiento y planeación de las ciudades. Actualmente, las ciudades necesitan de proyectos, no solo de crecimiento y expansión por la cada vez mayor densidad poblacional que se evidencia en ciertos centros urbanos principales, sino de la construcción de un ambiente sano, seguro, eficiente y limpio, siendo las áreas naturales importantes para alcanzar estos objetivos.
En este trabajo se abordarán principalmente 3 objetivos del desarrollo sostenible: El objetivo número 11 “Lograr que las ciudades y los asentamientos humanos sean inclusivos, seguros, resilientes y sostenibles”; el objetivo 13 “Adoptar medidas urgentes para combatir el cambio climático y sus efectos” y el objetivo 15 “Gestionar sosteniblemente los bosques, luchar contra la desertificación, detener e invertir la degradación de las tierras y detener la pérdida de la biodiversidad”. Estos objetivos serán evaluados en la Reserva Natural Thomas van der Hammen, la cual ha sido centro de debate en múltiples discusiones sobre urbanismo en la ciudad de Bogotá, Colombia. Este estudio busca hacer un análisis espacial de las coberturas de la reserva TVDH, comparando dos períodos de tiempo diferentes e identificar la transformación de las coberturas como evidencia para soportar la idea de que la urbanización de la reserva conllevará a la pérdida de diferentes servicios ecosistémicos valiosos para la ciudad de Bogotá.
Datos y métodos
Área de estudio
La reserva natural Thomás Van der Hammen fue declarada elemento fundamental de la Estructura Ecológica Principal de Bogotá, declarada y regulada por los Acuerdos CAR 011 de 2011 y 021 de 2014 (Reserva forestal del norte de Bogotá, 2017). Cuenta con los últimos relictos de bosques andinos de la sabana de Bogotá y su área de influencia funciona como un elemento importante de conectividad entre La Sabana, el río Bogotá, los cerros Orientales y los cerros de Cota. La reserva se ubica al norte de la ciudad de Bogotá, entre las localidades de Usaquén y Suba, compuesto por 1395,16 hectáreas donde se pueden encontrar el bosque de Las Mercedes, parte del humedal La Conejera, la quebrada La Salitrosa y varios canales que alimentan al río Bogotá y al humedal Guaymaral (Puentes, 2016). A pesar de su declaración como reserva, en diferentes períodos de alcaldías en el distrito se ha planteado la propuesta de utilizar los predios para la construcción de viviendas y ampliación de vías que atraviesen sus terrenos. Estos proyectos urbanos son impulsados bajo el supuesto de que la reserva no tiene las características para regular los procesos ecológicos que se le atribuyen, llamando de manera despectiva esta zona de importancia ambiental como potreros sin ningún valor ecológico, debido a la intervención de industrias, zonas residenciales, agricultura y pastos.
Ya se han realizado múltiples investigaciones que demuestran el elevado potencial ecológico de la reserva, encontrando alrededor de 28 especies de aves migratorias, 153 especies de aves residentes nativas (Fundación humedales Bogotá, 2016) y cuatro especies de mamíferos nativos (Mendoza & Sánchez, 2014), además de representar el último pulmón de la ciudad con la mejor representatividad de flora y fauna de la sabana de Bogotá.
Métodos
Se utilizaron dos imágenes satelitales, una Landsat 5 tomada el 22 de marzo de 1988 y una Sentinel 2B tomada el 24 de enero de 2018, con el fin de evidenciar los cambios que han ocurrido en la reserva en los últimos 30 años. Las imágenes fueron corregidas cambiando los Niveles Digitales por Reflectancia Top pf Atmosphere (TOA) y no fue necesario hacer corrección de nubes, puesto que ambas imágenes utilizadas tenían poca nubosidad en el área de estudio. La imagen Landsat fue descargada de la página earthexplorer y la sentinel se descargó la página copernicus. Se realizó para cada una de las imágenes la clasificación supervisada y no supervisada para la identificación de seis coberturas: 1) Zonas naturales, 2) Zonas artificiales, 3) Pastos, 4) Cultivos, 5) Cuerpos de agua y 6) vías. No hubo ningún inconveniente con la imagen Sentinel, pero la resolución espacial de la imagen Landsat no permitió delimitar correctamente las vías destapadas que se encuentran en los límites de la reserva por lo que esta categoría fue excluida para los análisis con Landsat 5.
La clasificación no supervisada se realizó con base en tres índices de vegetación: el de diferencia normalizada (NDVI), el de suelo ajustado (SAVI) y el perpendicular (PVI), los cuales fueron los parámetros de ajustes de clasificación. La clasificación supervisada utilizó datos de entrenamiento obtenidos por polígonos seleccionados manualmente en las escenas que abarcaban el área de la reserva, realizándose este paso con ayuda del software ArcGIS 10.5. Para cada uno de los períodos estudiados, debido a la diferencia de resolución espacial, se utilizaron puntos de entrenamiento diferentes, abarcando las coberturas mencionadas anteriormente.
Resultados
Cargue de archivos del área de interés.
Se buscó la imagen satelital más reciente que abarcara el área de la reserva y que el porcentaje de nubosidad fuera lo más baja posible, encontrando una imágen Sentinel 2B tomada el 24 de enero de 2018, con un porcentaje de nubosidad de 5.12% y una imagen antigua para determinar el cambio de coberturas en dos instantes de tiempo, encontrando una imagen Landsat 5 tomada el 22 de marzo de 1988.
# Definición de ruta de archivos, bandas S2B de 10 m de
# resolución espacial
lista.tif_10 <- paste0("Sentinel2/S2B_MSIL1C_20180124T152629_N0206_R025_T18NWL_20180124T202241.SAFE/GRANULE/L1C_T18NWL_A004627_20180124T152632/RT_T18NWL_20180124T152629_B0",
c(2, 3, 4, 8), ".tif")
# Cargue de bandas
bog_S2_10 <- stack(lista.tif_10)
# Cambio de nombre de las bandas espectrales
names(bog_S2_10) <- c("Blue", "Green", "Red", "NIR")
bog_S2_10## class : RasterStack
## dimensions : 10980, 10980, 120560400, 4 (nrow, ncol, ncell, nlayers)
## resolution : 10, 10 (x, y)
## extent : 499980, 609780, 490200, 6e+05 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : Blue, Green, Red, NIR# Ahora para la imagen Landsat 5
lista.tif_30 <- paste0("Landsat/LT04_L1TP_008057_19880322_20170209_01_T1_B",
c(1, 2, 3, 4, 5, 7), ".tif")
bog_land_30 <- stack(lista.tif_30)
names(bog_land_30) <- c("Azul", "Verde", "Rojo", "NIR", "SWIR1",
"SWIR2")
bog_land_30## class : RasterStack
## dimensions : 6901, 7731, 53351631, 6 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 441585, 673515, 374685, 581715 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : Azul, Verde, Rojo, NIR, SWIR1, SWIR2
## min values : 0, 0, 0, 0, 0, 0
## max values : 255, 255, 255, 255, 255, 255
Como se evidencia en las salidas anteriores, las imágenes landsat están con la información cruda captada por el sensor en niveles digitales, por lo que debe realizarse una transformación de estos valores para trabajar con información real de la superficie de estudio. Aunque no es la información de mayor calidad (que en sensores remotos correspondería a reflectancia del fondo de la atmósfera, BOA), se trabajará con la reflectancia del tope de la atmósfera, TOA. Para la imagen sentinel se realizó la transformación por medio la herramienta Semiauthomatic Classification Plugin de QGIS 3.4, mientras que las imágenes landsat se transaformarion mediante la función de radCor de la librería RStoolbox.
mtlfile <- "Landsat/LT04_L1TP_008057_19880322_20170209_01_T1_MTL.txt"
metadata <- readMeta(mtlfile)
lsat <- stackMeta(metadata)
lsat_cal <- radCor(lsat, metaData = metadata, method = "apref")
lsat_cal <- dropLayer(lsat_cal, 6)
writeRaster(lsat_cal, "lsat_cal_1988.tif", overwrite = TRUE)
names(lsat_cal) <- c("Blue", "Green", "Red", "NIR", "SWIR1",
"SWIR2")
lsat_cal## class : RasterStack
## dimensions : 6901, 7731, 53351631, 6 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 441585, 673515, 374685, 581715 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## names : Blue, Green, Red, NIR, SWIR1, SWIR2
## min values : 0, 0, 0, 0, 0, 0
## max values : 0.3423814, 0.7435552, 0.6653971, 0.8323059, 0.5509866, 0.7758228Ploteo de bandas color natural y falso color de la imagen
Aquí se muestra una primera observación de las imágenes descargadas, con dos combinaciones diferentes de bandas: Color natural, la cual combina las bandas 4 (Red), 3 (Green) y 2 (Blue) para Landsat y Sentinel 2B y la composición de falso color con las bandas 8 (Near Infra Red, NIR), 4 y 3 para Sentinel 2B y las bandas 5 (NIR), 4 y 3 para Landsat. En la imagen puede apreciarse una cobertura parcial del departamento de Cundinamarca y en la esquina inferior derecha se evidencia gran parte del casco urbano de Bogotá, Distrito Capital.
plotRGB(bog_S2_10, r = 3, g = 2, b = 1, axes = T, stretch = "lin",
main = "Composición Color natural Sentinel 2B 2018")
plotRGB(bog_S2_10, r = 4, g = 3, b = 2, axes = T, stretch = "lin",
main = "Composición Falso Color Sentinel 2B 2018")
plotRGB(lsat_cal, r = 3, g = 2, b = 1, axes = T, stretch = "lin",
main = "Composición Color natural Landsat 1988")
plotRGB(lsat_cal, r = 4, g = 3, b = 2, axes = T, stretch = "lin",
main = "Composición Color natural Landsat 1988")
Relación entre bandas
Una matriz scatterplot puede ser útil para explorar la relación entre las capas estudiadas. Con esta información podemos determinar la redundancia de la información y el porcentaje de variabilidad explicada por cada unas de las bandas. Se observa como las bandas tanto de Sentinel como de Landsat están muy correlacionadas, sobre todo con aquellas que se distribuyen en el espectro visible. La banda NIR es la que demuestra tener mayor variabilidad en los datos, por lo que es ampliamente utilizada al momento de calcular índices. En los datos landsat, al incluir las bandas de onda corta ubicadas en valores mayores del espectro electromagnético, varía bastante con respecto a las bandas del espectro visible y el NIR.
pairs(bog_S2_10, main = "Análisis multi-banda Sentinel 2B")
lsat88 <- stack("D:/Magister en Geomática/2 Semestre/Percepción Remota Avanzada Teoría/Primera Entrega/lsat_cal_1988.tif")
names(lsat88) <- c("Azul", "Verde", "Rojo", "NIR", "SWIR1", "SWIR2")
pairs(lsat88, main = "Análisis multi-banda Landsat")
Recorte del área de estudio
Se realizó un recorte cudrado en la zona de Bogotá para que los procesos de clasificación de coberturas demanden menos tiempo de procesamiento, para que al final de cada proceso se haga un nuevo recorte con el shapefile que delimita la reserva TVDH. Se optó por esta medida porque al utilizar desde un principio el Shapefile de la reserva se generaban errores con el procesamiento de la imagen, los cuales se solucionaban al trabajar con el siguiente recorte cuadrado de la zona.
` ## Recorte del área de estudio
Se realizó un recorte cudrado en la zona de Bogotá para que los procesos de clasificación de coberturas demanden menos tiempo de procesamiento, para que al final de cada proceso se haga un nuevo recorte con el shapefile que delimita la reserva TVDH. Se optó por esta medida porque al utilizar desde un principio el Shapefile de la reserva se generaban errores con el procesamiento de la imagen, los cuales se solucionaban al trabajar con el siguiente recorte cuadrado de la zona.
## Llamada del Shapefile
reserva_lim <- shapefile("Shapefiles/Reserva")
## Reproyección para que las coordenadas de ambos archivos
## coincidan para el recorte
reserva_proj <- spTransform(reserva_lim, crs(bog_S2_10))
## Raster recortado
reserva_10_clip <- crop(bog_S2_10, extent(reserva_proj))
## Raster recortado con shapefile como máscara
reserva_10_mask <- mask(reserva_10_clip, reserva_proj)
plotRGB(reserva_10_mask, main = "Imagen Sentinel 2B 2018 con máscara de shapefile = 10 m",
r = 3, g = 2, b = 1, stretch = "lin", axes = TRUE, ylab = "Norte",
xlab = "Este")
## Ahora con la imagen Landsat 5
reserva_30_clip <- crop(lsat_cal, extent(reserva_proj))
reserva_30_mask <- mask(reserva_30_clip, reserva_proj)
plotRGB(reserva_30_mask, main = "Imagen Landsat 5 1988 con máscara de shapefile = 30 m",
r = 3, g = 2, b = 1, stretch = "lin", axes = TRUE, ylab = "Norte",
xlab = "Este")
Clasificación No Supervisada
Mediante el método de Lloyd, se realizó la clasificación no supervisada de la zona de estudio. Con este método se particiona un conjunto de n observaciones en k clusters, tomando en cuenta la varianza interna de cada uno de los clusters. Los objetos se clasifican de manera que la suma de las varianzas internas en los clusters sean lo menor posible. En este caso, se tomaron en cuenta 7 clusters, de acuerdo con las coberturas que ofrecen mayor información sobre la zona de estudio. Se intentaron utilizar los 3 índices de vegetación mencionados anteriormente, pero al momento de la generación de clusters se generaron errores que no pudieron solucionarse, por lo que solo se utilizó el índice de vegetación diferencial normalizado.
ndvi <- (reserva_10_mask[["NIR"]] - reserva_10_mask[["Red"]])/(reserva_10_mask[["NIR"]] +
reserva_10_mask[["Red"]])
nr <- getValues(ndvi) ## Conversión de raster a vector/matriz
str(nr)## num [1:980100] 0.715 0.713 0.697 0.679 0.701 ...set.seed(99)
# 10 clusters, 500 iteraciones, empezar con 5 sets aleatorios
# utilizando el método Lloyd
kmcluster <- kmeans(na.omit(nr), centers = 12, iter.max = 500,
nstart = 5, algorithm = "Lloyd")
str(kmcluster)## List of 9
## $ cluster : int [1:980100] 11 11 11 11 11 11 11 12 12 11 ...
## $ centers : num [1:12, 1] 0.1519 0.0774 0.4875 0.4055 0.2372 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:12] "1" "2" "3" "4" ...
## .. ..$ : NULL
## $ totss : num 72257
## $ withinss : num [1:12] 37.5 39.1 37.5 39 39.5 ...
## $ tot.withinss: num 484
## $ betweenss : num 71774
## $ size : int [1:12] 70026 95335 69382 67553 69045 80042 66696 79412 63419 113225 ...
## $ iter : int 305
## $ ifault : NULL
## - attr(*, "class")= chr "kmeans"# objetos NDVI para establecer los valores de cluster en un
# nuevo raster
knr <- raster(ndvi)
values(knr) <- kmcluster$cluster
knr## class : RasterLayer
## dimensions : 990, 990, 980100 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 598020, 607920, 522600, 532500 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 1, 12 (min, max)mycolor <- c("#006633", "#333333", "#00ff00", "#00ffff", "#0000ff",
"#ff3399", "#ff9900") #(1-Zonas naturales, 2-vías, 3-pastos, 4-agua artificial, 5-cuerpos de agua, 6-cultivos, 7-Zonas artificiales)
plot(ndvi, col = rev(terrain.colors(7)), main = "Sentinel 2B-NDVI")
plot(knr, main = "Clasificación No Supervisada Sentinel 2b",
col = mycolor)
ndvi88 <- (reserva_30_mask[["NIR"]] - reserva_30_mask[["Red"]])/(reserva_30_mask[["NIR"]] +
reserva_30_mask[["Red"]])
nr88 <- getValues(ndvi88) ## Conversión de raster a vector/matriz
str(nr88)## num [1:108900] 0.503 0.441 0.434 0.496 0.503 ...set.seed(99)
# 10 clusters, 500 iteraciones, empezar con 5 sets aleatorios
# utilizando el método Lloyd
kmcluster88 <- kmeans(na.omit(nr88), centers = 12, iter.max = 500,
nstart = 5, algorithm = "Lloyd")
str(kmcluster88)## List of 9
## $ cluster : int [1:108900] 3 4 4 3 3 4 4 4 11 4 ...
## $ centers : num [1:12, 1] 0.367 0.73 0.518 0.42 0.612 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:12] "1" "2" "3" "4" ...
## .. ..$ : NULL
## $ totss : num 2650
## $ withinss : num [1:12] 2.08 3.45 2.4 2.26 2.82 ...
## $ tot.withinss: num 32.6
## $ betweenss : num 2618
## $ size : int [1:12] 8310 6351 13380 10242 13154 6050 10291 6968 14241 3305 ...
## $ iter : int 76
## $ ifault : NULL
## - attr(*, "class")= chr "kmeans"# objetos NDVI para establecer los valores de cluster en un
# nuevo raster
knr88 <- raster(ndvi88)
values(knr88) <- kmcluster88$cluster
knr88## class : RasterLayer
## dimensions : 330, 330, 108900 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 598005, 607905, 522615, 532515 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## source : memory
## names : layer
## values : 1, 12 (min, max)mycolor <- c("#006633", "#333333", "#00ff00", "#00ffff", "#0000ff",
"#ff3399", "#ff9900") #(1-Zonas naturales, 2-vías, 3-pastos, 4-agua artificial, 5-cuerpos de agua, 6-cultivos, 7-Zonas artificiales)
plot(ndvi88, col = rev(terrain.colors(7)), main = "Landsat 5-NDVI")
plot(knr88, main = "Clasificación No Supervisada Landsat 5",
col = mycolor)
Reclasificación
Al momento de realizar la clasificación no supervisada se definieron 12 clusters los cuales fueron reclasificados en ArcGIS con la herramienta Raster calculator dejando las seis coberturas de estudio para las imágenes Sentinel y cinco coberturas (excluyendo vías) para Landsat.
mycolor1 <- c("#00ff00", "#ff3399", "#006633", "#ff9900", "#0000ff") # Pastos, Cultivos, Z. naturales, Z. artificiales, Agua
mycolor2 <- c("#0000ff", "#ff9900", "#ff3399", "#00ff00", "#333333",
"#006633")
reclas2018 <- stack("CNSR_S2b.tif")
reclas1988 <- stack("CNSR_lsat.tif")
TVDH_lim <- shapefile("Shapefiles/TVDH")
TVDH_proj <- spTransform(TVDH_lim, crs(reclas2018))
TVDH_2018_clip <- crop(reclas2018, extent(TVDH_proj))
TVDH_2018_mask <- mask(TVDH_2018_clip, TVDH_proj)
TVDH_proj <- spTransform(TVDH_lim, crs(reclas1988))
TVDH_1988_clip <- crop(reclas1988, extent(TVDH_proj))
TVDH_1988_mask <- mask(TVDH_1988_clip, TVDH_proj)
plot(TVDH_1988_mask, main = "Reclasificación Landsat 5 1988",
col = mycolor1)
plot(TVDH_2018_mask, main = "Reclasificación Sentinel 2b 2018",
col = mycolor2)
Como se aprecia en la clasificación anterior, una de las tendencias de cambio al interior de la reserva es la del crecimiento de áreas artificiales (color naranja) y una pérdida de las áreas clasificadas como zonas naturales (verde oscuro) y cuerpos de agua (azul). Se ha incrementado las áreas de pastos (verde claro) y cultivos (rosa).
Clasificación Supervisada CART para Sentinel
Esta clasificación se desarrolló a partir de polígonos de entrenamiento generados en ArcGIS. Cada imagen cuenta con sus propios polìgonos de entrenamiento y se presentan de manera independientemente para no poner demasiada cantidad de líneas en los chunk.
Parámetros iniciales de clasificación Sentinel
nlcd <- brick("D:/Magister en Geomática/2 Semestre/Percepción Remota Avanzada Teoría/Primera Entrega/CSReserva2.tif")
names(nlcd) <- "Clasificacion2018"
nlcdclass <- c("Cuerpos de agua", "Pastos", "Zonas artificiales",
"Cultivos", "Zonas naturales", "Agua artifical", "Vías")
classdf <- data.frame(classvalue1 = c(1, 3, 6, 8, 14, 16, 26),
classnames1 = nlcdclass)
classcolor <- c("#0000ff", "#00ff00", "#ff9900", "#ff3399", "#006633",
"#0000ff", "#333333")
# Now we ratify (RAT = 'Raster Attribute Table') the ncld2011
# (define RasterLayer as a categorical variable). This is
# helpful for plotting.
cobTVDH <- nlcd[[1]]
cobTVDH <- ratify(cobTVDH)
rat <- levels(cobTVDH)[[1]]
rat$landcover <- nlcdclass
levels(cobTVDH) <- ratGeneración de sitios de muestreo Sentinel
set.seed(99)
# Sampling
samp2018 <- sampleStratified(cobTVDH, size = 200, na.rm = TRUE,
sp = TRUE)
table(samp2018$Clasificacion2018)##
## 1 3 6 8 14 16 26
## 200 200 200 200 200 200 200plt <- levelplot(cobTVDH, col.regions = classcolor, main = "Distribución de sitios de entrenamiento 2018")
print(plt + layer(sp.points(samp2018, pch = 3, cex = 0.5, col = 1)))
Parámetros iniciales de clasificación Landsat
nlcd88 <- brick("D:/Magister en Geomática/2 Semestre/Percepción Remota Avanzada Teoría/Primera Entrega/CS88.tif")
names(nlcd88) <- "Clasificacion1998"
nlcd88class <- c("Zonas naturales", "Zonas artificiales", "Cultivos",
"Pastos", "Cuerpos de agua")
class88df <- data.frame(classvalue1 = c(17, 21, 24, 29, 40),
classnames1 = nlcd88class)
classcolor1 <- c("#006633", "#ff9900", "#ff3399", "#00ff00",
"#0000ff")
# Now we ratify (RAT = 'Raster Attribute Table') the ncld2011
# (define RasterLayer as a categorical variable). This is
# helpful for plotting.
cobTVDH88 <- nlcd88[[1]]
cobTVDH88 <- ratify(cobTVDH88)
rat88 <- levels(cobTVDH88)[[1]]
rat88$landcover <- nlcd88class
levels(cobTVDH88) <- rat88Generación de sitios de muestreo Landsat
set.seed(99)
# Muestreo
samp1988 <- sampleStratified(cobTVDH88, size = 200, na.rm = TRUE,
sp = TRUE)
table(samp1988$Clasificacion1998)##
## 17 21 24 29 40
## 200 200 200 200 200plt <- levelplot(cobTVDH88, col.regions = classcolor1, main = "Distribución de sitios de entrenamiento 1988")
print(plt + layer(sp.points(samp1988, pch = 3, cex = 0.5, col = 1)))
Extraer valores para sitio de muestreo
A partir de aquí se muestran conjuntamente Landsat y Sentinel
# Extract the layer values for the locations
sampvals <- extract(reserva_10_mask, samp2018, df = TRUE)
# sampvals no longer has the spatial information. To keep the
# spatial information you use `sp=TRUE` argument in the
# `extract` function. drop the ID column
sampvals <- sampvals[, -1]
# combine the class information with extracted values
sampdata <- data.frame(classvalue = samp2018@data$Clasificacion2018,
sampvals)
sampvals88 <- extract(reserva_30_mask, samp1988, df = TRUE)
sampvals88 <- sampvals88[, -1]
sampdata88 <- data.frame(classvalue = samp1988@data$Clasificacion1998,
sampvals88)Entrenamiento del clasificador
# Train the model
cart <- rpart(as.factor(classvalue) ~ ., data = sampdata, method = "class",
minsplit = 5)
# print(model.class) Plot the trained classification tree
plot(cart, uniform = TRUE, main = "Árbol de clasificación 2018")
text(cart, cex = 0.8)
# Train the model
cart88 <- rpart(as.factor(classvalue) ~ ., data = sampdata88,
method = "class", minsplit = 5)
plot(cart88, uniform = TRUE, main = "Árbol de clasificación 1988")
text(cart88, cex = 0.8)
Clasificación
pr2019 <- predict(reserva_10_mask, cart, type = "class")
pr2019 <- ratify(pr2019)
rat <- levels(pr2019)[[1]]
rat$legend <- classdf$classnames1
levels(pr2019) <- rat
levelplot(pr2019, maxpixels = 1e+06, col.regions = classcolor,
scales = list(draw = FALSE), main = "Decisión de árbol de clasificación para sentinel 2B 2018")
mycolor3 <- c("#006633", "#ff9900", "#ff3399", "#00ff00", "#0000ff")
pr2019_88 <- predict(reserva_30_mask, cart88, type = "class")
pr2019_88 <- ratify(pr2019_88)
rat88 <- levels(pr2019_88)[[1]]
rat88$legend <- class88df$classnames1
levels(pr2019_88) <- rat88
levelplot(pr2019_88, maxpixels = 1e+06, col.regions = mycolor3,
scales = list(draw = FALSE), main = "Decisión de árbol de clasificación para Landsat 1988")
Comparación de resultados
TVDH_lim <- shapefile("Shapefiles/TVDH")
TVDH_proj_2 <- spTransform(TVDH_lim, crs(reserva_10_mask))
TVDH_pr2019_clip <- crop(pr2019, extent(TVDH_proj_2))
TVDH_pr2019_mask <- mask(TVDH_pr2019_clip, TVDH_proj_2)
TVDH_proj_3 <- spTransform(TVDH_lim, crs(reserva_30_mask))
TVDH_pr2019_88_clip <- crop(pr2019_88, extent(TVDH_proj_3))
TVDH_pr2019_88_mask <- mask(TVDH_pr2019_88_clip, TVDH_proj_3)
par(mfrow = c(1, 2))
plot(TVDH_pr2019_mask, main = "Clasificación Supervisada 2018",
col = mycolor1)
plot(TVDH_pr2019_88_mask, main = "Clasificación Supervisada 1988",
col = mycolor2)
En ambas clasificaciones se pudo apreciar como la imagen Landsat tiende a sobrestimar las coberturas de zonas naturales, especialmente en la clasificación no supervisada. Dado que la clasificación supervisada se basa en datos en los que hay certeza de la identidad de las coberturas, se optará por elegir este método para realizar las conclusiones finales.
Evaluación del Modelo
## K-fold cross validation
set.seed(99)
j <- kfold(sampdata, k = 5, by = sampdata$classvalue)
table(j)## j
## 1 2 3 4 5
## 280 280 280 280 280x <- list()
for (k in 1:5) {
train <- sampdata[j != k, ]
test <- sampdata[j == k, ]
cart <- rpart(as.factor(classvalue) ~ ., data = train, method = "class",
minsplit = 5)
pclass <- predict(cart, test, type = "class")
# crear data.frame usando la referencia y predicción
x[[k]] <- cbind(test$classvalue, as.integer(pclass))
}
# Combinación de las cinco listas de elementos en un
# data.frame
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c("observed", "predicted")
conmat <- table(y)
# change the name of the classes
colnames(conmat) <- classdf$classnames
rownames(conmat) <- classdf$classnames
conmat## predicted
## observed Cuerpos de agua Pastos Zonas artificiales Cultivos
## Cuerpos de agua 184 0 2 6
## Pastos 0 159 0 30
## Zonas artificiales 3 4 151 10
## Cultivos 12 61 14 60
## Zonas naturales 0 1 0 7
## Agua artifical 10 0 8 0
## Vías 0 3 19 14
## predicted
## observed Zonas naturales Agua artifical Vías
## Cuerpos de agua 0 8 0
## Pastos 0 0 11
## Zonas artificiales 7 10 15
## Cultivos 30 1 22
## Zonas naturales 192 0 0
## Agua artifical 0 182 0
## Vías 0 0 164Exactitud general del modelo (Estadístico Kappa)
n <- sum(conmat)
n## [1] 1400# número de casos correctamente clasificados por clases
diag <- diag(conmat)
# Overall Accuracy
OA <- sum(diag)/n
OA## [1] 0.78## Kappa
# Casos de verdad observada por clases
rowsums <- apply(conmat, 1, sum)
p <- rowsums/n
# Casos predichos por clase
colsums <- apply(conmat, 2, sum)
q <- colsums/n
expAccuracy <- sum(p * q)
kappa <- (OA - expAccuracy)/(1 - expAccuracy)
kappa## [1] 0.7433333# Exactitud del productor
PA <- diag/colsums
# User accuracy
UA <- diag/rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc## producerAccuracy userAccuracy
## Cuerpos de agua 0.8803828 0.920
## Pastos 0.6973684 0.795
## Zonas artificiales 0.7783505 0.755
## Cultivos 0.4724409 0.300
## Zonas naturales 0.8384279 0.960
## Agua artifical 0.9054726 0.910
## Vías 0.7735849 0.820
Los métodos utilizados en el trabajo presentan en general una buena respuesta ante las coberturas que se encuentran presentes en la reserva TVDH exceptuando los cultivos. Estas coberturas, al ser parte de mozaicos con pastos, tiende a confundirse con estos, sobre todo aquellos que en su respuesta espectral demuestran tener alto contenido de humedad. Esto se evidencia en la matriz de confusiòn en donde se indican 61 pixeles de cultivo clasficados como pasto. Las vías fueron muy difíciles de delimitar puesto que en la reserva el ancho máximo encontrado fue de aproximadamente 8 metros, lo que no alcanza a ser cubierto por las resoluciones espaciales utilizadas (10m y 30m).
Las diferencias encontradas en ambas clasificaciones también pueden verse maximizadas por las diferencias entre las resoluciones espaciales. Habría sido ideal realizar esta actividad con dos imágenes sentinel para obtener un mejor nivel de detalle, pero la misión en que fue lanzado este satélite es muy reciente como para evidenciar cambios. Se pudo haber utilizado dos imágenes Landsat, pero se quiso evidenciar la utilidad que tiene el programa RStudio para el manejo de datos de diferentes sensores remotos y las diferencias y problemàticas que pueden arrojar cada uno de ellos.
Conclusiones
A pesar de que la reserva efectivamente tiene una proporción considerable de terrenos ocupados por mozaicos de pastos y cultivos, el principal cambio observado durante los últimos 30 años se evidencia en la expansión de los territorios artificializados como agrupaciones de viviendas, colegios y construcciones para la práctica agrícola. Pero este cambio no solo se evidencia dentro de la reserva, sino de una manera más profunda en los alrededores, sobretodo en la parte sur de la reserva.
Las transformaciones en la periferia de la reserva son de igual manera importantes de evaluar, puesto que la reserva representa un área potencial para recuperar los procesos de conectividad ecológica entre los Cerros Orientales y el valle aluvial del río Bogotá que puede evidenciarse en las imágenes utilizadas en este trabajo (Torres & Niño, 2018). En este mismo trabajo de Torres & Niño se menciona lo que es evidente de estos análisis y es la profunda transformación antrópica que tiene la reserva. Son muy pocos los espacios naturales que aún se conservan y cada vez más se va viendo afectada por la agricultura y la ganadería. No obstante, esto NO debe considerarse como una excusa para seguir transformando la reserva y continuar con el patrón de urbanización que la ha consumido en los últimos 30 años.
Las zonas que actualmente se clasificaron como pastos y cultivos pueden ser objetos de restauraciòn y uso sostenible para darle un mayor valor ecológico a la reserva. Lograr este objetivo nos acerca más a la realización del objetivo número 11 del desarrollo sostenible, “Lograr que las ciudades y los asentamientos humanos sean inclusivos, seguros, resilientes y sostenibles”. De igual forma, al ser una zona potencial para la regulaciòn de la contaminación ambiental, regulación de la temperatura de la ciudad y en general de la calidad de vida de los bogotanos, una adecuada gestión de la reserva TVDH conduce al Distrito Capital hacia el objetivo 13 de desarrollo sostenible, “Adoptar medidas urgentes para combatir el cambio climático y sus efectos”.
Finalmente, lograr transformar las zonas que actualmente están destinadas a la agricultura y los pastos limpios requiere de la participación activa de la población en programas de reforestación, cuyos resultados sean protegidos y respaldados por el gobierno distrital. LOgrar este objetivo acercaría a la ciudad hacia el objetivo 15 “Gestionar sosteniblemente los bosques, luchar contra la desertificación, detener e invertir la degradación de las tierras y detener la pérdida de la biodiversidad”. Es este último objetivo el que más se ve obstaculizado por el gobierno distrital, no solo por su constante insistencia de urbanizar la reserva, sino por el entorpecimiento de las actividades elaboradas por diferentes actores locales que buscan la protecciòn de la reserva.
Los resultados presentados en este informe son de carácter académico y aún requieren de otras consideraciones para aproximarse a la realidad de la reserva TVDH. El uso de otros algoritmos de clasificación como el Random Forest o Redes Neuronales, complementados con herramientas más sólidas de series de tiempo y el uso de la reflectancia del fondo de la atmòsfera (BOA), pueden representar un mejor escenario de la reserva para delinear las estrategias que debe seguir la ciudad para la protección de este espacio natural, estratégico e importante para el ambiente de la ciudad de Bogotá y la calidad de vida de sus habitantes. Protejan la Reserva Thomas van der Hammen
Bibliografía
Fundación Humedales de Bogotá (2016). Reserva Thomás Van der Hammen, hogar de aves y naturaleza que un alcalde llama potrero. Tomado de: ´` = Problemática van der Hammen
Pérez P. A. (2000). La estructura ecológica principal de la sabana de Bogotá. Sociedad Geográfica de Colombia. Academia de Ciencias Geográficas. 1 – 37.
Reserva Forestal del Norte de Bogotá (2017). Descripción de la Reserva. Tomado de: Reserva TVDH
Torres N., Angélica & Niño M. Freddy (2018). Los servicios ecosistémicos provistos por la reserva. Análisis de los tres escenarios. En: Bogotá y la reserva Thomas van der Hammen. Aportes desde la Economía Ecológica (Villamil et al. 2018). Instituto de Estudios Ambientales y Universidad Nacional de Colombia. 80 pp.