Source code for grass.temporal.temporal_raster_algebra

"""!@package grass.temporal

Temporal raster algebra

(C) 2013 by the GRASS Development Team
This program is free software under the GNU General Public
License (>=v2). Read the file COPYING that comes with GRASS
for details.

:authors: Thomas Leppelt and Soeren Gebbert

.. code-block:: python

    >>> p = TemporalRasterAlgebraLexer()
    >>> p.build()
    >>> p.debug = True
    >>> expression =  'R = A[0,1,0] / B[0,0,1] * 20 + C[0,1,1] - 2.45'
    >>> p.test(expression)
    R = A[0,1,0] / B[0,0,1] * 20 + C[0,1,1] - 2.45
    LexToken(NAME,'R',1,0)
    LexToken(EQUALS,'=',1,2)
    LexToken(NAME,'A',1,4)
    LexToken(L_SPAREN,'[',1,5)
    LexToken(INT,0,1,6)
    LexToken(COMMA,',',1,7)
    LexToken(INT,1,1,8)
    LexToken(COMMA,',',1,9)
    LexToken(INT,0,1,10)
    LexToken(R_SPAREN,']',1,11)
    LexToken(DIV,'/',1,13)
    LexToken(NAME,'B',1,15)
    LexToken(L_SPAREN,'[',1,16)
    LexToken(INT,0,1,17)
    LexToken(COMMA,',',1,18)
    LexToken(INT,0,1,19)
    LexToken(COMMA,',',1,20)
    LexToken(INT,1,1,21)
    LexToken(R_SPAREN,']',1,22)
    LexToken(MULT,'*',1,24)
    LexToken(INT,20,1,26)
    LexToken(ADD,'+',1,29)
    LexToken(NAME,'C',1,31)
    LexToken(L_SPAREN,'[',1,32)
    LexToken(INT,0,1,33)
    LexToken(COMMA,',',1,34)
    LexToken(INT,1,1,35)
    LexToken(COMMA,',',1,36)
    LexToken(INT,1,1,37)
    LexToken(R_SPAREN,']',1,38)
    LexToken(SUB,'-',1,40)
    LexToken(FLOAT,2.45,1,42)

"""

from __future__ import annotations

try:
    from ply import yacc
except ImportError:
    pass

import grass.pygrass.modules as pymod

from .space_time_datasets import RasterDataset
from .temporal_raster_base_algebra import (
    TemporalRasterAlgebraLexer,
    TemporalRasterBaseAlgebraParser,
)


[docs]class TemporalRasterAlgebraParser(TemporalRasterBaseAlgebraParser): """The temporal raster algebra class""" def __init__( self, pid: int | None = None, run: bool = False, debug: bool = True, spatial: bool = False, register_null: bool = False, dry_run: bool = False, nprocs: int = 1, time_suffix=None, ) -> None: TemporalRasterBaseAlgebraParser.__init__( self, pid=pid, run=run, debug=debug, spatial=spatial, register_null=register_null, dry_run=dry_run, nprocs=nprocs, time_suffix=time_suffix, ) if spatial is True: self.m_mapcalc = pymod.Module("r.mapcalc", region="union", run_=False) else: self.m_mapcalc = pymod.Module("r.mapcalc") self.m_mremove = pymod.Module("g.remove")
[docs] def parse(self, expression, basename=None, overwrite: bool = False): # Check for space time dataset type definitions from temporal algebra lx = TemporalRasterAlgebraLexer() lx.build() lx.lexer.input(expression) while True: tok = lx.lexer.token() if not tok: break if tok.type in {"STVDS", "STRDS", "STR3DS"}: raise SyntaxError("Syntax error near '%s'" % (tok.type)) self.lexer = TemporalRasterAlgebraLexer() self.lexer.build() self.parser = yacc.yacc(module=self, debug=self.debug, write_tables=False) self.overwrite = overwrite self.count = 0 self.stdstype = "strds" self.maptype = "raster" self.mapclass = RasterDataset self.basename = basename self.expression = expression self.parser.parse(expression) return self.process_chain_dict
[docs] def p_statement_assign(self, t) -> None: # The expression should always return a list of maps. """ statement : stds EQUALS expr """ TemporalRasterBaseAlgebraParser.p_statement_assign(self, t)
[docs] def p_ts_neighbour_operation(self, t) -> None: # Spatial and temporal neighbour operations via indexing # Examples: # A[1,0] # B[-2] # C[-2,1,3] """ expr : stds L_SPAREN number COMMA number R_SPAREN | stds L_SPAREN number R_SPAREN | stds L_SPAREN number COMMA number COMMA number R_SPAREN """ # Check input stds. maplist = self.check_stds(t[1]) row_neigbour = None col_neigbour = None if len(t) == 5: t_neighbour = t[3] elif len(t) == 7: t_neighbour = 0 row_neigbour = t[3] col_neigbour = t[5] elif len(t) == 9: t_neighbour = t[7] row_neigbour = t[3] col_neigbour = t[5] if self.run: resultlist = [] max_index = len(maplist) for map_i in maplist: # Get map index and temporal extent. map_index = maplist.index(map_i) new_index = map_index + t_neighbour if new_index < max_index and new_index >= 0: map_i_t_extent = map_i.get_temporal_extent() # Get neighbouring map and set temporal extent. map_n = maplist[new_index] # Generate an intermediate map for the result map list. map_new = self.generate_new_map(map_n, bool_op="and", copy=True) map_new.set_temporal_extent(map_i_t_extent) # Create r.mapcalc expression string for the operation. if "cmd_list" in dir(map_new) and len(t) == 5: cmdstring = "%s" % (map_new.cmd_list) elif "cmd_list" not in dir(map_new) and len(t) == 5: cmdstring = "%s" % (map_n.get_id()) elif "cmd_list" in dir(map_new) and len(t) in {7, 9}: cmdstring = "%s[%s,%s]" % ( map_new.cmd_list, row_neigbour, col_neigbour, ) elif "cmd_list" not in dir(map_new) and len(t) in {7, 9}: cmdstring = "%s[%s,%s]" % ( map_n.get_id(), row_neigbour, col_neigbour, ) # Set new command list for map. map_new.cmd_list = cmdstring # Append map with updated command list to result list. resultlist.append(map_new) t[0] = resultlist
if __name__ == "__main__": import doctest doctest.testmod()