A continuación, se detallan los pasos esenciales para la construcción de una clasificación de coberturas mediante el uso de una escena producto del satélite Landsat 8, proceso para el cual se implementan algoritmos supervisados y No Supervisados, además de construir como paso intermedio un índice de vegetación NDVI. Su propósito es apoyar la caracterización de municipios priorizados sobre la zona occidental de Bogotá, siendo algunos partes de la sabana de Bogotá. La aplicación de técnicas de percepción remota pretende en esta oportunidad apoyar la clasificación general de coberturas, valiéndose para esto del contraste entre dos métodos de clasificación, a saber: No supervisada y supervisada.

de manera adicional se ha buscado apoyar esta caracterización mediante el cálculo de un índice fundamental, de amplia difusión, el Índice de vegetación de diferencia normalizada (NDVI) cuyo uso en este caso apoya la caracterización de coberturas.

El flujo de trabajo que continuación se detalla involucra diferentes fases: * Introducción y reconocimiento a características de la escena * Corte de la escena a zona priorizada y exploración numérica general} * Cálculo del NDVI, exploración de características y análisis de componentes principales * Clasificación no supervisada * Clasificación supervisada mediante arboles de decisión * Conclusiones preliminares

Exploración de imágenes

En primer lugar, se ubican y descomprimen los archivos obtenidos a través del portal Earth Explorer, esta escena fue capturada en 2018-04-03 y agrupa con una baja cantidad de nubes los municipios conforman el área de interés.

A continuación, se realiza una exploración general al producto obtenido, escena con identificador LC80080572018076LGN00. En primer lugar se realizará una indagación a las bandas que operan en el rango visible del espectro electromagnético.

División del espectro electromagnético

División del espectro electromagnético

# Llamar basnda satelitaes
# azul
b2 <- raster(raslist[2])

# Verde
b3 <- raster(raslist[3])
# Rojo
b4 <- raster(raslist[4])
# Infrarojo cercano (NIR)
b5 <- raster(raslist[5])

la exploración de las bandas puede hacerse al llamar el objetos uno por uno

#Exploración del objeto
b2
## class      : RasterLayer 
## dimensions : 7741, 7581, 58684521  (nrow, ncol, ncell)
## resolution : 30, 30  (x, y)
## extent     : 447285, 674715, 362985, 595215  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : D:/PercepcionRemota_2020/IMG/L82018/LC08_L1TP_008057_20180317_20180403_01_T1_B2.tif 
## names      : LC08_L1TP_008057_20180317_20180403_01_T1_B2 
## values     : 0, 65535  (min, max)

las características impresas muestran información general de la banda:

Comparación de valores puntuales en dos raster Realizado esto se procede a comparar valores que comparten posición en el espacio, pero son leídos en bandas diferentes

# ¿qué sistema de referencia coordenado tiene?
crs(b2)
## CRS arguments:
##  +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
# ¿Cuantas celdas tiene?
ncell(b2)
## [1] 58684521
# ¿cuales son las dimenciones de la banda?
dim(b2)
## [1] 7741 7581    1
# ¿cuál es su resolución espacial?
res(b2)
## [1] 30 30
# ¿cuál es el número de bandas del archivo?
nlayers(b2)
## [1] 1
# Identificar si dos raster comparten iguales caractersiticas de extensión, npumero de columnas, crs, resolución, etc
compareRaster(b2,b3)
## [1] TRUE

una vez con las bandas es posible generar un archivo raster combinando bandas, con lo cual se crea un raster multiespectral

# Combinación de bandas
b543<-stack(b5,b4,b3)
# Exploración de caracteristicas de "b543"
b543
## class      : RasterStack 
## dimensions : 7741, 7581, 58684521, 3  (nrow, ncol, ncell, nlayers)
## resolution : 30, 30  (x, y)
## extent     : 447285, 674715, 362985, 595215  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      : LC08_L1TP_008057_20180317_20180403_01_T1_B5, LC08_L1TP_008057_20180317_20180403_01_T1_B4, LC08_L1TP_008057_20180317_20180403_01_T1_B3 
## min values :                                           0,                                           0,                                           0 
## max values :                                       65535,                                       65535,                                       65535

