###Detección de cuerpos de agua a partir de imágenes de Sentinel-2 con NDWI y MNDWI utilizando la banda SWIR modificada a una resolución espacial de 10 m producida mediante pan-sharpening

El uso de imágenes satelitales en percepción remota para el estudio de las dinámicas de la superficie terrestre se ha convertido en un elemento fundamental debido a su extensión de aplicaciones. Las imágenes Sentinel-2 son imágenes ópticas de alta resolución espacial y espectral respecto a otros satélites y debido a la disponibilidad y acceso a los datos que produce este satélite se ha incrementado su uso en las aplicaciones e investigaciones de percepción remota. De igual manera, en percepción remota también existen limitaciones como la compensación en las características de los sensores por lo que fusión de datos de múltiples fuentes en percepción remota ha obtenido una gran atención e implementación en la actualidad. En el presente informe se calcularon los índices espectrales de agua NDWI y MNDWI, con una limitación en la resolución espacial de debido a que la banda SWIR en Sentinel-2 es de 20 m, por lo que se ejecutó una metodología que modifica la resolución espacial de la banda SWIR con fusión de datos con enfoque pan-sharpening y se calculan índices MNDWI con la banda SWIR de 10 m.

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.2.1, released 2020/12/29
## Path to GDAL shared files: D:/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: D:/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## 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.
## Overwritten PROJ_LIB was D:/Documents/R/win-library/4.0/rgdal/proj
## 
## 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
#DIRECTORIO DE TRABAJO
setwd('D:/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)

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.

#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: "D:\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: "D:\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: "D:\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: "D:\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: "D:\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: "D:\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")

#/////////////////////////////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)