t.rast.climatologies: new addon for GRASS 8

src/temporal/t.rast.climatologies/Makefile

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

src/temporal/t.rast.climatologies/t.rast.climatologies.html

@@ -0,0 +1,40 @@
+<h2>DESCRIPTION</h2>
+
+<b>t.rast.climatologies</b> is a module to calculate climatologies, i.e., long term stats, for single days
+or months in a time series. If the <em>s</em> flag is used, the module outputs a new space time
+raster dataset of relative temporal type with the aggregated maps.
+
+<h2>EXAMPLE</h2>
+
+<h3>Daily climatologies</h3>
+
+Starting from a space time raster dataset of daily granularity (or granularity lower than one day),
+the module will create a new space time raster dataset of relative temporal type containing the
+long term average for each day along years.
+
+<div class="code">
+    <pre>
+        t.rast.climatologies input=myinput output=dailyoutput
+    </pre>
+</div>
+
+<h3>Monthly climatologies</h3>
+
+Starting from a space time raster dataset of monthly granularity (or lower than one month),
+the module will create two new space time raster datasets containing the long term minimum
+and maximum for each month along years.
+<div class="code">
+    <pre>
+        t.rast.climatologies input=myinput granularity=month method=min,max output=monthlyoutputmin,monthlyoutputmax
+    </pre>
+</div>
+
+<h2>SEE ALSO</h2>
+
+<a href="https://grass.osgeo.org/grass-stable/manuals/r.series.html">r.series</a>,
+<a href="https://grass.osgeo.org/grass-stable/manuals/t.rast.series.html">t.rast.series</a>
+
+
+<h2>AUTHOR</h2>
+
+Luca Delucchi, <i>Fondazione Edmund Mach</i>

src/temporal/t.rast.climatologies/t.rast.climatologies.py