Para reunir todas las bandas en un solo archivo es posible usar a función RasterStack, este objeto puede reunir todas las bandas sin que ello refiera una combinación de estas. Sin embargo, vale la pena aclarar que para ello se han seleccionado todas las bandas que comparten características, incluidas en ellas extensión y resolución espacial.

Landsat 8 cuenta con una banda pancromática cuya resolución espacial es de 15 metros, ello genera que su extensión sea diferente a la del resto de bandas, por lo que para la compilación en un solo stack esta banda ha sido descartada.

Bandas y composiciones
Las composiciones de color mediante el uso de tres bandas permiten identificar o resaltar características sobre terreno, un ejemplo de ello está dado por la combinación de falso color natural obtenida a partir de las tres bandas que registran información en el rango visible del espectro. Adicionalmente resulta de gran utilidad la combinación de falso color – vegetación obtenida al involucrar la banda que registra información ene l rango de infrarrojo cercano del espectro; estaba banda en particular resulta de gran utilidad para la identificación y caracterización de coberturas vegetales.

Esats combinaciones y sus resultados son expuestas a continuación

# Ploteo individual de bandas
par(mfrow=c(2,2))
plot(landsat@layers[[2]], main= "Azul", col= gray(0:100/100)) # Banda Azul (2)
plot(landsat@layers[[3]], main= "Verde", col= gray(0:100/100)) # Banda Verde (3)
plot(landsat@layers[[4]], main= "Rojo", col= gray(0:100/100)) # Banda Rojo (4)
plot(landsat@layers[[5]], main= "NIR", col= gray(0:100/100)) # Banda NIR (5)

b543<-stack(landsat@layers[[5]],landsat@layers[[4]],landsat@layers[[3]])
b432<-stack(landsat@layers[[4]],landsat@layers[[3]],landsat@layers[[2]])
# Plotear composición Color verdadero
par(mfrow=c(1,2))

plotRGB(b432, axes=TRUE,stretch="hist",main="Landsat Falso color verdadero")
plotRGB(b543, axes=TRUE,stretch="hist",main="Landsat Falso Color - Vegetación NIR")

Ahora bien, la exploración de escenas es un paso imprescindible para reconocer las características que pueden ser aprovechadas en diferentes análisis, sin embargo, resulta usual que estos se realicen de manera particularizada, es decir, para áreas específicas sobre el espacio geográfico. Ante ello resulta importante conocer las posibilidades para seleccionar bandas específicas útiles al trabajo y extraer secciones o porciones de área de la escena útiles al trabajo delimitado. Operaciones que a continuación se ilustran

Extracción de bandas especificas

Los ejercicios que aquí se realizan solo aprovechan las primeras 7 bandas, ante lo cual no es necesario contar con todos los archivos, una forma sencilla de hacer esta extracción es:

# Crear RasterStack para bandas 1 a 7
landsat<-subset(landsat,1:7)

El número de bandas y rango del espectro que cada una cubre varía en diferentes sensores, por lo cual resulta mucho más eficiente nombrar las bandas de acuerdo a su característica principal, en el caso de Landsat 8 se utilizará la clasificación que para este sensor opera.\

Caracterización por bandas, Landsat 8

Caracterización por bandas, Landsat 8

# Renombras bandas de acuerdo a caracterisca espectral
names(landsat)<-c("Costa-Aerosol","Azul","Verde","Rojo","NIR","SWIR1","SWIR2")
names(landsat)
## [1] "Costa.Aerosol" "Azul"          "Verde"         "Rojo"         
## [5] "NIR"           "SWIR1"         "SWIR2"

Corte y enmascaramiento de imagen

para generar un correcto corte de la imagen es necesario superponer un objeto geométrico, en este caso han sido seleccionados algunos municipios cercanos a Bogotá

una mirada rápida sobre la identificación de la zona de estudio

ggplot(data = L8) + # llamar dato
    geom_sf(size = .4, color = "black", fill = "gray95") + # caracteristicas del poligono
    geom_sf(data = Muns, aes(fill = Municipio)) +
    xlab("Longitud") + ylab("Latitud") +
    ggtitle("Escena Landsat 8", subtitle = "(PathRow: 08-57)")

