# 8. Cloud Masking and Mosaicking

# Project set up
# Only if needed:
# install.packages(c("terra", ggplot2))

#Sets the working directory to the folder you have your data in
setwd("D:/ACADEMIA/TRAINING/ECOCap4Africa/QGIS")

library(terra)


# Get list of all .jp2 band files with cloud
all_sentinel10mbands <- list.files("S2B_MSIL2A_20250110T102319_N0511_R065_T30NYN_20250110T141704.SAFE/GRANULE/L2A_T30NYN_A040989_20250110T103111/IMG_DATA/R10m",
                             pattern = glob2rx("*B*.jp2$"), full.names = T)
# Load the SWIR band
Band11 <- rast("S2B_MSIL2A_20250110T102319_N0511_R065_T30NYN_20250110T141704.SAFE/GRANULE/L2A_T30NYN_A040989_20250110T103111/IMG_DATA/R20m/T30NYN_20250110T102319_B11_20m.jp2")

# Open sentinel bands (Now you have a list of all  the Sentinel 20m bands in your folder. You could choose to open file individually)
# load blue band
blue_band <-rast(all_sentinel10mbands[2])
plot(blue_band, main= "Sentinel 20m Blue Band")

# It's more efficient to open all of the bands together as a stack 
# stack all 10m sentinel bands
sentinel10m <- rast(all_sentinel10mbands)

#select needed bands
#bandselect <- all_sentinel10mbands[2:4]
#bandselect <- all_sentinel20mbands[c(2:5,10)]

#Resample the SWIR band_11 using the 10m raster file
Rs_band11 <- resample(Band11, sentinel10m, method="near")

#Combine the band 11 to the four 10m Sentinel-2 bands
sentinel10_add <- c(sentinel10m,Rs_band11)

#Plot and visualize your data
plotRGB(sentinel10_add, r=4, g=3, b=2, axes=T, stretch='lin')


### Get list of all .jp2 files in QI_data folder
maskbands <- list.files("S2B_MSIL2A_20250110T102319_N0511_R065_T30NYN_20250110T141704.SAFE/GRANULE/L2A_T30NYN_A040989_20250110T103111/QI_DATA",
                                   pattern = glob2rx("*B*.jp2$"), full.names = T)

#Select the 20m cloud mask layer
maskband_20m <- rast(maskbands[2])
maskband_20m

#Resample the cloud mask band to 10m
maskband_rs <- resample(maskband_20m,sentinel10_add, method="near")
maskband_rs
plot(maskband_rs)

#?mask
#Mask the couds from the 10m sentinel composite bands using the resampled 10m cloud mask layer
sentinel_masked <- mask(sentinel10_add,maskband_rs,maskvalues=100)
#Plot and visualize your data
plotRGB(sentinel_masked, r=4, g=3, b=2, axes=T, stretch='lin')

# Write the raster to your working directory and open it in GQIS for better visualization and inspection
writeRaster(sentinel_masked, filename = "MaskedSentinel.tif", datatype='INT4S', overwrite =TRUE)


# Get list of all .jp2 band files without cloud
T30NYN <- list.files("S2A_MSIL2A_20250105T102421_N0511_R065_T30NYN_20250105T143159.SAFE/GRANULE/L2A_T30NYN_A049826_20250105T102729/IMG_DATA/R10m",
                                   pattern = glob2rx("*B*.jp2$"), full.names = T)
T30NYN_raster <- rast(T30NYN)

# Load the SWIR band
b11 <- rast("S2A_MSIL2A_20250105T102421_N0511_R065_T30NYN_20250105T143159.SAFE/GRANULE/L2A_T30NYN_A049826_20250105T102729/IMG_DATA/R20m/T30NYN_20250105T102421_B11_20m.jp2")

#Resample b11 to 10m
b11_re <- resample(b11,T30NYN_raster, method="near")

#Add resmpled band 11 to sentinel 10m bands
T30NYN_add <- c(T30NYN_raster,b11_re)

plotRGB(T30NYN_add, r=4, g=3, b=2, axes=T, stretch='lin')

# Mosaicking two files having the same extent to fill in holes in Masked layer
?mosaic
b2mosaic <- mosaic(T30NYN_add,sentinel_masked,fun=mean)
plotRGB(b2mosaic, r=4, g=3, b=2, axes=T, stretch='lin')



# Get list of all .jp2 band files from an adjecent tile
T30NXN <- list.files("S2A_MSIL2A_20250105T102421_N0511_R065_T30NXN_20250105T143159.SAFE/GRANULE/L2A_T30NXN_A049826_20250105T102729/IMG_DATA/R10m",
                     pattern = glob2rx("*B*.jp2$"), full.names = T)
T30NXN
T30NXN_raster <- rast(T30NXN)

# Load the SWIR band
b11 <- rast("S2A_MSIL2A_20250105T102421_N0511_R065_T30NXN_20250105T143159.SAFE/GRANULE/L2A_T30NXN_A049826_20250105T102729/IMG_DATA/R20m/T30NXN_20250105T102421_B11_20m.jp2")

#Resample b11 to 10m
b11_re <- resample(b11,T30NXN_raster, method="near")

#Add resmpled band 11 to sentinel 10m bands
T30NXN_add <- c(T30NXN_raster,b11_re)

#Mosaick the two different tiles for a wider coverage
all_tiles <- mosaic(T30NXN_add,b2mosaic, fun=mean)
plotRGB(all_tiles, r=4, g=3, b=2, axes=T, stretch='lin')
