Parallelized implementation of mapcalc expressions in space time dataset extract and mapcalc modules.

New space time raster3d dataset mapcalc module.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@51375 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/python/temporal/Makefile

@@ -8,7 +8,7 @@ PYDIR = $(ETC)/python
 GDIR = $(PYDIR)/grass
 DSTDIR = $(GDIR)/temporal
 
-MODULES = base core abstract_dataset abstract_map_dataset abstract_space_time_dataset space_time_datasets space_time_datasets_tools metadata spatial_extent temporal_extent datetime_math temporal_granularity temporal_relationships unit_tests aggregation
+MODULES = base core abstract_dataset abstract_map_dataset abstract_space_time_dataset space_time_datasets space_time_datasets_tools metadata spatial_extent temporal_extent datetime_math temporal_granularity temporal_relationships unit_tests aggregation extract mapcalc
 
 PYFILES := $(patsubst %,$(DSTDIR)/%.py,$(MODULES) __init__)
 PYCFILES := $(patsubst %,$(DSTDIR)/%.pyc,$(MODULES) __init__)

lib/python/temporal/__init__.py

@@ -13,3 +13,5 @@ from temporal_granularity import *
 from temporal_relationships import *
 from unit_tests import *
 from aggregation import *
+from extract import *
+from mapcalc import *

temporal/Makefile

@@ -23,6 +23,7 @@ SUBDIRS = \
 	t.rast.extract \
 	t.rast3d.list \
 	t.rast3d.extract \
+	t.rast3d.mapcalc \
 	t.vect.list \
 	t.vect.what.strds \
 	t.vect.observe.strds

temporal/t.info/t.info.py

@@ -15,7 +15,7 @@
 #############################################################################
 
 #%module
-#% description: List informtion about space time and map datasets
+#% description: List information about space time datasets and maps
 #% keywords: temporal
 #% keywords: info
 #%end

temporal/t.rast.extract/t.rast.extract.py

@@ -74,160 +74,7 @@ def main():
     # Make sure the temporal database exists
     tgis.create_temporal_database()
     