#plot(L8)
#plot(Muns, add = TRUE, col= "gray" )
 #mapview(L8, zcol = "PathRow")+ mapview(Muns,zcol="MPIO_CNMBR")
#rm(L8)

los sistemas de referencia entre la escena y el archivo shapefile son diferentes es necesario homogenizar los CRS antes de ejecutar el corte y enmascaramiento

La diferencia entre estos dos métodos radica en que el corte asume los puntos extremos de extensión de la capa de corte para generar un recuerdo homogéneo en la imagen de salida, mientras que en enmascaramiento se ajusta a los limites geométricos del área a estudiar.

#Ejemplos en ploteo de las dos imagenes
b432Crop<-stack(landsatCrop$Rojo,landsatCrop$Verde,landsatCrop$Azul)
b432Mask<-stack(landsatMask$Rojo,landsatMask$Verde,landsatMask$Azul)

par(mfrow=c(1,2))

plotRGB(b432Crop, axes=TRUE,stretch="hist",main="Color verdadero - Corte")
plotRGB(b432Mask, axes=TRUE,stretch="hist",main="Color verdadero - Mascara")

Estos resultados pueden ser alojados en disco para facilitar su manipulación en un SIG convencional de escritorio

Una alternativa al formato GeoTif es el formato Grid. La notable ventaja de este formato es que permite guardar los nombres de las bandas, los cuales ya previamente han sido definidos, sin embargo, no son muchos los SIG de escritorio que pueden leerlo, frente a esto QGIS es una opción que permite su manipulación.

Relación entre bandas

Una matriz de correlación puede ayudar a identificar posibles relaciones entre el comportamiento registrado por las bandas, información útil para utilizar al máximo los datos recolectados, adquiriendo así la mayor cantidad de información posible con el menor número de insumos necesarios para ello.

#Matriz de correlación bandas 1 y 2
pairs(landsatMask[[1:2]], main = "Costa aerosol(1) - Azul (2)")

el diagrama anterior evidencia una fuerte correlación entre los valores de la banda Costa aerosol y Azul, con lo cual es posible afirmar que al usar una sola de estas bandas no se perderá una cantidad considerable de información que podría aportar la banda con la que se confronta, escenario opuesto al que muestra la relación entre las bandas Roja y NIR, en donde su baja correlación evidencia las diferencias sustanciales en la información que cada una captura, siendo entonces deseable generar contrastes entre estos dos registros.

#Matriz de correlación bandas 4 y 5
pairs(landsatMask[[4:5]], main = "Rojo (4) - NIR (5)")

Valor en píxel y perfiles espectrales

La exploración de relaciones por banda suponen un paso más en la apuesta por caracterizar la zona de estudio mediante la relación de informacióon que el sensor proporciona, esta da una idea de cómo utilizar estos insumos mediante la combinación de bandas, sin embargo, un análisis particularizado demandará conocer el comportamiento de los elementos que allí aparecen (cuerpos de agua, cultivos, bosques, construcción, suelos desnudos, afloramiento rocoso, etc.). Para ello se utilizarán muestras interpretadas sobre la imagen, buscando generalizar el comportamiento espectral de estos elementos en cada una de las bandas. Estas áreas serán utilizadas para generar un muestreo aleatorio que permita acercarse al valor que cada píxel que se tiene dentro de las zonas previamente interpretadas.

Con las categorías definidas en los polígonos es posible crear puntos de muestreo al interior de los mismos, esto con el fin de ligar el valor de pixel a las respectivas categorías

Los perfiles espectrales permiten conocer el comportamiento de ciertas coberturas en cada una de las bandas, es decir, en diferentes rangos del espectro electromagnético, su conocimiento permite identificar comportamientos típicos los cuales resultan de gran utilidad para la posterior distinción y agrupación de estas entidades tomando como criterio su comportamiento espectral

#Contruir perfil espectral
#Crear vector con colores de ploteo

