Sentinel-2 repaso

Sentinel2-repaso.R

## Working with Sentinel-2 bands in R


# Crear proyecto de RStudio


## Cargar paquetes necesarios

library(raster)


## Cargar capas en R (Octubre y Julio)

bandas.Oct <- list.files("C:/Users/FRS/Dropbox/AGRESTA/Rprojects/Remote-sensing-with-R-intro/S2B_MSIL2A_20191013T115219_N0213_R123_T28RDS_20191013T153823.SAFE/GRANULE/L2A_T28RDS_A013591_20191013T115219/IMG_DATA/R20m/", full.names = TRUE)

octubre <- stack(bandas.Oct)
octubre
plot(octubre, 1)
plotRGB(octubre, r = 12, g = 13, b = 14)

julio <- stack("C:/Users/FRS/Dropbox/AGRESTA/CursoTeledeteccion/Dia2/imagenes/pila_20m_20190725.tif")
julio

# Renombrar capas Julio

names(julio) <- c("B2", "B3", "B4", "B5", "B6", "B7", "B8", "B11", "B12")
names(julio)



## Recortar zona del incendio

octubre <- crop(octubre, c(426166, 446166, 3097108, 3115000))
plot(octubre, 2)

julio <- crop(julio, c(426166, 446166, 3097108, 3115000))
plot(julio, 1)




## Normalised Burnt Ratio (NBR)

# $$
# NBR = \frac {NIR (B8) - SWIR(B12)} {NIR(B8) + SWIR(B12)}
# $$



## Calcular NBR

names(octubre)

nbr.octubre <- (octubre[[10]] - octubre[[9]]) / (octubre[[10]] + octubre[[9]])

plot(nbr.octubre)

names(julio)

nbr.julio <- (julio[[7]] - julio[[9]]) / (julio[[7]] + julio[[9]])

plot(nbr.julio)


# Histogramas

hist(nbr.julio, main = "NBR Julio")

hist(nbr.octubre, main = "NBR Octubre")


# Calcular diferencia NBR
nbr.dif <- nbr.julio - nbr.octubre

hist(nbr.dif)

plot(nbr.dif)

plot(nbr.dif, col = tinter::tinter("Gray 10"))


## Guardar raster para visualizar en QGIS

writeRaster(nbr.dif, "nbrdif.tif")