Mapeo de cuerpos de agua utilizando índices espectrales NDWI y MNDWI con la banda SWIR fusionada de imágenes Sentinel-2 a partir de algoritmos pan-sharpening

La siguiente metodología de procesamiento y análisis de imágenes satelitales en R se compone de:

  1. Proceso 1: Preliminar
  2. Proceso 2: Implementación de algoritmos pan-sharpening
  3. Proceso 3: Cálculo de índices espectrales de agua
  4. Proceso 4: Evaluación de la fusión de imágenes
  5. Proceso 5: Clasificación de imágenes implementando CART
  6. Proceso 6: Evaluación de la clasificación de imágenes implementando CART

1. Proceso 1: Preliminar

library (terra)
## terra version 1.1.4
library(raster)
## Loading required package: sp
library(rgdal)
## rgdal: version: 1.5-23, (SVN revision 1121)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.1.4, released 2020/10/20
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
## 
## Attaching package: 'rgdal'
## The following object is masked from 'package:terra':
## 
##     project
library(ggplot2)
library(knitr)
## 
## Attaching package: 'knitr'
## The following object is masked from 'package:terra':
## 
##     spin
library(rpart)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:raster':
## 
##     intersect, select, union
## The following objects are masked from 'package:terra':
## 
##     collapse, desc, intersect, near, select, union
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(RColorBrewer)

#DIRECTORIO DE TRABAJO
setwd('/Users/lauranatalia/Desktop/IMG_DATA')

#BANDAS 10M S2
B4 <- rast('R10m/T18NXL_20200104T152639_B04_10m.jp2')
B3 <- rast('R10m/T18NXL_20200104T152639_B03_10m.jp2')
B2 <- rast('R10m/T18NXL_20200104T152639_B02_10m.jp2')
B8 <- rast('R10m/T18NXL_20200104T152639_B08_10m.jp2')

#BANDAS 20M S2
B4_1 <- rast('R20m/T18NXL_20200104T152639_B04_20m.jp2')
B3_1 <- rast('R20m/T18NXL_20200104T152639_B03_20m.jp2')
B2_1 <- rast('R20m/T18NXL_20200104T152639_B02_20m.jp2')
B8_1 <- rast('R20m/T18NXL_20200104T152639_B8A_20m.jp2')
B11 <- rast('R20m/T18NXL_20200104T152639_B11_20m.jp2')

Salidas gráficas de la Banda 2, Banda 3, Banda 4, Banda 8 y Banda 11.

par(mfrow = c(3,2))
plot(B2, main = "B2 - AZUL", col = gray(0:100/100))
plot(B3, main = "B3 - VERDE", col = gray(0:100/100))
plot(B4, main = "B4 - ROJO", col = gray(0:100/100))
plot(B8, main = "B8 - VNIR", col = gray(0:100/100))
plot(B11, main = "B11 - SWIR", col = gray(0:100/100))

Coeficiente de correlación de Pearson entre las bandas antes visualizadas

#COEFICIENTE CORRELACIÓN DE PEARSON
S2 <- c(B4_1, B3_1, B2_1, B11)
pairs(S2[[1:4]], hist=TRUE, cor=TRUE, use="pairwise.complete.obs", maxpixels=100000)

2. Proceso 2: Implementación de algoritmos pan-sharpening

La fusión de datos con enfoque pan-sharpening se realiza en el software Matlab, utilizando el siguiente de código libre, desarrollado en G. Vivone, L. Alparone, J. Chanussot, M. Dalla Mura, A. Garzelli, G. Licciardi, R. Restaino, and L. Wald, “A Critical Comparison Among Pansharpening Algorithms”, IEEE Transactions on Geoscience and Remote Sensing, vol. 53, no. 5, pp. 2565-2586, Mayo 2015.

3. Proceso 3: Cálculo de índices espectrales de agua

