temporal_raster_algebra.py 6.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183
  1. """!@package grass.temporal
  2. Temporal raster algebra
  3. (C) 2013 by the GRASS Development Team
  4. This program is free software under the GNU General Public
  5. License (>=v2). Read the file COPYING that comes with GRASS
  6. for details.
  7. :authors: Thomas Leppelt and Soeren Gebbert
  8. .. code-block:: python
  9. >>> p = TemporalRasterAlgebraLexer()
  10. >>> p.build()
  11. >>> p.debug = True
  12. >>> expression = 'R = A[0,1,0] / B[0,0,1] * 20 + C[0,1,1] - 2.45'
  13. >>> p.test(expression)
  14. R = A[0,1,0] / B[0,0,1] * 20 + C[0,1,1] - 2.45
  15. LexToken(NAME,'R',1,0)
  16. LexToken(EQUALS,'=',1,2)
  17. LexToken(NAME,'A',1,4)
  18. LexToken(L_SPAREN,'[',1,5)
  19. LexToken(INT,0,1,6)
  20. LexToken(COMMA,',',1,7)
  21. LexToken(INT,1,1,8)
  22. LexToken(COMMA,',',1,9)
  23. LexToken(INT,0,1,10)
  24. LexToken(R_SPAREN,']',1,11)
  25. LexToken(DIV,'/',1,13)
  26. LexToken(NAME,'B',1,15)
  27. LexToken(L_SPAREN,'[',1,16)
  28. LexToken(INT,0,1,17)
  29. LexToken(COMMA,',',1,18)
  30. LexToken(INT,0,1,19)
  31. LexToken(COMMA,',',1,20)
  32. LexToken(INT,1,1,21)
  33. LexToken(R_SPAREN,']',1,22)
  34. LexToken(MULT,'*',1,24)
  35. LexToken(INT,20,1,26)
  36. LexToken(ADD,'+',1,29)
  37. LexToken(NAME,'C',1,31)
  38. LexToken(L_SPAREN,'[',1,32)
  39. LexToken(INT,0,1,33)
  40. LexToken(COMMA,',',1,34)
  41. LexToken(INT,1,1,35)
  42. LexToken(COMMA,',',1,36)
  43. LexToken(INT,1,1,37)
  44. LexToken(R_SPAREN,']',1,38)
  45. LexToken(SUB,'-',1,40)
  46. LexToken(FLOAT,2.45,1,42)
  47. """
  48. from __future__ import print_function
  49. try:
  50. import ply.lex as lex
  51. import ply.yacc as yacc
  52. except:
  53. pass
  54. # i18N
  55. import gettext
  56. from .temporal_raster_base_algebra import TemporalRasterBaseAlgebraParser,\
  57. TemporalRasterAlgebraLexer
  58. import grass.pygrass.modules as pymod
  59. from .space_time_datasets import RasterDataset
  60. ###############################################################################
  61. class TemporalRasterAlgebraParser(TemporalRasterBaseAlgebraParser):
  62. """The temporal raster algebra class"""
  63. def __init__(self, pid=None, run=False, debug=True, spatial=False,
  64. register_null=False, dry_run=False, nprocs=1):
  65. TemporalRasterBaseAlgebraParser.__init__(self, pid=pid, run=run, debug=debug,
  66. spatial=spatial, register_null=register_null,
  67. dry_run=dry_run, nprocs=nprocs)
  68. self.m_mapcalc = pymod.Module('r.mapcalc')
  69. self.m_mremove = pymod.Module('g.remove')
  70. def parse(self, expression, basename = None, overwrite=False):
  71. # Check for space time dataset type definitions from temporal algebra
  72. l = TemporalRasterAlgebraLexer()
  73. l.build()
  74. l.lexer.input(expression)
  75. while True:
  76. tok = l.lexer.token()
  77. if not tok: break
  78. if tok.type == "STVDS" or tok.type == "STRDS" or tok.type == "STR3DS":
  79. raise SyntaxError("Syntax error near '%s'" %(tok.type))
  80. self.lexer = TemporalRasterAlgebraLexer()
  81. self.lexer.build()
  82. self.parser = yacc.yacc(module=self, debug=self.debug)
  83. self.overwrite = overwrite
  84. self.count = 0
  85. self.stdstype = "strds"
  86. self.maptype = "raster"
  87. self.mapclass = RasterDataset
  88. self.basename = basename
  89. self.expression = expression
  90. self.parser.parse(expression)
  91. return self.process_chain_dict
  92. ######################### Temporal functions ##############################
  93. def p_statement_assign(self, t):
  94. # The expression should always return a list of maps.
  95. """
  96. statement : stds EQUALS expr
  97. """
  98. TemporalRasterBaseAlgebraParser.p_statement_assign(self, t)
  99. def p_ts_neighbour_operation(self, t):
  100. # Spatial and temporal neighbour operations via indexing
  101. # Examples:
  102. # A[1,0]
  103. # B[-2]
  104. # C[-2,1,3]
  105. """
  106. expr : stds L_SPAREN number COMMA number R_SPAREN
  107. | stds L_SPAREN number R_SPAREN
  108. | stds L_SPAREN number COMMA number COMMA number R_SPAREN
  109. """
  110. # Check input stds.
  111. maplist = self.check_stds(t[1])
  112. row_neigbour = None
  113. col_neigbour = None
  114. if len(t) == 5:
  115. t_neighbour = t[3]
  116. elif len(t) == 7:
  117. t_neighbour = 0
  118. row_neigbour = t[3]
  119. col_neigbour = t[5]
  120. elif len(t) == 9:
  121. t_neighbour = t[7]
  122. row_neigbour = t[3]
  123. col_neigbour = t[5]
  124. if self.run:
  125. resultlist = []
  126. max_index = len(maplist)
  127. for map_i in maplist:
  128. # Get map index and temporal extent.
  129. map_index = maplist.index(map_i)
  130. new_index = map_index + t_neighbour
  131. if new_index < max_index and new_index >= 0:
  132. map_i_t_extent = map_i.get_temporal_extent()
  133. # Get neighbouring map and set temporal extent.
  134. map_n = maplist[new_index]
  135. # Generate an intermediate map for the result map list.
  136. map_new = self.generate_new_map(map_n, bool_op = 'and', copy = True)
  137. map_new.set_temporal_extent(map_i_t_extent)
  138. # Create r.mapcalc expression string for the operation.
  139. if "cmd_list" in dir(map_new) and len(t) == 5:
  140. cmdstring = "%s" %(map_new.cmd_list)
  141. elif "cmd_list" not in dir(map_new) and len(t) == 5:
  142. cmdstring = "%s" %(map_n.get_id())
  143. elif "cmd_list" in dir(map_new) and len(t) in (7, 9):
  144. cmdstring = "%s[%s,%s]" %(map_new.cmd_list, row_neigbour, col_neigbour)
  145. elif "cmd_list" not in dir(map_new) and len(t) in (7, 9):
  146. cmdstring = "%s[%s,%s]" %(map_n.get_id(), row_neigbour, col_neigbour)
  147. # Set new command list for map.
  148. map_new.cmd_list = cmdstring
  149. # Append map with updated command list to result list.
  150. resultlist.append(map_new)
  151. t[0] = resultlist
  152. ###############################################################################
  153. if __name__ == "__main__":
  154. import doctest
  155. doctest.testmod()