-    mapset =  grass.gisenv()["MAPSET"]
-
-    if input.find("@") >= 0:
-        id = input
-    else:
-        id = input + "@" + mapset
-
-    sp = tgis.space_time_raster_dataset(id)
-    
-    dbif = tgis.sql_database_interface_connection()
-    dbif.connect()
-    
-    if sp.is_in_db(dbif) == False:
-	dbif.close()
-        grass.fatal(_("Space time %s dataset <%s> not found") % (sp.get_new_map_instance(None).get_type(), id))
-
-    if expression and not base:
-	dbif.close()
-        grass.fatal(_("Please specify base="))
-
-    sp.select(dbif)
-
-    if output.find("@") >= 0:
-        out_id = output
-    else:
-        out_id = output + "@" + mapset
-
-    # The new space time raster dataset
-    new_sp = tgis.space_time_raster_dataset(out_id)
-    if new_sp.is_in_db():
-        if grass.overwrite() == False:
-	    dbif.close()
-            grass.fatal(_("Space time raster dataset <%s> is already in database, use overwrite flag to overwrite") % out_id)
-            
-    rows = sp.get_registered_maps("id", where, "start_time", dbif)
-
-    new_maps = {}
-    if rows:
-	num_rows = len(rows)
-	
-	grass.percent(0, num_rows, 1)
-	
-	# Run the r.mapcalc expression
-        if expression:
-	    count = 0
-	    proc_count = 0
-	    proc_list = []
-	    for row in rows:
-		count += 1
-		
-		grass.percent(count, num_rows, 1)
-		
-		map_name = "%s_%i" % (base, count)
-
-		expr = "%s = %s" % (map_name, expression.replace(sp.get_id(), row["id"]))
-		expr = expr.replace(sp.base.get_name(), row["id"])
-
-		map_id = map_name + "@" + mapset
-
-		new_map = sp.get_new_map_instance(map_id)
-
-		# Check if new map is in the temporal database
-		if new_map.is_in_db(dbif):
-		    if grass.overwrite() == True:
-			# Remove the existing temporal database entry
-			new_map.delete(dbif)
-			new_map = sp.get_new_map_instance(map_id)
-		    else:
-			grass.error(_("Raster map <%s> is already in temporal database, use overwrite flag to overwrite"))
-			continue
-
-		grass.verbose(_("Apply r.mapcalc expression: \"%s\"") % expr)
-		
-		proc_list.append(Process(target=run_mapcalc, args=(expr,)))
-		proc_list[proc_count].start()
-		proc_count += 1
-		
-		if proc_count == nprocs:
-		    proc_count = 0
-		    exitcodes = 0
-		    for proc in proc_list:
-			proc.join()
-			exitcodes += proc.exitcode
-			
-		    if exitcodes != 0:
-			dbif.close()
-			grass.fatal(_("Error while r.mapcalc computation"))
-			
-		    # Empty proc list
-		    proc_list = []
-		new_maps[row["id"]] = new_map
-	
-	grass.percent(0, num_rows, 1)
-	
-	# Insert the new space time raster dataset
-	if new_sp.is_in_db():
-	    if grass.overwrite() == True:
-		new_sp.delete(dbif)
-		new_sp = tgis.space_time_raster_dataset(out_id)
-		
-	temporal_type, semantic_type, title, description = sp.get_initial_values()
-	new_sp.set_initial_values(temporal_type, semantic_type, title, description)
-	new_sp.insert(dbif)
-	
-	# Register the maps in the database
-        count = 0
-        for row in rows:
-            count += 1
-	    
-	    grass.percent(count, num_rows, 1)
-
-            old_map = sp.get_new_map_instance(row["id"])
-            old_map.select(dbif)
-            
-            if expression:
-		
-		if new_maps.has_key(row["id"]):
-		    new_map = new_maps[row["id"]]
-
-                # Read the raster map data
-                new_map.load()
-                
-                # In case of a null map continue, do not register null maps
-                if new_map.metadata.get_min() == None and new_map.metadata.get_max() == None:
-                    if not register_null:
-                        continue
-
-                # Set the time stamp
-                if old_map.is_time_absolute():
-                    start, end, tz = old_map.get_absolute_time()
-                    new_map.set_absolute_time(start, end, tz)
-                else:
-                    start, end = old_map.get_relative_time()
-                    new_map.set_relative_time(start, end)
-
-                # Insert map in temporal database
-                new_map.insert(dbif)
-
-                new_sp.register_map(new_map, dbif)
-            else:
-                new_sp.register_map(old_map, dbif)          
-                
-        # Update the spatio-temporal extent and the raster metadata table entries
-        new_sp.update_from_registered_maps(dbif)
-	
-	grass.percent(num_rows, num_rows, 1)
-        
-    dbif.close()
-
-###############################################################################
-
-def run_mapcalc(expr):
-    """Helper function to run r.mapcalc in parallel"""
-    return grass.run_command("r.mapcalc", expression=expr, overwrite=grass.overwrite(), quiet=True)
+    tgis.extract_dataset(input, output, "raster", where, expression, base, nprocs, register_null)
     
 ###############################################################################
 

temporal/t.rast.extract/test.t.rast.extract.sh

@@ -2,8 +2,7 @@
 
 # We need to set a specific region in the
 # @preprocess step of this test. We generate
-# raster with r.mapcalc and create two space time raster inputs
-# with relative and absolute time
+# raster with r.mapcalc 
 # The region setting should work for UTM and LL test locations
 g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
 
@@ -17,15 +16,13 @@ r.mapcalc --o expr="prec_6 = rand(0, 650)"
 t.create --o type=strds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
 t.register -i type=rast input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="3 months"
 
-# The first @test
-# We create the space time raster inputs and register the raster maps with absolute time interval
-
+# The @test
 t.rast.extract --o --v input=precip_abs1 output=precip_abs2 where="start_time > '2001-06-01'" \
            expression=" if(precip_abs1 > 400, precip_abs1, null())" base=new_prec nprocs=2
 t.info type=strds input=precip_abs2
 
 t.rast.extract --o --v -n input=precip_abs1 output=precip_abs3 where="start_time > '2001-06-01'" \
