t.rast.seasons: new temporal module to get seasonally aggregated data

src/temporal/t.rast.seasons/Makefile

@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = t.rast.aggregate.seasons
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

src/temporal/t.rast.seasons/t.rast.seasons.html

@@ -0,0 +1,29 @@
+<h2>DESCRIPTION</h2>
+<em><b>t.rast.aggregate.seasons</b></em> calculate a new space time
+raster dataset with statics at astronomical seasonal level from an input space time
+raster dataset. It use <a href="https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.ds.html">t.rast.aggregate.ds</a>
+to calculate aggragated data.
+Astronomical seasons are defined as:
+<ul>
+    <li><em>Spring</em> 20 March - 21 June</li>
+    <li><em>Summer</em> 21 June - 20 September</li>
+    <li><em>Autumn</em> 20 September - 21 December</li>
+    <li><em>Winter</em> 21 December - 20 March (following year)</li>
+</ul>
+<p>
+
+<h2>EXAMPLES</h2>
+Calculate seasonal data from
+<div class="code"><pre>
+    t.rast.aggregate.seasons input=mystrds basename=mystrds_seasons
+</pre></div>
+
+<h2>SEE ALSO</h2>
+<em>
+    <a href="https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.ds.html">t.rast.aggregate.ds</a>
+    <a href="https://grass.osgeo.org/grass-stable/manuals/r.series.html">r.null</a>
+</em>
+
+<h2>AUTHOR</h2>
+
+Luca Delucchi, Fondazione Edmund Mach

src/temporal/t.rast.seasons/t.rast.seasons.py

@@ -0,0 +1,203 @@
+#!/usr/bin/env python
+
+
+################################################
+#
+# MODULE:       t.rast.aggregate.seasons
+# AUTHOR(S):    Luca Delucchi
+# PURPOSE:      Create new .
+#
+# COPYRIGHT:    (C) 2018 by Luca Delucchi
+#
+#               This program is free software under the GNU General Public
+#               License (>=v2). Read the file COPYING that comes with GRASS
+#               for details.
+#
+################################################
+
+# %module
+# % description: Calculates seasonal data according to astronomical seasons.
+# % keyword: temporal
+# % keyword: raster
+# % keyword: aggregation
+# % keyword: series
+# %end
+
+# %option G_OPT_STRDS_INPUT
+# %end
+
+# %option
+# % key: years
+# % type: string
+# % label: List of years, separator could be comma (list) or minus (range)
+# % multiple: yes
+# % required: no
+# %end
+
+# %option
+# % key: basename
+# % type: string
+# % label: Basename of the new generated output maps and space time raster dataset
+# % description: A numerical suffix separated by an underscore will be attached to create a unique identifier
+# % required: yes
+# % multiple: no
+# % gisprompt:
+# %end
+
+# %option
+# % key: method
+# % type: string
+# % description: Aggregate operation to be performed on the raster maps
+# % required: yes
+# % multiple: no
+# % options: average,count,median,mode,minimum,min_raster,maximum,max_raster,stddev,range,sum,variance,diversity,slope,offset,detcoeff,quart1,quart3,perc90,quantile,skewness,kurtosis
+# % answer: average
+# %end
+
+# %option
+# % key: nprocs
+# % type: integer
+# % description: Number of r.null processes to run in parallel
+# % required: no
+# % multiple: no
+# % answer: 1
+# %end
+
+# %flag
+# % key: n
+# % description: Register Null maps
+# %end
+
+import copy
+import sys
+from datetime import datetime
+
+import grass.temporal as tgis
+import grass.script as gscript
+import grass.pygrass.modules as pymod
+from grass.pygrass.vector import VectorTopo
+from grass.pygrass.vector.geometry import Point
+
+
+def main():
+    strds = options["input"]
+
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+    mapset = tgis.core.get_current_mapset()
+
+    try:
+        vals = options["years"].split("-")
+        years = range(vals)
+    except:
+        try:
+            years = options["years"].split(",")
+        except:
+            if strds.find("@") >= 0:
+                id_ = strds
+            else:
+                id_ = strds + "@" + gscript.gisenv()["MAPSET"]
+            dataset = tgis.dataset_factory("strds", id_)
+            dataset.select(dbif)
+            ext = dataset.get_temporal_extent()
+            years = range(ext.start_time.year, ext.end_time.year)
+
+    method = options["method"]
+    basename = options["basename"]
+    nprocs = int(options["nprocs"])
+    register_null = flags["n"]
+
+    seasons_name = ["spring", "summer", "autumn", "winter"]
+
+    # create new space time vector datasets one for each year to be used as sampler
+    for year in years:
+        season_vect = []
+        for seas in seasons_name:
+            name = f"sample_{seas}_{year}"
+            vect = VectorTopo(name)
+            vect.open("w")
+            point = Point(0, 0)
+            vect.write(point, cat=1)
+            vect.close()
+            map_layer = tgis.space_time_datasets.VectorDataset(f"{name}@{mapset}")
+            if seas == "spring":
+                extent = tgis.AbsoluteTemporalExtent(
+                    start_time=datetime(int(year), 3, 20),
+                    end_time=datetime(int(year), 6, 21),
+                )
+            elif seas == "summer":
+                extent = tgis.AbsoluteTemporalExtent(
+                    start_time=datetime(int(year), 6, 21),
+                    end_time=datetime(int(year), 9, 20),
+                )
+            elif seas == "autumn":
+                extent = tgis.AbsoluteTemporalExtent(
+                    start_time=datetime(int(year), 9, 20),
+                    end_time=datetime(int(year), 12, 21),
+                )
+            elif seas == "winter":
+                extent = tgis.AbsoluteTemporalExtent(
+                    start_time=datetime(int(year), 12, 21),
+                    end_time=datetime(int(year) + 1, 3, 20),
+                )
+            map_layer.set_temporal_extent(extent=extent)
+            season_vect.append(map_layer)
+        outsp = tgis.open_new_stds(
+            f"sample_seasons_{year}",
+            "stvds",
+            "absolute",
+            f"Season vector year {year}",
+            f"Season vector for the year {year}",
+            "mean",
+            dbif,
+            gscript.overwrite(),
+        )
+        tgis.register_map_object_list(
+            "vector",
+            season_vect,
+            outsp,
+            False,
+            None,
+            dbif,
+        )
+
+    process_queue = pymod.ParallelModuleQueue(int(nprocs))
+
+    # create t.rast.aggregate.ds module to be copied
+    mod = pymod.Module("t.rast.aggregate.ds")
+    mod.inputs.input = strds
+    mod.inputs.method = method
+    mod.inputs.basename = basename
+    mod.inputs.type = "stvds"
+    mod.flags.quiet = True
+    mod.flags.n = register_null
+
+    count = 0
+
+    # for each year calculate seasonal aggregation
+    for year in years:
+        mymod = copy.deepcopy(mod)
+        mymod.inputs.sample = f"sample_seasons_{year}@{mapset}"
+        mymod.outputs.output = f"{basename}_{year}"
+        process_queue.put(mymod)
+
+        if count % 10 == 0:
+            gscript.percent(count, len(years), 1)
+
+    # Wait for unfinished processes
+    process_queue.wait()
+
+    # remove space time vector datasets
+    for year in years:
+        remod = pymod.Module("t.register")
+        remod.inputs.inputs = f"sample_seasons_{year}@{mapset}"
+        remod.flags.d = True
+        remod.flags.f = True
+        remod.flags.quiet = True
+
+
+if __name__ == "__main__":
+    options, flags = gscript.parser()
+    sys.exit(main())

