diff --git a/src/temporal/t.rast.aggregate.seasons/Makefile b/src/temporal/t.rast.aggregate.seasons/Makefile
new file mode 100644
index 0000000000..26a7530f29
--- /dev/null
+++ b/src/temporal/t.rast.aggregate.seasons/Makefile
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = t.rast.aggregate.seasons
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
diff --git a/src/temporal/t.rast.aggregate.seasons/t.rast.aggregate.seasons.html b/src/temporal/t.rast.aggregate.seasons/t.rast.aggregate.seasons.html
new file mode 100644
index 0000000000..5b2b315868
--- /dev/null
+++ b/src/temporal/t.rast.aggregate.seasons/t.rast.aggregate.seasons.html
@@ -0,0 +1,46 @@
+<h2>DESCRIPTION</h2>
+<em><b>t.rast.aggregate.seasons</b></em> aggregates an input space time
+raster dataset at astronomical seasons level using a statistical method selected by the user.
+ It uses <a href="https://grass.osgeo.org/grass-stable/manuals/t.rast.aggregate.ds.html">t.rast.aggregate.ds</a>
+to detect and copy the astronomical seasons granularity.
+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>
+
+Using the <em>output</em> option it will create a unified space time raster dataset with the name set in this option,
+otherwise it will create a space time raster dataset for each year in the input space time raster dataset.
+
+<h2>EXAMPLES</h2>
+Calculate seasonal data from an input space time raster dataset with a unified output space time raster dataset
+<div class="code"><pre>
+    t.rast.aggregate.seasons input=mystrds basename=mystrds_seasons output=outstrds
+</pre></div>
+
+Calculate seasonal data from an input space time raster dataset with yearly output space time raster datasets
+<div class="code"><pre>
+    t.rast.aggregate.seasons input=mystrds basename=mystrds_seasons
+</pre></div>
+
+<h2>NOTES</h2>
+If one of the season is not fully covered by the input space time
+raster dataset that season will not be created by the module. For
+example if you have a daily space time raster dataset for just one year
+it will create only raster maps for 3 seasons (Spring, Summer, Autumn),
+the Winter one will be avoided since your input space time raster
+dataset as data only until the end of the year but it miss the data
+from January to March.
+
+<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
diff --git a/src/temporal/t.rast.aggregate.seasons/t.rast.aggregate.seasons.py b/src/temporal/t.rast.aggregate.seasons/t.rast.aggregate.seasons.py
new file mode 100644
index 0000000000..bb9e65e45f
--- /dev/null
+++ b/src/temporal/t.rast.aggregate.seasons/t.rast.aggregate.seasons.py
@@ -0,0 +1,303 @@
+#!/usr/bin/env python
+
+
+################################################
+#
+# MODULE:       t.rast.aggregate.seasons
+# AUTHOR(S):    Luca Delucchi
+# PURPOSE:      Aggregates an input strds with astronomical seasons granularity.
+#
+# COPYRIGHT:    (C) 2023 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 G_OPT_STRDS_OUTPUT
+# % label: The name of a singular space time raster dataset
+# % description: Using this option all the yearly space time raster datasets will be merged in a singular space time raster dataset
+# % required: no
+# %end
+
+# %option
+# % key: years
+# % type: string
+# % label: List of years, separator may 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 datasets if output option is not used
+# % 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 processes to run in parallel
+# % required: no
+# % multiple: no
+# % answer: 1
+# %end
+
+# %flag
+# % key: n
+# % description: Register Null maps
+# %end
+
+# %flag
+# % key: m
+# % description: Use meteorological seasons
+# %end
+
+
+import copy
+import atexit
+from datetime import datetime
+
+import grass.temporal as tgis
+import grass.script as gs
+import grass.pygrass.modules as pymod
+from grass.pygrass.vector import VectorTopo
+from grass.pygrass.vector.geometry import Point
+
+remove_dataset = {"stvds": [], "strds": []}
+
+
+def cleanup():
+    """
+    Clean up temporary maps
+    """
+    # remove space time vector datasets
+    for typ, maps in remove_dataset.items():
+        for map in maps:
+            remod = pymod.Module("t.remove", run_=False)
+            remod.inputs.inputs = map
+            remod.inputs.type = typ
+            remod.flags.f = True
+            remod.flags.d = True
+            remod.flags.quiet = True
+            remod.run()
+
+
+def leap_year(year):
+    """Function to calculate if it is a leap year or not
+
+    Args:
+        year (int): the year to check
+
+    Returns:
+        bool: True if it is a leap year otherwise Falses
+    """
+    # divided by 100 means century year (ending with 00)
+    # century year divided by 400 is leap year
+    if (year % 400 == 0) and (year % 100 == 0):
+        return True
+
+    # not divided by 100 means not a century year
+    # year divided by 4 is a leap year
+    elif (year % 4 == 0) and (year % 100 != 0):
+        return True
+
+    # if not divided by both 400 (century year) and 4 (not century year)
+    # year is not leap year
+    else:
+        return False
+
+
+def main():
+    options, flags = gs.parser()
+    strds = options["input"]
+
+    output_name = options["output"]
+
+    meteorological = flags["m"]
+
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+    mapset = tgis.core.get_current_mapset()
+
+    if options["years"] != "":
+        try:
+            vals = options["years"].split("-")
+            years = range(vals)
+        except:
+            try:
+                years = options["years"].split(",")
+            except:
+                gs.fatal(_("Invalid option years"))
+    else:
+        if strds.find("@") >= 0:
+            id_ = strds
+        else:
+            id_ = f'{strds}@{gs.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 = ["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 = 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}")
+            year_int = int(year)
+            # check if it is meteorological or astronomic season for each season
+            if seas == "spring":
+                if meteorological:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 3, 1),
+                        end_time=datetime(year_int, 5, 31),
+                    )
+                else:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 3, 20),
+                        end_time=datetime(year_int, 6, 21),
+                    )
+            elif seas == "summer":
+                if meteorological:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 6, 1),
+                        end_time=datetime(year_int, 8, 30),
+                    )
+                else:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 6, 21),
+                        end_time=datetime(year_int, 9, 20),
+                    )
+            elif seas == "autumn":
+                if meteorological:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 9, 1),
+                        end_time=datetime(year_int, 11, 30),
+                    )
+                else:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 9, 20),
+                        end_time=datetime(year_int, 12, 21),
+                    )
+            elif seas == "winter":
+                if meteorological:
+                    if leap_year(year):
+                        extent = tgis.AbsoluteTemporalExtent(
+                            start_time=datetime(year_int, 12, 1),
+                            end_time=datetime(year_int + 1, 2, 29),
+                        )
+                    else:
+                        extent = tgis.AbsoluteTemporalExtent(
+                            start_time=datetime(year_int, 12, 1),
+                            end_time=datetime(year_int + 1, 2, 28),
+                        )
+                else:
+                    extent = tgis.AbsoluteTemporalExtent(
+                        start_time=datetime(year_int, 12, 21),
+                        end_time=datetime(year_int + 1, 3, 20),
+                    )
+            map_layer.set_temporal_extent(extent=extent)
+            season_vect.append(map_layer)
+        temp_season = f"sample_seasons_{year}"
+        outsp = tgis.open_new_stds(
+            temp_season,
+            "stvds",
+            "absolute",
+            f"Season vector year {year}",
+            f"Season vector for the year {year}",
+            "mean",
+            dbif,
+        )
+        tgis.register_map_object_list(
+            "vector",
+            season_vect,
+            outsp,
+            False,
+            None,
+            dbif,
+        )
+        remove_dataset["stvds"].append(temp_season)
+
+    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 = False
+    mod.flags.n = register_null
+
+    count = 0
+
+    outputs = []
+    # for each year calculate seasonal aggregation
+    for year in years:
+        print(year)
+        mymod = copy.deepcopy(mod)
+        mymod.inputs.sample = f"sample_seasons_{year}@{mapset}"
+        if output_name:
+            myout = f"{output_name}_{year}"
+            remove_dataset["strds"].append(myout)
+            outputs.append(myout)
+            mymod.outputs.output = myout
+        else:
+            mymod.outputs.output = f"{basename}_{year}"
+        print(mymod.get_bash())
+        process_queue.put(mymod)
+
+        if count % 10 == 0:
+            gs.percent(count, len(years), 1)
+
+    # Wait for unfinished processes
+    process_queue.wait()
+
+    if len(outputs) > 1:
+        pymod.Module("t.merge", inputs=",".join(outputs), output=output_name)
+
+
+if __name__ == "__main__":
+    atexit.register(cleanup)
+    main()
diff --git a/src/temporal/t.rast.aggregate.seasons/testsuite/test_seasons.py b/src/temporal/t.rast.aggregate.seasons/testsuite/test_seasons.py
new file mode 100644
index 0000000000..9dd52aaf1b
--- /dev/null
+++ b/src/temporal/t.rast.aggregate.seasons/testsuite/test_seasons.py
@@ -0,0 +1,130 @@
+"""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
+from grass.gunittest.case import TestCase
+from grass.gunittest.main import test
+
+
+class TestClimatologies(TestCase):
+    @classmethod
+    def setUpClass(cls):
+        """Initiate the temporal GIS and set the region"""
+        tgis.init(raise_fatal_error=True)
+        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="monthly")
+
+    def test_no_years(self):
+        """Test on all years"""
+        self.assertModule(
+            "t.rast.aggregate.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(), 7)
+
+    def test_one_year(self):
+        """Test just one year"""
+        self.assertModule(
+            "t.rast.aggregate.seasons",
+            input="monthly",
+            basename="oneseason",
+            years=2002,
+            overwrite=True,
+            verbose=True,
+        )
+        out = tgis.open_old_stds("oneseason_2002", type="strds")
+        self.assertEqual(out.metadata.get_number_of_maps(), 3)
+
+    def test_error_missing_basename(self):
+        """Test if basename is missing"""
+        self.assertModuleFail(
+            "t.rast.aggregate.seasons",
+            input="monthly",
+        )
+
+    def test_meteorological(self):
+        self.assertModule(
+            "t.rast.aggregate.seasons",
+            input="monthly",
+            output="meteo_seasons",
+            basename="meteoseasons",
+            overwrite=True,
+            verbose=True,
+            flags=["m"],
+        )
+        out = tgis.open_old_stds("meteo_seasons", type="strds")
+        self.assertEqual(out.metadata.get_number_of_maps(), 8)
+
+
+if __name__ == "__main__":
+    test()