-           expression=" if(precip_abs1 > 400, precip_abs1, null())" base=new_prec nprocs=4
+           expression=" if(precip_abs1@PERMANENT > 400, precip_abs1@PERMANENT, null())" base=new_prec nprocs=4
 t.info type=strds input=precip_abs3
 
 t.unregister type=rast maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6

temporal/t.rast.mapcalc/t.rast.mapcalc.py

@@ -87,273 +87,8 @@ def main():
     # Make sure the temporal database exists
     tgis.create_temporal_database()
     
-    # We need a database interface for fast computation
-    dbif = tgis.sql_database_interface_connection()
-    dbif.connect()
-
-    mapset =  grass.gisenv()["MAPSET"]
-    
-    input_name_list = inputs.split(",")
-    
-    # Process the first input
-    if input_name_list[0].find("@") >= 0:
-	id = input_name_list[0]
-    else:
-	id = input_name_list[0] + "@" + mapset
-	
-    first_input = tgis.space_time_raster_dataset(id)
-    
-    if first_input.is_in_db(dbif) == False:
-	dbif.close()
-        grass.fatal(_("Space time %s dataset <%s> not found") % (sp.get_new_map_instance(None).get_type(), id))
-
-    # Fill the object with data from the temporal database
-    first_input.select(dbif)
+    tgis.dataset_mapcalculator(inputs, output, "raster", expression, base, method, nprocs, register_null, spatial)
     
