1"""
2Extract functions for space time raster, 3d raster and vector datasets
3
4(C) 2012-2013 by the GRASS Development Team
5This program is free software under the GNU General Public
6License (>=v2). Read the file COPYING that comes with GRASS
7for details.
8
9:authors: Soeren Gebbert
10"""
11from .core import get_tgis_message_interface, get_current_mapset, SQLDatabaseInterfaceConnection
12from .abstract_map_dataset import AbstractMapDataset
13from .open_stds import open_old_stds, check_new_stds, open_new_stds
14from .datetime_math import create_suffix_from_datetime
15from .datetime_math import create_time_suffix
16from .datetime_math import create_numeric_suffix
17from multiprocessing import Process
18import grass.script as gscript
19from grass.exceptions import CalledModuleError
20
21############################################################################
22
23
24def extract_dataset(input, output, type, where, expression, base, time_suffix,
25                    nprocs=1, register_null=False, layer=1,
26                    vtype="point,line,boundary,centroid,area,face", ):
27    """Extract a subset of a space time raster, raster3d or vector dataset
28
29       A mapcalc expression can be provided to process the temporal extracted
30       maps.
31       Mapcalc expressions are supported for raster and raster3d maps.
32
33       :param input: The name of the input space time raster/raster3d dataset
34       :param output: The name of the extracted new space time raster/raster3d
35                     dataset
36       :param type: The type of the dataset: "raster", "raster3d" or vector
37       :param where: The temporal SQL WHERE statement for subset extraction
38       :param expression: The r(3).mapcalc expression or the v.extract where
39                         statement
40       :param base: The base name of the new created maps in case a mapclac
41                   expression is provided
42       :param time_suffix: string to choose which suffix to use: gran, time, num%*
43                          (where * are digits)
44       :param nprocs: The number of parallel processes to be used for mapcalc
45                     processing
46       :param register_null: Set this number True to register empty maps
47                            (only raster and raster3d maps)
48       :param layer: The vector layer number to be used when no timestamped
49              layer is present, default is 1
50       :param vtype: The feature type to be extracted for vector maps, default
51              is point,line,boundary,centroid,area and face
52    """
53
54    # Check the parameters
55    msgr = get_tgis_message_interface()
56
57    if expression and not base:
58        msgr.fatal(_("You need to specify the base name of new created maps"))
59
60    mapset = get_current_mapset()
61
62    dbif = SQLDatabaseInterfaceConnection()
63    dbif.connect()
64
65    sp = open_old_stds(input, type, dbif)
66    # Check the new stds
67    new_sp = check_new_stds(output, type, dbif, gscript.overwrite())
68    if type == "vector":
69        rows = sp.get_registered_maps(
70            "id,name,mapset,layer", where, "start_time", dbif)
71    else:
72        rows = sp.get_registered_maps("id", where, "start_time", dbif)
73
74    new_maps = {}
75    if rows:
76        num_rows = len(rows)
77
78        msgr.percent(0, num_rows, 1)
79
80        # Run the mapcalc expression
81        if expression:
82            count = 0
83            proc_count = 0
84            proc_list = []
85
86            for row in rows:
87                count += 1
88
89                if count % 10 == 0:
90                    msgr.percent(count, num_rows, 1)
91
92                if sp.get_temporal_type() == 'absolute' and time_suffix == 'gran':
93                    old_map = sp.get_new_map_instance(row["id"])
94                    old_map.select(dbif)
95                    suffix = create_suffix_from_datetime(old_map.temporal_extent.get_start_time(),
96                                                         sp.get_granularity())
97                    map_name = "{ba}_{su}".format(ba=base, su=suffix)
98                elif sp.get_temporal_type() == 'absolute' and time_suffix == 'time':
99                    old_map = sp.get_new_map_instance(row["id"])
100                    old_map.select(dbif)
101                    suffix = create_time_suffix(old_map)
102                    map_name = "{ba}_{su}".format(ba=base, su=suffix)
103                else:
104                    map_name = create_numeric_suffix(base, count, time_suffix)
105
106                # We need to modify the r(3).mapcalc expression
107                if type != "vector":
108                    expr = expression
109                    expr = expr.replace(sp.base.get_map_id(), row["id"])
110                    expr = expr.replace(sp.base.get_name(), row["id"])
111
112                    expr = "%s = %s" % (map_name, expr)
113
114                    # We need to build the id
115                    map_id = AbstractMapDataset.build_id(map_name, mapset)
116                else:
117                    map_id = AbstractMapDataset.build_id(map_name, mapset,
118                                                         row["layer"])
119
120                new_map = sp.get_new_map_instance(map_id)
121
122                # Check if new map is in the temporal database
123                if new_map.is_in_db(dbif):
124                    if gscript.overwrite():
125                        # Remove the existing temporal database entry
126                        new_map.delete(dbif)
127                        new_map = sp.get_new_map_instance(map_id)
128                    else:
129                        msgr.error(_("Map <%s> is already in temporal database"
130                                     ", use overwrite flag to overwrite") %
131                                    (new_map.get_map_id()))
132                        continue
133
134                # Add process to the process list
135                if type == "raster":
136                    msgr.verbose(_("Applying r.mapcalc expression: \"%s\"")
137                                 % expr)
138                    proc_list.append(Process(target=run_mapcalc2d,
139                                             args=(expr,)))
140                elif type == "raster3d":
141                    msgr.verbose(_("Applying r3.mapcalc expression: \"%s\"")
142                                 % expr)
143                    proc_list.append(Process(target=run_mapcalc3d,
144                                             args=(expr,)))
145                elif type == "vector":
146                    msgr.verbose(_("Applying v.extract where statement: \"%s\"")
147                                 % expression)
148                    if row["layer"]:
149                        proc_list.append(Process(target=run_vector_extraction,
150                                                 args=(row["name"] + "@" +
151                                                       row["mapset"], map_name,
152                                                       row["layer"], vtype,
153                                                       expression)))
154                    else:
155                        proc_list.append(Process(target=run_vector_extraction,
156                                                 args=(row["name"] + "@" +
157                                                       row["mapset"], map_name,
158                                                       layer, vtype,
159                                                       expression)))
160
161                proc_list[proc_count].start()
162                proc_count += 1
163
164                # Join processes if the maximum number of processes are
165                # reached or the end of the loop is reached
166                if proc_count == nprocs or count == num_rows:
167                    proc_count = 0
168                    exitcodes = 0
169                    for proc in proc_list:
170                        proc.join()
171                        exitcodes += proc.exitcode
172                    if exitcodes != 0:
173                        dbif.close()
174                        msgr.fatal(_("Error in computation process"))
175
176                    # Empty process list
177                    proc_list = []
178
179                # Store the new maps
180                new_maps[row["id"]] = new_map
181
182        msgr.percent(0, num_rows, 1)
183
184        temporal_type, semantic_type, title, description = sp.get_initial_values()
185        new_sp = open_new_stds(output, type, sp.get_temporal_type(), title,
186                               description, semantic_type, dbif,
187                               gscript.overwrite())
188
189        # collect empty maps to remove them
190        empty_maps = []
191
192        # Register the maps in the database
193        count = 0
194        for row in rows:
195            count += 1
196
197            if count % 10 == 0:
198                msgr.percent(count, num_rows, 1)
199
200            old_map = sp.get_new_map_instance(row["id"])
201            old_map.select(dbif)
202
203            if expression:
204                # Register the new maps
205                if row["id"] in new_maps:
206                    new_map = new_maps[row["id"]]
207
208                    # Read the raster map data
209                    new_map.load()
210
211                    # In case of a empty map continue, do not register empty
212                    # maps
213                    if type == "raster" or type == "raster3d":
214                        if new_map.metadata.get_min() is None and \
215                            new_map.metadata.get_max() is None:
216                            if not register_null:
217                                empty_maps.append(new_map)
218                                continue
219                    elif type == "vector":
220                        if new_map.metadata.get_number_of_primitives() == 0 or \
221                           new_map.metadata.get_number_of_primitives() is None:
222                            if not register_null:
223                                empty_maps.append(new_map)
224                                continue
225
226                    # Set the time stamp
227                    new_map.set_temporal_extent(old_map.get_temporal_extent())
228
229                    # Insert map in temporal database
230                    new_map.insert(dbif)
231
232                    new_sp.register_map(new_map, dbif)
233            else:
234                new_sp.register_map(old_map, dbif)
235
236        # Update the spatio-temporal extent and the metadata table entries
237        new_sp.update_from_registered_maps(dbif)
238
239        msgr.percent(num_rows, num_rows, 1)
240
241        # Remove empty maps
242        if len(empty_maps) > 0:
243            names = ""
244            count = 0
245            for map in empty_maps:
246                if count == 0:
247                    names += "%s" % (map.get_name())
248                else:
249                    names += ",%s" % (map.get_name())
250                count += 1
251            if type == "raster":
252                gscript.run_command("g.remove", flags='f', type='raster',
253                                    name=names, quiet=True)
254            elif type == "raster3d":
255                gscript.run_command("g.remove", flags='f', type='raster_3d',
256                                    name=names, quiet=True)
257            elif type == "vector":
258                gscript.run_command("g.remove", flags='f', type='vector',
259                                    name=names, quiet=True)
260
261    dbif.close()
262
263###############################################################################
264
265
266def run_mapcalc2d(expr):
267    """Helper function to run r.mapcalc in parallel"""
268    try:
269        gscript.run_command("r.mapcalc", expression=expr,
270                            overwrite=gscript.overwrite(), quiet=True)
271    except CalledModuleError:
272        exit(1)
273
274
275def run_mapcalc3d(expr):
276    """Helper function to run r3.mapcalc in parallel"""
277    try:
278        gscript.run_command("r3.mapcalc", expression=expr,
279                            overwrite=gscript.overwrite(), quiet=True)
280    except CalledModuleError:
281        exit(1)
282
283
284def run_vector_extraction(input, output, layer, type, where):
285    """Helper function to run r.mapcalc in parallel"""
286    try:
287        gscript.run_command("v.extract", input=input, output=output,
288                            layer=layer, type=type, where=where,
289                            overwrite=gscript.overwrite(), quiet=True)
290    except CalledModuleError:
291        exit(1)
292