TGIS: add spatial filter to STDS

[ninsbl]

This PR implements spatial filtering in the temporal library, in particular the get_registered_maps and get_registered_maps_as_objects functions.

The API should be unchanged, so that module code using those functions should not be affected.

If the library implementation is acceptable, it can be included in relevant module and library code.

Currently, only maps that are spatially not overlapping / intersecting with a given spatial_extent (e.g. from the computational region) are removed from the results of the functions named above.

If there is a need to support other spatial relations (like within or contains) the implementation needs to be modified a bit, and subsequent modificatins of module code may have to use an option (with the allowed spatial relation to the relevant spatial extent (e.g. computational region)) instead of a flag.

Any feedback is welcome. Module code will be worked on in separate PR(s).

See: https://lists.osgeo.org/pipermail/grass-dev/2022-December/095837.html
for discussion on ML.

[metzm]

One small question, no change request. For me ready to merge as it is.

[metzm]

This adds 3 spatial filters to the where statement, should work.

Maybe it would be more efficient to shift the spatial filter to the actual extent of the STDS to reduce the where statements?

[veroandreo]

Please go ahead with merging, I'll test once the functionality is included in t.* modules. Relationships like contain and within seem indeed desirable in the future 🤓

[ninsbl]

Thanks for your feedback, @metzm and @veroandreo.

@metzm : I am interested in improving the performance, esp. for larger STDS. Unfortunately, I don` think I can follow your idea of "shifting the spatial filter to the actual extent of the STDS to reduce the where statements"? Can you explain?
In locations that are not latlong, only one spatial filter is added to the where clause, consisting of 4 (STRDS) to 6 (STDS with 3D) conditions. In latlong locations two additional such filters will be added to the where clause to capture wrapping around the 0 and 180/-180. Talking of which, for the new relations (is_contained + contains) I have to think through the LL conditions once more, I think...

Now, I implemented handling of different spatial relations and added that to more relevant functions in the abstract_space_time_dataset.py. Spatial relations currently supported are overlaps (returning only maps that spatially overlap / intersect with a given spatial extent like the computational region), is_contained (returning only maps that are fully within a given spatial extent like the computational region), and contains (returning only maps that are fully cover / contain a given spatial extent like the computational region).

I am fully open to renaming the spatial relations, if you think other terms are more appropriate.

In relevant modules (e.g. t.rast.univar t.rast.list, ...), I suggest adding an optional option to process / return only maps with the given spatial relation to the computational region. If the option is not given/used, all maps are processed / returned. That way, the current behavior of the modules remains the default, to not break the current API.

[ninsbl]

OK, I tested now with some artificial data in WGS84 and data / extents that wrap around the -180/18+ meridian and it seems to work as expected.

import grass.script as gs
import grass.temporal as tgis
from grass.pygrass.modules.interface import Module

tgis.init()

spatial_extent = gs.parse_command("g.region", flags="3g", w=10, e=190, t=1, b=0)
Module("r.mapcalc", expression="map_a=1")
Module("t.create", output="strds1", title="strds1", description="strds1")
Module("t.register", input="strds1", maps="map_a", start="2022-12-14", end="2022-12-15")

strds = tgis.open_old_stds("strds1", type="strds")

spatial_extent = gs.parse_command("g.region", flags="3g", w=-160, e=20, t=1, b=0)

trds.get_registered_maps(where="", spatial_extent=spatial_extent, spatial_relation="overlaps")[0]["id"]
# Out[45]: 'map_a@PERMANENT'


spatial_extent = gs.parse_command("g.region", flags="3g", w=5, e=200, t=1, b=0)

strds.get_registered_maps(where="", spatial_extent=spatial_extent, spatial_relation="is_contained")[0]["id"]
# Out[50]: 'map_a@PERMANENT'

In [51]: strds.get_registered_maps(where="", spatial_extent=spatial_extent, spatial_relation="contains")[0]["id"]
# IndexError: list index out of range

strds.get_registered_maps(where="", spatial_extent=spatial_extent, spatial_relation="overlaps")[0]["id"]
# Out[52]: 'map_a@PERMANENT'

If there are no objections I will merge this rather soon and then add a new PR to make use of the modified function in a first module (likely t.rast.univar). Then it is also more easy to test, and implementation for regular users becomes clear...

I would like to add a doctest to the function, but doctests seem deactivated for this module... In unittests I am not sure how to include a WGS84 test... Hints are welcome...