-    # All additional inputs in reverse sorted order to avoid wrong name substitution
-    input_name_list = input_name_list[1:]
-    input_name_list.sort()
-    input_name_list.reverse()
-    input_list = []
-        
-    for input in input_name_list:
-
-	if input.find("@") >= 0:
-	    id = input
-	else:
-	    id = input + "@" + mapset
-	    
-	sp = tgis.space_time_raster_dataset(id)
-	
-	if sp.is_in_db(dbif) == False:
-	    dbif.close()
-	    grass.fatal(_("Space time raster dataset <%s> not found in temporal database") % (id))
-
-	sp.select(dbif)
-	
-	input_list.append(copy.copy(sp))
-
-    # Create the new space time raster dataset
-    if output.find("@") >= 0:
-        out_id = output
-    else:
-        out_id = output + "@" + mapset
-        
-    new_sp = tgis.space_time_raster_dataset(out_id)
-    
-    # Check if in database
-    if new_sp.is_in_db(dbif):
-        if grass.overwrite() == False:
-	    dbif.close()
-            grass.fatal(_("Space time raster dataset <%s> is already in database, use overwrite flag to overwrite") % out_id)
- 
-    # Sample all inputs by the first input and create a sample matrix
-    if spatial:
-        grass.message(_("Start spatio-temporal sampling"))
-    else:
-        grass.message(_("Start temporal sampling"))
-    map_matrix = []
-    id_list = []
-    sample_map_list = []
-    # First entry is the first dataset id
-    id_list.append(first_input.get_name())
-    
-    if len(input_list) > 0:
-	has_samples = False
-	for dataset in input_list:
-	    list = dataset.sample_by_dataset(stds=first_input, method=method, spatial=spatial, dbif=dbif)
-	    
-	    # In case samples are not found
-	    if not list and len(list) == 0:
-		dbif.close()
-		grass.message(_("No samples found for map calculation"))
-		return 0
-	    
-	    # The fist entries are the samples
-	    map_name_list = []
-	    if has_samples == False:
-		for entry in list:
-		    granule = entry["granule"]
-		    # Do not consider gaps
-		    if granule.get_id() == None:
-			continue
-		    sample_map_list.append(granule)
-		    map_name_list.append(granule.get_name())
-		# Attach the map names
-		map_matrix.append(copy.copy(map_name_list))
-		has_samples = True
-		
-	    map_name_list = []
-	    for entry in list:
-		maplist = entry["samples"]
-		granule = entry["granule"]
-		
-		# Do not consider gaps in the sampler
-		if granule.get_id() == None:
-		    continue
-		
-		if len(maplist) > 1:
-		    grass.warning(_("Found more than a single map in a sample granule. "\
-		    "Only the first map is used for computation. "\
-		    "Use t.rast.aggregate.ds to create synchronous datasets."))
-		
-		# Store all maps! This includes non existent maps, identified by id == None 
-		map_name_list.append(maplist[0].get_name())
-	    
-	    # Attach the map names
-	    map_matrix.append(copy.copy(map_name_list))
-
-	    id_list.append(dataset.get_name())
-    else:
-	list = first_input.get_registered_maps_as_objects(dbif=dbif)
-	
-	if list == None:
-	    dbif.close()
-            grass.message(_("No maps in input dataset"))
-            return 0
-	
-	map_name_list = []
-	for map in list:
-	    map_name_list.append(map.get_name())
-	    sample_map_list.append(map)
-	
-	# Attach the map names
-	map_matrix.append(copy.copy(map_name_list))
-    
-   
-    # Needed for map registration
-    map_list = []
-	
-    if len(map_matrix) > 0:
-	
-	grass.message(_("Start r.mapcalc computation"))
-	    
-	count = 0
-	# Get the number of samples
-	num = len(map_matrix[0])
-	
-	# Parallel processing
-        proc_list = []
-        proc_count = 0
-	
-	# For all samples
-        for i in range(num):
-            
-            count += 1
-	    grass.percent(count, num, 1)
-
-	    # Create the r.mapcalc statement for the current time step
-	    map_name = "%s_%i" % (base, count)   
-	    expr = "%s = %s" % (map_name, expression)
-            
-            # Check that all maps are in the sample
-            valid_maps = True
-            # Replace all dataset names with their map names of the current time step
-            for j in range(len(map_matrix)):
-		if map_matrix[j][i] == None:
-		    valid_maps = False
-		    break
-		# Substitute the dataset name with the map name
-		expr = expr.replace(id_list[j], map_matrix[j][i])
-
-	    # Proceed with the next sample
-	    if valid_maps == False:
-		continue
-		
-	    # Create the new map id and check if the map is already in the database
-	    map_id = map_name + "@" + mapset
-
-	    new_map = first_input.get_new_map_instance(map_id)
-
-	    # Check if new map is in the temporal database
-	    if new_map.is_in_db(dbif):
-		if grass.overwrite() == True:
-		    # Remove the existing temporal database entry
-		    new_map.delete(dbif)
-		    new_map = first_input.get_new_map_instance(map_id)
-		else:
-		    grass.error(_("Raster map <%s> is already in temporal database, use overwrite flag to overwrite"))
-		    continue
-
-	    # Set the time stamp
-	    if sample_map_list[i].is_time_absolute():
-		start, end, tz = sample_map_list[i].get_absolute_time()
-		new_map.set_absolute_time(start, end, tz)
-	    else:
-		start, end = sample_map_list[i].get_relative_time()
-		new_map.set_relative_time(start, end)
-	    
-	    map_list.append(new_map)
-	    
-	    # Start the parallel r.mapcalc computation
-	    grass.verbose(_("Apply r.mapcalc expression: \"%s\"") % expr)
-	    
-	    proc_list.append(Process(target=run_mapcalc, args=(expr,)))
-	    proc_list[proc_count].start()
-	    proc_count += 1
-	    
-	    if proc_count == nprocs:
-		proc_count = 0
-		exitcodes = 0
-		for proc in proc_list:
-		    proc.join()
-		    exitcodes += proc.exitcode
-		    
-		if exitcodes != 0:
-		    dbif.close()
-		    grass.fatal(_("Error while r.mapcalc computation"))
-		    
-		# Empty process list
-		proc_list = []
-		
-	# Register the new maps in the output space time dataset
-	grass.message(_("Start map registration in temporal database"))
-	    
-	# Overwrite an existing dataset if requested
-	if new_sp.is_in_db(dbif):
-	    if grass.overwrite() == True:
-		new_sp.delete(dbif)
-		new_sp = tgis.space_time_raster_dataset(out_id)
-		
-	# Copy the ids from the first input
-	temporal_type, semantic_type, title, description = first_input.get_initial_values()
-	new_sp.set_initial_values(temporal_type, semantic_type, title, description)
-	# Insert the dataset in the temporal database
-	new_sp.insert(dbif)
-    
-	count = 0
-	
-	for new_map in map_list:
-
-            count += 1
-	    grass.percent(count, num, 1)
-	    
-	    # Read the raster map data
-	    new_map.load()
-	    
-	    # In case of a null map continue, do not register null maps
-	    if new_map.metadata.get_min() == None and new_map.metadata.get_max() == None:
-		if not register_null:
-		    continue
-
-	    # Insert map in temporal database
-	    new_map.insert(dbif)
-
-	    new_sp.register_map(new_map, dbif)
-
-        # Update the spatio-temporal extent and the raster metadata table entries
-        new_sp.update_from_registered_maps(dbif)
-	
-	grass.percent(1, 1, 1)
-        
-    dbif.close()
-
-###############################################################################
-
-def run_mapcalc(expr):
-    """Helper function to run r.mapcalc in parallel"""
-    return grass.run_command("r.mapcalc", expression=expr, overwrite=grass.overwrite(), quiet=True)
-
 ###############################################################################
 
 if __name__ == "__main__":

