1"""
2Aggregation methods for space time raster datasets
3
4Usage:
5
6.. code-block:: python
7
8    import grass.temporal as tgis
9
10    tgis.aggregate_raster_maps(dataset, mapset, inputs, base, start, end, count, method, register_null, dbif)
11
12(C) 2012-2013 by the GRASS Development Team
13This program is free software under the GNU General Public
14License (>=v2). Read the file COPYING that comes with GRASS
15for details.
16
17:author: Soeren Gebbert
18"""
19
20import grass.script as gscript
21from grass.exceptions import CalledModuleError
22from .space_time_datasets import RasterDataset
23from .datetime_math import create_suffix_from_datetime
24from .datetime_math import create_time_suffix
25from .datetime_math import create_numeric_suffix
26from .core import get_current_mapset, get_tgis_message_interface, init_dbif
27from .spatio_temporal_relationships import SpatioTemporalTopologyBuilder, \
28    create_temporal_relation_sql_where_statement
29###############################################################################
30
31
32def collect_map_names(sp, dbif, start, end, sampling):
33    """Gather all maps from dataset using a specific sample method
34
35       :param sp: The space time raster dataset to select aps from
36       :param dbif: The temporal database interface to use
37       :param start: The start time of the sample interval, may be relative or
38              absolute
39       :param end: The end time of the sample interval, may be relative or
40              absolute
41       :param sampling: The sampling methods to use
42    """
43
44    use_start = False
45    use_during = False
46    use_overlap = False
47    use_contain = False
48    use_equal = False
49    use_follows = False
50    use_precedes = False
51
52    # Initialize the methods
53    if sampling:
54        for name in sampling.split(","):
55            if name == "start":
56                use_start = True
57            if name == "during":
58                use_during = True
59            if name == "overlap":
60                use_overlap = True
61            if name == "contain":
62                use_contain = True
63            if name == "equal":
64                use_equal = True
65            if name == "follows":
66                use_follows = True
67            if name == "precedes":
68                use_precedes = True
69
70    else:
71        use_start = True
72
73    if sp.get_map_time() != "interval":
74        use_start = True
75        use_during = False
76        use_overlap = False
77        use_contain = False
78        use_equal = False
79        use_follows = False
80        use_precedes = False
81
82    where = create_temporal_relation_sql_where_statement(start, end,
83                                                         use_start,
84                                                         use_during,
85                                                         use_overlap,
86                                                         use_contain,
87                                                         use_equal,
88                                                         use_follows,
89                                                         use_precedes)
90
91    rows = sp.get_registered_maps("id", where, "start_time", dbif)
92
93    if not rows:
94        return None
95
96    names = []
97    for row in rows:
98        names.append(row["id"])
99
100    return names
101
102###############################################################################
103
104
105def aggregate_raster_maps(inputs, base, start, end, count, method,
106                          register_null, dbif,  offset=0):
107    """Aggregate a list of raster input maps with r.series
108
109       :param inputs: The names of the raster maps to be aggregated
110       :param base: The basename of the new created raster maps
111       :param start: The start time of the sample interval, may be relative or
112                    absolute
113       :param end: The end time of the sample interval, may be relative or
114                  absolute
115       :param count: The number to be attached to the basename of the new
116                    created raster map
117       :param method: The aggreation method to be used by r.series
118       :param register_null: If true null maps will be registered in the space
119                            time raster dataset, if false not
120       :param dbif: The temporal database interface to use
121       :param offset: Offset to be added to the map counter to create the map ids
122    """
123
124    msgr = get_tgis_message_interface()
125
126    msgr.verbose(_("Aggregating %s raster maps") % (len(inputs)))
127    output = "%s_%i" % (base, int(offset) + count)
128
129    mapset = get_current_mapset()
130    map_id = output + "@" + mapset
131    new_map = RasterDataset(map_id)
132
133    # Check if new map is in the temporal database
134    if new_map.is_in_db(dbif):
135        if gscript.overwrite() is True:
136            # Remove the existing temporal database entry
137            new_map.delete(dbif)
138            new_map = RasterDataset(map_id)
139        else:
140            msgr.error(_("Raster map <%(name)s> is already in temporal "
141                         "database, use overwrite flag to overwrite" %
142                       ({"name": new_map.get_name()})))
143            return
144
145    msgr.verbose(_("Computing aggregation of maps between %(st)s - %(end)s" % {
146                   'st': str(start), 'end': str(end)}))
147
148    # Create the r.series input file
149    filename = gscript.tempfile(True)
150    file = open(filename, 'w')
151
152    for name in inputs:
153        string = "%s\n" % (name)
154        file.write(string)
155
156    file.close()
157    # Run r.series
158    try:
159        if len(inputs) > 1000:
160            gscript.run_command("r.series", flags="z", file=filename,
161                                output=output, overwrite=gscript.overwrite(),
162                                method=method)
163        else:
164            gscript.run_command("r.series", file=filename,
165                                output=output, overwrite=gscript.overwrite(),
166                                method=method)
167
168    except CalledModuleError:
169        dbif.close()
170        msgr.fatal(_("Error occurred in r.series computation"))
171
172    # Read the raster map data
173    new_map.load()
174
175    # In case of a null map continue, do not register null maps
176    if new_map.metadata.get_min() is None and \
177       new_map.metadata.get_max() is None:
178        if not register_null:
179            gscript.run_command("g.remove", flags='f', type='raster',
180                                name=output)
181            return None
182
183    return new_map
184
185##############################################################################
186
187
188def aggregate_by_topology(granularity_list, granularity, map_list, topo_list,
189                          basename, time_suffix, offset=0, method="average",
190                          nprocs=1, spatial=None, dbif=None, overwrite=False,
191                          file_limit=1000):
192    """Aggregate a list of raster input maps with r.series
193
194       :param granularity_list: A list of AbstractMapDataset objects.
195                                The temporal extents of the objects are used
196                                to build the spatio-temporal topology with the
197                                map list objects
198       :param granularity: The granularity of the granularity list
199       :param map_list: A list of RasterDataset objects that contain the raster
200                        maps that should be aggregated
201       :param topo_list: A list of strings of topological relations that are
202                         used to select the raster maps for aggregation
203       :param basename: The basename of the new generated raster maps
204       :param time_suffix: Use the granularity truncated start time of the
205                           actual granule to create the suffix for the basename
206       :param offset: Use a numerical offset for suffix generation
207                      (overwritten by time_suffix)
208       :param method: The aggregation method of r.series (average,min,max, ...)
209       :param nprocs: The number of processes used for parallel computation
210       :param spatial: This indicates if the spatial topology is created as
211                       well: spatial can be None (no spatial topology), "2D"
212                       using west, east, south, north or "3D" using west,
213                       east, south, north, bottom, top
214       :param dbif: The database interface to be used
215       :param overwrite: Overwrite existing raster maps
216       :param file_limit: The maximum number of raster map layers that
217                          should be opened at once by r.series
218       :return: A list of RasterDataset objects that contain the new map names
219                and the temporal extent for map registration
220    """
221    import grass.pygrass.modules as pymod
222    import copy
223
224    msgr = get_tgis_message_interface()
225
226    dbif, connected = init_dbif(dbif)
227
228    topo_builder = SpatioTemporalTopologyBuilder()
229    topo_builder.build(mapsA=granularity_list, mapsB=map_list, spatial=spatial)
230
231    # The module queue for parallel execution
232    process_queue = pymod.ParallelModuleQueue(int(nprocs))
233
234    # Dummy process object that will be deep copied
235    # and be put into the process queue
236    r_series = pymod.Module("r.series", output="spam", method=[method],
237                            overwrite=overwrite, quiet=True, run_=False,
238                            finish_=False)
239    g_copy = pymod.Module("g.copy", raster=['spam', 'spamspam'],
240                          quiet=True, run_=False, finish_=False)
241    output_list = []
242    count = 0
243
244    for granule in granularity_list:
245        msgr.percent(count, len(granularity_list), 1)
246        count += 1
247
248        aggregation_list = []
249
250        if "equal" in topo_list and granule.equal:
251            for map_layer in granule.equal:
252                aggregation_list.append(map_layer.get_name())
253        if "contains" in topo_list and granule.contains:
254            for map_layer in granule.contains:
255                aggregation_list.append(map_layer.get_name())
256        if "during" in topo_list and granule.during:
257            for map_layer in granule.during:
258                aggregation_list.append(map_layer.get_name())
259        if "starts" in topo_list and granule.starts:
260            for map_layer in granule.starts:
261                aggregation_list.append(map_layer.get_name())
262        if "started" in topo_list and granule.started:
263            for map_layer in granule.started:
264                aggregation_list.append(map_layer.get_name())
265        if "finishes" in topo_list and granule.finishes:
266            for map_layer in granule.finishes:
267                aggregation_list.append(map_layer.get_name())
268        if "finished" in topo_list and granule.finished:
269            for map_layer in granule.finished:
270                aggregation_list.append(map_layer.get_name())
271        if "overlaps" in topo_list and granule.overlaps:
272            for map_layer in granule.overlaps:
273                aggregation_list.append(map_layer.get_name())
274        if "overlapped" in topo_list and granule.overlapped:
275            for map_layer in granule.overlapped:
276                aggregation_list.append(map_layer.get_name())
277
278        if aggregation_list:
279            msgr.verbose(_("Aggregating %(len)i raster maps from %(start)s to"
280                           " %(end)s")  %({"len": len(aggregation_list),
281                           "start": str(granule.temporal_extent.get_start_time()),
282                           "end": str(granule.temporal_extent.get_end_time())}))
283
284            if granule.is_time_absolute() is True and time_suffix == 'gran':
285                suffix = create_suffix_from_datetime(granule.temporal_extent.get_start_time(),
286                                                     granularity)
287                output_name = "{ba}_{su}".format(ba=basename, su=suffix)
288            elif granule.is_time_absolute() is True and time_suffix == 'time':
289                suffix = create_time_suffix(granule)
290                output_name = "{ba}_{su}".format(ba=basename, su=suffix)
291            else:
292                output_name = create_numeric_suffix(basename, count + int(offset),
293                                                    time_suffix)
294
295            map_layer = RasterDataset("%s@%s" % (output_name,
296                                                 get_current_mapset()))
297            map_layer.set_temporal_extent(granule.get_temporal_extent())
298
299            if map_layer.map_exists() is True and overwrite is False:
300                msgr.fatal(_("Unable to perform aggregation. Output raster "
301                             "map <%(name)s> exists and overwrite flag was "
302                             "not set" % ({"name": output_name})))
303
304            output_list.append(map_layer)
305
306            if len(aggregation_list) > 1:
307                # Create the r.series input file
308                filename = gscript.tempfile(True)
309                file = open(filename, 'w')
310                for name in aggregation_list:
311                    string = "%s\n" % (name)
312                    file.write(string)
313                file.close()
314
315                mod = copy.deepcopy(r_series)
316                mod(file=filename, output=output_name)
317                if len(aggregation_list) > int(file_limit):
318                    msgr.warning(_("The limit of open files (%i) was "\
319                                   "reached (%i). The module r.series will "\
320                                   "be run with flag z, to avoid open "\
321                                   "files limit exceeding."%(int(file_limit),
322                                                             len(aggregation_list))))
323                    mod(flags="z")
324                process_queue.put(mod)
325            else:
326                mod = copy.deepcopy(g_copy)
327                mod(raster=[aggregation_list[0],  output_name])
328                process_queue.put(mod)
329
330    process_queue.wait()
331
332    if connected:
333        dbif.close()
334
335    msgr.percent(1, 1, 1)
336
337    return output_list
338