#CALCULO DE INDICE NDWI DE 20M
NDWI <- ((B3-B8)/(B3+B8))
#plot(NDWI, main = "NDWI")
#writeRaster(NDWI, filename="NDWI.tif", overwrite=TRUE)

#CALCULO DE INDICE MNDWI DE 20M
MNDWI <- ((B3_1-B11)/(B3_1+B11))
#plot(MNDWI, main = "MNDWI")
#writeRaster(MNDWI, filename="MNDWI.tif", overwrite=TRUE)

B11_IHS <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_231.tif')
B11_IHS_1 <- B11_IHS[[2]]
res(B11_IHS_1)
## [1] 10 10
B11_B <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_241.tif')
B11_B_1 <- B11_B[[2]]
res(B11_B_1)
## [1] 10 10
B11_GS <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_251.tif')
B11_GS_1 <- B11_GS[[2]]
res(B11_GS_1)
## [1] 10 10
B11_PCA <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_261.tif')
B11_PCA_1 <- B11_PCA[[2]]
res(B11_PCA_1)
## [1] 10 10
#CALCULO DE INDICE MNDWI DE 10M IHS
MNDWI1 <- ((B3-B11_IHS_1)/(B3+B11_IHS_1))
#plot(MNDWI1, main = "MNDWI")
#writeRaster(MNDWI1, filename="MNDWI1.tif", overwrite=TRUE)

#CALCULO DE INDICE MNDWI DE 10M BROVEY
MNDWI2 <- ((B3-B11_B_1)/(B3+B11_B_1))
#plot(MNDWI2, main = "MNDWI - BROVEY")
#writeRaster(MNDWI2, filename="MNDWI2.tif", overwrite=TRUE)

#CALCULO DE INDICE MNDWI DE 10M GS
MNDWI3 <- ((B3-B11_GS_1)/(B3+B11_GS_1))
#plot(MNDWI3, main = "MNDWI - GS")
#writeRaster(MNDWI3, filename="MNDWI3.tif", overwrite=TRUE)

#CALCULO DE INDICE MNDWI DE 10M PCA
MNDWI4 <- ((B3-B11_PCA_1)/(B3+B11_PCA_1))
#plot(MNDWI4, main = "MNDWI - PCA")
#writeRaster
par(mfrow = c(3,2))
plot(NDWI, main = "NDWI 10M")
plot(MNDWI, main = "MNDWI 20M")
plot(MNDWI1, main = "MNDWI 10 M - IHS")
plot(MNDWI2, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3, main = "MNDWI 10 M - GS")
plot(MNDWI4, main = "MNDWI 10 M - PCA")

#AREAE <- readOGR(dsn = "AREA DE ESTUDIO", layer = "AREA_ESTUDIO")
A1 <- readOGR(dsn = "AREA DE ESTUDIO", layer = "A1")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/lauranatalia/Desktop/IMG_DATA/AREA DE ESTUDIO", layer: "A1"
## with 1 features
## It has 1 fields
A2 <- readOGR(dsn = "AREA DE ESTUDIO", layer = "A2")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/lauranatalia/Desktop/IMG_DATA/AREA DE ESTUDIO", layer: "A2"
## with 1 features
## It has 1 fields
A3 <- readOGR(dsn = "AREA DE ESTUDIO", layer = "A3")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/lauranatalia/Desktop/IMG_DATA/AREA DE ESTUDIO", layer: "A3"
## with 1 features
## It has 1 fields
A4 <- readOGR(dsn = "AREA DE ESTUDIO", layer = "A4")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/lauranatalia/Desktop/IMG_DATA/AREA DE ESTUDIO", layer: "A4"
## with 1 features
## It has 1 fields
A5 <- readOGR(dsn = "AREA DE ESTUDIO", layer = "A5")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/lauranatalia/Desktop/IMG_DATA/AREA DE ESTUDIO", layer: "A5"
## with 1 features
## It has 1 fields
A6 <- readOGR(dsn = "AREA DE ESTUDIO", layer = "A6")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/lauranatalia/Desktop/IMG_DATA/AREA DE ESTUDIO", layer: "A6"
## with 1 features
## It has 1 fields
#NDWI
NDWI1 <- crop(NDWI, A1)
NDWI2 <- crop(NDWI, A2)
NDWI3 <- crop(NDWI, A3)
NDWI4 <- crop(NDWI, A4)
NDWI5 <- crop(NDWI, A5)
NDWI6 <- crop(NDWI, A6)

