SENTINEL

Sentinel 2

Sentinel-2 is an Earth observation mission developed by ESA as part of the Copernicus Programme to perform terrestrial observations in support of services such as forest monitoring, land cover changes detection, and natural disaster management. It consists of two identical satellites, Sentinel-2A and Sentinel-2B. Details about the mission, user and technical guides as well as other documentation can be found at the ESA dedicated website [1] and a number of other sites.

Overview of Sentinel-2 mission capabilities

Sentinel-2 mission has the following features:

• Multi-spectral data with 13 bands in the visible, near infrared (NIR), and short wave infrared (SWIR).
• Systematic global coverage of land surfaces from 56° S to 84° N, coastal waters, and all of the Mediterranean Sea
• Revisiting every 5 days under the same viewing angles. At high latitudes, Sentinel-2 swath overlaps and some regions are observed twice or more every 5 days, but with different viewing angles.
• Spatial resolutions of 10 m (Bands 2, 3, 4 and 8 - Blue, Green, Red and NIR), 20 m (Bands 5, 6 and 7 - Red Edge bands, Band 8A - Narrow infrared and Bands 11 and 12 - SWIR) and 60 m (Bands 1 - Coastal aerosol, Band 9 - Water vapour and Band 10 - SWIR/Cirrus)
• 290 km field of view
• Free and open data policy

Sentinel-2 products

• Level-1B product provides radiometrically corrected imagery in Top-Of-Atmosphere (TOA) radiance values and in sensor geometry. Level-1B product is composed of an ensemble of granules that are 25 km across track (AC) by 23 km along track (AL). Additionally, this product includes the refined geometry which is used to generate the Level-1C product. Level-1B pixel coordinates refer to the centre of each pixel.

• Level-1C product is composed of 100 km² single tiles (ortho-images in UTM/WGS84 projection). The Level-1C product results from using a Digital Elevation Model (DEM) to project the image in cartographic coordinates. Per-pixel radiometric measurements are provided in TOA reflectances with all parameters to transform them into radiances. Level-1C products are resampled with a constant Ground Sampling Distance (GSD) of 10, 20 and 60 m depending on the native resolution of the different spectral bands. In Level-1C products, pixel coordinates refer to the upper left corner of the pixel. Before September 2016, Level-1C was delivered as multi-tile product.

• Level-2A product provides Bottom Of Atmosphere (BOA) reflectance images derived from the associated Level-1C products. Therefore, each Level-2A product is also a 100 km² tile in cartographic geometry (UTM/WGS84 projection). BOA reflectances in Level-2A products are scaled by a factor of 10000.

Data download

Sentinel-2 data might be obtained from the Copernicus Open Access Hub, either by using the online interactive interface and browsing to your area of interest or by means of the API Hub. There are several other sites from which it is possible to download Sentinel-2 imagery as well. Here, some of them:

• Sentinel-app
• Amazon S3
• Copernicus Data and Exploitation Platform – Deutschland (CODE-DE)
• mundialis
• RemotePixel
• Land Viewer

For those preferring or needing to download data via the API, there also exist several options:

• sentinelsat: Python utilities to access the API of Copernicus Sentinels Scientific Data Hub (useful for all Sentinel satellite imagery)
• Sentinel-download: Automated download of Sentinel-2 L1C data from ESA (through wget)

GRASS GIS has a dedicated set of AddOns, r.sentinel.download and r.sentinel.import, to download and import Sentinel data into GRASS GIS database, respectively. This is probably the most convenient and straight-forward way if working with GRASS. Note that r.sentinel.download is build on top of the sentinelsat library and it also requires python-pandas.

To install these addons, just run:

g.extension extension=r.sentinel

Data pre-processing

As from May 2017, ESA started to provide Level 2A products for Europe and they will reprocess the archive during this year (2017). However, it is (still) possible to download Level-1C products and customize the atmospheric correction by using different algorithms (i.e.: 6S, 6SV, arcsi) and/or a different DEM (Sen2Cor uses STRM 90m by default).

Sentinel-2 in GRASS GIS

To exemplify the handling of Sentinel-2 data with GRASS GIS, we will use S2 L1C scene from North Carolina (NC) from November 2016. Here are the details:

S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502
830f2a5a-9914-427c-b4e0-69c6c14c71d4
2016-11-10T15:05:02Z