mycolor <- c('darkred','forestgreen','darkgreen','darkblue','yellow','deeppink','chartreuse','gold','cornsilk4')
#transform ms from a data.frame to a matrix
perf <- as.matrix(perf)
# First create an empty plot
plot(0, ylim=c(6000,25000), xlim = c(1,7), type='n', xlab="Bandas", ylab = "Nivel digital (ND)")
# add the different classes
for (i in 1:nrow(perf)){
  lines(perf[i,], type = "l", lwd = 3, lty = 1, col = mycolor[i])
}
# Title
title(main="Perfil espectral imagen landsat", font.main = 2)
# Legend
legend("topleft", rownames(perf),
       cex=0.8, col=mycolor, lty = 1, lwd =3, bty = "n")

El perfil expuesto permite identificar las características que los objetos diferenciados en la imagen asumen en cada una de las bandas, esto en términos de los valores espectrales que permiten singularizar un elemento frente a otro. Los datos usados para su construcción son el insumo base para la clasificación supervisada, pues su labor es servir de datos base para el entrenamiento del modelo con el cual se construirá la interpretación final.

Operaciones matemáticas

Realizados los ejercicios de exploración de imágenes y bandas, se ha seleccionado una zona en la sabana de Bogotá para realizar la exploración numérica respectiva, en esta ocasión se ha seleccionado una imagen del sensor Landsat 8 Nivel 1 LC80080572018076LGN00.

La exploración que a continuación se llevará a cabo tiene como fin identificar las características estadísticas que la información capturada tiene a la hora de identificar relaciones espaciales en el estudio de objetos específicos de la cobertura terrestre, como insumo base se utilizan los archivos previamente manipulados en la sección de Exploración de imágenes

Cálculo Indicie de vegetación

De acuerdo con Jensen (2016) el Normalized Difference Vegetation Index -NDVI es el resultado de una sencilla operación entre las bandas Rojo (4) y NIR (5), cuya relación para su obtención es el resultado de:
\(NDVI=\frac{p_{nir}-p_{rojo}}{p_{nir}+p_{rojo}}\) 

en donde \(p\) corresponde a la banda con las características específicas de acuerdo al sensor, para el caso de Landsat 8. En este caso corresponden a las bandas 5 y 4 respectivamente.

# NIR = 5, rojo = 4.

ndvi <- (b5-b4)/(b5+b4)
plot(ndvi, col = rev(terrain.colors(10)), main = "Landsat 8 - NDVI")
plot(Muns,add=TRUE)

Histograma

El índice permite una identificación rápida de elementos, cuya prioridad son áreas con cobertura vegetal. Ademas de ello permite diferenciar desde cuerpos de agua y construcciones, hasta la identificación de particularidades en la vegetación presente en la zona. El histograma resulta de gran utilidad a la hora de identificar valores anómalos o outliers, siempre tomando en cuenta que los valores más cercanos a 1 indicaran una mayor presencia y desarrollo de coberturas vegetales.

#Distribución de valor ndvi para área seleccionada
ndvi2<-as.data.frame(getValues(ndvi))
ndvi2<-as.data.frame(ndvi2[!is.na(ndvi2$`getValues(ndvi)`),])
names(ndvi2)<-"Valores"
ggplot(ndvi2, aes(ndvi2$Valores, fill = cut(ndvi2$Valores, 50))) +
  geom_histogram(show.legend = FALSE) +
 scale_fill_discrete(h = c(0, 160))+
  labs(x = "NDVI", y = "Frecuencia") +
  ggtitle("Distribución general valores NDVI")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Umbral de análisis

El NDVI permite diferenciar algunas coberturas expuestas en la escena, este índice oscila entre -1 y 1, siendo los valores más altos indicativos de mayor cobertura vegetal. Un ejemplo aplicado de ello puede realizarse identificando los pixeles con valores por encima de 0.4, definiéndoles como valores característicos de coberturas vegetales. A continuación, se asumirá este intervalo para enmascarar la imagen.

veg <- reclassify(ndvi, cbind(-Inf,0.4, NA))
plot(veg, main='Vegetación')

y a continuación se privilegian solo aquella áreas que están entre 0.25 y 0.3

Un contraste entre la imagen combinada en Color verdadero y la clasificación permite hacer algunas inferencias

viewRGB(landsatMask, 4, 3, 2)+ mapview(land)
## Warning in rasterCheckSize(x, maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '
## Warning: 'mapview::addHomeButton' is deprecated.
## Use 'leafem::addHomeButton' instead.
## See help("Deprecated") and help("mapview-deprecated").
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '

El mismo método puede ser usado para crear diferentes tipos de clasificación y así lograr diferenciar tipos de cobertura a partir del comportamiento del índice NDVI. Sin embargo, no es claro si esta clasificación privilegia en la imagen zonas con suelos desnudos, coberturas tipo arbusto o pastos, pues no logra tomar en este caso un elemento particular.

vegc <- reclassify(ndvi, c(-Inf,0.25,1, 0.25,0.3,2, 0.3,0.4,3, 0.4,0.5,4, 0.5,Inf, 5))
plot(vegc,col = rev(terrain.colors(4)), main = 'NDVI basado en umbral')

Análisis de componentes principales ACP

set.seed(1)
sr <- sampleRandom(landsatMask, 10000)
plot(sr[,c(4,5)], main = "NIR-Red plot")

pca <- prcomp(sr, scale = TRUE)
pca
## Standard deviations (1, .., p=7):
## [1] 2.17890210 1.04142328 0.99069002 0.35038190 0.22394855 0.09985967 0.05885665
## 
## Rotation (n x k) = (7 x 7):
##                      PC1         PC2         PC3         PC4         PC5
## Costa.Aerosol 0.41217472 -0.21483440  0.33120398 -0.52056977  0.13006151
## Azul          0.42940697 -0.17402851  0.28373104 -0.29289677  0.03269374
## Verde         0.44109644 -0.19924847  0.09212357  0.39330239 -0.05132774
## Rojo          0.44418176  0.01848075  0.09636280  0.63334095 -0.12948199
## NIR           0.07156128 -0.71661099 -0.65097696 -0.08266589 -0.19128840
## SWIR1         0.35639749  0.31801819 -0.51658595 -0.02806030  0.70969036
## SWIR2         0.34888102  0.51850930 -0.31842744 -0.28247130 -0.64989487
##                        PC6         PC7
## Costa.Aerosol -0.200529886  0.58833748
## Azul           0.013193182 -0.78599553
## Verde          0.752792131  0.18229079
## Rojo          -0.612142555  0.02168817
## NIR           -0.117830982 -0.01635969
## SWIR1          0.006958914 -0.02209061
## SWIR2          0.065319502  0.04017580
screeplot(pca)

El análisis de componentes principales ACP permite identificar qué banda agrupa exitosamente la mayor cantidad de información, esto en contraste con el total de las bandas involucradas en el análisis, de su rápida lectura destaca que la banda Costa aerosol agrupa una buena cantidad de características presentes en las demás bandas. Si se tiene en cuenta que su correlación con la banda Azul es bastante alta se puede inferir que el uso de cualquiera de estas dos bandas es de gran importancia a la hora de caracterizar los elementos a singularizar en la imagen.

Adicionalmente la dispersión hallada entre la banda Roja y NIR indican una gran diferencia entre ambas, es decir, cada una de ellas representa con elementos que no son de fácil interpretación en la otra. Ante estos resultados y apoyados en el ACP se infiere que una combinación de bandas con características fuertes en este análisis será NIR-Rojo-Azul, características de gran utilidad en la diferenciación y caracterización de unidades espaciales sobre la zona priorizada.

pci <- predict(landsatMask, pca, index = 1:2)
plot(pci[[1]])

Comparación entre la combinación mencionada y los resultados de ACP

pc2 <- reclassify(pci[[2]], c(-Inf,0,1,0,Inf,NA))
par(mfrow = c(1,2))
# Comparar imagenes
viewRGB(landsatMask, 5, 4, 3)+ mapview(pc2)
## Warning in rasterCheckSize(x, maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '
## Warning: 'mapview::addHomeButton' is deprecated.
## Use 'leafem::addHomeButton' instead.
## See help("Deprecated") and help("mapview-deprecated").
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '

Clasificación no supervisada

Los algoritmos de clasificación para imagenes de satélite son variados, la practica se desarrolla fundamentalmente utilizando el K-Means, con el fin de ilustrar de manera general sus posibilidades y funcionamiento.