#MNDWI 20 M
MNDWI_1 <- crop(MNDWI, A1)
MNDWI_2 <- crop(MNDWI, A2)
MNDWI_3 <- crop(MNDWI, A3)
MNDWI_4 <- crop(MNDWI, A4)
MNDWI_5 <- crop(MNDWI, A5)
MNDWI_6 <- crop(MNDWI, A6)

#MNDWI1 10 M
MNDWI1_1 <- crop(MNDWI1, A1)
MNDWI1_2 <- crop(MNDWI1, A2)
MNDWI1_3 <- crop(MNDWI1, A3)
MNDWI1_4 <- crop(MNDWI1, A4)
MNDWI1_5 <- crop(MNDWI1, A5)
MNDWI1_6 <- crop(MNDWI1, A6)

#MNDWI2 10 M
MNDWI2_1 <- crop(MNDWI2, A1)
MNDWI2_2 <- crop(MNDWI2, A2)
MNDWI2_3 <- crop(MNDWI2, A3)
MNDWI2_4 <- crop(MNDWI2, A4)
MNDWI2_5 <- crop(MNDWI2, A5)
MNDWI2_6 <- crop(MNDWI2, A6)

#MNDWI3 10 M
MNDWI3_1 <- crop(MNDWI3, A1)
MNDWI3_2 <- crop(MNDWI3, A2)
MNDWI3_3 <- crop(MNDWI3, A3)
MNDWI3_4 <- crop(MNDWI3, A4)
MNDWI3_5 <- crop(MNDWI3, A5)
MNDWI3_6 <- crop(MNDWI3, A6)

#MNDWI4 10 M
MNDWI4_1 <- crop(MNDWI4, A1)
MNDWI4_2 <- crop(MNDWI4, A2)
MNDWI4_3 <- crop(MNDWI4, A3)
MNDWI4_4 <- crop(MNDWI4, A4)
MNDWI4_5 <- crop(MNDWI4, A5)
MNDWI4_6 <- crop(MNDWI4, A6)
#SUBAREA 1
par(mfrow = c(3,2))
plot(NDWI1, main = "NDWI 10 M")
plot(MNDWI_1, main = "MNDWI 20 M")
plot(MNDWI1_1, main = "MNDWI 10 M - IHS")
plot(MNDWI2_1, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3_1, main = "MNDWI 10 M - GS")
plot(MNDWI4_1, main = "MNDWI 10 M - PCA")

#SUBAREA 2
par(mfrow = c(3,2))
plot(NDWI2, main = "NDWI 10 M")
plot(MNDWI_2, main = "MNDWI 20 M")
plot(MNDWI1_2, main = "MNDWI 10 M - IHS")
plot(MNDWI2_2, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3_2, main = "MNDWI 10 M - GS")
plot(MNDWI4_2, main = "MNDWI 10 M - PCA")

#SUBAREA 3
par(mfrow = c(3,2))
plot(NDWI3, main = "NDWI 10 M")
plot(MNDWI_3, main = "MNDWI 20 M")
plot(MNDWI1_3, main = "MNDWI 10 M - IHS")
plot(MNDWI2_3, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3_3, main = "MNDWI 10 M - GS")
plot(MNDWI4_3, main = "MNDWI 10 M - PCA")