In the zipped folder delivered by ESA, the data itself is hidden under:

S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502.SAFE/GRANULE/S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_N02.04/IMG_DATA/

S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B01.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B02.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B03.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B04.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B05.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B06.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B07.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B08.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B09.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B10.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B11.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B12.jp2
S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B8A.jp2

Importing data

To import Sentinel-2 data into GRASS GIS, the simplest is to use the r.sentinel.import addon. This nice module allows to list all files in a folder and import all the bands at once (with re-projection on the fly if needed), only import some of them by a pattern or just link the data through r.external. Here, are some examples on how to use r.sentinel.import:

# list all sentinel bands in a folder
r.sentinel.import -p input=S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502.SAFE
# import all sentinel bands within the folder
r.sentinel.import input=S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502.SAFE
# import only band 4 and 8
r.sentinel.import input=S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502.SAFE pattern='B0[4|8]'
# link files through r.external (-l flag)
r.sentinel.import -l imput=S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502.SAFE
# import all bands with re-projection on the fly
r.sentinel.import -r imput=S2A_OPER_PRD_MSIL1C_PDMC_20161111T123605_R097_V20161110T160502_20161110T160502.SAFE

Sentinel-2 data can be also easily imported into GRASS GIS by means of r.in.gdal (in a UTM location) or using r.import that provides re-projection on the fly. In this example, we will use r.import to re-project and import into the NC sample location delivered with GRASS GIS.

# list all jp2 files and import them
for map in `ls *.jp2` ; do
 outmap=`basename $map .jp2`
 r.import input=$map output=$outmap
done

If you prefer shorter map names, then it is possible to use g.rename or the addon g.rename.many. The other option is editing the file names before importing with sed for example. Let us keep the S2A, date and band number parts only:

for map in `ls *.jp2` ; do
 outmap1=`basename $map .jp2`
 outmap2=`echo $outmap1 | sed s/OPER_MSI_L1C_TL_MTI__//g | sed s/T210803_A007244_T17SPV//g`
 r.import input=$map output=$outmap2
done

The output resolution estimated by r.import is slightly different than the nominal resolution of Sentinel-2 bands.

r.info S2A_20161110_B08 -g | grep res
nsres=9.80178836938829
ewres=9.80224757766693

We need to set the proper resolution at import by means of the resolution and resolution_value options in r.import.

for i in 02 03 04 08 ; do
 g.message message="importing band ${i} with R 10m"
 r.import S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B${i}.jp2 out=S2A_20161110_B${i} resolution=value resolution_value=10 extent=input
done

for i in 05 06 07 8A 11 12 ; do 
  g.message message="importing band ${i} with R 20m"
  r.import S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B${i}.jp2 out=S2A_20161110_B${i} resolution=value resolution_value=20 extent=input
done

for i in 01 09 10 ; do 
 g.message message="importing band ${i} with R 60m"
 r.import S2A_OPER_MSI_L1C_TL_MTI__20161110T210803_A007244_T17SPV_B${i}.jp2 out=S2A_20161110_B${i} resolution=value resolution_value=60 extent=input
done

Unfortunately, nodata are not encoded in the Sentinel data and need to be set manually by replacing 0 with NULL. This should be done in the original projection, before re-projecting the data.

In L2A products, there are artifacts along edges of swath paths, cell values are too low, caused by a bug in sen2cor. All cells that are bordering a NULL cell must be set to NULL. This step might need to be repeated twice.

Because neighboring Sentinel tiles are overlapping each other, file sizes can be reduced by patching together neighboring tiles with the same datatake sensing start time.

Before the actual re-projection, the computational region must be set properly. In the original projection, four 10 meter cells fall exactly into one 20 meter cell, and nine 20 meter cells fall exactly into one 60 meter cell. This can be achieved by first setting the region to the input extents (obtained with r.proj -g), then aligning the region to a resolution of 60 meter, after that setting the region's resolution to the actual resolution of the input band, this time without alignment. Alternatively, a reference grid can be created and the region for the current input can then be aligned to that reference grid.

Now the data can be re-projected, preferably with method=bicubic_f or method=lanczos_f to avoid the sawtooth patterns of linear features which can appear with method=nearest (default).

Atmospheric correction with i.atcorr

