Introduction of temporal operators in t.rast.mapcalc

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

lib/python/temporal/mapcalc.py

@@ -24,9 +24,41 @@ def dataset_mapcalculator(inputs, output, type, expression, base, method,
        raster/raster3d datasets, using a specific sampling method 
        to select temporal related maps.
 
-       A mapcalc expression can be provided to process the temporal 
-       extracted maps.
-       Mapcalc expressions are supported for raster and raster3d maps.
+       A mapcalc expression must be provided to process the temporal 
+       selected maps. Temporal operators are available in addition to 
+       the r.mapcalc operators:
+       
+       Supported operators for relative and absolute time:
+       * td() - the time delta of the current interval in days 
+         and fractions of days or the unit in case of relative time
+       * start_time() - The start time of the interval from the begin of the time 
+         series in days and fractions of days or the unit in case of relative time
+       * end_time() - The end time of the current interval from the begin of the time 
+         series in days and fractions of days or the unit in case of relative time
+         
+       Supported operators for absolute time:
+       * start_doy() - Day of year (doy) from the start time [1 - 366]
+       * start_dow() - Day of week (dow) from the start time [1 - 7], 
+         the start of the week is monday == 1
+       * start_year() - The year of the start time [0 - 9999]
+       * start_month() - The month of the start time [1 - 12]
+       * start_week() - Week of year of the start time [1 - 54]
+       * start_day() - Day of month from the start time [1 - 31]
+       * start_hour() - The hour of the start time [0 - 23]
+       * start_minute() - The minute of the start time [0 - 59]
+       * start_second() - The second of the start time [0 - 59]
+       
+       * end_doy() - Day of year (doy) from the start time [1 - 366]
+       * end_dow() - Day of week (dow) from the end time [1 - 7], 
+         end is monday == 1
+       * end_year() - The year of the end time [0 - 9999]
+       * end_month() - The month of the end time [1 - 12]
+       * end_woy() - Week of year (woy) of the end time [1 - 54]
+       * end_day() - Day of month from the start time [1 - 31]
+       * end_hour() - The hour of the end time [0 - 23]
+       * end_minute() - The minute of the end time [0 - 59]
+       * end_second() - The second of the end time [0 - 59]
+       
 
        @param input The name of the input space time raster/raster3d dataset
        @param output The name of the extracted new space time raster(3d) dataset
@@ -207,8 +239,10 @@ def dataset_mapcalculator(inputs, output, type, expression, base, method,
 
             # Create the r.mapcalc statement for the current time step
             map_name = "%s_%i" % (base, count)
-            expr = "%s = %s" % (map_name, expression)
-
+            expr = "%s=%s" % (map_name, expression)
+            # Remove spaces and new lines
+            expr = expr.replace(" ", "")
+            
             # Check that all maps are in the sample
             valid_maps = True
             # Replace all dataset names with their map names of the 
@@ -248,6 +282,9 @@ def dataset_mapcalculator(inputs, output, type, expression, base, method,
             else:
                 start, end, unit = sample_map_list[i].get_relative_time()
                 new_map.set_relative_time(start, end, unit)
+            
+            # Parse the temporal expressions
+            expr = _operator_parser(expr,  sample_map_list[0],  sample_map_list[i])
 
             map_list.append(new_map)
 
@@ -255,9 +292,9 @@ def dataset_mapcalculator(inputs, output, type, expression, base, method,
             core.verbose(_("Apply mapcalc expression: \"%s\"") % expr)
 
             if type == "raster":
-                proc_list.append(Process(target=run_mapcalc2d, args=(expr,)))
+                proc_list.append(Process(target=_run_mapcalc2d, args=(expr,)))
             else:
-                proc_list.append(Process(target=run_mapcalc3d, args=(expr,)))
+                proc_list.append(Process(target=_run_mapcalc3d, args=(expr,)))
             proc_list[proc_count].start()
             proc_count += 1
 
@@ -342,15 +379,274 @@ def dataset_mapcalculator(inputs, output, type, expression, base, method,
 
 ###############################################################################
 
-def run_mapcalc2d(expr):
+def _run_mapcalc2d(expr):
     """Helper function to run r.mapcalc in parallel"""
     return core.run_command("r.mapcalc", expression=expr, 
                             overwrite=core.overwrite(), quiet=True)
 
 ###############################################################################
 
-
-def run_mapcalc3d(expr):
+def _run_mapcalc3d(expr):
     """Helper function to run r3.mapcalc in parallel"""
     return core.run_command("r3.mapcalc", expression=expr, 
                             overwrite=core.overwrite(), quiet=True)
+
+###############################################################################
+
+def _operator_parser(expr, first, current):
+    """This method parses the expression string and substitutes 
+       the temporal operators with numerical values.
+       
+       Supported operators for relative and absolute time are:
+       * td() - the time delta of the current interval in days 
+         and fractions of days or the unit in case of relative time
+       * start_time() - The start time of the interval from the begin of the time 
+         series in days and fractions of days or the unit in case of relative time
+       * end_time() - The end time of the current interval from the begin of the time 
+         series in days and fractions of days or the unit in case of relative time
+         
+       Supported operators for absolute time:
+       * start_doy() - Day of year (doy) from the start time [1 - 366]
+       * start_dow() - Day of week (dow) from the start time [1 - 7], 
+         the start of the week is monday == 1
+       * start_year() - The year of the start time [0 - 9999]
+       * start_month() - The month of the start time [1 - 12]
+       * start_week() - Week of year of the start time [1 - 54]
+       * start_day() - Day of month from the start time [1 - 31]
+       * start_hour() - The hour of the start time [0 - 23]
+       * start_minute() - The minute of the start time [0 - 59]
+       * start_second() - The second of the start time [0 - 59]
+       
+       * end_doy() - Day of year (doy) from the start time [1 - 366]
+       * end_dow() - Day of week (dow) from the end time [1 - 7], 
+         end is monday == 1
+       * end_year() - The year of the end time [0 - 9999]
+       * end_month() - The month of the end time [1 - 12]
+       * end_woy() - Week of year (woy) of the end time [1 - 54]
+       * end_day() - Day of month from the start time [1 - 31]
+       * end_hour() - The hour of the end time [0 - 23]
+       * end_minute() - The minute of the end time [0 - 59]
+       * end_second() - The second of the end time [0 - 59]
+       
+       The modified expression is returned.
+       
+    """
+    is_time_absolute = first.is_time_absolute()
+    
+    expr = _parse_td_operator(expr, is_time_absolute, first, current)
+    expr = _parse_start_time_operator(expr, is_time_absolute, first, current)
+    expr = _parse_end_time_operator(expr, is_time_absolute, first, current)
+    expr = _parse_start_operators(expr, is_time_absolute, current)
+    expr = _parse_end_operators(expr, is_time_absolute, current)
+    
+    return expr
+
+###############################################################################
+
+def _parse_start_operators(expr, is_time_absolute, current):
+    """
+       Supported operators for absolute time:
+       * start_doy() - Day of year (doy) from the start time [1 - 366]
+       * start_dow() - Day of week (dow) from the start time [1 - 7], 
+         the start of the week is monday == 1
+       * start_year() - The year of the start time [0 - 9999]
+       * start_month() - The month of the start time [1 - 12]
+       * start_week() - Week of year of the start time [1 - 54]
+       * start_day() - Day of month from the start time [1 - 31]
+       * start_hour() - The hour of the start time [0 - 23]
+       * start_minute() - The minute of the start time [0 - 59]
+       * start_second() - The second of the start time [0 - 59]
+    """
+    
+    if not is_time_absolute:
+        core.fatal(_("The temporal operators <%s> supports only absolute time."%("start_*")))
+    
+    start, end, tz = current.get_absolute_time()
+        
+    if expr.find("start_year()") >= 0:
+        expr = expr.replace("start_year()", str(start.year))
+        
+    if expr.find("start_month()") >= 0:
+        expr = expr.replace("start_month()", str(start.month))
+        
+    if expr.find("start_week()") >= 0:
+        expr = expr.replace("start_week()", str(start.isocalendar()[1]))
+        
+    if expr.find("start_day()") >= 0:
+        expr = expr.replace("start_day()", str(start.day))
+        
+    if expr.find("start_hour()") >= 0:
+        expr = expr.replace("start_hour()", str(start.hour))
+        
+    if expr.find("start_minute()") >= 0:
+        expr = expr.replace("start_minute()", str(start.minute))
+        
+    if expr.find("start_second()") >= 0:
+        expr = expr.replace("start_second()", str(start.second))
+        
+    if expr.find("start_dow()") >= 0:
+        expr = expr.replace("start_dow()", str(start.isoweekday()))
+        
+    if expr.find("start_doy()") >= 0:        
+        year = datetime(start.year, 1, 1)
+        delta = start - year
+        
+        expr = expr.replace("start_doy()", str(delta.days + 1))
+        
+    return expr
+
+###############################################################################
+
+def _parse_end_operators(expr, is_time_absolute, current):
+    """
+       Supported operators for absolute time:
+       * end_doy() - Day of year (doy) from the end time [1 - 366]
+       * end_dow() - Day of week (dow) from the end time [1 - 7], 
+         the end of the week is monday == 1
+       * end_year() - The year of the end time [0 - 9999]
+       * end_month() - The month of the end time [1 - 12]
+       * end_week() - Week of year of the end time [1 - 54]
+       * end_day() - Day of month from the end time [1 - 31]
+       * end_hour() - The hour of the end time [0 - 23]
+       * end_minute() - The minute of the end time [0 - 59]
+       * end_second() - The minute of the end time [0 - 59]
+       
+       In case of time instances the end_* expression will be replaced by null()
+    """
+    
+    if not is_time_absolute:
+        core.fatal(_("The temporal operators <%s> supports only absolute time."%("end_*")))
+    
+    start, end, tz = current.get_absolute_time()
+        
+    if expr.find("end_year()") >= 0:
+        if not end:
+            expr = expr.replace("end_year()", "null()")
+        else:
+            expr = expr.replace("end_year()", str(end.year))
+        
+    if expr.find("end_month()") >= 0:
+        if not end:
+            expr = expr.replace("end_month()", "null()")
+        else:
+            expr = expr.replace("end_month()", str(end.month))
+        
+    if expr.find("end_week()") >= 0:
+        if not end:
+            expr = expr.replace("end_week()", "null()")
+        else:
+            expr = expr.replace("end_week()", str(end.isocalendar()[1]))
+        
+    if expr.find("end_day()") >= 0:
+        if not end:
+            expr = expr.replace("end_day()", "null()")
+        else:
+            expr = expr.replace("end_day()", str(end.day))
+        
+    if expr.find("end_hour()") >= 0:
+        if not end:
+            expr = expr.replace("end_hour()", "null()")
+        else:
+            expr = expr.replace("end_hour()", str(end.hour))
+        
+    if expr.find("end_minute()") >= 0:
+        if not end:
+            expr = expr.replace("end_minute()", "null()")
+        else:
+            expr = expr.replace("end_minute()", str(end.minute))
+        
+    if expr.find("end_second()") >= 0:
+        if not end:
+            expr = expr.replace("end_second()", "null()")
+        else:
+            expr = expr.replace("end_second()", str(end.second))
+        
+    if expr.find("end_dow()") >= 0:
+        if not end:
+            expr = expr.replace("end_dow()", "null()")
+        else:
+            expr = expr.replace("end_dow()", str(end.isoweekday()))
+        
+    if expr.find("end_doy()") >= 0:
+        if not end:
+            expr = expr.replace("end_doy()", "null()")
+        else:
+            year = datetime(end.year, 1, 1)
+            delta = end - year
+            
+            expr = expr.replace("end_doy()", str(delta.days + 1))
+        
+    return expr
+
+###############################################################################
+
+def _parse_td_operator(expr, is_time_absolute, first, current):
+    """Parse the time delta operator td(). This operator
+       represents the size of the current time interval
+       in days and fraction of days for absolute time,
+       and in relative units in case of relative time.
+       
+       In case of time instances, the td() operator will be of type null(). 
+    """
+    if expr.find("td()") >= 0:
+        td = "null()"
+        if is_time_absolute:
+            start, end, tz = current.get_absolute_time()
+            if end != None:
+                td = time_delta_to_relative_time(end - start)
+        else:
+            start, end, unit = current.get_relative_time()
+            if end != None:
+                td = end - start
+        expr = expr.replace("td()", str(td))
+        
+    return expr
+
+###############################################################################
+
+def _parse_start_time_operator(expr, is_time_absolute, first, current):
+    """Parse the start_time() operator. This operator represent 
+    the time difference between the start time of the sample space time
+    raster dataset and the start time of the current interval or instance. The time
+    is measured  in days and fraction of days for absolute time,
+    and in relative units in case of relative time."""
+    if expr.find("start_time()") >= 0:
+        if is_time_absolute:
+            start1, end, tz = first.get_absolute_time()
+            start, end, tz = current.get_absolute_time()
+            x = time_delta_to_relative_time(start - start1)
+        else:
+            start1, end, unit = first.get_relative_time()
+            start, end, unit = current.get_relative_time()
+            x = start - start1
+        expr = expr.replace("start_time()", str(x))
+        
+    return expr
+
+###############################################################################
+
+def _parse_end_time_operator(expr, is_time_absolute, first, current):
+    """Parse the end_time() operator. This operator represent 
+    the time difference between the start time of the sample space time
+    raster dataset and the start time of the current interval. The time
+    is measured  in days and fraction of days for absolute time,
+    and in relative units in case of relative time.
+    
+    The end_time() will be represented by null() in case of a time instance.
+    """
+    if expr.find("end_time()") >= 0:
+        x = "null()"
+        if is_time_absolute:
+            start1, end, tz = first.get_absolute_time()
+            start, end, tz = current.get_absolute_time()
+            if end:
+                x = time_delta_to_relative_time(end - start1)
+        else:
+            start1, end, unit = first.get_relative_time()
+            start, end, unit = current.get_relative_time()
+            if end:
+                x = end - start1
+        expr = expr.replace("end_time()", str(x))
+        
+    return expr
+

lib/python/temporal/space_time_datasets_tools.py

@@ -666,7 +666,7 @@ def sample_stds_by_stds_topology(intype, sampletype, inputs, sampler, header,
 
     for st in sts:
         if st.is_in_db(dbif) == False:
-            core.fatal(_("Dataset <%s> not found in temporal database") % (id))
+            core.fatal(_("Dataset <%s> not found in temporal database") % (st.get_id()))
         st.select(dbif)
 
     if sst.is_in_db(dbif) == False:
@@ -868,6 +868,9 @@ def create_space_time_dataset(name, type, temporaltype, title, descr, semantic,
 
     mapset = core.gisenv()["MAPSET"]
     id = name + "@" + mapset
+    
+    print id
+    print overwrite
 
     sp = dataset_factory(type, id)
 

temporal/t.rast.mapcalc/test.t.rast.mapcalc.operators.sh

@@ -0,0 +1,57 @@
+#!/bin/sh
+# 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
+
+r.mapcalc --o expr="prec_1 = rand(0, 550)"
+r.mapcalc --o expr="prec_2 = rand(0, 450)"
+r.mapcalc --o expr="prec_3 = rand(0, 320)"
+r.mapcalc --o expr="prec_4 = rand(0, 510)"
+r.mapcalc --o expr="prec_5 = rand(0, 300)"
+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.create --o type=strds temporaltype=absolute output=precip_abs2 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"
+t.register --o type=rast input=precip_abs2 maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+
+t.info precip_abs1
+t.info precip_abs2
+
+# The first @test
+t.rast.mapcalc --o --v -sn inputs=precip_abs1,precip_abs2 output=precip_abs3 \
+    expression="if(start_time() >= 0 && end_time() >= 0, (precip_abs1*td() + precip_abs2) / td(), null()) " base=new_prec \
+           method=equal nprocs=5
+t.info type=strds input=precip_abs3
+
+t.rast.mapcalc --o --v -sn inputs=precip_abs1,precip_abs2 output=precip_abs3 \
+    expression="if(start_year()>=0&&start_month()>=0&&start_day()>=0&&start_hour()>=0&&start_minute()>=0&&start_second()>=0, (precip_abs1*td() + precip_abs2) / td(), null()) " base=new_prec \
+           method=equal nprocs=5
+t.info type=strds input=precip_abs3
+
+t.rast.mapcalc --o --v -sn inputs=precip_abs1,precip_abs2 output=precip_abs3 \
+    expression="if(end_year()>=0&&end_month()>=0&&end_day()>=0&&end_hour()>=0&&end_minute()>=0&&end_second()>=0, (precip_abs1*td() + precip_abs2) / td(), null()) " base=new_prec \
+           method=equal nprocs=5
+t.info type=strds input=precip_abs3
+
+# The first @test
+t.rast.mapcalc --o --v -sn inputs=precip_abs1,precip_abs2 output=precip_abs3 \
+    expression="if(start_doy() >= 0 && start_dow() >= 0, (precip_abs1*td() + precip_abs2) / td(), null()) " base=new_prec \
+           method=equal nprocs=5
+t.info type=strds input=precip_abs3
+
+# The first @test
+t.rast.mapcalc --o --v -sn inputs=precip_abs1,precip_abs2 output=precip_abs3 \
+    expression="if(end_doy() >= 0 && end_dow() >= 0, (precip_abs1*td() + precip_abs2) / td(), null()) " base=new_prec \
+           method=equal nprocs=5
+t.info type=strds input=precip_abs3
+
+
+# @postprocess
+t.unregister type=rast maps=prec_1,prec_2,prec_3,prec_4,prec_5,prec_6
+t.unregister type=rast maps=new_prec_1,new_prec_2,new_prec_3,new_prec_4,new_prec_5,new_prec_6
+t.remove type=strds input=precip_abs1,precip_abs2,precip_abs3