t.rast.univar, t.rast3d.univar: Add support for zones

python/grass/temporal/univar_statistics.py

@@ -29,7 +29,15 @@
 
 
 def print_gridded_dataset_univar_statistics(
-    type, input, output, where, extended, no_header=False, fs="|", rast_region=False
+    type,
+    input,
+    output,
+    where,
+    extended,
+    no_header=False,
+    fs="|",
+    rast_region=False,
+    zones=None,
 ):
     """Print univariate statistics for a space time raster or raster3d dataset
 
@@ -43,8 +51,14 @@ def print_gridded_dataset_univar_statistics(
     :param rast_region: If set True ignore the current region settings
            and use the raster map regions for univar statistical calculation.
            Only available for strds.
+    :param zones: raster map with zones to calculate statistics for
     """
 
+    stats_module = {
+        "strds": "r.univar",
+        "str3ds": "r3.univar",
+    }[type]
+
     # We need a database interface
     dbif = SQLDatabaseInterfaceConnection()
     dbif.connect()
@@ -54,9 +68,14 @@ def print_gridded_dataset_univar_statistics(
     if output is not None:
         out_file = open(output, "w")
 
-    rows = sp.get_registered_maps("id,start_time,end_time", where, "start_time", dbif)
+    strds_cols = (
+        "id,start_time,end_time,semantic_label"
+        if type == "strds"
+        else "id,start_time,end_time"
+    )
+    rows = sp.get_registered_maps(strds_cols, where, "start_time", dbif)
 
-    if not rows:
+    if not rows and rows != [""]:
         dbif.close()
         err = "Space time %(sp)s dataset <%(i)s> is empty"
         if where:
@@ -66,40 +85,60 @@ def print_gridded_dataset_univar_statistics(
         )
 
     if no_header is False:
-        string = ""
-        string += "id" + fs + "start" + fs + "end" + fs + "mean" + fs
-        string += "min" + fs + "max" + fs
-        string += "mean_of_abs" + fs + "stddev" + fs + "variance" + fs
-        string += "coeff_var" + fs + "sum" + fs + "null_cells" + fs + "cells"
-        string += fs + "non_null_cells"
+        cols = (
+            ["id", "semantic_label", "start", "end"]
+            if type == "strds"
+            else ["id", "start", "end"]
+        )
+        if zones:
+            cols.append("zone")
+        cols.extend(
+            [
+                "mean",
+                "min",
+                "max",
+                "mean_of_abs",
+                "stddev",
+                "variance",
+                "coeff_var",
+                "sum",
+                "null_cells",
+                "cells",
+                "non_null_cells",
+            ]
+        )
         if extended is True:
-            string += fs + "first_quartile" + fs + "median" + fs
-            string += "third_quartile" + fs + "percentile_90"
+            cols.extend(["first_quartile", "median", "third_quartile", "percentile_90"])
+        string = fs.join(cols)
 
         if output is None:
             print(string)
         else:
             out_file.write(string + "\n")
 
+    flag = "g"
+
+    if extended is True:
+        flag += "e"
+    if type == "strds" and rast_region is True:
+        flag += "r"
+
     for row in rows:
         string = ""
         id = row["id"]
         start = row["start_time"]
         end = row["end_time"]
+        semantic_label = (
+            ""
+            if type != "strds" or not row["semantic_label"]
+            else row["semantic_label"]
+        )
 
-        flag = "g"
-
-        if extended is True:
-            flag += "e"
-        if type == "strds" and rast_region is True:
-            flag += "r"
-
-        if type == "strds":
-            stats = gscript.parse_command("r.univar", map=id, flags=flag)
-        elif type == "str3ds":
-            stats = gscript.parse_command("r3.univar", map=id, flags=flag)
+        univar_stats = gscript.read_command(
+            stats_module, map=id, flags=flag, zones=zones
+        ).rstrip()
 
-        if not stats:
+        if not univar_stats:
             if type == "strds":
                 gscript.warning(
                     _("Unable to get statistics for raster map " "<%s>") % id
@@ -109,19 +148,36 @@ def print_gridded_dataset_univar_statistics(
                     _("Unable to get statistics for 3d raster map" " <%s>") % id
                 )
             continue
+        eol = ""
 
-        string += str(id) + fs + str(start) + fs + str(end)
-        string += fs + str(stats["mean"]) + fs + str(stats["min"])
-        string += fs + str(stats["max"]) + fs + str(stats["mean_of_abs"])
-        string += fs + str(stats["stddev"]) + fs + str(stats["variance"])
-        string += fs + str(stats["coeff_var"]) + fs + str(stats["sum"])
-        string += fs + str(stats["null_cells"]) + fs + str(stats["cells"])
-        string += fs + str(int(stats["cells"]) - int(stats["null_cells"]))
-        if extended is True:
-            string += fs + str(stats["first_quartile"]) + fs + str(stats["median"])
+        for idx, stats_kv in enumerate(univar_stats.split(";")):
+            stats = gscript.utils.parse_key_val(stats_kv)
             string += (
-                fs + str(stats["third_quartile"]) + fs + str(stats["percentile_90"])
+                f"{id}{fs}{semantic_label}{fs}{start}{fs}{end}"
+                if type == "strds"
+                else f"{id}{fs}{start}{fs}{end}"
             )
+            if zones:
+                if idx == 0:
+                    zone = str(stats["zone"])
+                    string = ""
+                    continue
+                string += f"{fs}{zone}"
+                if "zone" in stats:
+                    zone = str(stats["zone"])
+                    eol = "\n"
+                else:
+                    eol = ""
+            string += f'{fs}{stats["mean"]}{fs}{stats["min"]}'
+            string += f'{fs}{stats["max"]}{fs}{stats["mean_of_abs"]}'
+            string += f'{fs}{stats["stddev"]}{fs}{stats["variance"]}'
+            string += f'{fs}{stats["coeff_var"]}{fs}{stats["sum"]}'
+            string += f'{fs}{stats["null_cells"]}{fs}{stats["n"]}'
+            string += f'{fs}{stats["n"]}'
+            if extended is True:
+                string += f'{fs}{stats["first_quartile"]}{fs}{stats["median"]}'
+                string += f'{fs}{stats["third_quartile"]}{fs}{stats["percentile_90"]}'
+            string += eol
 
         if output is None:
             print(string)

temporal/t.rast.univar/t.rast.univar.html

@@ -7,10 +7,15 @@ <h2>DESCRIPTION</h2>
 By default it returns the name of the map, the start and end date of
 dataset and the following values: mean, minimum and maximum vale,
 mean_of_abs, standard deviation, variance, coeff_var, number of null
-cells, total number of cell.
+cells, total number of cells.
 <p>
 Using the <em>e</em> flag it can calculate also extended statistics:
 first quartile, median value, third quartile and percentile 90.
+<p>
+If a <em>zones</em> raster map is provided, statistics are computed for
+each zone (category) in that input raster map. The <em>zones</em> option
+does not support Spatio-Temporal-Raster-Datasets (STRDS) but only a single,
+static raster map.
 
 <h2>EXAMPLE</h2>
 
@@ -32,8 +37,10 @@ <h2>SEE ALSO</h2>
 <em>
 <a href="t.create.html">t.create</a>,
 <a href="t.info.html">t.info</a>
+<a href="r.univar.html">t.create</a>,
 </em>
 
 <h2>AUTHOR</h2>
 
 S&ouml;ren Gebbert, Th&uuml;nen Institute of Climate-Smart Agriculture
+Stefan Blumentrath, (Support for zones)

temporal/t.rast.univar/t.rast.univar.py

@@ -31,6 +31,12 @@
 # %option G_OPT_STRDS_INPUT
 # %end
 
+# %option G_OPT_R_INPUT
+# % key: zones
+# % description: Raster map used for zoning, must be of type CELL
+# % required: no
+# %end
+
 # %option G_OPT_F_OUTPUT
 # % required: no
 # %end
@@ -60,8 +66,7 @@
 # % guisection: Formatting
 # %end
 
-import grass.script as grass
-
+import grass.script as gs
 
 ############################################################################
 
@@ -72,12 +77,13 @@ def main():
 
     # Get the options
     input = options["input"]
+    zones = options["zones"]
     output = options["output"]
     where = options["where"]
     extended = flags["e"]
     no_header = flags["u"]
     rast_region = bool(flags["r"])
-    separator = grass.separator(options["separator"])
+    separator = gs.separator(options["separator"])
 
     # Make sure the temporal database exists
     tgis.init()
@@ -87,11 +93,24 @@ def main():
     if output == "-":
         output = None
 
+    # Check if zones map exists and is of type CELL
+    if zones:
+        if gs.raster.raster_info(zones)["datatype"] != "CELL":
+            gs.fatal(_("Zoning raster must be of type CELL"))
+
     tgis.print_gridded_dataset_univar_statistics(
-        "strds", input, output, where, extended, no_header, separator, rast_region
+        "strds",
+        input,
+        output,
+        where,
+        extended,
+        no_header,
+        separator,
+        rast_region,
+        zones,
     )
 
 
 if __name__ == "__main__":
-    options, flags = grass.parser()
+    options, flags = gs.parser()
     main()