#SUBAREA 4
par(mfrow = c(3,2))
plot(NDWI4, main = "NDWI 10 M")
plot(MNDWI_4, main = "MNDWI 20 M")
plot(MNDWI1_4, main = "MNDWI 10 M - IHS")
plot(MNDWI2_4, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3_4, main = "MNDWI 10 M - GS")
plot(MNDWI4_4, main = "MNDWI 10 M - PCA")

#SUBAREA 5
par(mfrow = c(3,2))
plot(NDWI5, main = "NDWI 10 M")
plot(MNDWI_5, main = "MNDWI 20 M")
plot(MNDWI1_5, main = "MNDWI 10 M - IHS")
plot(MNDWI2_5, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3_5, main = "MNDWI 10 M - GS")
plot(MNDWI4_5, main = "MNDWI 10 M - PCA")

#SUBAREA 6
par(mfrow = c(3,2))
plot(NDWI6, main = "NDWI 10 M")
plot(MNDWI_6, main = "MNDWI 20 M")
plot(MNDWI1_6, main = "MNDWI 10 M - IHS")
plot(MNDWI2_6, main = "MNDWI 10 M - BROVEY")
plot(MNDWI3_6, main = "MNDWI 10 M - GS")
plot(MNDWI4_6, main = "MNDWI 10 M - PCA")

4. Proceso 4: Evaluación de la fusión de imágenes

#/////////////////////////////COEFICIENTE DE CORRELATION/////////
RGB = rast('PANSHARPENING/RGB4112_10.tif')
PAN1_1 <- rast('R10m/T18NXL_20200104T152639_B03_10m.jp2')
IHS1 <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_231.tif')
B1 <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_241.tif')
GS1 <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_251.tif')
PCA1 <- rast('PANSHARPENING/B11/T18NXL_20200104T152639_B04_261.tif')
S3 <- c(PAN1_1, IHS1, B1, GS1, PCA1)
S4 <- c(RGB, IHS1, B1, GS1, PCA1)
pairs(S3[[1:5]], hist=TRUE, cor=TRUE, use="pairwise.complete.obs", maxpixels=100000)

pairs(S4[[1:7]], hist=TRUE, cor=TRUE, use="pairwise.complete.obs", maxpixels=100000)

5. Proceso 5: Clasificación de imágenes implementando CART

La implementación de CART se realiza a partir del desarrollo computacional en Lizarazo, I, 2021. A tutorial on land cover classification from satellite imagery using R. Available at https://rpubs.com/ials2un/cart_landsat8.

#COMPOSICI?N DE IM?GENES
#NDWI 10 METROS
NDWIC <- c(B3,B8,NDWI)
names(NDWIC) <- c("GREEN", "NIR", "NDWI")   

#MNDWI 20 M
MNDWIC <- c(B3_1,B11,MNDWI)
names(MNDWIC) <- c("GREEN", "SWIR", "MNDWI")  

#MNDWI 10M
MNDWIBC <- c(B3,B11_B_1,MNDWI2)
names(MNDWIBC) <- c("GREEN", "SWIR10M", "MNDWIB")

#ENTRENAMIENTO DE DATOS

#DATOS DE REFERENCIA