src/temporal/t.rast.seasons/testsuite/test_seasons.py

@@ -0,0 +1,114 @@
+"""Test for t.rast.seasons
+(C) 2023 by the GRASS Development Team
+This program is free software under the GNU General Public
+License (>=v2). Read the file COPYING that comes with GRASS
+for details.
+@author Luca Delucchi
+"""
+
+import grass.temporal as tgis
+import datetime
+from grass.gunittest.case import TestCase
+from grass.gunittest.gmodules import SimpleModule
+
+
+class TestClimatologies(TestCase):
+    @classmethod
+    def setUpClass(cls):
+        """Initiate the temporal GIS and set the region"""
+        tgis.init(True)  # Raise on error instead of exit(1)
+        cls.use_temp_region()
+        cls.runModule("g.region", s=0, n=80, w=0, e=120, b=0, t=50, res=10, res3=10)
+
+        cls.runModule("r.mapcalc", expression="a_1 = 5", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_2 = 10", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_3 = 15", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_4 = 20", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_5 = 25", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_6 = 30", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_7 = 35", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_8 = 40", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_9 = 45", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_10 = 50", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_11 = 55", overwrite=True)
+        cls.runModule("r.mapcalc", expression="a_12 = 60", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_1 = 5", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_2 = 10", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_3 = 15", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_4 = 20", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_5 = 25", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_6 = 30", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_7 = 35", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_8 = 40", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_9 = 45", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_10 = 50", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_11 = 55", overwrite=True)
+        cls.runModule("r.mapcalc", expression="b_12 = 60", overwrite=True)
+
+        cls.runModule(
+            "t.create",
+            type="strds",
+            temporaltype="absolute",
+            output="monthly",
+            title="Monthly test",
+            description="Monthly test",
+            overwrite=True,
+        )
+        cls.runModule(
+            "t.register",
+            flags="i",
+            type="raster",
+            input="monthly",
+            maps="a_1,a_2,a_3,a_4,a_5,a_6,a_7,a_8,a_9,a_10,a_11,a_12",
+            start="2001-01-01",
+            increment="1 month",
+            overwrite=True,
+        )
+        cls.runModule(
+            "t.register",
+            flags="i",
+            type="raster",
+            input="monthly",
+            maps="b_1,b_2,b_3,b_4,b_5,b_6,b_7,b_8,b_9,b_10,b_11,b_12",
+            start="2002-01-01",
+            increment="1 month",
+            overwrite=True,
+        )
+
+    @classmethod
+    def tearDownClass(cls):
+        """Remove the time series"""
+        cls.runModule("t.remove", flags="rf", type="strds", inputs="daily")
+
+    def test_no_years(self):
+        """Test on all years"""
+        self.assertModule(
+            "t.rast.seasons",
+            input="monthly",
+            output="monthly_seasons",
+            basename="seasons",
+            overwrite=True,
+            verbose=True,
+        )
+        out = tgis.open_old_stds("monthly_seasons", type="strds")
+        self.assertEqual(out.metadata.get_number_of_maps(), 8)
+
+    def test_one_year(self):
+        """Test just one year"""
+        self.assertModule(
+            "t.rast.seasons",
+            input="monthly",
+            basename="season_2002",
+            years=2002,
+            overwrite=True,
+            verbose=True,
+        )
+        out = tgis.open_old_stds("monthly_seasons", type="strds")
+        self.assertEqual(out.metadata.get_number_of_maps(), 4)
+
+    def test_error_missing_basename(self):
+        """Test if basename is missing"""
+        self.assertModuleFail(
+            "t.rast.seasons",
+            input="monthly",
+        )