Sentinel2-intro.R

Sentinel2-intro.R

## Working with Sentinel-2 bands in R


## Primeros pasos

# - Crear proyecto de RStudio
#
# - Guardar capas (ZIP) en la misma carpeta


## Cargar paquetes necesarios

library(raster)

# La primera vez habrá que instalarlo

## Descomprimir ZIP

#unzip("S2B_MSIL2A_20191013T115219_N0213_R123_T28RDS_20191013T153823.zip")

## Listar capas disponibles

azul <- raster("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/T28RDS_20191013T115219_B02_20m.jp2")

azul
plot(azul)








bandas <- 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)



## Cargar capas en R

bandas <- stack(bandas)



## Cargar capas en R

bandas



## Si esto no funciona, probablemente te falte el driver 'JP2OpenJPEG'

# - Instala GDAL usando el [instalador de OSGEO](https://trac.osgeo.org/osgeo4w/)
#
# - Usa `gdalUtils::gdal_chooseInstallation("JP2OpenJPEG")`
#
# - Usa `gdalUtils::gdal_translate()` para convertir de jp2 a GeoTiff
#
# - Cargar GeoTiff


## Mapa rápido

plot(bandas, 2)

names(bandas)


## Zoom rápido

zoom(bandas)


## Imagen visible

plotRGB(bandas, r = 12, g = 13, b = 14)



## Recortar zona del incendio

incendio <- crop(bandas, c(426166, 446166, 3097108, 3115000))
incendio


## Mapa

plot(incendio, 8)




# Máscara de nubes y sombras


## Cargar raster de nubes

nubes <- raster("C:/Users/FRS/Dropbox/AGRESTA/Rprojects/Remote-sensing-with-R-intro/S2B_MSIL2A_20191013T115219_N0213_R123_T28RDS_20191013T153823.SAFE/GRANULE/L2A_T28RDS_A013591_20191013T115219/QI_DATA/MSK_CLDPRB_20m.jp2")


## Cargar raster de nubes

nubes



## Mapa de nubes

plot(nubes)



## Recortar a zona del incendio

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



## Eliminar zonas donde Prob. Nubes >20%

nubes[nubes > 20] <- NA
incendio <- mask(incendio, nubes)

writeRaster(nubes, "nubes.tif")

## Eliminar zonas donde Prob. Nubes >20%

plot(incendio, 1)



## Sentinel 2a incluye máscara de nubes y sombras

names(incendio)

#Banda 11: SCL




## Sentinel 2a incluye máscara de nubes y sombras

#knitr::include_graphics("images/SC.png")

#https://earth.esa.int/web/sentinel/technical-guides/sentinel-2-msi/level-2a/algorithm


## Sentinel 2a incluye máscara de nubes y sombras

freq(incendio[[11]])




## Eliminar nubes, sombras...

incendio[incendio[[11]] < 4 | incendio[[11]] > 6] <- NA



## Eliminar nubes, sombras...

plot(incendio, 1)






# Calcular índices

# Normalised Difference Vegetation Index (NDVI)

## Calcular NDVI

# $$
# NDVI = \frac {NIR - Red} {NIR + Red}
# $$

#Valores de -1 (agua) a +1 (bosques)

#https://www.sentinel-hub.com/eoproducts/ndvi-normalized-difference-vegetation-index


## Bandas Sentinel-2

#knitr::include_graphics("images/bands10m.jpg")



## Bandas Sentinel-2

names(incendio)



## Calcular NDVI

# $$
# NDVI = \frac {NIR (B8) - Red(B4)} {NIR(B8) + Red(B4)}
# $$

## Calcular NDVI

ndvi <- (incendio[[10]] - incendio[[4]]) / (incendio[[10]] + incendio[[4]])


## NDVI

ndvi



## NDVI

hist(ndvi, main = "NDVI")



## NDVI

plot(ndvi)



## NDVI

plot(ndvi, col = viridis::viridis(9))




# Normalised Burnt Ratio (NBR)


## Normalised Burnt Ratio (NBR)

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

## Normalised Burnt Ratio (NBR)

#knitr::include_graphics("images/bands20m.jpg")

#https://www.earthdatascience.org/courses/earth-analytics/multispectral-remote-sensing-modis/normalized-burn-index-dNBR/


## Calcular NBR

names(incendio)



## Calcular NBR

nbr <- (incendio[[10]] - incendio[[9]]) / (incendio[[10]] + incendio[[9]])



## NBR

nbr



## NBR

plot(nbr)


## NBR

hist(nbr, main = "NBR")




## Guardar rasters

writeRaster(ndvi, "ndvi.tif")
writeRaster(nbr, "nbr.tif")




# Ejercicios

## Calcular NBR pre-incendio y diferencia pre-post incendio

# - Valores negativos (< -0.1): Buen crecimiento
#
# - Valores positivos: severidad creciente



#################################################


# 1. Cargar capas de Julio usando stack

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

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

names(julio)

names(julio) <- c("B2", "B3", "")

# 2. Calcular NBR Julio
# 3. Calcular la diferencia en NBR.Oct - NBR.Julio
# 4. Plotear el raster de diferencia