RDatos <- vect("Cobertura_tierra_2010_2012/SHP/Cobertura_tierra_2010_2012_C1.shp")
unique(RDatos$CODIGO)
##  [1]    233    323   3131   3132   3222    111    142    231  32111    243
## [11]    244  31121    112    125    131    242    315    514 321111  31221
## [21]   3221    122    241    225    121    314    512   2151    141    333
## [31] 321113    411    211    245  31111    511   3232 321112 321121  32222
## [41]    413    513  32122   1315   2222    313    332   5141  31211
unique(RDatos$LEYENDA3N)
##  [1] "2.3.3. Pastos enmalezados"                              
##  [2] "3.2.3. Vegetacion secundaria o en transicion"           
##  [3] "3.1.3. Bosque fragmentado"                              
##  [4] "3.2.2. Arbustal"                                        
##  [5] "1.1.1. Tejido urbano continuo"                          
##  [6] "1.4.2. Instalaciones recreativas"                       
##  [7] "2.3.1. Pastos limpios"                                  
##  [8] "3.2.1. Herbazal"                                        
##  [9] "2.4.3. Mosaico de cultivos, pastos y espacios naturales"
## [10] "2.4.4. Mosaico de pastos con espacios naturales"        
## [11] "3.1.1. Bosque denso"                                    
## [12] "1.1.2. Tejido urbano discontinuo"                       
## [13] "1.2.5. Obras hidraulicas"                               
## [14] "1.3.1. Zonas de extraccion minera"                      
## [15] "2.4.2. Mosaico de pastos y cultivos"                    
## [16] "3.1.5. Plantacion forestal"                             
## [17] "5.1.4. Cuerpos de agua artificiales"                    
## [18] "3.1.2. Bosque abierto"                                  
## [19] "1.2.2. Red vial, ferroviaria y terrenos asociados"      
## [20] "2.4.1. Mosaico de cultivos"                             
## [21] "2.2.5. Cultivos confinados"                             
## [22] "1.2.1. Zonas industriales o comerciales"                
## [23] "3.1.4. Bosque de galeria y ripario"                     
## [24] "5.1.2. Lagunas, lagos y cienagas naturales"             
## [25] "2.1.5. Tuberculos"                                      
## [26] "1.4.1. Zonas verdes urbanas"                            
## [27] "3.3.3. Tierras desnudas y degradadas"                   
## [28] "4.1.1. Zonas Pantanosas"                                
## [29] "2.1.1. Otros cultivos transitorios"                     
## [30] "2.4.5. Mosaico de cultivos con espacios naturales"      
## [31] "5.1.1. Rios (50 m)"                                     
## [32] "4.1.3. Vegetacion acuatica sobre cuerpos de agua"       
## [33] "5.1.3. Canales"                                         
## [34] "2.2.2. Cultivos permanentes arbustivos"                 
## [35] "3.3.2. Afloramientos rocosos"
unique(RDatos$CCLC)
## [1] 2 3 1 5 4
#Guardo los puntos para siempre trabajar con el mismo árbol
#PT <- spatSample(RDatos, size=350, method="random", strata = RDatos$CCLC, chess="")
#writeVector(PT, "PT", filetype="ESRI Shapefile", overwrite=FALSE)

PT<- vect("Cobertura_tierra_2010_2012/SHP/PT.shp")
unique(PT$CCLC)
## [1] 2 3 1 5 4
NPT <- as(PT, "Spatial")
NRDatos<- as(RDatos, "Spatial")
PT$class <- over(NPT, NRDatos)$CCLC
#CLASIFICACION CON NDWI 10 M

XY <- as.matrix(geom(PT)[,c('x','y')])
DataFNDWI <- extract(NDWIC, XY)
#head(DataF)
SDatos<- data.frame(class = as.factor(PT$class), DataFNDWI)
#head(SDatos)
#tail(SDatos)
summary(SDatos$class)
##   1   2   3   4   5 
## 500 500 500 500 500
ModeloCART <- rpart(as.factor(class)~., data = SDatos, method = 'class', minsplit = 3)
plot(ModeloCART, uniform=TRUE, main="?rbol de Clasificaci?n")
text(ModeloCART, cex=0.8 ,digits=3, col="blue")

ClasificacionNDWI <- predict(NDWIC, ModeloCART, na.rm = TRUE)
plot(ClasificacionNDWI)

