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