Unit 21 - Sentinel spatio-temporal

Create a new mapset in jena-region location, eg. sentinel-ndvi (Settings ‣ GRASS working enviroment ‣ Create new mapset).

Let’s download Sentinel L2A products for spring/summer 2019. There are 7 products available as shown below.

i.sentinel.download -l settings=sentinel.txt map=jena_boundary area_relation=Contains \
start=2019-04-01 end=2019-10-31 producttype=S2MSI2A clouds=10

7 Sentinel product(s) found
a0ae6f58-4890-4382-bbd8-571874bfc65e 2019-06-26T10:20:31Z  1% S2MSI2A
caa11e7b-454d-4301-86b9-4c11659cc8a1 2019-04-17T10:20:31Z  3% S2MSI2A
31ad53f4-146a-41a8-bce6-d9e99dfd7f66 2019-04-22T10:20:29Z  3% S2MSI2A
e6ecc89f-0a06-498b-9e8d-354d1e8ead4e 2019-10-14T10:20:31Z  5% S2MSI2A
3fe2df76-019d-4611-8ffb-560683ae4500 2019-08-25T10:20:31Z  6% S2MSI2A
be91b224-0d6c-45d0-95a5-5cf5743a349e 2019-09-04T10:20:21Z  7% S2MSI2A
60149bc8-5ab0-4a79-b964-24672739d5b9 2019-04-07T10:20:21Z  9% S2MSI2A

Download selected Sentinel scenes:

i.sentinel.download settings=sentinel.txt map=jena_boundary area_relation=Contains \
start=2019-04-01 end=2019-10-31 producttype=S2MSI2A clouds=10 \
output=geodata/sentinel/2019

Note

Pre-downloaded Sentinel scenes are available in the sample dataset geodata/sentinel/2019.

Data can be imported by i.sentinel.import similarly as done in Unit 20 - Sentinel downloader. At fisrt check list of bands to be imported by -p flag. Since NDVI is going to be computed only 4th and 8th band are selected.

i.sentinel.import -p input=/home/user/geodata/sentinel/2019 pattern="B0(4|8)_10m"

i.sentinel.import -l -c input=/home/user/geodata/sentinel/2019 pattern="B0(4|8)_10m"

Note

In GRASS GIS 7.9 use register_output in order to create timestamps directly.

i.sentinel.import -l -c input=/home/user/geodata/sentinel/2019 pattern="B0(4|8)_10m" \
register_output=/home/user/sentinel-timestamps.txt

Example

T32UPB_20190407T102021_B04_10m|2019-04-07 10:26:45.035007|S2_4
T32UPB_20190407T102021_B08_10m|2019-04-07 10:26:45.035007|S2_8
T32UPB_20190422T102029_B04_10m|2019-04-22 10:26:50.312683|S2_4
...

Otherwise follow cookbook below. Unfortunately there are no timestamps for imported maps required by t.register. For this purpose a simple Python script has been designed.

#!/usr/bin/env python3

#%module
#% description: Timestamps Sentinel bands from current mapset.
#%end
#%option G_OPT_F_OUTPUT
#%end

import sys
import os
from datetime import datetime, timedelta

import grass.script as gs

from grass.pygrass.gis import Mapset

def main():
    mapset = Mapset()
    mapset.current()

    with open(options['output'], 'w') as fd:
        for rast in mapset.glist('raster'):
            items = rast.split('_')
            d = datetime.strptime(items[1], '%Y%m%dT%H%M%S')
            #fd.write("{0}|{1}{2}".format(rast, iso_date, os.linesep))
            ## workaround
            dd = d + timedelta(seconds=1)
            fd.write("{0}|{1}|{2}{3}".format(
                rast,
                d.strftime('%Y-%m-%d %H:%M:%S'),
                dd.strftime('%Y-%m-%d %H:%M:%S'),
                os.linesep))
        
    return 0

if __name__ == "__main__":
    options, flags = gs.parser()
    
    sys.exit(main())

Timestamps can be easily determined from raster map name, for example L2A_T32UPB_A021941_20190904T102045 raster map will be timestamped by 2019-09-04 10:20:45, see line 30.

Sample script to download: sentinel-timestamp.py

By running this script a timestamp file will be produced.

sentinel-timestamp.py output=sentinel-timestamps.txt

At this moment a new space time dataset can be created by t.create and all imported Sentinel bands registered by t.register.

t.create output=sentinel title="Sentinel L2A 2019" desc="Jena region"
t.register input=sentinel file=/home/user/sentinel-timestamps.txt

Let’s check basic metadata (t.info) and list of registered maps (t.rast.list).

t.info input=sentinel

...
| Start time:................. 2019-04-07 10:26:45.035007
| End time:................... 2019-10-14 10:26:46.742599
...
| Number of registered maps:.. 14

t.rast.list input=sentinel