El algortimo realiza una clasificación automatica basada en el comportameinto similar en los valores de pixel, así lo que el resultado arrojará será una categorzación por agrupación en valores de pixel. Una vez realizada la grupación por valores el algoritmo busca el centroide de estas nuevas clasificaciones, luego agrupa de forma iterativa los píxeles en la clase más cercana utilizando para ello la distancia mínima.

Texto alternativo

Texto alternativo

A continuación se realizará una prueba de este algortimo a los datos obtenidos del calculo de NDVI

Con el fin de ilustrar de amnera puntualizada los resultados de este indice se ha priorizado el municipio de Bojacá, municipio en el cuál el NDVI se comporta de la siquiente manera.

#Operación soreb NDVI
# NIR = 5, rojo = 4.
ndvi <- (b5-b4)/(b5+b4)

Bj<-Muns[Muns$MPIO_CCDGO==25099,]
ndvi<-crop(ndvi,Bj)
#ndvi<-mask(ndvi,Bj)


nr<- getValues(ndvi) # Estracción de valores ndvi a matriz numérica independiente

nr1<-as.data.frame(!is.na(nr)) # Data frame para creación de histograma

#Distribución de valor ndvi para área seleccionada
ggplot(nr1, aes(nr, fill = cut(nr, 100))) +
  geom_histogram(show.legend = FALSE) +
 scale_fill_discrete(h = c(0, 180))+
  labs(x = "NDVI", y = "Frecuencia") +
  ggtitle("Histograma NDVI")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

K Means

Los valores del NDVI han sido almacanado en una amtriz diferente, sobre la cual se realizara la agrupación de valores similares bajo la aplicación del algoritmo K-means

la creación del objeto kmncluster arroja como resultado una lista con 9 elementos, cuya extensión es de 267124, emanados del objeto original, procedente de los valores calculados en el NDVI del área priorizada. Esta característica es de vital importancia pues si la longitud de los elementos varía, será imposible generar un raster para el resultado de la aplicación de K-means, Es por ello que resulta recomendable trabajar con imágenes cortadas crop y no con imágenes enmascaradas mask pues en esta última operación se crean valores NA los cuales serán eliminados en el agrupamiento de valores, en consecuencia las longitudes de los objetos serán incompatibles.

Con los resultados obtenidos será necesario crear un raster con la misma extensión geográfica que los datos de origen

el objeto knr agrupa el resultado en formato Raster, es decir, aquí se almacenan los resultados. Vale la pena recordar que la iteración realizada buscó definir 7 clases, por lo cual el resultado defnirá esos valores agrupados en objettos temáticos con colores difereniales

#Enmascarar resultados para ploteo
ndvi2<-mask(ndvi,Bj)
knr2<-mask(knr,Bj)
# Define a color vector for 10 clusters (learn more about setting the color later)
mycolor <- c("#fef65b","#ff0000", "#daa520","#0000ff","#00ff00","#cbbeb5",
             "#c3ff5b")
#, "#ff7373", "#00ff00", "#808080"
par(mfrow = c(1,2))
plot(ndvi2, col = rev(terrain.colors(10)), main = 'Landsat-NDVI')
plot(Bj,add=TRUE)
plot(knr2, main = 'Clasificación No Supervisada', col = mycolor )
plot(Bj,add=TRUE)

# Comparación visual a clasificación no supervisada
m1<-viewRGB(landsatMask,5,4,3)
## Warning in rasterCheckSize(x, maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '
## Warning: 'mapview::addHomeButton' is deprecated.
## Use 'leafem::addHomeButton' instead.
## See help("Deprecated") and help("mapview-deprecated").
m2<-mapview(knr2, col = mycolor) # Clasificación No Supervisada

sync(m1,m2)
## Warning: 'mapview::sync' is deprecated.
## Use 'leafsync::sync' instead.
## See help("Deprecated") and help("leafsync-deprecated").
## Warning: 'mapview::latticeView' is deprecated.
## Use 'leafsync::latticeView' instead.
## See help("Deprecated") and help("mapview-deprecated").
rm(m1,m2)

Clasificación supervisada