#ASIGNACIÓN DE UNA CLASE PARA VALORES INFERIORES A 1
ClasificacionNDWI$X6 <- ClasificacionNDWI$X5
ClasificacionNDWI$X6[is.na(ClasificacionNDWI)] <- 1.0
#CLASIFICACION CON MNDWI 20 CM
DataFMNDWI <- extract(MNDWIC, XY)
#head(DataF)
SDatosMNDWI<- data.frame(class = as.factor(PT$class), DataFMNDWI)
#head(SDatosMNDWI)
#tail(SDatosMNDWI)
summary(SDatosMNDWI$class)
##   1   2   3   4   5 
## 500 500 500 500 500
ModeloCART2 <- rpart(as.factor(class)~., data = SDatosMNDWI, method = 'class', minsplit = 3)
plot(ModeloCART2, uniform=TRUE, main="Árbol de Clasificación")
text(ModeloCART2, cex=0.8 ,digits=3, col="blue")

ClasificacionMNDWI <- predict(MNDWIC, ModeloCART2, na.rm = TRUE)
plot(ClasificacionMNDWI)

#ASIGNACIÓN DE UNA CLASE PARA VALORES INFERIORES A 1
ClasificacionMNDWI$X6 <- ClasificacionMNDWI$X5
ClasificacionMNDWI$X6[is.na(ClasificacionMNDWI)] <- 1.0
#CLASIFICACION CON MNDWI 10 M - SWIR BROVEY
DataFMNDWIB <- extract(MNDWIBC, XY)
#head(DataF)
SDatosMNDWIB<- data.frame(class = as.factor(PT$class), DataFMNDWIB)
#head(SDatosMNDWI)
#tail(SDatosMNDWI)
summary(SDatosMNDWIB$class)
##   1   2   3   4   5 
## 500 500 500 500 500
ModeloCART3 <- rpart(as.factor(class)~., data = SDatosMNDWIB, method = 'class', minsplit = 3)
plot(ModeloCART3, uniform=TRUE, main="Árbol de Clasificación")
text(ModeloCART3, cex=0.8 ,digits=3, col="blue")

ClasificacionMNDWIB <- predict(MNDWIBC, ModeloCART3, na.rm = TRUE)
plot(ClasificacionMNDWIB)

#ASIGNACIÓN DE UNA CLASE PARA VALORES INFERIORES A 1
ClasificacionMNDWIB$X6 <- ClasificacionMNDWIB$X5
ClasificacionMNDWIB$X6[is.na(ClasificacionMNDWIB)] <- 1.0

6. Proceso 6: Evaluación de la clasificación de imágenes implementando CART

Evaluación Cuantitativa

#NDWI
set.seed(63)
# number of folds
k <- 5
j <- sample(rep(1:k, each = round(nrow(SDatos))/k))
table(j)
## j
##   1   2   3   4   5 
## 500 500 500 500 500
x <- list()
for (k in 1:5) {
  train <- SDatos[j!= k, ]
  test <- SDatos[j == k, ]
  cart <- rpart(as.factor(class)~., data=train, method = 'class',
                minsplit = 5)
  pclass <- predict(cart, test, na.rm = TRUE)
  # assign class to maximum probablity
  pclass <- apply(pclass, 1, which.max)
  # create a data.frame using the reference and prediction
  x[[k]] <- cbind(test$class, as.integer(pclass))
}
class <- c("Territorios Artificializados", "Territorios Agrícolas", "Bosques y Áreas Seminaturales", "Áreas Húmedas", "Superficies de Agua")
classdf <- data.frame(classvalue = seq(1:5), classnames = class)

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
# confusion matrix
conmat <- table(y)
# 
print(conmat)
##         predicted
## observed   1   2   3   4   5
##        1 291 120  42  43   4
##        2  54 298 115  33   0
##        3  27  80 369  23   1
##        4  52 192  92 127  37
##        5   8  61  51  41 339
n <- sum(conmat)
diag <- diag(conmat)