To exemplify the atmospheric correction of Sentinel-2 data, we download the L1C scene S2A_OPER_PRD_MSIL1C_PDMC_20161029T092602_R054_V20161028T155402_20161028T155402 that overlaps completely with the elevation map in North Carolina location. This scene was imported with region cropping (see r.import) into PERMANENT mapset. The computational region was set to the extent of the elevation map. The atmospheric correction, as implemented by i.atcorr, is applied separately to each band of the scene and needs certain parameters to be set.

In the first step we create a file containing the 6S parameters for a particular scene and band. To create a 6S file, we need to obtain the following information:

• geometrical conditions,
• moth, day, decimal hours in GMT, decimal longitude and latitude of measurement,
• atmospheric model,
• aerosol model,
• visibility or aerosol optical depth,
• mean target elevation above sea level,
• sensor height and,
• sensor band.

Some of those details are encompassed in the metadata of the scene/bands and others are codified in different tables provided in the i.atcorr manual page. Let's go one by one:

Geometrical conditions

For Sentinel-2A, the geometrical conditions take the value 25 and for Sentinel-2B, the geometrical conditions value is 26 (See table A in i.atcorr manual page).

Day, time, longitude and latitude of measurement

Day and time of the measurement are hidden in the file name (i.e., the second datum in the file name with format YYYYMMDDTHHMMSS), and are also noted in the metadata file, which is included in the downloaded scene (file with .xml extension). Our sample scene was taken on October 28th (20161028) at 15:54:02 (155402). Note that the time has to be specified in decimal hours in Greenwich Mean Time (GMT). Luckily, the time in the scene name is in GMT and we can convert it to decimal hours as follows: 15 + 54/60 + 2/3600 = 15.901. Longitude and latitude refer to the centre of the computational region (which can be smaller than the scene), and must be in WGS84 decimal coordinates. To obtain the coordinates of the center, we can run:

g.region -bg

The longitude and latitude of the center are stored are ll_clon=-78.691 and ll_clat=35.749.

Atmospheric model

We can choose between various atmospheric models as defined at the beginning of this manual. For North Carolina, we can choose 2 - midlatitude summer.

Aerosol model

We can also choose between various aerosol models as defined at the beginning of this manual. For North Carolina, we can choose 1 - continental model.

Visibility or Aerosol Optical Depth

For Sentinel-2 scenes, the visibility is not measured, and therefore we have to estimate the aerosol optical depth instead, e.g. from AERONET. With a bit of luck, you can find a station nearby your location which measured the Aerosol Optical Depth at 500 nm at the same time as the scene was taken. In our case, on 28th October 2016, the EPA-Res_Triangle_Pk station measured AOD = 0.07 (approximately).

Mean target elevation above sea level

Mean target elevation above sea level refers to the mean elevation of the computational region. You can estimate it from the digital elevation model, e.g. by running:

r.univar -g elevation

In our case, mean=110. In the 6S file it will be displayed in [-km], i.e., -0.110.

Sensor height

Since the sensor is on board a satellite, the sensor height will be set to -1000.

Sensor band

The overview of satellite bands can be found in table F in i.atcorr manual page. For Sentinel-2A, the band numbers span from 166 to 178, and for Sentinel-2B, from 179 to 191.

Finally, here is what the 6S file would look like for Band 02 of our scene. In order to use it in the i.atcorr module, we can save it in a text file, for example params_B02.txt.

25
10 28 15.901 -78.691 35.749
2
1
0
0.07
-0.110
-1000
167

Now, we compute the atmospheric correction for the selected band B02 of our Sentinel 2 scene. We have to specify the following parameters:

• input = raster band to be processed,
• parameters = path to 6S file created in the previous step (we could also enter the values directly),
• output = name for the output corrected raster band,
• range = from 1 to the QUANTIFICATION_VALUE stored in the metadata file. It is 10000 for both Sentinel-2A and Sentinel-2B.
• rescale = the output range of values for the corrected bands. This is up to the user to choose, for example: 0-255, 0-1, 1-10000.

If the data is available, the following parameters can be specified as well:

• elevation = raster of digital elevation model,
• visibility = raster of visibility model.

Finally, this is the command would look like to apply atmospheric correction to band B02:

i.atcorr input=B02 parameters=params_B02.txt output=B02.atcorr range=1,10000 rescale=0,255 elevation=elevation

To apply atmospheric correction to the remaining bands, only the last line in the 6S parameters file (i.e., the sensor band) needs to be changed. The other parameters will remain the same.