@@ -0,0 +1,244 @@
+#!/usr/bin/env python
+
+################################################
+#
+# MODULE:       t.rast.climatologies
+# AUTHOR(S):    Luca Delucchi, Fondazione Edmund Mach
+# PURPOSE:      t.rast.climatologies calculates climatologies for space time raster datasets
+#
+# 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 climatologies from a space time raster dataset of absolute temporal type
+# % keyword: temporal
+# % keyword: raster
+# % keyword: aggregation
+# % keyword: series
+# % keyword: time
+# %end
+
+# %option G_OPT_STRDS_INPUT
+# %end
+
+# %option G_OPT_STRDS_OUTPUT
+# % required: no
+# %end
+
+# %option
+# % key: basename
+# % type: string
+# % label: Basename of the newly generated output maps
+# % description: Either a numerical suffix or the start time (s-flag) 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: yes
+# % 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: quantile
+# % type: double
+# % description: Quantile to calculate for method=quantile
+# % required: no
+# % multiple: yes
+# % options: 0.0-1.0
+# %end
+
+# %option
+# % key: granularity
+# % type: string
+# % label: Aggregate by day or month
+# % required: yes
+# % multiple: no
+# % options: day, month
+# % answer: day
+# %end
+
+# %option
+# % key: nprocs
+# % type: integer
+# % description: Number of processes to run in parallel
+# % required: no
+# % multiple: no
+# % answer: 1
+# %end
+
+# %flag
+# % key: s
+# % description: Do not create a space time raster dataset as output
+# %end
+
+import copy
+from datetime import datetime
+
+
+def main():
+    import grass.pygrass.modules as pymod
+    import grass.script as gscript
+    import grass.temporal as tgis
+
+    options, flags = gscript.parser()
+    strds = options["input"]
+    output = options["output"]
+    method = options["method"]
+    gran = options["granularity"]
+    basename = options["basename"]
+    nprocs = options["nprocs"]
+    quantile = options["quantile"]
+    space_time = flags["s"]
+
+    # check if quantile value is used correctly
+    if "quantile" in method and not quantile:
+        gscript.fatal(_("Number requested methods and output maps do not match."))
+    elif quantile and "quantile" not in method:
+        gscript.warning(
+            _("Quantile option set but quantile not selected in method option")
+        )
+
+    # Check if number of methods and output maps matches
+    if "quantile" in method:
+        len_method = len(method.split(",")) - 1
+    else:
+        len_method = len(method.split(","))
+
+    if not space_time:
+        if (len(list(filter(None, quantile.split(",")))) + len_method) != len(
+            output.split(",")
+        ):
+            gscript.fatal(_("Number requested methods and output maps do not match."))
+
+    tgis.init()
+    # We need a database interface
+    dbif = tgis.SQLDatabaseInterfaceConnection()
+    dbif.connect()
+
+    insp = tgis.open_old_stds(strds, "strds", dbif)
+    temporal_type, semantic_type, title, description = insp.get_initial_values()
+    if temporal_type != "absolute":
+        gscript.fatal(
+            _(
+                "Space time raster dataset temporal type is not absolute, this module requires absolute time"
+            )
+        )
+    maps = insp.get_registered_maps_as_objects(None, "start_time", dbif)
+    if maps is None:
+        gscript.fatal(
+            _(
+                "No maps selected in the space time raster dataset {};"
+                " it might be empty or the where option returns no data".format(strds)
+            )
+        )
+    # start the r.series module to be used in a ParallelModuleQueue
+    mod = pymod.Module("r.series")
+    mod.inputs.method = method
+    mod.flags.quiet = True
+    if quantile:
+        mod.inputs.quantile = quantile
+    process_queue = pymod.ParallelModuleQueue(int(nprocs))
+    mapset = tgis.core.get_current_mapset()
+    # depending on granularity it calculates daily or monthly climatologies
+    outmaps = []
+    if gran == "day":
+        outunit = "days"
+        # for each day
+        for doy in range(1, 367):
+            doystr = datetime.strptime(f"2000 {doy}", "%Y %j").strftime("%m_%d")
+            thiswhere = f"strftime('%m_%d', start_time) == '{doystr}'"
+            selemaps = insp.get_registered_maps_as_objects(
+                thiswhere, "start_time", dbif
+            )
+            maps_name = [sam.get_id() for sam in selemaps]
+            # check if there are maps for that day
+            if len(maps_name) > 0:
+                outname = f"{basename}_{doystr}"
+                runmod = copy.deepcopy(mod)
+                runmod.inputs.input = ",".join(maps_name)
+                runmod.outputs.output = outname
+                process_queue.put(runmod)
+                map_layer = tgis.space_time_datasets.RasterDataset(
+                    f"{outname}@{mapset}"
+                )
+                extent = tgis.RelativeTemporalExtent(
+                    start_time=doy - 1,
+                    end_time=doy,
+                    unit=outunit,
+                )
+                map_layer.set_temporal_extent(extent=extent)
+                outmaps.append(map_layer)
+
+            if doy % 10 == 0:
+                gscript.percent(doy, 366, 1)
+    else:
+        outunit = "months"
+        for month in range(1, 13):
+            monthstr = "{:02d}".format(month)
+            thiswhere = f"strftime('%m', start_time) == '{monthstr}'"
+            selemaps = insp.get_registered_maps_as_objects(
+                thiswhere, "start_time", None
+            )
+            maps_name = [sam.get_id() for sam in selemaps]
+            if len(maps_name) > 0:
+                outname = f"{basename}_{monthstr}"
+                runmod = copy.deepcopy(mod)
+                runmod.inputs.input = ",".join(maps_name)
+                runmod.outputs.output = outname
+                process_queue.put(runmod)
+                map_layer = tgis.space_time_datasets.RasterDataset(
+                    f"{outname}@{mapset}"
+                )
+                extent = tgis.RelativeTemporalExtent(
+                    start_time=month - 1,
+                    end_time=month,
+                    unit=outunit,
+                )
+                map_layer.set_temporal_extent(extent=extent)
+                outmaps.append(map_layer)
+            gscript.percent(month, 12, 1)
+
+    if not space_time:
+        # create new space time raster dataset
+        output_strds = tgis.open_new_stds(
+            output,
+            "strds",
+            "relative",
+            f"{gran} {method} climatologies",
+            f"Climatologies created from {strds}, {gran} {method} maps",
+            semantic_type,
+            dbif,
+            gscript.overwrite(),
+        )
+        register_null = False
+        # register maps into space time raster dataset
+        tgis.register_map_object_list(
+            "rast",
+            outmaps,
+            output_strds,
+            register_null,
+            outunit,
+            dbif,
+        )
+
+        # Update the raster metadata table entries
+        output_strds.metadata.update(dbif)
+
+    dbif.close()
+    return True
+
+
+if __name__ == "__main__":
+    main()