"""
Aggregation methods for space time raster datasets
Usage:
.. code-block:: python
import grass.temporal as tgis
tgis.aggregate_raster_maps(dataset, mapset, inputs, base, start, end, count, method, register_null, dbif)
(C) 2012-2013 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: Soeren Gebbert
"""
import grass.script as gscript
from grass.exceptions import CalledModuleError
from .space_time_datasets import RasterDataset
from .datetime_math import create_suffix_from_datetime
from .datetime_math import create_time_suffix
from .datetime_math import create_numeric_suffix
from .core import get_current_mapset, get_tgis_message_interface, init_dbif
from .spatio_temporal_relationships import (
SpatioTemporalTopologyBuilder,
create_temporal_relation_sql_where_statement,
)
###############################################################################
[docs]def collect_map_names(sp, dbif, start, end, sampling):
"""Gather all maps from dataset using a specific sample method
:param sp: The space time raster dataset to select aps from
:param dbif: The temporal database interface to use
:param start: The start time of the sample interval, may be relative or
absolute
:param end: The end time of the sample interval, may be relative or
absolute
:param sampling: The sampling methods to use
"""
use_start = False
use_during = False
use_overlap = False
use_contain = False
use_equal = False
use_follows = False
use_precedes = False
# Initialize the methods
if sampling:
for name in sampling.split(","):
if name == "start":
use_start = True
if name == "during":
use_during = True
if name == "overlap":
use_overlap = True
if name == "contain":
use_contain = True
if name == "equal":
use_equal = True
if name == "follows":
use_follows = True
if name == "precedes":
use_precedes = True
else:
use_start = True
if sp.get_map_time() != "interval":
use_start = True
use_during = False
use_overlap = False
use_contain = False
use_equal = False
use_follows = False
use_precedes = False
where = create_temporal_relation_sql_where_statement(
start,
end,
use_start,
use_during,
use_overlap,
use_contain,
use_equal,
use_follows,
use_precedes,
)
rows = sp.get_registered_maps("id", where, "start_time", dbif)
if not rows:
return None
names = []
for row in rows:
names.append(row["id"])
return names
###############################################################################
[docs]def aggregate_raster_maps(
inputs, base, start, end, count, method, register_null, dbif, offset=0
):
"""Aggregate a list of raster input maps with r.series
:param inputs: The names of the raster maps to be aggregated
:param base: The basename of the new created raster maps
:param start: The start time of the sample interval, may be relative or
absolute
:param end: The end time of the sample interval, may be relative or
absolute
:param count: The number to be attached to the basename of the new
created raster map
:param method: The aggreation method to be used by r.series
:param register_null: If true null maps will be registered in the space
time raster dataset, if false not
:param dbif: The temporal database interface to use
:param offset: Offset to be added to the map counter to create the map ids
"""
msgr = get_tgis_message_interface()
msgr.verbose(_("Aggregating %s raster maps") % (len(inputs)))
output = "%s_%i" % (base, int(offset) + count)
mapset = get_current_mapset()
map_id = output + "@" + mapset
new_map = RasterDataset(map_id)
# Check if new map is in the temporal database
if new_map.is_in_db(dbif):
if gscript.overwrite() is True:
# Remove the existing temporal database entry
new_map.delete(dbif)
new_map = RasterDataset(map_id)
else:
msgr.error(
_(
"Raster map <%(name)s> is already in temporal "
"database, use overwrite flag to overwrite"
% ({"name": new_map.get_name()})
)
)
return
msgr.verbose(
_(
"Computing aggregation of maps between %(st)s - %(end)s"
% {"st": str(start), "end": str(end)}
)
)
# Create the r.series input file
filename = gscript.tempfile(True)
file = open(filename, "w")
for name in inputs:
string = "%s\n" % (name)
file.write(string)
file.close()
# Run r.series
try:
if len(inputs) > 1000:
gscript.run_command(
"r.series",
flags="z",
file=filename,
output=output,
overwrite=gscript.overwrite(),
method=method,
)
else:
gscript.run_command(
"r.series",
file=filename,
output=output,
overwrite=gscript.overwrite(),
method=method,
)
except CalledModuleError:
dbif.close()
msgr.fatal(_("Error occurred in r.series computation"))
# Read the raster map data
new_map.load()
# In case of a null map continue, do not register null maps
if new_map.metadata.get_min() is None and new_map.metadata.get_max() is None:
if not register_null:
gscript.run_command("g.remove", flags="f", type="raster", name=output)
return None
return new_map
##############################################################################
[docs]def aggregate_by_topology(
granularity_list,
granularity,
map_list,
topo_list,
basename,
time_suffix,
offset=0,
method="average",
nprocs=1,
spatial=None,
dbif=None,
overwrite=False,
file_limit=1000,
):
"""Aggregate a list of raster input maps with r.series
:param granularity_list: A list of AbstractMapDataset objects.
The temporal extents of the objects are used
to build the spatio-temporal topology with the
map list objects
:param granularity: The granularity of the granularity list
:param map_list: A list of RasterDataset objects that contain the raster
maps that should be aggregated
:param topo_list: A list of strings of topological relations that are
used to select the raster maps for aggregation
:param basename: The basename of the new generated raster maps
:param time_suffix: Use the granularity truncated start time of the
actual granule to create the suffix for the basename
:param offset: Use a numerical offset for suffix generation
(overwritten by time_suffix)
:param method: The aggregation method of r.series (average,min,max, ...)
:param nprocs: The number of processes used for parallel computation
:param spatial: This indicates if the spatial topology is created as
well: spatial can be None (no spatial topology), "2D"
using west, east, south, north or "3D" using west,
east, south, north, bottom, top
:param dbif: The database interface to be used
:param overwrite: Overwrite existing raster maps
:param file_limit: The maximum number of raster map layers that
should be opened at once by r.series
:return: A list of RasterDataset objects that contain the new map names
and the temporal extent for map registration
"""
import grass.pygrass.modules as pymod
import copy
msgr = get_tgis_message_interface()
dbif, connection_state_changed = init_dbif(dbif)
topo_builder = SpatioTemporalTopologyBuilder()
topo_builder.build(mapsA=granularity_list, mapsB=map_list, spatial=spatial)
# The module queue for parallel execution
process_queue = pymod.ParallelModuleQueue(int(nprocs))
# Dummy process object that will be deep copied
# and be put into the process queue
r_series = pymod.Module(
"r.series",
output="spam",
method=[method],
overwrite=overwrite,
quiet=True,
run_=False,
finish_=False,
)
g_copy = pymod.Module(
"g.copy", raster=["spam", "spamspam"], quiet=True, run_=False, finish_=False
)
output_list = []
count = 0
for granule in granularity_list:
msgr.percent(count, len(granularity_list), 1)
count += 1
aggregation_list = []
if "equal" in topo_list and granule.equal:
for map_layer in granule.equal:
aggregation_list.append(map_layer.get_name())
if "contains" in topo_list and granule.contains:
for map_layer in granule.contains:
aggregation_list.append(map_layer.get_name())
if "during" in topo_list and granule.during:
for map_layer in granule.during:
aggregation_list.append(map_layer.get_name())
if "starts" in topo_list and granule.starts:
for map_layer in granule.starts:
aggregation_list.append(map_layer.get_name())
if "started" in topo_list and granule.started:
for map_layer in granule.started:
aggregation_list.append(map_layer.get_name())
if "finishes" in topo_list and granule.finishes:
for map_layer in granule.finishes:
aggregation_list.append(map_layer.get_name())
if "finished" in topo_list and granule.finished:
for map_layer in granule.finished:
aggregation_list.append(map_layer.get_name())
if "overlaps" in topo_list and granule.overlaps:
for map_layer in granule.overlaps:
aggregation_list.append(map_layer.get_name())
if "overlapped" in topo_list and granule.overlapped:
for map_layer in granule.overlapped:
aggregation_list.append(map_layer.get_name())
if aggregation_list:
msgr.verbose(
_("Aggregating %(len)i raster maps from %(start)s to" " %(end)s")
% (
{
"len": len(aggregation_list),
"start": str(granule.temporal_extent.get_start_time()),
"end": str(granule.temporal_extent.get_end_time()),
}
)
)
if granule.is_time_absolute() is True and time_suffix == "gran":
suffix = create_suffix_from_datetime(
granule.temporal_extent.get_start_time(), granularity
)
output_name = "{ba}_{su}".format(ba=basename, su=suffix)
elif granule.is_time_absolute() is True and time_suffix == "time":
suffix = create_time_suffix(granule)
output_name = "{ba}_{su}".format(ba=basename, su=suffix)
else:
output_name = create_numeric_suffix(
basename, count + int(offset), time_suffix
)
map_layer = RasterDataset("%s@%s" % (output_name, get_current_mapset()))
map_layer.set_temporal_extent(granule.get_temporal_extent())
if map_layer.map_exists() is True and overwrite is False:
msgr.fatal(
_(
"Unable to perform aggregation. Output raster "
"map <%(name)s> exists and overwrite flag was "
"not set" % ({"name": output_name})
)
)
output_list.append(map_layer)
if len(aggregation_list) > 1:
# Create the r.series input file
filename = gscript.tempfile(True)
file = open(filename, "w")
for name in aggregation_list:
string = "%s\n" % (name)
file.write(string)
file.close()
mod = copy.deepcopy(r_series)
mod(file=filename, output=output_name)
if len(aggregation_list) > int(file_limit):
msgr.warning(
_(
"The limit of open files (%i) was "
"reached (%i). The module r.series will "
"be run with flag z, to avoid open "
"files limit exceeding."
% (int(file_limit), len(aggregation_list))
)
)
mod(flags="z")
process_queue.put(mod)
else:
mod = copy.deepcopy(g_copy)
mod(raster=[aggregation_list[0], output_name])
process_queue.put(mod)
process_queue.wait()
if connection_state_changed:
dbif.close()
msgr.percent(1, 1, 1)
return output_list
