Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

t.rast.aggregate: Adjustments to temporally weighted aggregation #1799

Open
wants to merge 13 commits into
base: main
Choose a base branch
from
Open
69 changes: 59 additions & 10 deletions python/grass/temporal/aggregation.py
Original file line number Diff line number Diff line change
Expand Up @@ -221,6 +221,7 @@ def aggregate_by_topology(
nprocs=1,
spatial=None,
dbif=None,
weighting=False,
overwrite=False,
file_limit=1000,
):
Expand Down Expand Up @@ -288,7 +289,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())
Expand Down Expand Up @@ -316,7 +316,43 @@ 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:
for map_layer in granule.overlaps:
set_list.add(map_layer)
if granule.overlapped:
for map_layer in granule.overlapped:
set_list.add(map_layer)
if granule.contains:
for map_layer in granule.contains:
set_list.add(map_layer)
if granule.equal:
for map_layer in granule.equal:
set_list.add(map_layer)
if granule.during:
for map_layer in granule.during:
set_list.add(map_layer)
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")
Expand Down Expand Up @@ -359,12 +395,17 @@ def aggregate_by_topology(
if len(aggregation_list) > 1:
# Create the r.series input file
filename = gscript.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 i, name in enumerate(aggregation_list):
echoix marked this conversation as resolved.
Show resolved Hide resolved
string = name
echoix marked this conversation as resolved.
Show resolved Hide resolved
weight = aggregation_weights[i]
echoix marked this conversation as resolved.
Show resolved Hide resolved
file.write("%s|%f\n" % (string, weight))
echoix marked this conversation as resolved.
Show resolved Hide resolved
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):
Expand All @@ -380,8 +421,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()
Expand Down
35 changes: 35 additions & 0 deletions temporal/t.rast.aggregate/t.rast.aggregate.html
Original file line number Diff line number Diff line change
Expand Up @@ -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
Comment on lines +202 to +205
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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
for i in `seq 0 11` ; do
r.mapcalc expression="map_${i} = ${i}"
echo map_${i} >> map_list.txt

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This example makes it much easier to understand what the weight really does as values in input maps are known and constant integer values

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>

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Let's see the result:
<div class="code"><pre>
t.rast.list input=temperature_10daily_weighted columns=name,start_time,end_time,min,max
name|start_time|end_time|min|max
tendaily_temp_weighted_2021_05_01|2021-05-01 00:00:00|2021-05-11 00:00:00|0.3|0.3
tendaily_temp_weighted_2021_05_11|2021-05-11 00:00:00|2021-05-21 00:00:00|1.6|1.6
tendaily_temp_weighted_2021_05_21|2021-05-21 00:00:00|2021-05-31 00:00:00|3.1|3.1
tendaily_temp_weighted_2021_05_31|2021-05-31 00:00:00|2021-06-10 00:00:00|4.5|4.5
tendaily_temp_weighted_2021_06_10|2021-06-10 00:00:00|2021-06-20 00:00:00|5.9|5.9
tendaily_temp_weighted_2021_06_20|2021-06-20 00:00:00|2021-06-30 00:00:00|7.4|7.4
tendaily_temp_weighted_2021_06_30|2021-06-30 00:00:00|2021-07-10 00:00:00|8.7|8.7
tendaily_temp_weighted_2021_07_10|2021-07-10 00:00:00|2021-07-20 00:00:00|10.3|10.3
tendaily_temp_weighted_2021_07_20|2021-07-20 00:00:00|2021-07-30 00:00:00|11.0|11.0
<pre><div>

This is especially useful when harmonizing STRDS with granules that are not
integer multiplications of each other.

<h2>SEE ALSO</h2>

<em>
Expand Down
12 changes: 11 additions & 1 deletion temporal/t.rast.aggregate/t.rast.aggregate.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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 gcore


Expand All @@ -128,13 +133,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 not sampling == "related":
gcore.fatal(_("Weighting only works with sampling: 'related'"))
echoix marked this conversation as resolved.
Show resolved Hide resolved

topo_list = sampling.split(",")

tgis.init()
Expand Down Expand Up @@ -203,6 +212,7 @@ def main():
method=method,
nprocs=nprocs,
spatial=None,
weighting=weighting,
overwrite=gcore.overwrite(),
file_limit=file_limit,
)
Expand Down
41 changes: 39 additions & 2 deletions temporal/t.rast.aggregate/testsuite/test_aggregation_absolute.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,11 +55,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"""
Expand Down Expand Up @@ -236,6 +236,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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There is splitlines for these cases (which actually behaves as expected unlike os.linesep here). (For future reference, see #2258, hopefully more straightforward.)

if (len(item) > 0 and item.startswith("b_"))
echoix marked this conversation as resolved.
Show resolved Hide resolved
]
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
Expand Down