PROGYSAT#1: download S2 data

Jean-Baptiste Féret

2023-12-05

Source: vignettes/biodivMapR_11.Rmd

This tutorial illustrates how to download and preprocess Sentinel-2 data for the application of biodivMapR in order to map spectral diversity over a site in Suriname.

This first step is based on the R package preprocS2, which needs to be installed, along with additional packages such as sen2r

This includes:

  • S2 data download based on a spatial footprint

  • atmospheric correction of S2 images to produce L2A images with Sen2Cor

  • cropping, stacking and writing the reflectance raster corresponding to the spatial footprint, as well as the corresponding cloud mask

Define structure for data, R scripts and results

This tutorial will assume that the file structure is defined as follows:

  • 01_DATA: a directory including original vector and raster files: Sentinel-2 data products level L1C and L2A before further process (full tile)

  • 02_PROGRAMS: a directory including all r scripts should be located. All R scripts provided in this tutorial will assume that this is your working directory

  • 03_RESULTS: a directory including files resulting from the application of the R scripts on the data located in 01_DATA, or in this same 03_RESULTS directory

Download a vector file corresponding to the spatial footprint of the study area

You can use QGIS, Google Earth, or any relevant GIS or R package to define the spatial footprint of your area of interest. Save it as kml file or shapefile in 01_DATA.

Here, we drew a rectangular polygon over a forested area in the Eastern part of the Brokopondo reservoir in Suriname.

First, let’s download the kml file corresponding to our area of interest from my gitlab repository, and save it locally in the 01_DATA directory.

# define data directory
Path_Data <- '../01_DATA'
dir.create(Path_Data,showWarnings = F, recursive = T)
# name zip file including plots located on the tile
destkml <- file.path(Path_Data,'Suriname_Progysat_Workshop.kml')
# url for the zip file
url <- 'https://gitlab.com/jbferet/myshareddata/-/raw/master/PROGYSAT/01_DATA/Suriname_Progysat_Workshop.kml'
download.file(url = url, destfile = destkml)

The study area corresponding to the Eastern part of the Brokopondo reservoir is represented by the red square in the figure below.

 

 

Download Sentinel-2 data and perform atmospheric corrections with sen2Cor

Then, we use Sentinel-hub in order to identify an acquisition corresponding to our requirement.

Once the S2 tile and the date of acquisition are identified, we can proceed to the image download. Here, the acquisition is from October 16th, 2020 showed was one of the cleanest image available from the Sentinel-2 archive.

 

 

We use the R package preprocS2 to download the corresponding Sentinel-2 SAFE archive. Here, level-2A images (i.e. atmospherically corrected with Sen2Cor) are not available, and only Level-1C images can be downloaded from the Copernicus Hub.

The function get_S2_L2A_Image allows to automatically download S2 images (from Copernicus Hub, Google Cloud or Google Earth Engine when properly parameterized) and to perform atmospheric corrections with Sen2Cor if only L1C data is available.

################################################################################
## S2 download & atmospheric correction using preprocS2
################################################################################
library(preprocS2)
# define date of S2 acquisition
dateAcq <- '2020-10-16'
# define path for study area
path_vector <- destkml
# define output directory where SAFE L2C and L2A files are stored. L1C is 
# automatically deleted after production of L2A to save space
DirWrite <- '../01_DATA/S2_Images'
dir.create(DirWrite,showWarnings = F, recursive = T)

Path_S2 <- get_S2_L2A_Image(l2a_path = DirWrite, 
                            spatial_extent = path_vector, 
                            dateAcq = dateAcq,
                            DeleteL1C = TRUE, 
                            Sen2Cor = TRUE)

Save cropped & stacked BOA reflectance and corresponding cloud mask

Sentinel-2 products are made available in the SAFE format, including image data in JPEG2000 format, quality indicators (e.g. defective pixels mask), auxiliary data and metadata.

The preprocS2 function extract_from_S2_L2A aims at extracting all information required for further processing, including reflectance data, cloud mask and metadata, and rearrange it in a more user-friendly fashion (personal opinion). It also crops S2 data following a spatial extent defined by path_vector, and it harmonizes the spatial resolution to either 20 meters, or 10 meters using user defined interpolation (bilinear as default, change it with the input variable interpolation set to the method of your choice available in GDAL).

If the vector file and the raster data do not share the same projection, the vector file is also reprojected to avoid possible errors in future processings.

################################################################################
## write a stacked reflectance raster corresponding to the study area
################################################################################
# Result directory where data is stored 
result_path <- '../03_RESULTS'
dir.create(path = result_path,showWarnings = FALSE,recursive = TRUE)

##____________________________________________________________________________##
##      Extract, resample all bands to 10m spatial resolution & stack data    ##
##----------------------------------------------------------------------------##
# define spatial resolution
resolution <- 10
# define source of data
S2source <- 'SAFE'
S2obj <- preprocS2::extract_from_S2_L2A(Path_dir_S2 = Path_S2,
                                        path_vector = path_vector,
                                        S2source = S2source,
                                        resolution = resolution)

# update vector file if needed (reprojection and write as shapefile in the same 
# location as the original file)
path_vector <- S2obj$path_vector

A stars object is obtained as output from extract_from_S2_L2A. This is a proxy stars object, in order to handle large rasters. This stars object includes the reflectance data and corresponding cloud mask.

These will be written on your local space with the functions save_cloud_s2 and save_reflectance_s2

# create specific result directory corresponding to granule name
results_site_path <- file.path(result_path,basename(S2obj$S2_Bands$GRANULE))
dir.create(path = results_site_path,showWarnings = FALSE,recursive = TRUE)

## Write CLOUD MASK
cloudmasks <- preprocS2::save_cloud_s2(S2_stars = S2obj$S2_Stack,
                                       Cloud_path = results_site_path,
                                       S2source = S2source, SaveRaw = T)

## Write REFLECTANCE
# Save Reflectance file as ENVI image with BIL interleaves
# metadata files are also important to account for offset applied on S2 L2A products 
Refl_path <- file.path(results_site_path,paste(basename(S2obj$S2_Bands$GRANULE),'_Refl',sep = ''))
tile_S2 <- get_tile(S2obj$S2_Bands$GRANULE)
dateAcq_S2 <- get_date(S2obj$S2_Bands$GRANULE)
preprocS2::save_reflectance_s2(S2_stars = S2obj$S2_Stack, 
                               Refl_path = Refl_path,
                               S2Sat = NULL, 
                               tile_S2 = tile_S2, 
                               dateAcq_S2 = dateAcq_S2,
                               Format = 'ENVI', 
                               datatype = 'Int16', 
                               MTD = S2obj$S2_Bands$metadata, 
                               MTD_MSI = S2obj$S2_Bands$metadata_MSI)

The cropped image and binary mask obtained from this preprocessing should correspond to the figure below:

 

 

Once S2 acquisition and corresponding cloud mask are downloaded and written on your disk, you can then run biodivMapR using the modality of your choice:

  • you can run biodivMapR in standard mode using directly the raster corresponding to BOA reflectance data, and applying the original workflow: radiometric filtering -> spectral transformation & component selection -> spectral species mapping -> diversity mapping

  • you can run biodivMapR in custom mode using user-defined custom layers, such as spectral indices, biophysical variables, or a combination of relevant features derived from multiple sensors for a customized workflow: [optional radiometric filtering -> optional spectral transformation & component selection -> ] spectral species mapping -> diversity mapping. Then, the selection of the relevant features used to compute remotely sensed variability is at the user’s discretion.