Unit 24 - Spatio-temporal scripting

There is another GRASS Python library to be introduced - GRASS GIS Temporal Framework.

In a script below similar computation like in Unit 11 - PyGRASS scripting and Unit 12 - Script User Interface will be performed. But instead of processing single Sentinel scene, a spatio-temporal dataset of multiple Sentinel scenes will be processed.

Let’s define input and output parameters:

  • b4 - Name of input 4th band space time raster dataset (line 19)

  • b8 - Name of input 8th band space time raster dataset (line 23)

  • clouds - Name of the input mask (region+clouds) space time raster dataset (line 27)

  • clouds - Name for output stats file (line 30)

  • basename - Basename for output raster maps (line 33)

  • threshold - Threshold for removing small areas (line 38)

Focus on functions used from GRASS GIS Temporal Framework:

  • initialization must be done by init function, see line 138

  • space time datasets are open on lines 140-142 by open_old_stds

  • raster maps registered in reference dataset (b4) are listed on line 146 by get_registered_maps

  • related raster maps in two other datasets (b8, cl) are searched on lines 149-152 by get_registered_maps with where parameter

  1#!/usr/bin/env python
  2#
  3##############################################################################
  4#
  5# MODULE:       ndvi-tgrass-v1
  6#
  7# AUTHOR(S):    martin
  8#
  9# PURPOSE:      NDVI TGRASS version 1
 10#
 11# DATE:         Sat Feb  3 15:45:35 2018
 12#
 13##############################################################################
 14
 15#%module
 16#% description: NDVI TGRASS script version 1
 17#%end                
 18#%option G_OPT_STRDS_INPUT
 19#% key: b4
 20#% description: Name of input 4th band space time raster dataset
 21#%end
 22#%option G_OPT_STRDS_INPUT
 23#% key: b8
 24#% description: Name of input 8th band space time raster dataset
 25#%end
 26#%option G_OPT_STRDS_INPUT
 27#% key: mask
 28#% description: Name of input mask space time raster dataset
 29#%end
 30#%option G_OPT_F_OUTPUT
 31#%end
 32#%option
 33#% key: basename
 34#% description: Basename for output raster maps
 35#% required: yes
 36#%end
 37#%option
 38#% key: threshold
 39#% description: Threshold for removing small areas
 40#% answer: 1600
 41#%end
 42
 43import sys
 44import os
 45import atexit
 46
 47from grass.pygrass.modules import Module
 48from grass.script import parser
 49from grass.script.vector import vector_db_select
 50    
 51def cleanup():
 52    Module('g.remove', flags='f', name='ndvi', type='raster')
 53    Module('g.remove', flags='f', name='ndvi_class', type='raster')
 54    Module('g.remove', flags='f', name='ndvi_class', type='vector')
 55
 56def compute(b4, b8, msk, output):
 57
 58    Module("g.region",
 59           overwrite = True,
 60           vector = msk,
 61           align = b4)
 62
 63    Module("r.mask",
 64           overwrite = True,
 65           raster = msk)
 66
 67    Module("i.vi",
 68           overwrite = True,
 69           red = b4,
 70           output = "ndvi",
 71           nir = b8)
 72                
 73    recode_str="""-1:0.1:1
 740.1:0.5:2
 750.5:1:3"""
 76
 77    Module("r.recode",
 78           overwrite = True,
 79           input = "ndvi",
 80           output = "ndvi_class",
 81           rules = "-",
 82           stdin_ = recode_str)
 83
 84    colors_str="""1 grey
 852 255 255 0
 863 green"""
 87    Module("r.colors",
 88           map = "ndvi_class",
 89           rules = "-",
 90           stdin_ = colors_str)
 91    
 92    Module("r.to.vect",
 93           flags = 'sv',
 94           overwrite = True,
 95           input = "ndvi_class",
 96           output = "ndvi_class",
 97           type = "area")
 98
 99    Module("v.clean",
100           overwrite = True,
101           input = "ndvi_class",
102           output = output,
103           tool = "rmarea",
104           threshold = options['threshold'])
105
106def stats(output, date, fd):
107    fd.write('-' * 80)
108    fd.write(os.linesep)
109    fd.write('NDVI class statistics ({0}: {1})'.format(output, date))
110    fd.write(os.linesep)
111    fd.write('-' * 80)
112    fd.write(os.linesep)
113    from subprocess import PIPE
114    ret = Module('v.report', map=output, option='area',
115                 stdout_=PIPE)
116    for line in ret.outputs.stdout.splitlines()[1:]: # skip first line (cat|label|area)
117        # parse line (eg. 1||2712850)
118        data = line.split('|')
119        cat = data[0]
120        area = float(data[-1])
121        fd.write('NDVI class {0}: {1:.1f} ha'.format(cat, area/1e4))
122        fd.write(os.linesep)
123
124    # v.to.rast: use -c flag for updating statistics if exists
125    Module('v.rast.stats', flags='c', map=output, raster='ndvi',
126           column_prefix='ndvi', method=['minimum','maximum','average'])
127    
128    data = vector_db_select(output)
129    for vals in data['values'].itervalues():
130        # unfortunately we need to cast values by float
131        fd.write('NDVI class {0}: {1:.4f} (min) {2:.4f} (max) {3:.4f} (mean)'.format(
132            vals[0], float(vals[2]), float(vals[3]), float(vals[4])))
133        fd.write(os.linesep)
134        
135def main():
136    import grass.temporal as tgis
137
138    tgis.init()
139
140    sp4 = tgis.open_old_stds(options['b4'], 'raster')
141    sp8 = tgis.open_old_stds(options['b8'], 'raster')
142    msk = tgis.open_old_stds(options['mask'], 'raster')
143
144    idx = 1
145    fd = open(options['output'], 'w')
146    for item in sp4.get_registered_maps(columns='name,start_time'):
147        b4 = item[0]
148        date=item[1]
149        b8 = sp8.get_registered_maps(columns='name',
150                                     where="start_time = '{}'".format(date))[0][0]
151        ms = msk.get_registered_maps(columns='name',
152                                     where="start_time = '{}'".format(date))[0][0]
153        output = '{}_{}'.format(options['basename'], idx)
154        compute(b4, b8, ms, output)
155        stats(output, date, fd)
156        cleanup()
157        idx += 1
158
159    fd.close()
160    
161    return 0
162
163if __name__ == "__main__":
164    options, flags = parser()
165    sys.exit(main())

