Process Sentinel2 reflectance data from Theia

S2_reflectance.R

## CALCUL NDVI avec masque nuages

## Image couleur naturel Sentinel2


library(raster)
library(rgdal)
library(stringr)

  

# repertoire de base : contient les images 
d_base <- "F:/Travail/Sentinel2/theia_download-master/T31TEJ_MONTPELLIER"
# extraire repertoire des images
d_idprod <- list.dirs(d_base, full.names=T, recursive=F)
d_idprod <- str_subset(d_idprod, "SENTINEL2[AB]")
# extraire date des images
dates_prod <- str_extract(d_idprod, "[:digit:]{8}")
sat_prod <- str_extract(d_idprod, "SENTINEL2[AB]")
# extraire chemin fichiers TIF : nuages, B2, B3, B4, B8
fun_data_path <- function(depart, motif) {
  f <- list.files(depart, motif, full.names=T, recursive=T)
}
f_prod_clm <- sapply(d_idprod, fun_data_path, "_CLM_R1.tif$", USE.NAMES=F)
f_prod_b4 <- sapply(d_idprod, fun_data_path, "_FRE_B4.tif$", USE.NAMES=F)
f_prod_b3 <- sapply(d_idprod, fun_data_path, "_FRE_B3.tif$", USE.NAMES=F)
f_prod_b2 <- sapply(d_idprod, fun_data_path, "_FRE_B2.tif$", USE.NAMES=F)
f_prod_b8 <- sapply(d_idprod, fun_data_path, "_FRE_B8.tif$", USE.NAMES=F)

# systeme de coordonnees de l'image Sentinel2 : UTM 31
#pcs_utm31 <- "+init=EPSG:32631"
# etendue a extraire
extcrop <- extent(c(549000, 559000, 4829000, 4844100))


# #######################################################################
# partie 2 : calcul NDVI

# 
fun_image_crop <- function(f_prod_r, f_prod_b, f_prod_g, f_prod_nir, 
                           ext_clip, f_reflect) {
  lf <- list(B4=f_prod_r, B3=f_prod_b, B2=f_prod_g, B8=f_prod_nir)
  stack_reflect <- stack(lf, quick=T)
  stack_reflect2 <- crop(stack_reflect, ext_clip, filename=f_reflect, 
                        datatype="INT2S", NAflag=-9999, overwrite=T)
  return(stack_reflect2)
}

fun_cloudm_crop <- function(f_prod_clm, ext_clip, f_boolm) {
  rast_clm <- raster(f_prod_clm)
  rast_clm2 <- (rast_clm > 0)
  rast_clm2 <- crop(rast_clm2, ext_clip, filename=f_boolm, 
                         datatype="INT1U", NAflag=99, overwrite=T)
  return(rast_clm2)
}

# repertoire de sortie
d_output <- paste0(d_base, "/OUTPUT")
dir.create(d_output)
# chemin rasters sortie
f_mask <- paste0(d_output, "/MASK_", sat_prod, "_", dates_prod, ".tif")
f_ndvi_mask <- paste0(d_output, "/NDVI_MASK_", sat_prod, "_", dates_prod, ".tif")
f_reflectance <- paste0(d_output, "/S2R_", sat_prod, "_", dates_prod, ".tif")

for (X in 1:length(dates_prod)) {
  stack_reflect <- fun_image_crop(f_prod_b4[X], f_prod_b3[X], f_prod_b2[X], f_prod_b8[X], 
                                  extcrop, f_reflectance[X])
  rast_m <- fun_cloudm_crop(f_prod_clm[X], extcrop, f_mask[X])
  rast_ndvi <- calc(stack_reflect, 
                    function(v){(v[4]-v[1])/(v[4]+v[1])})
  rast_ndvi_mask <- mask(rast_ndvi, rast_m, 
                         filename=f_ndvi_mask[X], maskvalue=1, overwrite=T)
}



# #######################################################################
# partie 4 : plot couleurs naturelles (RGB)
fun_plot_natcol <- function(f_reflect, f_plot, titre, ext_plot) {
  brick_reflect <- brick(f_reflect)
  png(f_plot)
  plotRGB(brick_reflect, stretch="lin", ext=ext_plot, main=titre)
  dev.off()
}

f_plot <- paste0(d_output, "/Plot_COULNAT_", dates_prod, ".png")
titres <- paste("RGB",dates_prod)
for (X in 1:length(dates_prod)) {
  fun_plot_natcol(f_reflectance[X], f_plot[X], titres[X], extcrop)
}

# #######################################################################
# partie 4 : plot NDVI
library(RColorBrewer)
pal_YlGr <- brewer.pal(5,"YlGn")
lim <- c(-1, 0, 0.2, 0.5, 0.8, 1)
titres <- paste("NDVI",dates_prod)

func_plot_NDVI <- function(f_ndvi_mask, f_plot, titre, ext_plot, pal, brks) {
  rast_ndvi_mask_X <- raster(f_ndvi_mask)
  png(f_plot)
  plot(rast_ndvi_mask_X, ext=ext_plot, col=pal, colNA=8, breaks=brks, main=titre)
  dev.off()
}

# chemin rasters sortie
f_plot <- paste0(d_output, "/Plot_NDVI_", dates_prod, ".png")

for (X in 1:length(dates_prod)) {
  func_plot_NDVI(f_ndvi_mask[X], f_plot[X], titres[X], extcrop, pal_YlGr, lim)
}