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