temporal/t.rast3d.extract/t.rast3d.extract.py

@@ -45,6 +45,15 @@
 #% multiple: no
 #%end
 
+#%option
+#% key: nprocs
+#% type: integer
+#% description: The number of r3.mapcalc processes to run in parallel
+#% required: no
+#% multiple: no
+#% answer: 1
+#%end
+
 #%flag
 #% key: n
 #% description: Register Null maps
@@ -63,125 +72,13 @@ def main():
     where = options["where"]
     expression = options["expression"]
     base = options["base"]
+    nprocs = int(options["nprocs"])
     register_null = flags["n"]
 
-    if expression and not base:
-        grass.fatal(_("Please specify %s=")%"base")
-
     # Make sure the temporal database exists
     tgis.create_temporal_database()
     
-    mapset =  grass.gisenv()["MAPSET"]
-
-    if input.find("@") >= 0:
-        id = input
-    else:
-        id = input + "@" + mapset
-
-    sp = tgis.space_time_raster3d_dataset(id)
-    
-    if sp.is_in_db() == False:
-        grass.fatal(_("Space time %s dataset <%s> not found") % (sp.get_new_map_instance(None).get_type(), id))
-
-    dbif = tgis.sql_database_interface_connection()
-    dbif.connect()
-
-    sp.select(dbif)
-
-    if output.find("@") >= 0:
-        out_id = output
-    else:
-        out_id = output + "@" + mapset
-
-    # The new space time raster3d dataset
-    new_sp = tgis.space_time_raster3d_dataset(out_id)
-    if new_sp.is_in_db():
-        if grass.overwrite() == True:
-            new_sp.delete(dbif)
-	    new_sp = tgis.space_time_raster3d_dataset(out_id)
-        else:
-            grass.fatal(_("Space time raster3d dataset <%s> is already in database, use overwrite flag to overwrite") % out_id)
-
-    temporal_type, semantic_type, title, description = sp.get_initial_values()
-    new_sp.set_initial_values(temporal_type, semantic_type, title, description)
-    new_sp.insert(dbif)
-
-    rows = sp.get_registered_maps("id", where, "start_time", dbif)
-
-    if rows:
-	num_rows = len(rows)
-	
-	grass.percent(0, num_rows, 1)
-	
-        count = 0
-        for row in rows:
-            count += 1
-	    
-	    grass.percent(count, num_rows, 1)
-
-            old_map = sp.get_new_map_instance(row["id"])
-            old_map.select(dbif)
-            
-            if expression:
-
-                map_name = "%s_%i" % (base, count)
-
-                expr = "%s = %s" % (map_name, expression.replace(sp.get_id(), row["id"]))
-                expr = expr.replace(sp.base.get_name(), row["id"])
-
-                map_id = map_name + "@" + mapset
-
-                new_map = sp.get_new_map_instance(map_id)
-
-                # Check if new map is in the temporal database
-                if new_map.is_in_db(dbif):
-                    if grass.overwrite() == True:
-                        # Remove the existing temporal database entry
-                        new_map.delete(dbif)
-                        new_map = sp.get_new_map_instance(map_id)
-                    else:
-                        grass.error(_("Raster3d map <%s> is already in temporal database, use overwrite flag to overwrite"))
-                        continue
-
-                grass.verbose(_("Apply r3.mapcalc expression: \"%s\"") % expr)
-
-                ret = grass.run_command("r3.mapcalc", expression=expr, overwrite=grass.overwrite(), quiet=True)
-
-                if ret != 0:
-                    grass.error(_("Error while r3.mapcalc computation, continue with next map"))
-                    break
-
-                # Read the raster3d map data
-                new_map.load()
-                
-                # In case of a null map continue, do not register null maps
-                print new_map.metadata.get_min(), new_map.metadata.get_max()
-                if new_map.metadata.get_min() == None and new_map.metadata.get_max() == None:
-                    print("Found null map")
-                    if not register_null:
-                        continue
-
-                # Set the time stamp
-                if old_map.is_time_absolute():
-                    start, end, tz = old_map.get_absolute_time()
-                    new_map.set_absolute_time(start, end, tz)
-                else:
-                    start, end = old_map.get_relative_time()
-                    new_map.set_relative_time(start, end)
-
-                # Insert map in temporal database
-                new_map.insert(dbif)
-
-                new_sp.register_map(new_map, dbif)
-            else:
-                new_sp.register_map(old_map, dbif)
-
-        # Update the spatio-temporal extent and the raster3d metadata table entries
-        new_sp.update_from_registered_maps(dbif)
-	
-	grass.percent(num_rows, num_rows, 1)
-        
-    dbif.close()
+    tgis.extract_dataset(input, output, "raster3d", where, expression, base, nprocs, register_null)
 
 if __name__ == "__main__":
     options, flags = grass.parser()