Example of usage:

ndvi-tgrass.py b4=b4 b8=b8 clouds=clouds basename=ndvi_t output=stats.txt

--------------------------------------------------------------------------------
NDVI class statistics (ndvi_t_1: 2017-05-06 10:50:31)
--------------------------------------------------------------------------------
NDVI class 1: 6714.4 ha
NDVI class 2: 11480.8 ha
NDVI class 3: 29852.3 ha
NDVI class 1: -0.9929 (min) 0.9928 (max) -0.3230 (mean)
NDVI class 2: -0.9512 (min) 0.9993 (max) 0.3327 (mean)
NDVI class 3: -0.9512 (min) 0.9988 (max) 0.7312 (mean)
...
--------------------------------------------------------------------------------
NDVI class statistics (ndvi_t_4: 2017-07-05 10:50:31)
--------------------------------------------------------------------------------
NDVI class 1: 24122.4 ha
NDVI class 2: 7243.5 ha
NDVI class 3: 16682.1 ha
NDVI class 1: -0.9914 (min) 0.9710 (max) -0.0292 (mean)
NDVI class 2: -0.9286 (min) 0.9487 (max) 0.2506 (mean)
NDVI class 3: -0.9545 (min) 0.9982 (max) 0.8190 (mean)

Note

Script will run for a while in the case of Oslo region.

Sample script to download: ndvi-tgrass-v1.py