t.rast.aggregate: Adjustments to temporally weighted aggregation

python/grass/temporal/aggregation.py

@@ -224,6 +224,7 @@ def aggregate_by_topology(
     nprocs=1,
     spatial=None,
     dbif=None,
+    weighting=False,
     overwrite=False,
     file_limit=1000,
 ):
@@ -292,7 +293,6 @@ def aggregate_by_topology(
         count += 1
 
         aggregation_list = []
-
         if "equal" in topo_list and granule.equal:
             for map_layer in granule.equal:
                 aggregation_list.append(map_layer.get_name())
@@ -320,7 +320,38 @@ def aggregate_by_topology(
         if "overlapped" in topo_list and granule.overlapped:
             for map_layer in granule.overlapped:
                 aggregation_list.append(map_layer.get_name())
-
+        if "related" in topo_list:
+            aggregation_weights = []
+            set_list = set()
+            if granule.overlaps:
+                set_list.update(granule.overlaps)
+            if granule.overlapped:
+                set_list.update(granule.overlapped)
+            if granule.contains:
+                set_list.update(granule.contains)
+            if granule.equal:
+                set_list.update(granule.equal)
+            if granule.during:
+                set_list.update(granule.during)
+            if len(set_list) > 0:
+                for map_layer in set_list:
+                    aggregation_list.append(map_layer.get_name())
+                    t_granule_contained = map_layer.get_absolute_time()
+                    t_granule = granule.get_absolute_time()
+                    if None in t_granule_contained or None in t_granule:
+                        # no weight for this map_layer because no
+                        # overlap whatsoever
+                        aggregation_weights.append(0)
+                    else:
+                        # calculate the absolute temporal overlap between the
+                        # new granule and the map_layer
+                        overlap_abs = min(t_granule[1], t_granule_contained[1]) - max(
+                            t_granule[0], t_granule_contained[0]
+                        )
+                        # calculate the relative percentage of the overlap
+                        # with respect to the total granule duration
+                        overlap_rel = overlap_abs / (t_granule[1] - t_granule[0])
+                        aggregation_weights.append(overlap_rel)
         if aggregation_list:
             msgr.verbose(
                 _("Aggregating %(len)i raster maps from %(start)s to %(end)s")
@@ -364,12 +395,15 @@ def aggregate_by_topology(
             if len(aggregation_list) > 1:
                 # Create the r.series input file
                 filename = gs.tempfile(True)
-                file = open(filename, "w")
-                for name in aggregation_list:
-                    string = "%s\n" % (name)
-                    file.write(string)
-                file.close()
-
+                with open(filename, "w") as file:
+                    if weighting:
+                        for name, weight in zip(aggregation_list, aggregation_weights):
+                            file.write(f"{name}|{weight}\n")
+                    else:
+                        for name in aggregation_list:
+                            string = "%s\n" % (name)
+                            file.write(string)
+                # Perform aggregation
                 mod = copy.deepcopy(r_series)
                 mod(file=filename, output=output_name)
                 if len(aggregation_list) > int(file_limit):
@@ -385,8 +419,16 @@ def aggregate_by_topology(
                     mod(flags="z")
                 process_queue.put(mod)
             else:
-                mod = copy.deepcopy(g_copy)
-                mod(raster=[aggregation_list[0], output_name])
+                if weighting:
+                    mod = copy.deepcopy(r_series)
+                    mod(
+                        input=aggregation_list[0],
+                        output=output_name,
+                        weights=aggregation_weights[0],
+                    )
+                else:
+                    mod = copy.deepcopy(g_copy)
+                    mod(raster=[aggregation_list[0], output_name])
                 process_queue.put(mod)
 
     process_queue.wait()

temporal/t.rast.aggregate/t.rast.aggregate.html

@@ -193,6 +193,41 @@ <h4>Yearly aggregation</h4>
 yearly_avg_temp_2004|climate|2004-01-01 00:00:00|2005-01-01 00:00:00
 </pre></div>
 
+<h4>Weighted aggregation</h4>
+
+In this example, we create a STRDS of fictional temperature values for
+a weekly interval:
+
+<div class="code"><pre>
+MAPS="map_1 map_2 map_3 map_4 map_5 map_6 map_7 map_8 map_9 map_10"
+for map in ${MAPS} ; do
+    r.surf.random output=${map} min=278 max=298
+    echo ${map} >> map_list.txt
+done
+
+t.create type=strds temporaltype=absolute \
+         output=temperature_weekly \
+         title="Weekly Temperature" \
+         description="Test dataset with weekly temperature"
+
+t.register -i type=raster input=temperature_weekly \
+          file=map_list.txt start="2021-05-01" increment="1 weeks"
+</pre></div>
+
+We can now use the <b>-w</b> flag and the <b>sampling=related</b> option
+to calculate the 10-daily average temperature. The values of each 10-days
+granule are calculated weighted by relative temporal coverage of the input
+granules.
+
+<div class="code"><pre>
+t.rast.aggregate input=temperature_weekly output=temperature_10daily_weighted \
+  basename=tendaily_temp_weighted method=average granularity="10 days" \
+  sampling=related -w
+</pre></div>
+
+This is especially useful when harmonizing STRDS with granules that are not
+integer multiplications of each other.
+
 <h2>SEE ALSO</h2>
 
 <em>

temporal/t.rast.aggregate/t.rast.aggregate.py

@@ -99,7 +99,7 @@
 # %end
 
 # %option G_OPT_T_SAMPLE
-# % options: equal,overlaps,overlapped,starts,started,finishes,finished,during,contains
+# % options: equal,overlaps,overlapped,starts,started,finishes,finished,during,contains,related
 # % answer: contains
 # %end
 
@@ -111,6 +111,11 @@
 # % description: Register Null maps
 # %end
 
+# %flag
+# % key: w
+# % description: Aggregation weighted by temporal overlap between input rasters and rasters of defined granularity
+# %end
+
 import grass.script as gs
 
 ############################################################################
@@ -127,13 +132,17 @@ def main():
     gran = options["granularity"]
     base = options["basename"]
     register_null = flags["n"]
+    weighting = flags["w"]
     method = options["method"]
     sampling = options["sampling"]
     offset = options["offset"]
     nprocs = options["nprocs"]
     file_limit = options["file_limit"]
     time_suffix = options["suffix"]
 
+    if weighting and sampling != "related":
+        gs.fatal(_("Weighting only works with sampling: 'related'"))
+
     topo_list = sampling.split(",")
 
     tgis.init()
@@ -202,6 +211,7 @@ def main():
         method=method,
         nprocs=nprocs,
         spatial=None,
+        weighting=weighting,
         overwrite=gs.overwrite(),
         file_limit=file_limit,
     )

temporal/t.rast.aggregate/testsuite/test_aggregation_absolute.py

@@ -57,11 +57,11 @@ def setUpClass(cls):
     def tearDownClass(cls):
         """Remove the temporary region"""
         cls.del_temp_region()
-        cls.runModule("t.remove", flags="df", type="strds", inputs="A")
+        cls.runModule("t.remove", flags="rf", type="strds", inputs="A")
 
     def tearDown(self):
         """Remove generated data"""
-        self.runModule("t.remove", flags="df", type="strds", inputs="B")
+        self.runModule("t.remove", flags="rf", type="strds", inputs="B")
 
     def test_disaggregation(self):
         """Disaggregation with empty maps"""
@@ -238,6 +238,43 @@ def test_aggregation_3months(self):
         maps = "b_101" + os.linesep
         self.assertEqual(maps, lister.outputs.stdout)
 
+    def test_weighted_aggregation(self):
+        """Weighted aggregation to 3 weeks"""
+        self.assertModule(
+            "t.rast.aggregate",
+            input="A",
+            output="B",
+            basename="b",
+            granularity="3 weeks",
+            method="average",
+            sampling=["related"],
+            file_limit=0,
+            suffix="num%03",
+            flags="w",
+        )
+        # assert that there are 5 rasters
+        t_rast_list = SimpleModule("t.rast.list", input="B", flags="u").run()
+        rasters = [
+            item
+            for item in t_rast_list.outputs.stdout.split(os.linesep)
+            if item.startswith("b_")
+        ]
+        self.assertEqual(
+            5,
+            len(rasters),
+            ("Output STRDS does not contain the correct number of rasters."),
+        )
+
+        # assert that the granularity is correct and the whole time range is
+        # covered
+        info = SimpleModule("t.info", flags="g", input="B")
+        ref_dict = {
+            "start_time": "'2001-01-15 00:00:00'",
+            "end_time": "'2001-04-30 00:00:00'",
+            "granularity": "'21 days'",
+        }
+        self.assertModuleKeyValue(module=info, reference=ref_dict, precision=2, sep="=")
+
 
 if __name__ == "__main__":
     from grass.gunittest.main import test