temporal/t.rast3d.extract/test.t.rast3d.extract.sh

@@ -1,9 +1,8 @@
-# Test the extraction of a subset of a space time raster input
+# Test the extraction of a subset of a space time raster3d input
 
 # We need to set a specific region in the
 # @preprocess step of this test. We generate
-# raster with r.mapcalc and create two space time raster inputs
-# with relative and absolute time
+# raster with r.mapcalc 
 # The region setting should work for UTM and LL test locations
 g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
 
@@ -17,9 +16,7 @@ r3.mapcalc --o expr="prec_6 = rand(0, 650)"
 t.create --o type=str3ds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
 t.register -i type=rast3d input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="3 months"
 
-# The first @test
-# We create the space time raster inputs and register the raster maps with absolute time interval
-
+# The @test
 t.rast3d.extract --o input=precip_abs1 output=precip_abs2 where="start_time > '2001-06-01'" expression=" if(precip_abs1 > 400, precip_abs1, null())" base=new_prec
 
 t.info type=str3ds input=precip_abs2

temporal/t.rast3d.mapcalc/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../../
+
+PGM = t.rast3d.mapcalc
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script $(TEST_DST)

temporal/t.rast3d.mapcalc/t.rast3d.mapcalc.py

@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+############################################################################
+#
+# MODULE:	t.rast3d.mapcalc
+# AUTHOR(S):	Soeren Gebbert
+#
+# PURPOSE:	Perform r3.mapcalc expressions on maps of sampled space time raster3d datasets
+# COPYRIGHT:	(C) 2012 by the GRASS Development Team
+#
+#		This program is free software under the GNU General Public
+#		License (version 2). Read the file COPYING that comes with GRASS
+#		for details.
+#
+#############################################################################
+
+#%module
+#% description: Perform r3.mapcalc expressions on maps of sampled space time raster3d datasets
+#% keywords: temporal
+#% keywords: mapcalc
+#%end
+
+#%option G_OPT_STR3DS_INPUTS
+#%end
+
+#%option
+#% key: expression
+#% type: string
+#% description: The r3.mapcalc expression applied to each time step of the sampled data
+#% required: yes
+#% multiple: no
+#%end
+
+#%option G_OPT_T_SAMPLE
+#% key: method
+#% answer: during,overlap,contain,equal
+#%end
+
+#%option G_OPT_STR3DS_OUTPUT
+#%end
+
+#%option G_OPT_R_BASE
+#% required: yes
+#% description: Name of the base raster3d map
+#% gisprompt: old,grid3,3d-raster
+#%end
+
+#%option
+#% key: nprocs
+#% type: integer
+#% description: The number of r3.mapcalc processes to run in parallel
+#% required: no
+#% multiple: no
+#% answer: 1
+#%end
+
+#%flag
+#% key: n
+#% description: Register Null maps
+#%end
+
+#%flag
+#% key: s
+#% description: Check spatial overlap
+#%end
+
+from multiprocessing import Process
+import copy
+import grass.script as grass
+import grass.temporal as tgis
+
+############################################################################
+
+def main():
+
+    # Get the options
+    inputs = options["inputs"]
+    output = options["output"]
+    expression = options["expression"]
+    base = options["base"]
+    method = options["method"]
+    nprocs = int(options["nprocs"])
+    register_null = flags["n"]
+    spatial = flags["s"]
+    
+    # Create the method list
+    method = method.split(",")
+        
+    # Make sure the temporal database exists
+    tgis.create_temporal_database()
+    
+    tgis.dataset_mapcalculator(inputs, output, "raster3d", expression, base, method, nprocs, register_null, spatial)
+    
+###############################################################################
+
+if __name__ == "__main__":
+    options, flags = grass.parser()
+    main()
+