La clasificación supervisada es un método que se vale de muestras que identifican un tipo de cobertura sobre el espacio, tomando los valores que esta categoría asumen en las diferentes bandas de la imagen, a partir de estas muestras se generan un modelo con el cuál se agrupan los valores que mejor representan este comportamiento en la matriz numérica de la imagen. El paso inicial que se dio en pro de generar esta clasificación estuvo dado por la construcción del perfil espectral, las muestras que hicieron posible identificar y diferenciar el comportamiento de estos elementos sobre la imagen, ahora serán de utilidad en la creación de un modelo que permita delimitar en una imagen estas áreas de cobertura.

En este caso se privilegia la construcción de un árbol de clasificación, método que servirá para entrenar al modelo con le fin de identificar y diferenciar aquellos pixeles que respondan al comportamiento numérico que las categorías han señalado previamente.

PtMuestra2<-PtMuestra@data
PtMuestra2<-PtMuestra2[,-1]

# Entrenar el modelo
cart <- rpart(as.factor(Clase)~., data=PtMuestra2, method = 'class', minsplit = 5)
# print(model.class)
# Plot the trained classification tree
plot(cart, uniform=TRUE, main="Arbol de clasificación")
text(cart, cex = 0.55)

A partir de este árbol de clasificación para la imagen seleccionada se agrupan los valores de pixel que mejor representan este comportamiento.

Clasificación basada en arbol

pr2018 <- ratify(pr2018)
rat <- levels(pr2018)[[1]]
rat$legend <- Clases$classnames1
levels(pr2018) <- rat
pr2018<-crop(pr2018,Muns)
pr2018<-mask(pr2018,Muns)
# Escala de color para ploteo
classcolor <- c("#C49871", "#07F332", "#046C17", "#DE4CC3", "#A9DE4C", "#D82020", "#FBF65D", "#C0C0C0")
levelplot(pr2018, maxpixels = 1e6,
          col.regions = classcolor,
          scales=list(draw=FALSE),
          main = "Clasificación supervisada basada en árbol de decisiones")

Una visualización dinámica de los resultados ayuda a corroborar la calidad de estos, en este caso el resultado obtenido mediante la clasificación supervisada es adecuado, sin embargo, es necesario generar más muestra para diferenciar zonas de Bosque frente a Cultivos y Pastos pues las categorías en este resultado no logran diferencia de manera adecuada estas coberturas.

Llama poderosamente la atención la clasificación realizada sobre la zona occidental del área de análisis en donde, de acuerdo con la clasificación, predominan poderosamente los arbustos. Este resultado puede deberse al cambio en el gradiente altitudinal y con él, el cambio en las características de la vegetación asociada, pues es una zona que se encuentra en el limite de la altiplanicie, en donde el descenso altitudinal está acompañado de una variación en las características físicas del espacio. Esta es apenas una inferencia general que debe ser estudiada con mayor detalle.

#Visualizar resultados 
viewRGB(landsatMask,4,3,2)+mapview(pr2018, col = classcolor)
## Warning in rasterCheckSize(x, maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '
## Warning: 'mapview::addHomeButton' is deprecated.
## Use 'leafem::addHomeButton' instead.
## See help("Deprecated") and help("mapview-deprecated").
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 2049138 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  2049138 '

Conclusión

La exploración de una imagen de satélite debe realizarse en varios niveles, tanto de manera visual como estadística, el conocimiento detallado de sus componentes permite generar análisis basados en decisiones y pruebas, con lo cuál los resultados de los mismos tendrán un mayor soporte y por ende una mayor comprensión sobre los resultados generados para la indagación del fenómeno de interés.

Las clasificaciones no supervisadas son un paso necesario para la comprensión de los posibles agrupamientos que sobre una imagen se pueden realizar, sin embargo a mayor diferenciación sobre la zona y sus características, menos confiables serán los resultados de este procedimiento, se recomienda su uso solo como paso intermedio.

Por otro lado, las clasificaciones supervisadas minimizan la incertidumbre sobre los resultados generados, los cuales dependen en gran medida del proceso de muestreo que se realice, se recomienda que este sea lo más detallado posible. En este caso se usó un modelo sencillo de clasificación, pero más allá de este, existen modelos de mayor potencia, su estudio e incorporación para la interpretación de coberturas son un campo rico de análisis en la implementación y comprensión de la percepción remota como técnica de análisis espacial en diferentes disciplinas.