# OVERALL
OANDWI <- sum(diag)/n
OANDWI 
## [1] 0.5696
#KAPPA
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
# predicted cases per class
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
KNDWI <- (OANDWI - expAccuracy) / (1 - expAccuracy)
KNDWI
## [1] 0.462
#MNDWI 20 M
k <- 5
j <- sample(rep(1:k, each = round(nrow(SDatosMNDWI))/k))
table(j)
## j
##   1   2   3   4   5 
## 500 500 500 500 500
x <- list()
for (k in 1:5) {
  train <- SDatosMNDWI[j!= k, ]
  test <- SDatosMNDWI[j == k, ]
  cart <- rpart(as.factor(class)~., data=train, method = 'class',
                minsplit = 5)
  pclass <- predict(cart, test, na.rm = TRUE)
  # assign class to maximum probablity
  pclass <- apply(pclass, 1, which.max)
  # create a data.frame using the reference and prediction
  x[[k]] <- cbind(test$class, as.integer(pclass))
}
class <- c("Territorios Artificializados", "Territorios Agrícolas", "Bosques y Áreas Seminaturales", "Áreas Húmedas", "Superficies de Agua")
classdf <- data.frame(classvalue = seq(1:5), classnames = class)

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
# confusion matrix
conmat <- table(y)
# 
print(conmat)
##         predicted
## observed   1   2   3   4   5
##        1 309  49  40  99   3
##        2  46 264 121  68   1
##        3  27  73 369  31   0
##        4  39 121  88 215  37
##        5   5  52  56  51 336
n <- sum(conmat)
diag <- diag(conmat)

# OVERALL
OAMNDWI <- sum(diag)/n
OAMNDWI 
## [1] 0.5972
#KAPPA
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
# predicted cases per class
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
KMNDWI <- (OAMNDWI - expAccuracy) / (1 - expAccuracy)
KMNDWI
## [1] 0.4965
#MNDWI 10 M
set.seed(63)
# number of folds
k <- 5
j <- sample(rep(1:k, each = round(nrow(SDatosMNDWIB))/k))
table(j)
## j
##   1   2   3   4   5 
## 500 500 500 500 500
x <- list()
for (k in 1:5) {
  train <- SDatosMNDWIB[j!= k, ]
  test <- SDatosMNDWIB[j == k, ]
  cart <- rpart(as.factor(class)~., data=train, method = 'class',
                minsplit = 5)
  pclass <- predict(cart, test, na.rm = TRUE)
  # assign class to maximum probablity
  pclass <- apply(pclass, 1, which.max)
  # create a data.frame using the reference and prediction
  x[[k]] <- cbind(test$class, as.integer(pclass))
}
class <- c("Territorios Artificializados", "Territorios Agrícolas", "Bosques y Áreas Seminaturales", "Áreas Húmedas", "Superficies de Agua")
classdf <- data.frame(classvalue = seq(1:5), classnames = class)

y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
# confusion matrix
conmat <- table(y)
# 
print(conmat)
##         predicted
## observed   1   2   3   4   5
##        1 281  80  45  91   3
##        2  25 293 120  62   0
##        3  27  57 369  43   4
##        4  53 143  86 192  26
##        5  19  52  47  52 330
n <- sum(conmat)
diag <- diag(conmat)

# OVERALL
OAMNDWIB <- sum(diag)/n
OAMNDWIB 
## [1] 0.586
#KAPPA
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
# predicted cases per class
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
KMNDWIB <- (OAMNDWIB - expAccuracy) / (1 - expAccuracy)
KMNDWIB
## [1] 0.4825

Evaluación Cualitativa

#NDWI
ClasificacionNDWIA1 <- crop(ClasificacionNDWI, A1)
ClasificacionNDWIA2 <- crop(ClasificacionNDWI, A2)
ClasificacionNDWIA3 <- crop(ClasificacionNDWI, A3)
ClasificacionNDWIA4 <- crop(ClasificacionNDWI, A4)
ClasificacionNDWIA5 <- crop(ClasificacionNDWI, A5)
ClasificacionNDWIA6 <- crop(ClasificacionNDWI, A6)