temporal/t.rast3d.mapcalc/test.t.rast3d.mapcalc.sh

@@ -0,0 +1,37 @@
+# Test for t.rast.mapcalc 
+
+# We need to set a specific region in the
+# @preprocess step of this test. We generate
+# raster with r.mapcalc and create several space time raster datasets
+# The region setting should work for UTM and LL test locations
+g.region s=0 n=80 w=0 e=120 b=0 t=50 res=10 res3=10 -p3
+
+r3.mapcalc --o expr="prec_1 = rand(0, 550)"
+r3.mapcalc --o expr="prec_2 = rand(0, 450)"
+r3.mapcalc --o expr="prec_3 = rand(0, 320)"
+r3.mapcalc --o expr="prec_4 = rand(0, 510)"
+r3.mapcalc --o expr="prec_5 = rand(0, 300)"
+r3.mapcalc --o expr="prec_6 = rand(0, 650)"
+
+t.create --o type=str3ds temporaltype=absolute output=precip_abs1 title="A test" descr="A test"
+t.create --o type=str3ds temporaltype=absolute output=precip_abs2 title="A test" descr="A test"
+t.register -i type=rast3d input=precip_abs1 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6 start="2001-01-01" increment="3 months"
+t.register type=rast3d input=precip_abs2 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+
+t.info type=str3ds input=precip_abs1
+t.info type=str3ds input=precip_abs2
+
+# The first @test
+t.rast3d.mapcalc --o --v -n inputs=precip_abs1,precip_abs2 output=precip_abs3 \
+           expression=" precip_abs1 + precip_abs2" base=new_prec \
+           method=equal nprocs=6
+t.info type=str3ds input=precip_abs3
+
+t.rast3d.mapcalc --o --v -s inputs=precip_abs1,precip_abs2,precip_abs3 output=precip_abs4 \
+           expression=" (precip_abs1 + precip_abs2) / precip_abs3@PERMANENT " base=new_prec \
+           method=equal nprocs=6
+t.info type=str3ds input=precip_abs4
+
+t.unregister type=rast3d maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.unregister type=rast3d maps=new_prec_1,new_prec_2,new_prec_3,new_prec_4,new_prec_5,new_prec_6
+t.remove type=str3ds input=precip_abs1,precip_abs2,precip_abs3,precip_abs4