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

Manual import of Sentinel data is described below:

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 reprojecting the data.

In L2A products, there are artefacts 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 reprojection, 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 reprojected, preferably with method=bicubic_f or method=lanczos_f to avoid sawtooth patterns of linear features which can appear with method=nearest (default).

Atmospheric correction with i.atcorr

Example here

Alternative: Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI)

Estimate diverse 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.

• EVI = 2.5*(B08 - B04) / (B08 + 6*B04 - 7.5*B02 + 1)
• ARVI = (B08 - (2*B04 - B02)) / (B08 + (2*B04 - B02))
• NDII = (B08 - B11) / (B08 + B11)
• MSI = B11 / B08
• PSRI-NIR = (B04 - B02) / B08
• NDVI-GREEN = B03 * (B08 - B04) / (B08 + B04)
• MTCI = (B06 - B05) / (B05 - B04)
• IRECI = (B07 - B04) * B06 / B05
• GRVI1 = (B04 - B03) / (B04 + B03)
• GNDVI = (B08 - B03) / (B08 + B03)
• EVI2 = 2.5 * (B08 - B04) / (B08 + 2.4 * B04 + 1)
• ChlRE = B05 / B08

Classification

Example here

Further reading

IMAGICO.DE blog

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