name|mapset|start_time|end_time
T32UPB_20190407T102021_B04_10m|sentinel|2019-04-07 10:26:45.035007|None
T32UPB_20190407T102021_B08_10m|sentinel|2019-04-07 10:26:45.035007|None
T32UPB_20190417T102031_B04_10m|sentinel|2019-04-17 10:26:46.415802|None
...

NDVI ST computation

For NDVI computation 4th and 8th bands are required (Unit 05 - Simple computation). Map algebra is performed in the case of spatio-temporal data by t.rast.mapcalc which requires data separated into spatio-temporal datasets (see example in Unit 22 - Spatio-temporal scripting). Such datasets can be prepared by t.rast.extract.

t.rast.extract input=sentinel where="name like '%B04%'" output=b4
t.rast.extract input=sentinel where="name like '%B08%'" output=b8

Let’s check content of the new datasets by t.rast.list.

t.rast.list input=b4
t.rast.list input=b8

Set computational region by g.region including mask for area of interest by r.mask.

g.region vector=jena_boundary align=T32UPB_20190407T102021_B04_10m
r.mask vector=jena_boundary

NDVI (see Unit 05 - Simple computation) computation on spatio-temporal datasets can be performed in parallel (nproc).

t.rast.mapcalc input=b4,b8 output=ndvi \
expression="float(b8 - b4) / ( b8 + b4 )" \
basename=ndvi nproc=3

Note

In GRASS GIS 7.9 due to band reference support the computation is much more straighforward.

t.rast.mapcalc inputs=sentinel.S2_8,sentinel.S2_4 output=ndvi basename=ndvi \
expression="float(sentinel.S2_8 - sentinel.S2_4) / (sentinel.S2_8 + sentinel.S2_4)"

When computation is finished ndvi color table can be by t.rast.colors.

t.rast.colors input=ndvi color=ndvi

Fig. 114 Simple NDVI animation (no clouds mask applied) created by g.gui.animation.

Note

Load data as multiple raster maps instead of space time dataset. There is problem with sampling related to trick with endtime mentioned above.

Cloud mask

The result is not perfect enough. At least cloud mask should be applied before NDVI computation. This operation can be easily solved by applying new space time dataset containing computed raster masks. A sample Python script has been designed for this purpose below. The script also produces a timestamp file similarly to sentinel-timestamp.py. Such mask can be created by r.mask, see line 30. But in this case a mask should be kept for further usage. Note that r.mask module produces normal raster map with unique name MASK. To disable mask it is enough to rename MASK map by g.rename, see line 43.

#!/usr/bin/env python3

#%module
#% description: Creates raster mask maps based on clouds mask features.
#%end
#%option G_OPT_V_MAP
#% description: Name of AOI vector map
#%end
#%option G_OPT_F_OUTPUT
#%end

import sys
import os
from datetime import datetime, timedelta

import grass.script as gs

from grass.pygrass.gis import Mapset
from grass.pygrass.modules import Module
from grass.pygrass.vector import Vector
from grass.pygrass.utils import copy

def main():
    mapset = Mapset()
    mapset.current()

    with open(options['output'], 'w') as fd:
        for rast in mapset.glist('raster', pattern='*_B04_10m'):
            items = rast.split('_')
            d = datetime.strptime(items[1], '%Y%m%dT%H%M%S')
            ## workaround
            dd = d + timedelta(seconds=1)

            vect = '{}_{}_MSK_CLOUDS'.format(items[0], items[1])
            mask_vect = '{}_{}'.format(vect, options['map'].split('@')[0])
            if Vector(vect).exist():
                Module('v.overlay', ainput=options['map'], binput=vect, operator='not',
                       output=mask_vect, overwrite=True)
            else:
                copy(options['map'], mask_vect, 'vector')
            Module('r.mask', vector=mask_vect, overwrite=True)
            Module('g.remove', flags='f', type='vector', name=mask_vect)
            Module('g.rename', raster=['MASK', mask_vect])
            fd.write("{0}|{1}|{2}{3}".format(
                mask_vect,
                d.strftime('%Y-%m-%d %H:%M:%S'),
                dd.strftime('%Y-%m-%d %H:%M:%S'),
                os.linesep))
        
    return 0

if __name__ == "__main__":
    options, flags = gs.parser()
    
    sys.exit(main())

Sample script to download: sentinel-cloud-mask.py

sentinel-cloud-mask.py map=jena_boundary output=cloud-timestamps.txt

Now we can create a new space time dataset with raster cloud masks registered.

t.create output=clouds title="Sentinel L2A 2019 (clouds)" desc="Jena region"
t.register input=clouds file=cloud-timestamps.txt

And apply modified expression for map algebra.

t.rast.mapcalc input=b4,b8,clouds output=ndvi \
expression="ndvi = if(isnull(clouds), null(), float(b8 - b4) / ( b8 + b4 ))" \
basename=ndvi nproc=3

t.rast.colors in=ndvi color=ndvi

Fig. 115 Simple NDVI animation with clouds mask applied.