Unit 22 - Spatio-temporal scripting

GRASS offers a Python API for space-time processing. The usage is presented in the script described below.

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

UI

Let’s define input and output parameters of the script:

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

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

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

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

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

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

Python Temporal API

Used functions from GRASS GIS Temporal Framework:

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

  • space time datasets are open on lines 144-146 by open_old_stds

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

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

  1#!/usr/bin/env python3
  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: min_area
 39# % description: Threshold for removing small areas in m2
 40# % answer: 1600
 41# %end
 42
 43import sys
 44import os
 45import atexit
 46from subprocess import PIPE
 47
 48from grass.script import parser, parse_key_val
 49from grass.pygrass.modules import Module
 50    
 51def cleanup():
 52    Module('g.remove', flags='f', name='region_mask', type='vector')
 53    Module('g.remove', flags='f', name='ndvi', type='raster')
 54    Module('g.remove', flags='f', name='ndvi_class', type='raster')
 55    Module('g.remove', flags='f', name='ndvi_class_area', type='raster')
 56    Module('g.remove', flags='f', name='ndvi_class_filled_i', type='raster')
 57    Module('g.remove', flags='f', name='ndvi_class_filled', type='vector')
 58
 59def compute(b4, b8, msk, min_area_ha, output):
 60    Module("g.region",
 61           overwrite=True,
 62           raster=msk,
 63           align=b4)
 64
 65    Module("r.mask",
 66           overwrite=True,
 67           maskcats="*",
 68           raster=msk,
 69           layer="1")
 70
 71    Module("i.vi",
 72           overwrite=True,
 73           output="ndvi",
 74           viname="ndvi",
 75           red=b4,
 76           nir=b8,
 77           storage_bit=8)
 78
 79    Module("r.recode",
 80           overwrite=True,
 81           input="ndvi",
 82           output="ndvi_class",
 83           rules="-",
 84           stdin_="-1:0.1:1\n0.1:0.5:2\n0.5:1:3")
 85
 86    Module("r.reclass.area",
 87           overwrite=True,
 88           input="ndvi_class",
 89           output="ndvi_class_area",
 90           value=min_area_ha,
 91           mode="greater",
 92           method="reclass")
 93
 94    Module("r.grow.distance",
 95           overwrite=True,
 96           input="ndvi_class_area",
 97           value=output,
 98           metric="euclidean")
 99
100    Module("r.colors",
101           map=output,
102           rules="-",
103           stdin_="1 grey\n2 255 255 0\n3 green",
104           offset=0,
105           scale=1)
106    
107def stats(output, date, fd):
108    fd.write('-' * 80)
109    fd.write('\n')
110    fd.write('NDVI class statistics ({0}: {1})'.format(output, date))
111    fd.write('\n')
112    fd.write('-' * 80)
113    fd.write('\n')
114
115    ret = Module('r.stats', input=output, flags='ian', stdout_=PIPE)
116    for line in ret.outputs.stdout.splitlines():
117        # parse line (eg. 1 2737300.000000)
118        data = line.split(' ')
119        cat = data[0]
120        area = float(data[-1])
121        fd.write('NDVI class {0}: {1:.1f} ha\n'.format(cat, area/1e4)) 
122
123    fd.write('-' * 80)
124    fd.write('\n')
125    # we need integer map
126    Module('r.mapcalc', expression='ndvi_class_filled_i = int({})'.format(output), overwrite=True)
127    Module('r.to.vect', flags='v', input='ndvi_class_filled_i', output='ndvi_class_filled', type='area', overwrite=True)
128
129    Module('v.rast.stats', flags='c', map='ndvi_class_filled', raster='ndvi',
130           column_prefix='ndvi', method=['minimum','maximum','average'])
131    # v.db.select: don't print column names (-c)
132    ret = Module('v.db.select', flags='c', map='ndvi_class_filled', separator='comma', stdout_=PIPE)
133    for line in ret.outputs.stdout.splitlines():
134        # parse line (eg. 1,,-0.433962264150943,0.740350877192983,0.051388909449992)
135        cat,label,min,max,mean = line.split(',')
136        fd.write('NDVI class {0}: {1:.4f} (min) {2:.4f} (max) {3:.4f} (mean)\n'.format(
137        cat, float(min), float(max), float(mean)))
138        
139def main():
140    import grass.temporal as tgis
141
142    tgis.init()
143
144    sp4=tgis.open_old_stds(options['b4'], 'raster')
145    sp8=tgis.open_old_stds(options['b8'], 'raster')
146    msk=tgis.open_old_stds(options['mask'], 'raster')
147
148    min_area = int(options['min_area']) / 1e4
149    idx=1
150    fd=open(options['output'], 'w')
151    for item in sp4.get_registered_maps(columns='name,start_time'):
152        b4=item[0]
153        date=item[1]
154        b8=sp8.get_registered_maps(columns='name',
155                                     where="start_time='{}'".format(date))[0][0]
156        ms=msk.get_registered_maps(columns='name',
157                                     where="start_time='{}'".format(date))[0][0]
158        output='{}_{}'.format(options['basename'], idx)
159        compute(b4, b8, ms, min_area, output)
160        stats(output, date, fd)
161        cleanup()
162        idx += 1
163
164    fd.close()
165    
166    return 0
167
168if __name__ == "__main__":
169    options, flags=parser()
170    sys.exit(main())

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

Example of usage:

ndvi-tgrass.py b4=b4 b8=b8 mask=clouds basename=ndvi out=stats.txt

Possible output:

--------------------------------------------------------------------------------
NDVI class statistics (ndvi_1: 2019-04-07 10:26:45.035007)
--------------------------------------------------------------------------------
NDVI class 1: 122.4 ha
NDVI class 2: 4873.0 ha
NDVI class 3: 6333.2 ha
--------------------------------------------------------------------------------
NDVI class 1: -0.2073 (min) 0.5353 (max) 0.0557 (mean)
NDVI class 2: -0.2380 (min) 0.9989 (max) 0.3754 (mean)
NDVI class 3: -0.4533 (min) 0.9988 (max) 0.6468 (mean)
...
--------------------------------------------------------------------------------
NDVI class statistics (ndvi_7: 2019-10-14 10:26:46.742599)
--------------------------------------------------------------------------------
NDVI class 1: 159.8 ha
NDVI class 2: 2669.3 ha
NDVI class 3: 8602.4 ha
--------------------------------------------------------------------------------
NDVI class 1: -0.2253 (min) 0.7452 (max) 0.0478 (mean)
NDVI class 2: -1.0000 (min) 0.9994 (max) 0.2999 (mean)
NDVI class 3: -0.9978 (min) 0.9994 (max) 0.6992 (mean)