PRÁCTICAS DE ANÁLISIS DE IMÁGENES DE SENSORES REMOTOS CON R
Diego Alejandro Orjuela Hurtado
Abril 22 de 2020
1. Introducción
El presente documento ha sido creado para el curso de Percepción remota, curso de la maestría en Geomática de la Universidad Nacional de colombia, sede Bogotá. Contempla un análisis básico de imagenes satelitales Landsat8 para el municipio de Cáqueza, en el departamento de Cundinamarca, Colombia. la métodología contemplada se fundamenta en los ejercicios “Remote Sensing Image Analysis” creado por el Consorcio de Investigación de Sistemas Geoespaciales y Agropecuarios, o por sus siglas en ingles GFC (Aniruddha y Hijmans, 2016).
2. Datos
2.1. Librerias R
Para la ejecución del presente análisis, nos soportamos en el lenguaje de programación R, apoyado en un conjunto de librerías especializadas para tratamiento y despliegue de datos geoespaciales, serán mencionadas a continuación.
2.1.1. Librería SP:
library (sp)Librería que permite definir para el lenguaje las clases y métodos especificos para datos espaciales, necesaria para reconocer los datos donde resida información geoespacial (Pebesma, 2020), algunas de las funciones de está librería van desde trazado de datos como mapas, selección espacial, métodos para recuperar coordenadas, ploteo, resumen de datos, entre otros.
2.1.2. Librería raster:
library (raster)Responsable de habilitar el análisis y modelamiento para datos geográficos, encargada de implementar funciones de lectura, escritura, manipulación y modelado de datos raster, y operaciones básicas para datos vectoriales (Hijmans, Perpinan, & Mattiuzzi, 2020).
2.1.3. Librería rgdal:
library (rgdal)## rgdal: version: 1.4-8, (SVN revision 845)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
## Path to GDAL shared files: E:/Documentos/R/win-library/3.6/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: E:/Documentos/R/win-library/3.6/rgdal/proj
## Linking to sp version: 1.4-1Proporiona accero a la Librería de abstracción de datos ‘geoespaciales’ (GDAL), dando posibilidad la importación de datos ráster y datos vectoriales, facilita funciones proyección / transformación de datos geoespaciales (Bivand, 2019).
2.1.4. Librería rpart:
library(rpart)Librería que da acceso al método “Particionamiento recursivo para árboles de clasificación”, corresponde a una herramienta creada a partir de funcion propuesta por Breiman, Friedman, Olshen y Stone en 1984 (Atkinson, 2019).
2.1.5. Librería rasterVis:
Métodos avanzados para visualización e interacción con datos ráster, encargado de ejecutar métodos de visualización para datos cuantitativos y categórias (Perpinan, 2019).
library(rasterVis)## Loading required package: lattice## Loading required package: latticeExtra2.1.6. Librería dismo:
library(dismo)Libreria encargada de proporcionar las funciones para el modelado y predicción de distribuciones geográficas.forma de ocurrencia y número de eventos detectados en un localización geográfica (Hijmans, 2017).
2.2. Imagenes satelitales multiespectrales
Los datos a procesar, se han recuperado de la base de datos de imagenes del satelite Landsat 8, que abarca una région de la coordillera oriental de los Andes, más especificamente corresponden a la zona WRS-path: 8 y wrs-row: 57.
# Blue
b2 <- raster('data/LC08_L1TP_008057_20150104_20170415_01_T1_B2.tif')
# Green
b3 <- raster('data/LC08_L1TP_008057_20150104_20170415_01_T1_B3.tif')
# Red
b4 <- raster('data/LC08_L1TP_008057_20150104_20170415_01_T1_B4.tif')
# NIR
b5 <- raster('data/LC08_L1TP_008057_20150104_20170415_01_T1_B5.tif')
landsatRGBOri1 <- stack(b4, b3, b2)
plotRGB(landsatRGBOri1, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 - Composición en color verdadero")
landsatRGBOri2 <- stack(b5, b4, b3)
plotRGB(landsatRGBOri2, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 - Composición en falso color")
#Creación de la lista de Raster corregido por valores TOA
filenames <- paste0('data/caqueza/caq_TOA_b', 1:7, ".tif")
#Selección de nombres para las bandas
landsat <- stack(filenames)
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')El conjunto de datos seleccionados, ha sido recuperado por medio de la plataforma virtual Earth Explorer (USGS, 2020). y cuentan con las siguientes caracteristicas:
Nivel de categoría: T1
Porcentaje de cobertura por nubosidad: 5.11%
Calidad de la imagen: OLI=9 y TIRS=9
Paramétros de proyección (Sistema coordenado): UTM 18 - WGS 84
crs(b2)## CRS arguments:
## +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0-Extensión de la imagen:
extent(b2)## class : Extent
## xmin : 446085
## xmax : 673515
## ymin : 362985
## ymax : 595215-Resolución espacial para las bandas OLI: 30[m] x 30[m]
res(b2)## [1] 30 302.3. Clasificación de usos del suelo
Para la clasificación de usos del suelo, nos apoyaremos en 2 clasificaciones: la primera es la clasificación recuperada del “Mapas de Suelos del Territorio Colombiano a escala 1:100.000” que ofrece el IGAC (IGAC, 2012). la segunda clasificación en basada en el “National Land Cover Database” es un producto de cobertura terrestre para los EE. UU. una cobertura del suelo basada en Landsat de 30 m.
2.3.1. Datos de referencia IGAC
Son los datos oficiales de clasificación para uso de suelo utilizados en colombia, aun que el abanico de clasificaciones es muy amplio, solo utilizaremos los usos de suelo correspondientes a nuestra zona de estudio: “Zonas industriales o comerciales”, “Mosaico de pastos y cultivos”, “Bosque fragmentado”, “Tejido urbano continuo”, “Pastos limpios”, “Pastos enmalezados”, “Mosaico de cultivos, pastos y espacios naturales”, “Mosaico de pastos con espacios naturales”, “Vegetacion secundaria o en transicion”, “Herbazal”, “Zonas de extraccion minera” y “Rios (50 m)”.
nlcd1 <- brick('data/caqueza/landuse2.tif')
names(nlcd1) <- c('nlcdcaq1')
# The class names and colors for plotting
nlcdclass1 <- c( "Zonas industriales o comerciales", "Mosaico de pastos y cultivos", "Bosque fragmentado", "Tejido urbano continuo", "Pastos limpios", "Pastos enmalezados", "Mosaico de cultivos, pastos y espacios naturales", "Mosaico de pastos con espacios naturales", "Vegetacion secundaria o en transicion", "Herbazal", "Zonas de extraccion minera", "Rios (50 m)")
classdf1 <- data.frame(classvalue1 = c(1,2,3,4,5,6,7,8,9,10,11,12), classnames1 = nlcdclass1)
# Hex codes of colors
classcolor1 <- c( "#B50000", "#FBF65D", "#38814E", "#BFBFBF", "#92D050", "#AF963C", "#548235", "#C6E0B4", "#D2CDC0", "#D1D182", "#BF8F00", "#5475A8")
# Now we ratify (RAT = "Raster Attribute Table") the ncld2011 (define RasterLayer as a categorical variable). This is helpful for plotting.
nlcdcaq1 <- nlcd1[[1]]
nlcdcaq1 <- ratify(nlcdcaq1)
rat1 <- levels(nlcdcaq1)[[1]]
rat1$landcover <- nlcdclass1
levels(nlcdcaq1) <- rat1plt11 <- levelplot(nlcdcaq1, col.regions = classcolor1, main = 'Clasificación usos del suelo - IGAC')
print(plt11)
2.3.2. Datos de referencia NLCD
Para la clasificación NLCD, nos interesa analizarlo para este ejercicio dada la mayor simplicidad de en su clasificación y por ende reconocimiento. para construir el conjunto de datos clasificacdos partimos de los datos originalas de IGAC y reclasificamos manualmente para hacer coincidir con la NLCD. los subconjuntos encontrados para esta aplicación son: “Urbanización” (corresponde a los valores de “Devoloped”), “Plantado / Cultivado”, “Bosque”, “Pastos”, “Matorral”, “Herbáceo”, “Estéril”, “Agua”.
nlcd <- brick('data/caqueza/landuse4.tif')
names(nlcd) <- c('nlcdcaq')
nlcdclass <- c("Urbanización", "Plantado / Cultivado", "Bosque", "Pastos", "Matorral", "Herbáceo", "Esteril", "Agua")
classdf <- data.frame(classvalue = c(1,2,3,4,5,6,7,8), classnames = nlcdclass)
classcolor <- c( "#B50000", "#FBF65D", "#38814E", "#C8E6F8", "#AF963C", "#D1D182", "#D2CDC0", "#5475A8")
nlcdcaq <- nlcd[[1]]
nlcdcaq <- ratify(nlcdcaq)
rat <- levels(nlcdcaq)[[1]]
rat$landcover <- nlcdclass
levels(nlcdcaq) <- ratplt12 <- levelplot(nlcdcaq, col.regions = classcolor, main = 'Clasificación usos del suelo - NLCD')
print(plt12)
2.4. Zona de estudio
Para efectos de este documento, se ha seleccionado el municipio de Cáqueza, en Cundinamarca como zona de estudio. Se localiza su centro urbano en coordenadas geográficas 4°24′19″ Norte - 73°56′52″ Oeste. Tiene una población aproximada de 7300 habitantes y una extensión de aproximadamente 11100Ha; segun el documento de esquema de ordenamiento territorial de cual se dispone (año 2000-2009), Su área urbana ocupa aproximadamente el 0.71% del área total del municipio, el resto de la extensión es dedicada para usos de Protección ambiental (18.6%) y Suelos de uso agropecuario (67.91%), un 12.78% de los suelos no han sido clasificados (Concejo Municipal de Cáqueza, 2000).
landsatRGB <- landsat[[c(4,3,2)]]
plotRGB(landsatRGB, axes = TRUE, stretch = "lin", main = "Imagen mun. Caqueza - Composición en color verdadero")
landsatFCC <- landsat[[c(5,4,3)]]
plotRGB(landsatFCC, axes=TRUE, stretch="lin", main = "Imagen mun. Caqueza - Composición en falso color")
La extención en coordenadas proyectadas UTM 18 que contiene la zona rural y zona urbana del municipio es presentado a continuación:
extent(landsat)## class : Extent
## xmin : 607905
## xmax : 631335
## ymin : 478425
## ymax : 4949252.5. Perfiles espectrales
La gráfica del valores espectrales de los píxeles, que representan características de la superficie terrestre; podemos determinar mayor inflencia de cuerpos de agua, especies vegetales, edificaciones e incluso la ausencia de estos elementos a raiz de dichas respuestas. Los valores espectrales se extraen del conjunto de datos de las imagenes satelitales y puedes compararse al los usos de suelos (pueden extraer de cualquier conjunto de datos multiespectrales utilizando la función de extracción. En el ejemplo anterior, extrajimos valores de datos de Landsat para la obtención de muestras que reflejen los valores de respuesta espectral para cada una de las coberturas (Aniruddha & Hijmans, 2016).
landuse <- readOGR('data/caqueza/landuse4.shp')## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\usuario\R\data\caqueza\landuse4.shp", layer: "landuse4"
## with 97 features
## It has 13 fields# genera 5000 puntos de muestras del shape de usos de suelo
ptlanduse <- spsample(landuse, 5000, type='regular')
ptlanduse$class <- over(ptlanduse, landuse)$LEYENDA3N
# extrae los valores de reflectacia en las bandas landsat utilizadas
df <- extract(landsat, ptlanduse)
ms <- aggregate(df, list(ptlanduse$class), mean)
rownames(ms) <- ms[,1]
ms <- ms[,-1]
# Asignación de caracteristicas de ploteo
mycolor <- c("#5475A8", "#38814E", "#D1D182", "#AF963C", "#D2CDC0", "#C8E6F8", "#FBF65D", "#B50000")
ms <- as.matrix(ms)
plot(0, ylim=c(0,0.6), xlim = c(1,7), type='n', xlab="Bands", ylab = "Reflectance")
for (i in 1:nrow(ms)){ lines(ms[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])}
title(main="Spectral Profile from Landsat", font.main = 2)
legend("topleft", rownames(ms),
cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")
3.Métodos
Se analizaran los datos de 3 formas diferentes los datos de las imagenes: por medio de los indices espectrales normalizados, un método de clasificación no supervisada (k-means) y un método supervisado (Particionamiento recursivo).
3.1. Índices normalizados derivados de imágenes satelitales
vi <- function(img, k, i) {#"k" e "i" corresponden a las bandas a utilizar dependiendo el indice ND que se desee utilizar
bk <- img[[k]]
bi <- img[[i]]
vi <- (bk - bi) / (bk + bi)
return(vi)
}3.1.1. Normalized Difference Vegetation Index (NDVI)
Es un valor residual de las respuestas espectrales entre el inflarojo cercano y la banda del color rojo visible, se apoya en el conocimiento respecto a las respuesta espectrales, la vegetación saludable absorbe el espectro electromagnético en el rango del rojo (banda 4) y el azul (banda 2) visible. mas sin embargo tambien tiene una alta reflectancia en infrarrojo cercano (NIR, banda 5). Es aprovechando estas diferencias de la reflectancia y absorción en espectros, a raiz de esto es que podemos determinar la localización de vegetación basandonos en una operación conocida como NDVI:
NDVI= (Banda 5 - Banda 4)/(Banda 5 + Banda 4)Los resultados de la operación podemos visualizarlos acontinuación.
ndvi <- vi(landsat, 5, 4)3.1.2. Normalized Difference Water Index (NDWI)
Los cuerpos de agua dependiendo su estado pueden reflejar distintos componentes de espectro electromagnetico visible, más sin embargi el agua presenta una gran reflejo en rangos de espectro posteriores al visible (más allá del infrarrojo cercano NIR). Este índice utiliza las bandas de infrarrojo cercano (NIR) y de infrarrojo de onda corta (SWIR) para intentar localizar el porcentaje de agua en terreno (Ceccato, Gobron, Flasse, Pinty, & Tarantola, 2002), cabe aclarar que es un indice muy sensible a la acumulación de tierra y da como resultado cuerpos de agua sobreestimados, pero que en conjunto con indices tipo NDVI sirven para para evaluar el contexto de las áreas en cuestión.
NDWI = (NIR - SWIR) / (NIR + SWIR)Para efectos de calculo con Landsat 8 corresponden a las bandas Banda 5 y Banda 6 respectivamente.
ndwi <- vi(landsat, 5, 6)3.1.3. Modified Normalized Difference Water Index (MNDBI)
Es un cambio planteado a las fórmula de NDWI, modificada por Xu (2006), plantea la utilización de las bandas verde y SWIR, es un indice que al operarse por el espectro visible verde deja de dar valores superiores en vegetación; con lo cual se puede proceder a localizar cuerpos de agua de manera más apropiada.
MNDWI = (Verde - SWIR) / (Verde + SWIR)Los cuerpos de agua poseen valores mayores a 0.5, la vegetación presenta valores considerablemente más pequeños que se encuentran en el rango de 0 a 0.2.
mndwi <- vi(landsat, 3, 6)3.1.4. Normalized Difference Built-up Index (NDBI)
Hace parte de uno de los indices existentes para el reconocimiento de zonas urbanizadas, este índices se basa en que los valores de zonas de acopio y el suelo desnudo reflejan más el espectro SWIR, y reflejan en menor proporción el NIR, es este patrón observado en urbes y lugares estériles, al cual puede interpolarse para áreas con edificaciones (Longley, Goodchild, Maguire, & Rhind, 2005).
MNDWI = (SWIR - NIR) / (SWIR + NIR)Los valores negativos de NDBI representa cuerpos de agua, el valor del NDBI para la vegetación es bajo y un valor más alto para zonas de urbanización.
ndbi <- vi(landsat, 6, 5)3.2. Clasificación No supervisada
El proceso de clasificación no supervisada que usaremos es el método de K-Means, es un algoritmo que de agrupación, divide las observaciones en un numero determinado de agrupaciones (k agrupaciones), donde dividimos los datos en grupos que pueden ser iguales o superiores al número de clases; el algoritmo tiene como objetivo la clasificacion de subgrupos homogéneos, basandose en las similitudes de los objetos entre sí.
Primeramente delimitaremos la extensión de la imagen, utilizamos una localización cerca al Rio Negro y en presencia de una zona urbana.
#Se restringe la zona de los raster
e <- extent(621000.0, 625899.3, 478440.3, 482000.0)
e2 <- extent(615000.0, 621700.0, 483000.0, 488500.0)
ndvi_e <- crop(ndvi, e)
ndwi_e <- crop(ndwi, e)
mndwi_e <- crop(mndwi, e)
ndbi_e <- crop(ndbi, e)
ndbi_e2 <- crop(ndbi, e2)
RGB_e <- crop(landsatRGB, e)
RGB_e2 <- crop(landsatRGB, e2)
#Conversión de raster a matriz
nrV <- getValues(ndvi_e)
nrW <- getValues(mndwi_e)
nrB <- getValues(ndbi_e)Se recuperan los valores de las celdas y se transforman en matrices para permitir el analisis del algoritmo; los valores recuperados son producto de tomas de muestras aleatorias (en este caso usamos 500 iteraciones).
set.seed(99)
kmnclusterV <- kmeans(na.omit(nrV), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")
set.seed(99)
kmnclusterW <- kmeans(na.omit(nrW), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")
set.seed(99)
kmnclusterB <- kmeans(na.omit(nrB), centers = 10, iter.max = 500, nstart = 5, algorithm="Lloyd")El algoritmo regresará un conjunto de datos tipo Raster Layer, que contemplaran un conjunto de 10 grupos con valores de 1 a 10, cada grupo es seleccionado en función a la similitud de los valores e interacciones encontradas, si bien no se tienen datos de a que clase de cobertura corresponde cada valor, permite identificar conjuntos de datos relacionables por las respuestas espectrales del las imagenes satelitales.
3.3. Clasificación supervisada
Es un método qué para una muestra de clasificación se tienen n observaciones sobre una variable de clase Y que toma valores de 1 hasta k y con p variables predictoras X1 hasta Xp. El objetivo es encontrar un modelo para predecir los valores de Y a partir de nuevos valores de X. La solución es una partición de X en k conjuntos disyuntos (A1, A2, … Ak) de modo que el valor predicho de Y es j si X pertenece a Aj, donde j= 1, 2, …, k. Si las variables X toman valores ordenados, existen dos soluciones clásicas: una es análisis discriminante lineal y clasificación vecina más cercana, donde se producen conjuntos Aj lineales y no lineales, que n son fáciles de interpretar si p es grande.
Los métodos del árbol de clasificación tienen la posibilidad de producir conjuntos rectangulares Aj al dividir el conjunto de datos una variable X a la vez. En consecuencia, esto hace que los conjuntos sean más fáciles de interpretar. Si se tiene una gráfica de puntos donde se denotan las particiones, es más difícil de interpretar que un árbol de decisión. La primera está limitada a dos variables, mientras que la segunda no tienen límite de número de variables.
3.3.1. Generación de sitios de muestra
Los datos de validacion que vamos a utilizar corresponden a los dos sistemas de clasificación de usos del suelo, es a partir de este conjunto de datos que los algoritmos direccionan su aprendizaje.
# Generador que produce resultados aleatorios
set.seed(99)
# creación de las muestras
samp <- sampleStratified(nlcdcaq, size = 116, na.rm = TRUE, sp = TRUE)
samp1 <- sampleStratified(nlcdcaq1, size = 116, na.rm = TRUE, sp = TRUE)El aprendizaje se da por la generación aleatoria de puntos, estos son seleccionados aleatoriamente entre los poligonos de cada conjunto (para cada clasificacion de uso del suelo), el número de puntos de muestra son determinados por el usuario, entre mayor sea el numero de variables utilizadas (mayor sea el numero de clasificaciones para el uso del suelo en este caso) es mas improbable es que todos los componentes tenga peso en la caracterización de los datos.
plt1 <- levelplot(nlcdcaq1, col.regions = classcolor, main = 'Distribucion muestras de entrenamiento - Clasificación IGAC')
print(plt1 + layer(sp.points(samp1, pch = 3, cex = 0, col = 1)))
plt <- levelplot(nlcdcaq, col.regions = classcolor, main = 'Distribucion muestras de entrenamiento - Clasificación NLCD')
print(plt + layer(sp.points(samp, pch = 3, cex = 0, col = 1)))
3.3.2. Extracción de valores de los sitios de muestra
Una vez obtenidas las distribuciones de los puntos para los valores de cada conjunto (cobertura de tierra), se extraen los valores de respuesta espectral para todas las bandas relacionadas con ese punto georeferenciado.Estos valores serán las variables que predigan los “valores” para cada clase.
filenam <- paste0('data/caqueza/caq_TOA_b', 2:7, ".tif")
landsatAux <- stack(filenam)
names(landsatAux) <- c('blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
sampvals <- extract(landsatAux, samp, df = TRUE)
sampvals <- sampvals[, -1]
sampvals1 <- extract(landsatAux, samp1, df = TRUE)
sampvals1 <- sampvals1[, -1]
sampdata <- data.frame(classvalue = samp@data$nlcdcaq, sampvals)
sampdata1 <- data.frame(classvalue = samp1@data$nlcdcaq1, sampvals1)3.3.3. Entrenamiento del Clasificador
Por medio del algoritmo de “particionamiento recursivo para arboles de clasificación” incluido en la librería rpart de R, podemos ejecutar la predicción de las clases de uso de suelo, en función a los valores recolectados durante la generación de muestras para entrenamiento.
# Train the model
cart1 <- rpart(as.factor(classvalue)~., data=sampdata1, method = 'class', minsplit = 5)
# Plotear el arbol de clasificación
plot(cart1, uniform=TRUE, main="Arbol de Clasificación - IGAC")
text(cart1, cex = 0.8)
# Train the model
cart <- rpart(as.factor(classvalue)~., data=sampdata, method = 'class', minsplit = 5)
# Plotear el arbol de clasificación
plot(cart, uniform=TRUE, main="Arbol de Clasificación - NLCD")
text(cart, cex = 0.8)
En el diagrama de árbol de clasificación, los valores de clase se imprimen en los nodos, y cada numeración corresponde a un tipo de cobertura.
4. Resultados
4.1. Índices normalizados derivados de imágenes satelitales
4.1.1. Normalized Difference Vegetation Index (NDVI)
En coincidencia con las caracteristicas de la clasificación del uso del suelo, los valores caracterizados del NVDI indican una alta presencia de vegetación, los valores picos de distribución de valores se enfocan alrededor del 0.7 (cobertura vegetal densa) y en menor medida para valores que indicen zonas de escasa covertura vegetal.
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat-NDVI")
# view histogram of data
hist(ndvi,
main = "Distribucion de valores NDVI",
xlab = "NDVI",
ylab= "Frequency",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
4.1.2. Normalized Difference Water Index (NDWI)
Se puede constatar por medio de las gráficas y el raster, que los valores elevados de NDWI no sobrepasan el valor de 0.5, evidenciandose la poca captación de cuerpos de agua, pero si reflejandose en la altas tasas de humedad de la cobertura que se encuentran en valores entre 0 y 0.45.
plot(ndwi, col = rev(terrain.colors(10)), main = "Landsat-NDWI")
# view histogram of data
hist(ndwi,
main = "Distribucion de valores NDWI",
xlab = "NDWI",
ylab= "Frequency",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
4.1.3. Modified Normalized Difference Water Index (MNDWI)
Los valores negativos, se ven explicados por la ausencia de valores correspondientes a cuerpos de agua (similar al indice MDWI), los números negativos se dan por la diferencia entre la poca reflectancia del espectro visible verde versus los altos niveles en el NIR, sinonimo de una ya comprobada presencia de coberturas vegetal sana.
# For Landsat NIR = 5, swir = 6.
plot(mndwi, col = rev(terrain.colors(10)), main = "Landsat-MNDWI")
# view histogram of data
hist(mndwi,
main = "Distribucion de valores MNDWI",
xlab = "MNDWI",
ylab= "Frequency",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
4.1.4. Normalized Difference Built Index (NDBI)
Los escasos valores de respuesta positivos, son el reflejo de una ausencia de coberturas destapadas/urbanizadas en los datos análizados, coincidiendo con el bajo porcentaje del área urbana descrita para el municipio de Cáqueza.
plot(ndbi, col = rev(terrain.colors(10)), main = "Landsat-NDBI")
# view histogram of data
hist(ndbi,
main = "Distribucion de valores NDBI",
xlab = "NDBI",
ylab= "Frequency",
col = "wheat",
xlim = c(-0.5, 1),
breaks = 30,
xaxt = 'n')
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))
4.1.5. THRESHOLDING
Otra forma como podemos visualizar las estimaciones previamente obtenidas, es utilizando los valores indice para reconocer elementos caracteristicos en superficie, apoyandonos en los valores los indices previamente obtenidos, la bibliografia aporta rangos de valores para los indices, que nos pueden indicar la presencia de determinadas coberturas (Ceccato, Gobron, Flasse, Pinty, & Tarantola, 2002), (Xu, 2006), (Longley, Goodchild, Maguire, & Rhind, 2005).
4.1.5.1. Detección por medio de la NDVI
- NDVI = {-1 a 0} Representa cuerpo de agua
AguaNDVI <- reclassify(ndvi_e, cbind(0, Inf, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(AguaNDVI, axes = TRUE, stretch = "lin", main= 'Cuerpos de Agua - NDVI')## Warning in plot.window(...): "stretch" is not a graphical parameter## Warning in plot.xy(xy, type, ...): "stretch" is not a graphical parameter## Warning in axis(side = side, at = at, labels = labels, ...): "stretch" is not a
## graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "stretch" is not a
## graphical parameter## Warning in box(...): "stretch" is not a graphical parameter## Warning in title(...): "stretch" is not a graphical parameter
- NDVI = {-0.1 a 0.1} Representa afloramientos rocosos o arena
ArenaNDVI <- reclassify(ndvi_e, c(-Inf, -0.1, NA, -0.1, 0.1, 1, 0.1, Inf, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(ArenaNDVI, main='Afloramientos rocosos y arena - NDVI')
- NDVI = {0.2 a 0.5} Representa arbustos y/o pastizales
ArbNDVI <- reclassify(ndvi_e, c(-Inf, 0.2, NA, 0.2, 0.5, 1, 0.5, Inf, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(ArbNDVI, main='Arbustos y/o pastizales - NDVI')
- NDVI = {0.4 a 1.0} Representa vegetación densa.
vegNDVI <- reclassify(ndvi_e, cbind(-Inf, 0.4, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(vegNDVI, main='Vegetation densa - NDVI')
4.1.5.2. Detección por medio de la NDWI y MNDWI
- NDWI = {0 a 1} Representa humedad en la cobertura existente
AguaNDwI <- reclassify(ndwi_e, cbind(-Inf, 0, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(AguaNDwI, main='Cuerpos de Agua - NDWI')
AguaMNDwI <- reclassify(mndwi_e, c(-Inf, 0, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(AguaMNDwI, main='Cuerpo de agua - MNDWI')
Podemos observar como la sombra producto de la nubosidad ha causado una captación erronea de los valores de reflectancia, que se ven reflejados en el indice, si bien la toma de los valores ha ayudado a detectar el cuerpo de agua, son considerablemente pocos los pixeles con el valor indicado como para poder reconstruir el cause entero del rio.
4.1.5.3. Detección por medio de la NDBI
- NDBI = {0 a 1} Representa Urbanización y lugares estériles
Para la presente detección ha decidido utilizarse dos escenarios distintos, en ambos se observa en las imagenes landsat, la presencia de zona urbana.
UrbNDBI1 <- reclassify(ndbi_e, cbind(-Inf, 0, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(UrbNDBI1, main='Urbanización - NDBI')
UrbNDBI2 <- reclassify(ndbi_e, cbind(-Inf, 0, NA))
par(mfrow = c(1,2))
plotRGB(RGB_e2, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(UrbNDBI2, main='Urbanización - NDBI')
4.2. Clasificación No supervisada
El proceso de análisis de la clasificación supervisada se hizo en un segmento cercano a una urbanización y al cauce principal del rio negro; Los valores resultantes son producto de una clasificación utilizando los valores del NDVI, MNDWI y NDBI los resultados son presentados acontinuación:
# Use the ndvi object to set the cluster values to a new raster
knrV <- setValues(ndvi_e, kmnclusterV$cluster)
knrW <- setValues(mndwi_e, kmnclusterW$cluster)
knrB <- setValues(ndbi_e, kmnclusterB$cluster)4.2.1. Clasificación No supervisada sobre banda NVDI
# Define a color vector for 10 clusters (learn more about setting the color later)
mycolor <- c( "#B50000", "#FBF65D", "#38814E", "#BFBFBF", "#92D050", "#AF963C", "#548235", "#C6E0B4", "#D2CDC0", "#D1D182", "#BF8F00", "#5475A8")
plot(ndvi_e, col = rev(terrain.colors(10)), main = 'NDVI')
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(knrV, main = 'Clasificación no supervisada - NDVI', col = mycolor )
compndvi <- stack(RGB_e, knrV)
set.seed(1)
sr <- sampleRandom(compndvi, 10000)
plot(sr[,c(1,2)], main = "Landsat8 vs Clasificación no supervisada NDVI")
4.2.2. Clasificación No supervisada sobre banda MNWDI
# Define a color vector for 10 clusters (learn more about setting the color later)
plot(mndwi_e, col = rev(terrain.colors(10)), main = 'MNDWI')
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(knrW, main = 'Clasificación no supervisada - MNDWI', col = terrain.colors(10) )
compmndwi <- stack(RGB_e, knrW)
set.seed(1)
sr <- sampleRandom(compmndwi, 10000)
plot(sr[,c(1,2)], main = "Landsat8 vs Clasificación no supervisada MNDWI")
4.2.3. Clasificación No supervisada sobre banda NVBI
# Define a color vector for 10 clusters (learn more about setting the color later)
plot(ndbi_e, col = rev(terrain.colors(10)), main = 'NDBI')
par(mfrow = c(1,2))
plotRGB(RGB_e, axes = TRUE, stretch = "lin", main = "Imagen Landsat8 RGB")
plot(knrB, main = 'Clasificación no supervisada - NDbI', col = terrain.colors(10) )
compmndwi <- stack(RGB_e, knrW)
set.seed(1)
sr <- sampleRandom(compmndwi, 10000)
plot(sr[,c(1,2)], main = "Landsat8 vs Clasificación no supervisada NDBI")
4.3. Clasificación supervisada
4.3.1. Clasificación
Tras el trenamiento del modelo, ahora es apto para hacer predicciones, por consiguiente ahora se procedera a clasificar todas las celdas en la imagen digital de Landsat 8.
# Se efectua la predicciòn
pr <- predict(landsatAux, cart, type='class')Los resultados de la clasificación resultante son presentados acontinuación:
classdfAux <- data.frame(classvalue1 = c(1,2,3,6,7,8), classnames1 = c("Urbanizacion", "Plantado / Cultivado", "Bosque", "Herbáceo", "Mineria", "Agua"))
classcolorAux <- c( "#B50000", "#FBF65D", "#38814E", "#D1D182", "#D2CDC0", "#5475A8")
pr <- ratify(pr)
rat <- levels(pr)[[1]]
rat$legend <- classdfAux$classnames
levels(pr) <- rat
levelplot(pr, maxpixels = 1e6, col.regions = classcolorAux, scales=list(draw=FALSE), main = "Arbol de clasificación supervisada")
4.3.2. Evaluación del modelo
Para análizar la precisión del modelo, utilizaremos dos medidas de precisión: precisión general y kappa. Puede realizar la evaluación de precisión utilizando las muestras independientes.
Para evaluar cualquier modelo, puede se utilizada la validación cruzada (k-fold), donde los datos utilizados para ajustar el modelo se dividen en k grupos (5 para este ejemplo), uno de los grupos se usará para la prueba del modelo y los otros datos se usarán para el entrenamiento del.
set.seed(99)
j <- kfold(sampdata, k = 5, by=sampdata$classvalue)Este modo de validación permite entrenar y probar el modelo k veces (5 veces).
x <- 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')
# create a data.frame using the reference and prediction
x[[k]] <- cbind(test$classvalue, as.integer(pclass))
}La unión de las 5 validaciones las juntamos en una unica matriz que permite tenerminar el número de eventos previstos vs observados.
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) <- c("Urbanizacion", "Plantado / Cultivado", "Bosque", "Pastos", "Matorral", "Herbáceo", "Esteril", "Agua")
rownames(conmat) <- classdf$classnames
conmat## predicted
## observed Urbanizacion Plantado / Cultivado Bosque Pastos Matorral
## Urbanización 33 22 8 1 11
## Plantado / Cultivado 20 21 29 4 13
## Bosque 4 7 83 0 3
## Pastos 14 16 26 1 6
## Matorral 21 19 30 2 14
## Herbáceo 9 12 23 2 8
## Esteril 2 0 1 0 0
## Agua 3 2 9 0 3
## predicted
## observed Herbáceo Esteril Agua
## Urbanización 13 14 14
## Plantado / Cultivado 18 2 9
## Bosque 3 2 14
## Pastos 22 4 27
## Matorral 14 2 14
## Herbáceo 23 2 37
## Esteril 14 86 13
## Agua 22 6 71Calculamos la precisión general
# numero de casos
n <- sum(conmat)
diag <- diag(conmat)
# Overall Accuracy
OA <- sum(diag) / nCalculamos el valor Kappa
# observaciòn de los casos (verdaderos) por clase
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)Evaluación de la precisión producida vs precisión del usuario
# Producer accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc## producerAccuracy userAccuracy
## Urbanizacion 0.3113208 0.28448276
## Plantado / Cultivado 0.2121212 0.18103448
## Bosque 0.3971292 0.71551724
## Pastos 0.1000000 0.00862069
## Matorral 0.2413793 0.12068966
## Herbáceo 0.1782946 0.19827586
## Esteril 0.7288136 0.74137931
## Agua 0.3567839 0.612068975. Conclusiones
El comportamiento general de los índices normalizados de imágenes satelitales, demuestra que si bien tiene una importante utilidad en el reconocimiento de componentes de terreno, se ven afectados por ruido procedente de otros componentes del terreno que llegan a tener respuestas similares a los datos deseados. Considerando prudente la intercepción de datos resultante de varios indices para mejorar el resultado obtenido durante este procedimiento.
El NDVI es el que mejor respuesta produce, dada la importante extensión de vegetación densa, reflejada tanto en la clasificación de cobertura del suelo, como en los ejercicios de clasificación supervisada y no supervisada. Caso contrario la capacidad de detección de las zonas urbanas y cuerpos de agua, se vieron importantemente límitados debido a su poca extensión; por más que se intento reducir el área de análisis, las respuestas espectrales (reflejados en los índices NDWI, MNDWI y NDBI) eran bajas, esto límito y altero considerablemente la capacidad de los métodos de clasificación.
El método de clasificación supervisada, fue el que tuvo mayor dificultad en el aporte a una previsión adecuada, partiendo de las diferencias temporales entre las imagenes satelitales y la clasificación de cobertura del suelo dispobible (2017 y 2012 respectivamente), e incluyendo la incertidumble de los poligonos; visible por ejemplo en la ronda de los cuerpos de agua.causaron que los valores de reflectancia hallados para algunas clases del conjunto de “landuse” o no tuvisen datos suficientes para ser distinguidos de las demas coberturas, o simplemente no se reflejaban en los poligonos de la clasificación de usos del suelo.
-Consecuencia de esto ultimo, se envidencia en los arboles de clasificación y regresión, se ven seriamente afectados por el número disponible de “samples” (muestras de calibración), a mayor sea el rango de muestras más influencia se verá reflejada en los arboles de decisión. Razón por la cual, habian elementos de la clasificación de suelos (para ambas clasificaciones) que al contar con una población poco significativa, estos dejaban de tener un peso relevante en las predicciones posteriores. viendose la situación de que el modelo predictivo no reflejaba todas las categorias de clasificación y por defecto, si se disminuian el número del muestreo, comenzaban a darse clasificaciones erroneas; visibles en el proceso de análisis que llevamos acabo en este documento, dado que solo se tomaron 116 muestras de entrenamiento, y observandose para coberturas de gran extensión como pastos y bosques, valores de confianza inferiores al 40%.
6. Referencias
Aniruddha, G., & Hijmans, R. (2016). Remote Sensing Image Analysis — R Spatial. Retrieved April 22, 2020, from https://rspatial.org/raster/rs/index.html
Atkinson, B. (2019). rpart package | R Documentation. Retrieved April 22, 2020, from https://www.rdocumentation.org/packages/rpart/versions/4.1-15
Bivand, R. (2019). rgdal package | R Documentation. Retrieved April 22, 2020, from https://www.rdocumentation.org/packages/rgdal/versions/1.4-8
Ceccato, P., Gobron, N., Flasse, S., Pinty, B., & Tarantola, S. (2002). Designing a spectral index to estimate vegetation water content from remote sensing data: Part 2. Validation and applications. Remote Sensing of Environment, 82(2–3), 198–207. https://doi.org/10.1016/S0034-4257(02)00036-6
Concejo Municipal de Cáqueza. (2000). Esquema de ordenamiento territorial Cáqueza-Cundinamarca (2000 - 2009). Retrieved from https://caquezacundinamarca.micolombiadigital.gov.co/sites/caquezacundinamarca/content/files/000051/2526_acuerdo-0062000-eot.pdf
Hijmans, R. (2017). dismo package | R Documentation. Retrieved April 22, 2020, from https://www.rdocumentation.org/packages/dismo/versions/1.1-4
Hijmans, R., Perpinan, O., & Mattiuzzi, M. (2020). raster package | R Documentation. Retrieved April 22, 2020, from https://www.rdocumentation.org/packages/raster/versions/3.1-5
IGAC. (2012). Mapas de clasificación de las tierras por su vocación de uso a escala 1:100.000. Retrieved April 22, 2020, from https://geoportal.igac.gov.co/contenido/datos-abiertos-agrologia
Longley, P. A., Goodchild, M. F., Maguire, D. J., & Rhind, D. W. (2005). Geographical Information Systems and Science. Retrieved from http://www.amazon.com/Geographic-Information-Systems-Science-Longley/dp/0470870001/ref=sr_1_2?ie=UTF8&qid=1430849641&sr=8-2&keywords=Geographical+Information+Systems+and+Science+longley+2005
Pebesma, E. (2020). sp package | R Documentation. Retrieved April 22, 2020, from https://www.rdocumentation.org/packages/sp/versions/1.4-0
Perpinan, O. (2019). rasterVis package | R Documentation. Retrieved April 22, 2020, from https://www.rdocumentation.org/packages/rasterVis/versions/0.47
USGS. (2020). EarthExplorer. Retrieved April 22, 2020, from https://earthexplorer.usgs.gov/
Xu, H. (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14), 3025–3033. https://doi.org/10.1080/01431160600589179