Atmospheric correction with Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI)

The Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI) software was developed at Aberystwyth University, under contract from the Norwegian Space Centre and allows automatic correction to Top Of Atmosphere (TOA) or Surface Reflectance (SREF). The conversion to Surface Reflectance uses the 6S code, accessed through the Py6S Python interface. It was originally developed for Landsat data but it also supports Sentinel 2 satellites.

Aerosol Optical Thickness (AOT) is an important parameter in atmospheric correction models such as 6S. In some regions of the world it has been found that a constant can be used. However, for most regions a constant is not viable as the atmosphere is very variable. Incorrect parameterization for AOT can lead to large errors in atmospheric correction. ARCSI provides a method of deriving AOT from an image by using a dark object subtraction (DOS) to estimate surface reflectance in the blue channel. Taking this estimated reflectance as input, 6S is then numerically inverted to identify an AOT value which provides a surface reflectance value as close to that estimated from DOS as possible.

Instructions to install the software can be found at: https://spectraldifferences.wordpress.com/software/, but you can also have it through a docker image from: https://github.com/mundialis/docker-arcsi.

Here are some examples on how to use arcsi for Sentinel data. The first one only applies simple DOS1 correction and the second one demonstrates the application of the more advanced correction "DOSAOTSGL" which masks for clouds, cloud shadows and topographic shadows. The latter also derives AOT automatically.

# define name of Sentinel-2 scene - Note: omit .SAFE
S2IMG=S2A_MSIL1C_20170327T105021_N0204_R051_T31UFT_20170327T105021
DEM=srtm_30m_myregion.tif
OUTDIR=arcsi_output_AOT_inv
TMPDIR=~/tmp/arcsi

cd ${S2IMG}.SAFE/
mkdir ${TMPDIR}

# simple DOS1 correction example
arcsi.py --sensor sen2 -i MTD_MSIL1C.xml -o ${OUTDIR} \
	 --tmpath ${TMPDIR} -f KEA --stats -p RAD DOSAOTSGL SREF \
	 --aeroimg ${CONDA_PREFIX}/share/arcsi/WorldAerosolParams.kea \
	 --atmosimg ${CONDA_PREFIX}/share/arcsi/WorldAtmosphereParams.kea \
	 --dem ${DEM} --minaot 0.05 --maxaot 0.6 --simpledos

# DOSAOTSGL correction example
arcsi.py -s sen2 --stats --format KEA \
  -p CLOUDS DOSAOTSGL STDSREF SATURATE TOPOSHADOW FOOTPRINT METADATA SHARP \
  -o ${OUTDIR} --dem ${DEM} --tmpath ${TMPDIR} \
  --k clouds.kea meta.json sat.kea toposhad.kea valid.kea stdsref.kea \
  -i ${S2IMG}.SAFE/MTD_MSIL1C.xml

Where:

• -s or --sensor sen2 specifies the sensor, in the examples Sentinel-2
• -i MTD_MSIL1C.xml is the header file of the input image
• -o ${OUTDIR} specifies the output directory where all output files will be saved to
• --tmpath ${TMPDIR} is a directory where temporary files can be put during the processing; they will be deleted afterwards
• -f KEA specifies that the output image file format should be KEA; any GDAL supported format can be used
• --stats specifies that the output images should be populated with statistics and pyramids; makes display much faster but it is only available for KEA and HFA output formats
• -p or --prodlist RAD TOA SREF specifies the output products to be generated (RAD: radiance, TOA: top-of-atmosphere, SREF: surface reflectance)
• --aeroimg ../WorldAerosolParams.kea is an image file which is sampled to identify the aerosol model to be used for the input image (this can be downloaded from ARCSI downloads)
• --atmosimg ../WorldAtmosphereParams.kea is an image file which is sampled to identify the atmosphere model to be used for the input image (this can be downloaded from ARCSI downloads)
• --dem ${DEM} is an elevation model covering the scene. The higher the resolution of the DEM available the better
• --minaot 0.05 is the lower limit of the AOT values tested when estimating the AOT value to be used to correct the scene
• --maxaot 0.6 is the upper limit of the AOT values tested when estimating the AOT value to be used to correct the scene
• --simpledos specifies that a simple single valued DOS method should be used to provide an estimate of the Blue SREF used in the inversion

Further Sentinel-2 examples (taken from mundialis github):

