diff --git a/python/grass/temporal/aggregation.py b/python/grass/temporal/aggregation.py index 354600855e0..f21611301bb 100644 --- a/python/grass/temporal/aggregation.py +++ b/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() diff --git a/temporal/t.rast.aggregate/t.rast.aggregate.html b/temporal/t.rast.aggregate/t.rast.aggregate.html index 55fbd963fa8..2fe33440d0a 100644 --- a/temporal/t.rast.aggregate/t.rast.aggregate.html +++ b/temporal/t.rast.aggregate/t.rast.aggregate.html @@ -193,6 +193,41 @@

Yearly aggregation

yearly_avg_temp_2004|climate|2004-01-01 00:00:00|2005-01-01 00:00:00 +

Weighted aggregation

+ +In this example, we create a STRDS of fictional temperature values for +a weekly interval: + +
+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"
+
+ +We can now use the -w flag and the sampling=related 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. + +
+t.rast.aggregate input=temperature_weekly output=temperature_10daily_weighted \
+  basename=tendaily_temp_weighted method=average granularity="10 days" \
+  sampling=related -w
+
+ +This is especially useful when harmonizing STRDS with granules that are not +integer multiplications of each other. +

SEE ALSO

diff --git a/temporal/t.rast.aggregate/t.rast.aggregate.py b/temporal/t.rast.aggregate/t.rast.aggregate.py index 0675407794a..af5561bc8bb 100755 --- a/temporal/t.rast.aggregate/t.rast.aggregate.py +++ b/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,6 +132,7 @@ def main(): gran = options["granularity"] base = options["basename"] register_null = flags["n"] + weighting = flags["w"] method = options["method"] sampling = options["sampling"] offset = options["offset"] @@ -134,6 +140,9 @@ def main(): 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, ) diff --git a/temporal/t.rast.aggregate/testsuite/test_aggregation_absolute.py b/temporal/t.rast.aggregate/testsuite/test_aggregation_absolute.py index 924e3950ca8..4791aba9e9a 100644 --- a/temporal/t.rast.aggregate/testsuite/test_aggregation_absolute.py +++ b/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