Source code for grass.temporal.aggregation

"""
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 aggregation 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