# SENTINEL CLOUDS MASKING (ONLY)
arcsi.py --sensor sen2 -i ${S2IMG}.SAFE/MTD_MSIL1C.xml -o ${S2IMG}.SAFE/Clouds \
   --tmpath ${TMPDIR} -f KEA --stats -p CLOUDS

# SENTINEL CLOUDS MASKING AND ATMCOR with 6S, with lookup of atmosphere profile
arcsi.py --sensor sen2 -i ${S2IMG}.SAFE/MTD_MSIL1C.xml -o ${S2IMG}.SAFE/OutputsAOTInvCL \
   --tmpath ${TMPDIR} -f KEA --stats -p CLOUDS RAD DOSAOTSGL SREF \
   --aeroimg Documents/arcsi-1.4.2/data/WorldAerosolParams.kea \
   --atmosimg Documents/arcsi-1.4.2/data/WorldAtmosphereParams.kea \
   --dem ${S2IMG}.SAFE/dem_VR_all --minaot 0.05 --maxaot 0.6 --simpledos

# SENTINEL CLOUDS MASKING AND ATMCOR with 6S but with fixed AOT (already known)
arcsi.py --sensor sen2 -i ${S2IMG}.SAFE/MTD_MSIL1C.xml -o ${S2IMG}.SAFE/OutputsAOTInvCL \
   --tmpath ${TMPDIR} -f KEA --stats -p CLOUDS SREF \
   --aeroimg Documents/arcsi-1.4.2/data/WorldAerosolParams.kea \
   --atmosimg Documents/arcsi-1.4.2/data/WorldAtmosphereParams.kea \
   --dem ${S2IMG}.SAFE/dem_VR_all --aot 0.3

# SENTINEL, OLD NAME style from 2016; simple DOS example
arcsi.py --sensor sen2 -i ${S2IMG}.SAFE/S2A_OPER_MTD_SAFL1C_PDMC_20170119T125545_R097_V20161120T160552_20161120T160552.xml \
   -o ${S2IMG}.SAFE/OutputsAOTInv --tmpath ${TMPDIR} -f KEA --stats -p RAD DOSAOTSGL SREF \
   --aeroimg Documents/arcsi-1.4.2/data/WorldAerosolParams.kea \
   --atmosimg Documents/arcsi-1.4.2/data/WorldAtmosphereParams.kea \
   --dem ${S2IMG}.SAFE/srtm_21_05_utm17 --minaot 0.05 --maxaot 0.6 --simpledos

Estimate indices

Given that Sentinel-2 MSI comes with 3 red edge bands, it allows to estimate several other vegetation indices, aside from the most classic ones. There is a quite exhaustive list at the Sentinel Hub website. Here, only some of them:

• Normalized Difference Vegetation Index: NDVI = (B08 - B04) / (B08 + B04)
• Enhanced Vegetation Index: EVI = 2.5*(B08 - B04) / (B08 + 6*B04 - 7.5*B02 + 1)
• Simple Ratio Index: SR = B8 / B4
• Atmospherically Resistant Vegetation Index: ARVI = (B08 - (2*B04 - B02)) / (B08 + (2*B04 - B02))
• Soil Adjusted Vegetation Index: SAVI = 1.5 * (B08 - B04) / (B08 + B04 + 0.5)
• Red Edge NDVI: NDVI-Re = (B08 - B06) / (B08 + B06)
• Red Edge Chlorophyll Index: Chl-Re = B05 / B08
• Inverted Red Edge Chlorophyll Index: IReCI = (B07 - B04) * B06 / B05
• Green Normalized Difference Vegetation Index: GNDVI = (B08 - B03) / (B08 + B03)

The following is an example of how to calculate NDVI with a simple mapcalc expression (assuming that bands were previously atmospherically corrected and they were rescaled to 1-10000):

# define bands
B4=S2A_20161110_B04
B8=S2A_20161110_B08

# NDVI
r.mapcalc --o \
  expression="NDVI_B08_B04 = \
  if(${B8} <= 10000 && ${B4} <= 10000, \
  float(${B8} - ${B4}) / float(${B8} + ${B4}), null())"

Note that a long list of vegetation and water indices can be obtained with i.vi and i.wi, respectively.

Classification

Example here

Further reading

IMAGICO.DE blog

• First look at Sentinel data
• Who is moving the cheese again?