t.rast.univar: Add region_relation option for spatial filtering STDS …

…by computational region (OSGeo#2793) * add spatial extent filter * pass only spatial relation * add spatial extent filter * fix 3D where * add tests for spatial filter * add more tests for spatial filter * fix failing test * add spatial relation doc * lower message severity * spatial_relation intro * Update temporal/t.rast.univar/t.rast.univar.html Co-authored-by: Veronica Andreo <[email protected]> * Update temporal/temporalintro.html Co-authored-by: Veronica Andreo <[email protected]> * rewrite user message * fix tests for empty strds * rename spatial filter option * rename spatial filter option * update r-flag description * rename spatial filter option * rename spatial filter option * improve spatial filter docs * fix failing test * Update python/grass/temporal/abstract_space_time_dataset.py Co-authored-by: Vaclav Petras <[email protected]> * Update python/grass/temporal/univar_statistics.py Co-authored-by: Vaclav Petras <[email protected]> * address code review * address code review * address code review, some linting * add back percentiles * black --------- Co-authored-by: Veronica Andreo <[email protected]> Co-authored-by: Vaclav Petras <[email protected]>
neteler · Sep 29, 2023 · b20af2c · b20af2c
1 parent b60356c
commit b20af2c
Show file tree

Hide file tree

Showing 7 changed files with 282 additions and 91 deletions.
diff --git a/python/grass/temporal/abstract_space_time_dataset.py b/python/grass/temporal/abstract_space_time_dataset.py
@@ -1367,15 +1367,22 @@ def get_registered_maps_as_objects_with_gaps(
         In case more map information are needed, use the select()
         method for each listed object.
 
+        The combination of the spatial_extent and spatial_relation parameters
+        can be used to return only map objects with the given spatial relation
+        to the provided spatial extent.
+
         :param where: The SQL where statement to select a
                      subset of the registered maps without "WHERE"
         :param dbif: The database interface to be used
-        :param spatial_extent: Return only maps with the provided spatial
-                     relation to the given spatial extent (requires
-                     spatial_relation parameter)
-        :param spatial_relation: Return only maps with the given spatial
-                     relation to the provided spatial extent (requires
-                     spatial_extent parameter)
+        :param spatial_extent: Spatial extent dict and projection information
+            e.g. from g.region -ug3 with GRASS GIS region keys
+            "n", "s", "e", "w", "b", "t", and  "projection".
+        :param spatial_relation: Spatial relation to the provided
+            spatial extent as a string with one of the following values:
+            "overlaps": maps that spatially overlap ("intersect")
+                        within the provided spatial extent
+            "is_contained": maps that are fully within the provided spatial extent
+            "contains": maps that contain (fully cover) the provided spatial extent
 
         :return: ordered object list, in case nothing found None is returned
         """
@@ -1436,17 +1443,24 @@ def get_registered_maps_as_objects_with_temporal_topology(
         In case more map information are needed, use the select()
         method for each listed object.
 
+        The combination of the spatial_extent and spatial_relation parameters
+        can be used to return only maps with the given spatial relation to
+        the provided spatial extent
+
         :param where: The SQL where statement to select a subset of
                      the registered maps without "WHERE"
         :param order: The SQL order statement to be used to order the
                      objects in the list without "ORDER BY"
         :param dbif: The database interface to be used
-        :param spatial_extent: Return only maps with the provided spatial
-                     relation to the given spatial extent (requires
-                     spatial_relation parameter)
-        :param spatial_relation: Return only maps with the given spatial
-                     relation to the provided spatial extent (requires
-                     spatial_extent parameter)
+        :param spatial_extent: Spatial extent dict and projection information
+            e.g. from g.region -ug3 with GRASS GIS region keys
+            "n", "s", "e", "w", "b", "t", and  "projection".
+        :param spatial_relation: Spatial relation to the provided
+            spatial extent as a string with one of the following values:
+            "overlaps": maps that spatially overlap ("intersect")
+                        within the provided spatial extent
+            "is_contained": maps that are fully within the provided spatial extent
+            "contains": maps that contain (fully cover) the provided spatial extent
 
         :return: The ordered map object list,
                 In case nothing found None is returned
@@ -1483,17 +1497,24 @@ def get_registered_maps_as_objects(
         In case more map information are needed, use the select()
         method for each listed object.
 
+        The combination of the spatial_extent and spatial_relation parameters
+        can be used to return only maps with the given spatial relation to
+        the provided spatial extent
+
         :param where: The SQL where statement to select a subset of
                       the registered maps without "WHERE"
         :param order: The SQL order statement to be used to order the
                       objects in the list without "ORDER BY"
         :param dbif: The database interface to be used
-        :param spatial_extent: Return only maps with the provided spatial
-                     relation to the given spatial extent (requires
-                     spatial_relation parameter)
-        :param spatial_relation: Return only maps with the given spatial
-                     relation to the provided spatial extent (requires
-                     spatial_extent parameter)
+        :param spatial_extent: Spatial extent dict and projection information
+            e.g. from g.region -ug3 with GRASS GIS region keys
+            "n", "s", "e", "w", "b", "t", and  "projection".
+        :param spatial_relation: Spatial relation to the provided
+            spatial extent as a string with one of the following values:
+            "overlaps": maps that spatially overlap ("intersect")
+                        within the provided spatial extent
+            "is_contained": maps that are fully within the provided spatial extent
+            "contains": maps that contain (fully cover) the provided spatial extent
 
         :return: The ordered map object list,
                 In case nothing found None is returned
@@ -1657,7 +1678,7 @@ def _update_where_statement_by_spatial_extent(
         :param str where: SQL WHERE statement to be updated
         :param dict spatial_extent: Spatial extent dict and projection information
             e.g. from g.region -ug3
-        :param dict spatial_relation: Spatial relation to the provided
+        :param str spatial_relation: Spatial relation to the provided
             spatial extent as a string with one of the following values:
             "overlaps": maps that spatially overlap ("intersect")
                         within the provided spatial extent
@@ -1715,10 +1736,10 @@ def _update_where_statement_by_spatial_extent(
         if self.get_type() == "str3ds":
             if spatial_relation == "overlaps":
                 spatial_where_template += " AND top > {b}" " AND bottom < {t}"
-        elif spatial_relation == "is_contained":
-            spatial_where_template += " AND top <= {t}" " AND bottom >= {b}"
-        elif spatial_relation == "contains":
-            spatial_where_template += " AND top >= {t}" " AND bottom <= {b}"
+            elif spatial_relation == "is_contained":
+                spatial_where_template += " AND top <= {t}" " AND bottom >= {b}"
+            elif spatial_relation == "contains":
+                spatial_where_template += " AND top >= {t}" " AND bottom <= {b}"
         spatial_where_template += ")"
 
         spatial_where_list = [spatial_where_template.format(**spatial_extent)]
@@ -1757,18 +1778,25 @@ def get_registered_maps(
         In case columns are not specified, each row includes all columns
         specified in the datatype specific view.
 
+        The combination of the spatial_extent and spatial_relation parameters
+        can be used to return only SQL rows of maps with the given spatial
+        relation to the provided spatial extent
+
         :param columns: Columns to be selected as SQL compliant string
         :param where: The SQL where statement to select a subset
                      of the registered maps without "WHERE"
         :param order: The SQL order statement to be used to order the
                      objects in the list without "ORDER BY"
         :param dbif: The database interface to be used
-        :param spatial_extent: Return only maps with the provided spatial
-                     relation to the given spatial extent (requires
-                     spatial_relation parameter)
-        :param spatial_relation: Return only maps with the given spatial
-                     relation to the provided spatial extent (requires
-                     spatial_extent parameter)
+        :param spatial_extent: Spatial extent dict and projection information
+            e.g. from g.region -ug3 with GRASS GIS region keys
+            "n", "s", "e", "w", "b", "t", and  "projection".
+        :param spatial_relation: Spatial relation to the provided
+            spatial extent as a string with one of the following values:
+            "overlaps": maps that spatially overlap ("intersect")
+                        within the provided spatial extent
+            "is_contained": maps that are fully within the provided spatial extent
+            "contains": maps that contain (fully cover) the provided spatial extent
 
         :return: SQL rows of all registered maps,
                 In case nothing found None is returned

diff --git a/python/grass/temporal/univar_statistics.py b/python/grass/temporal/univar_statistics.py
@@ -21,11 +21,12 @@
 from multiprocessing import Pool
 from subprocess import PIPE
 
+import grass.script as gs
+from grass.pygrass.modules import Module
+
 from .core import SQLDatabaseInterfaceConnection, get_current_mapset
 from .factory import dataset_factory
 from .open_stds import open_old_stds
-import grass.script as gs
-from grass.pygrass.modules import Module
 
 ###############################################################################
 
@@ -115,11 +116,14 @@ def print_gridded_dataset_univar_statistics(
     no_header=False,
     fs="|",
     rast_region=False,
+    region_relation=None,
     zones=None,
     percentile=None,
     nprocs=1,
 ):
-    """Print univariate statistics for a space time raster or raster3d dataset
+    """Print univariate statistics for a space time raster or raster3d dataset.
+    Returns None if the space time raster dataset is empty or if applied
+    filters (where, region_relation) do not return any maps to process.
 
     :param type: Type of Space-Time-Dataset, must be either strds or str3ds
     :param input: The name of the space time dataset
@@ -133,6 +137,12 @@ 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 region_relation: Process only maps with the given spatial relation
+           to the computational region. A string with one of the following values:
+           "overlaps": maps that spatially overlap ("intersect")
+                       within the provided spatial extent
+           "is_contained": maps that are fully within the provided spatial extent
+           "contains": maps that contain (fully cover) the provided spatial extent
     :param zones: raster map with zones to calculate statistics for
     """
     # We need a database interface
@@ -144,22 +154,40 @@ def print_gridded_dataset_univar_statistics(
     if output is not None:
         out_file = open(output, "w")
 
+    spatial_extent = None
+    if region_relation:
+        spatial_extent = gs.parse_command("g.region", flags="3gu")
+
     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)
+    rows = sp.get_registered_maps(
+        strds_cols,
+        where,
+        "start_time",
+        dbif,
+        spatial_extent=spatial_extent,
+        spatial_relation=region_relation,
+    )
 
     if not rows and rows != [""]:
         dbif.close()
-        err = "Space time %(sp)s dataset <%(i)s> is empty"
-        if where:
-            err += " or where condition does not return any maps"
-        gs.fatal(
-            _(err) % {"sp": sp.get_new_map_instance(None).get_type(), "i": sp.get_id()}
+        gs.verbose(
+            _(
+                "No maps found to process. "
+                "Space time {type} dataset <{id}> is either empty "
+                "or the where condition (if used) does not return any maps "
+                "or no maps with the requested spatial relation to the "
+                "computational region exist in the dataset."
+            ).format(type=sp.get_new_map_instance(None).get_type(), id=sp.get_id())
         )
 
+        if output is not None:
+            out_file.close()
+        return
+
     if no_header is False:
         cols = (
             ["id", "semantic_label", "start", "end"]
@@ -316,7 +344,7 @@ def print_vector_dataset_univar_statistics(
             + fs
         )
         string += "min" + fs + "max" + fs + "range"
-        if type == "point" or type == "centroid":
+        if type in ("point", "centroid"):
             string += (
                 fs
                 + "mean"

diff --git a/temporal/t.rast.univar/t.rast.univar.html b/temporal/t.rast.univar/t.rast.univar.html
@@ -4,18 +4,32 @@ <h2>DESCRIPTION</h2>
 non-null cells for each registered raster map of a space time raster
 dataset.
 <p>
-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 cells.
+By default it returns the name of the map, the semantic label of the
+map, the start and end date of the map and the following values:
+mean, minimum and maximum vale, mean_of_abs, standard deviation, variance,
+coeff_var, number of null 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,
+does not support space time raster datasets (STRDS) but only a single,
 static raster map.
+<p>
+Space time raster datasets may contain raster maps with varying spatial
+extent like for example series of scenes of satellite images. With the
+<em>region_relation</em> option, computations can be limited to
+maps of the space time raster dataset that have a given spatial relation
+to the current computational region. Supported spatial relations are:
+<ul>
+ <li>"overlaps": process only maps that spatially overlap ("intersect")
+with the current computational region</li>
+ <li>"is_contained": process only maps that are fully within the current
+computational region</li>
+ <li>"contains": process only maps that contain (fully cover) the current
+computational region</li>
+</ul>
 
 <h2>EXAMPLE</h2>
 
@@ -43,4 +57,4 @@ <h2>SEE ALSO</h2>
 <h2>AUTHOR</h2>
 
 S&ouml;ren Gebbert, Th&uuml;nen Institute of Climate-Smart Agriculture
-Stefan Blumentrath, (Support for zones)
+Stefan Blumentrath, (Support for zones, parallel processing, and spatial relations)
diff --git a/temporal/t.rast.univar/t.rast.univar.py b/temporal/t.rast.univar/t.rast.univar.py
@@ -59,6 +59,15 @@
 # % guisection: Selection
 # %end
 
+# %option
+# % key: region_relation
+# % description: Process only maps with this spatial relation to the current computational region
+# % guisection: Selection
+# % options: overlaps,contains,is_contained
+# % required: no
+# % multiple: no
+# %end
+
 # %option G_OPT_F_SEP
 # % label: Field separator character between the output columns
 # % guisection: Formatting
@@ -71,7 +80,7 @@
 
 # %flag
 # % key: r
-# % description: Ignore the current region settings and use the raster map regions for univar statistical calculation
+# % description: Use the raster map regions for univar statistical calculation instead of the current region
 # %end
 
 # %flag
@@ -103,6 +112,7 @@ def main():
     output = options["output"]
     nprocs = int(options["nprocs"])
     where = options["where"]
+    region_relation = options["region_relation"]
     extended = flags["e"]
     no_header = flags["u"]
     rast_region = bool(flags["r"])
@@ -136,11 +146,12 @@ def main():
         output,
         where,
         extended,
+        percentile=percentile,
         no_header=no_header,
         fs=separator,
-        rast_region=rast_region,
         zones=zones,
-        percentile=percentile,
+        rast_region=rast_region,
+        region_relation=region_relation,
         nprocs=nprocs,
     )