Unit 21 - Sentinel ST

This unit is dedicated to downloading and processing time-series Sentinel products.

Let’s download Sentinel products (still Sentinel L2A as we did in Unit 20 - Sentinel downloader) for year 2017. There are 7 products available as shown below. Is is enough for our purpose.

i.sentinel.download -l settings=sentinel.txt map=jena_boundary area_relation=Contains \
start=2017-01-01 end=2017-12-31 producttype=S2MSI2Ap
7 Sentinel product(s) found
e5df8b4f-a4c6-47bd-88f3-e16b7540cc7a 2017-05-27T10:20:31Z  2% S2MSI2Ap
b62afeda-a28d-475c-8220-91e24fc368ff 2017-05-17T10:20:31Z  2% S2MSI2Ap
9a6bc289-98e9-4489-84eb-1aac95aaa056 2017-08-15T10:20:21Z  3% S2MSI2Ap
35c72002-78a0-45f8-b39d-66c2d7b4ad87 2017-10-14T10:20:21Z  6% S2MSI2Ap
c0ae8085-c1bb-4a27-89f2-2138a0866586 2017-07-06T10:20:21Z 12% S2MSI2Ap
433ebfbc-5144-42f8-97dc-b9f4f1eb7b5a 2017-11-03T10:22:01Z 12% S2MSI2Ap
fc56a594-d9d8-4e93-8dec-af3a58b24080 2017-09-04T10:20:21Z 19% S2MSI2Ap

Note

For simplification you can find this data in sample dataset (download 7z-archive) in sentinel/2017 folder or as GRASS mapset as grassdata/st.tar.gz. You can use this mapset (extract to jena-utm location) and go to next part if you like.

Now we can start downloading data.

i.sentinel.download settings=sentinel.txt map=jena_boundary area_relation=Contains \
start=2017-01-01 end=2017-12-31 producttype=S2MSI2Ap output=/opt/geodata/sentinel-2017

When data are downloaded we can import them into our jena-utm location. It is a good idea to import data into newly created mapset (g.mapset, see Unit 16). In command below we force linking instead of importing data.

g.mapset -c mapset=st
i.sentinel.import -l -c input=/opt/geodata/sentinel-2017 pattern="B0[4|8]_10m"

In next step a new space time raster dataset will created and imported maps registered. Unfortunately we don’t have timestamps for imported maps required by t.register (as shown in Unit 25). For this purpose we will design a simple Python script as shown below.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
#!/usr/bin/env python

#%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[2], '%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_20170517T102031_B04_10m raster map will be timestamped by 2017-05-17 10:20:31, see line 30.

Note

To avoid error in computation we also set up endtime timestamp as starttime+1sec, see 27.

This workaround should be avoided. Determine how?

Sample script to download: sentinel-timestamp.py

Now we can run our script to produce timestamp file.

sentinel-timestamp.py output=sentinel.txt

At this moment we can create a new space time dataset (t.create) and register all imported Sentinel bands (t.register).

t.create output=sentinel title="Sentinel L2A 2017" desc="Jena region"
t.register input=sentinel file=sentinel.txt

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

t.info input=sentinel
t.rast.list input=sentinel

NDVI ST computation

For NDVI computation we will create new space time datasets for 4th and 8th bands respectively (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

Before computation let’s desired computation region by g.region.

g.region vector=jena_boundary align=L2A_T32UPB_20170517T102031_B04_10m

And perform NDVI computation based map algebra expression as described in Unit 05 - Perform computation in parallel (nproc) by t.rast.mapcalc.

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

When computation is finished we can set up ndvi color table in order to create fancy animation as shown below.

t.rast.colors input=ndvi color=ndvi
../_images/simple-animation.gif

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

Note

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

The result is not perfect as we can see. We should apply at least cloud mask before NDVI computation. This operation can be easily solved by applying new space time dataset containing computed raster masks. We will create a new Python script for this purpose. The script will also produce timestamp file similarly to sentinel-timestamp.py. The mask can be created as we know from Unit 05 - Perform computation by r.mask, see line 39. But in this case we need to keep mask for further usage. We will use simple trick, r.mask produces normal raster map with unique name MASK. It’s enough to rename this map by g.rename, see line 30.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
#!/usr/bin/env python

#%module
#% description: Creates raster mask maps based on clouds mask features.
#%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

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

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

            Module('r.mask', vector=vect, flags='i')
            Module('g.rename', raster=['MASK', vect])
            fd.write("{0}|{1}|{2}{3}".format(
                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-clouds-mask.py

sentinel-clouds-mask.py output=clouds.txt

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

t.create output=clouds title="Sentinel L2A 2017 (clouds)" desc="Jena region"
t.register input=clouds file=clouds.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 --overwrite

t.rast.colors in=ndvi color=ndvi

The result is a little bit better, provided cloud mask are not perfect as we already know.

../_images/simple-animation-clouds.gif

Fig. 110 Simple NDVI animation with clouds mask applied.