temporal_raster3d_algebra.py 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138
  1. """
  2. Temporal 3d 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. """
  9. from temporal_raster_base_algebra import *
  10. ###############################################################################
  11. class TemporalRaster3DAlgebraParser(TemporalRasterBaseAlgebraParser):
  12. """The temporal raster algebra class"""
  13. def __init__(self, pid=None, run=False, debug=True, spatial=False,
  14. nprocs=1, register_null=False):
  15. TemporalRasterBaseAlgebraParser.__init__(self, pid, run, debug,
  16. spatial, nprocs,
  17. register_null)
  18. self.m_mapcalc = pymod.Module('r3.mapcalc')
  19. self.m_remove = pymod.Module('g.remove')
  20. def parse(self, expression, basename=None, overwrite=False):
  21. self.lexer = TemporalRasterAlgebraLexer()
  22. self.lexer.build()
  23. self.parser = yacc.yacc(module=self, debug=self.debug)
  24. self.overwrite = overwrite
  25. self.count = 0
  26. self.stdstype = "str3ds"
  27. self.basename = basename
  28. self.expression = expression
  29. self.parser.parse(expression)
  30. def remove_empty_maps(self):
  31. """! Removes the intermediate vector maps.
  32. """
  33. if self.empty_maps:
  34. self.msgr.message(_("Removing empty 3D raster maps"))
  35. namelist = self.empty_maps.values()
  36. max = 100
  37. chunklist = [namelist[i:i + max] for i in range(0, len(namelist), max)]
  38. for chunk in chunklist:
  39. stringlist = ",".join(chunk)
  40. if self.run:
  41. m = copy.deepcopy(self.m_remove)
  42. m.inputs["type"].value = "rast3d"
  43. m.inputs["name"].value = stringlist
  44. m.flags["f"].value = True
  45. m.run()
  46. ######################### Temporal functions ##############################
  47. def p_statement_assign(self, t):
  48. # The expression should always return a list of maps.
  49. """
  50. statement : stds EQUALS expr
  51. """
  52. TemporalRasterBaseAlgebraParser.p_statement_assign(self, t)
  53. def p_ts_neighbor_operation(self, t):
  54. # Examples:
  55. # A[1,0,-1]
  56. # B[-2]
  57. # C[1,-2,1,3]
  58. """
  59. expr : stds L_SPAREN number COMMA number COMMA number R_SPAREN
  60. | stds L_SPAREN number R_SPAREN
  61. | stds L_SPAREN number COMMA number COMMA number COMMA number R_SPAREN
  62. """
  63. # Check input stds.
  64. maplist = self.check_stds(t[1])
  65. t_neighbor = 0
  66. row_neighbor = None
  67. col_neighbor = None
  68. depth_neighbor = None
  69. if len(t) == 5:
  70. t_neighbor = t[3]
  71. elif len(t) == 9:
  72. row_neighbor = t[3]
  73. col_neighbor = t[5]
  74. depth_neighbor = t[7]
  75. elif len(t) == 11:
  76. t_neighbor = t[9]
  77. row_neighbor = t[3]
  78. col_neighbor = t[5]
  79. depth_neighbor = t[7]
  80. if self.run:
  81. resultlist = []
  82. max_index = len(maplist)
  83. for map_i in maplist:
  84. # Get map index and temporal extent.
  85. map_index = maplist.index(map_i)
  86. new_index = map_index + t_neighbor
  87. if new_index < max_index and new_index >= 0:
  88. map_i_t_extent = map_i.get_temporal_extent()
  89. # Get neighboring map and set temporal extent.
  90. map_n = maplist[new_index]
  91. # Generate an intermediate map for the result map list.
  92. map_new = self.generate_new_map(map_n, bool_op='and',
  93. copy=True)
  94. map_new.set_temporal_extent(map_i_t_extent)
  95. # Create r.mapcalc expression string for the operation.
  96. if "cmd_list" in dir(map_new) and len(t) == 5:
  97. cmdstring = "%s" % (map_new.cmd_list)
  98. elif "cmd_list" not in dir(map_new) and len(t) == 5:
  99. cmdstring = "%s" % (map_n.get_id())
  100. elif "cmd_list" in dir(map_new) and len(t) in (9, 11):
  101. cmdstring = "%s[%s,%s,%s]" % (map_new.cmd_list,
  102. row_neighbor,
  103. col_neighbor,
  104. depth_neighbor)
  105. elif "cmd_list" not in dir(map_new) and len(t) in (9, 11):
  106. cmdstring = "%s[%s,%s,%s]" % (map_n.get_id(),
  107. row_neighbor,
  108. col_neighbor,
  109. depth_neighbor)
  110. # Set new command list for map.
  111. map_new.cmd_list = cmdstring
  112. # Append map with updated command list to result list.
  113. resultlist.append(map_new)
  114. t[0] = resultlist
  115. ###############################################################################
  116. if __name__ == "__main__":
  117. import doctest
  118. doctest.testmod()