CNDWI1 <- app(ClasificacionNDWIA1, fun = which.max)
CNDWI2 <- app(ClasificacionNDWIA2, fun = which.max)
CNDWI3 <- app(ClasificacionNDWIA3, fun = which.max)
CNDWI4 <- app(ClasificacionNDWIA4, fun = which.max)
CNDWI5 <- app(ClasificacionNDWIA5, fun = which.max)
CNDWI6 <- app(ClasificacionNDWIA6, fun = which.max)
#writeRaster(CNDWI, filename="./RESULTADOS/CNDWI.TIF", datatype="INT1U", overwrite=TRUE)
#unique(CNDWI$lyr.1)

#MNDWI 20 M
ClasificacionMNDWI1 <- crop(ClasificacionMNDWI, A1)
ClasificacionMNDWI2 <- crop(ClasificacionMNDWI, A2)
ClasificacionMNDWI3 <- crop(ClasificacionMNDWI, A3)
ClasificacionMNDWI4 <- crop(ClasificacionMNDWI, A4)
ClasificacionMNDWI5 <- crop(ClasificacionMNDWI, A5)
ClasificacionMNDWI6 <- crop(ClasificacionMNDWI, A6)

CMNDWI1 <- app(ClasificacionMNDWI1, fun = which.max)
CMNDWI2 <- app(ClasificacionMNDWI2, fun = which.max)
CMNDWI3 <- app(ClasificacionMNDWI3, fun = which.max)
CMNDWI4 <- app(ClasificacionMNDWI4, fun = which.max)
CMNDWI5 <- app(ClasificacionMNDWI5, fun = which.max)
CMNDWI6 <- app(ClasificacionMNDWI6, fun = which.max)
#writeRaster(CMNDWI1, filename="./RESULTADOS/CMNDWI1.TIF", datatype="INT1U", overwrite=TRUE)
#unique(CMNDWI1$lyr.1)

#MNDWI 20 M
ClasificacionMNDWIB1 <- crop(ClasificacionMNDWIB, A1)
ClasificacionMNDWIB2 <- crop(ClasificacionMNDWIB, A2)
ClasificacionMNDWIB3 <- crop(ClasificacionMNDWIB, A3)
ClasificacionMNDWIB4 <- crop(ClasificacionMNDWIB, A4)
ClasificacionMNDWIB5 <- crop(ClasificacionMNDWIB, A5)
ClasificacionMNDWIB6 <- crop(ClasificacionMNDWIB, A6)

CMNDWIB1 <- app(ClasificacionMNDWIB1, fun = which.max)
CMNDWIB2 <- app(ClasificacionMNDWIB2, fun = which.max)
CMNDWIB3 <- app(ClasificacionMNDWIB3, fun = which.max)
CMNDWIB4 <- app(ClasificacionMNDWIB4, fun = which.max)
CMNDWIB5 <- app(ClasificacionMNDWIB5, fun = which.max)
CMNDWIB6 <- app(ClasificacionMNDWIB6, fun = which.max)
#writeRaster(CMNDWIB1, filename="./RESULTADOS/CMNDWIB1.TIF", datatype="INT1U", overwrite=TRUE)
#unique(CMNDWIB1$lyr.1)
ColoresV <- c('#fc4e2a','#f7fcb9' , '#99d8c9','#006837','#41ae76', '#d9f0a3','#addd8e', '#2b8cbe','blue') %>%
colorRampPalette()

#NDWI
par(mfrow = c(3,2))
plot(CNDWI1,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CNDWI2,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CNDWI3,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CNDWI4,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CNDWI5,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CNDWI6,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))

#MNDWI 20 M
par(mfrow = c(3,2))
plot(CMNDWI1,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWI2,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWI3,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWI4,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWI5,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWI6,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))

#MNDWI 10 M
par(mfrow = c(3,2))
plot(CMNDWIB1,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWIB2,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWIB3,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWIB4,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWIB5,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))
plot(CMNDWIB6,  legend=TRUE, axes=TRUE,  maxcell=5000000, col=ColoresV(9))