1#!/usr/local/bin/python3.8 2# -*- coding: utf-8 -*- 3############################################################################ 4# 5# MODULE: t.rast.what 6# AUTHOR(S): Soeren Gebbert 7# 8# PURPOSE: Sample a space time raster dataset at specific vector point 9# coordinates and write the output to stdout using different 10# layouts 11# COPYRIGHT: (C) 2015-2017 by the GRASS Development Team 12# 13# This program is free software; you can redistribute it and/or modify 14# it under the terms of the GNU General Public License as published by 15# the Free Software Foundation; either version 2 of the License, or 16# (at your option) any later version. 17# 18# This program is distributed in the hope that it will be useful, 19# but WITHOUT ANY WARRANTY; without even the implied warranty of 20# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 21# GNU General Public License for more details. 22# 23############################################################################# 24 25#%module 26#% description: Sample a space time raster dataset at specific vector point coordinates and write the output to stdout using different layouts 27#% keyword: temporal 28#% keyword: sampling 29#% keyword: raster 30#% keyword: time 31#%end 32 33#%option G_OPT_V_INPUT 34#% key: points 35#% required: no 36#%end 37 38#%option G_OPT_M_COORDS 39#% required: no 40#% description: Comma separated list of coordinates 41#%end 42 43#%option G_OPT_STRDS_INPUT 44#% key: strds 45#%end 46 47#%option G_OPT_F_OUTPUT 48#% required: no 49#% description: Name for the output file or "-" in case stdout should be used 50#% answer: - 51#%end 52 53#%option G_OPT_T_WHERE 54#%end 55 56#%option G_OPT_M_NULL_VALUE 57#%end 58 59#%option G_OPT_F_SEP 60#%end 61 62#%option 63#% key: order 64#% type: string 65#% description: Sort the maps by category 66#% required: no 67#% multiple: yes 68#% options: id, name, creator, mapset, creation_time, modification_time, start_time, end_time, north, south, west, east, min, max 69#% answer: start_time 70#%end 71 72#%option 73#% key: layout 74#% type: string 75#% description: The layout of the output. One point per row (row), one point per column (col), all timsteps in one row (timerow) 76#% required: no 77#% multiple: no 78#% options: row, col, timerow 79#% answer: row 80#%end 81 82#%option 83#% key: nprocs 84#% type: integer 85#% description: Number of r.what processes to run in parallel 86#% required: no 87#% multiple: no 88#% answer: 1 89#%end 90 91#%flag 92#% key: n 93#% description: Output header row 94#%end 95 96#%flag 97#% key: i 98#% description: Use stdin as input and ignore coordinates and point option 99#%end 100 101## Temporary disabled the r.what flags due to test issues 102##%flag 103##% key: f 104##% description: Show the category labels of the grid cell(s) 105##%end 106 107##%flag 108##% key: r 109##% description: Output color values as RRR:GGG:BBB 110##%end 111 112##%flag 113##% key: i 114##% description: Output integer category values, not cell values 115##%end 116 117#%flag 118#% key: v 119#% description: Show the category for vector points map 120#%end 121 122import sys 123import copy 124import grass.script as gscript 125 126import pdb 127 128############################################################################ 129 130def main(options, flags): 131 # lazy imports 132 import grass.temporal as tgis 133 import grass.pygrass.modules as pymod 134 135 # Get the options 136 points = options["points"] 137 coordinates = options["coordinates"] 138 strds = options["strds"] 139 output = options["output"] 140 where = options["where"] 141 order = options["order"] 142 layout = options["layout"] 143 null_value = options["null_value"] 144 separator = gscript.separator(options["separator"]) 145 146 nprocs = int(options["nprocs"]) 147 write_header = flags["n"] 148 use_stdin = flags["i"] 149 vcat = flags["v"] 150 151 #output_cat_label = flags["f"] 152 #output_color = flags["r"] 153 #output_cat = flags["i"] 154 155 overwrite = gscript.overwrite() 156 157 if coordinates and points: 158 gscript.fatal(_("Options coordinates and points are mutually exclusive")) 159 160 if not coordinates and not points and not use_stdin: 161 gscript.fatal(_("Please specify the coordinates, the points option or use the 'i' flag to pipe coordinate positions to t.rast.what from stdin, to provide the sampling coordinates")) 162 163 if vcat and not points: 164 gscript.fatal(_("Flag 'v' required option 'points'")) 165 166 if use_stdin: 167 coordinates_stdin = str(sys.__stdin__.read()) 168 # Check if coordinates are given with site names or IDs 169 stdin_length = len(coordinates_stdin.split('\n')[0].split()) 170 if stdin_length <= 2: 171 site_input = False 172 elif stdin_length >= 3: 173 site_input = True 174 else: 175 site_input = False 176 177 # Make sure the temporal database exists 178 tgis.init() 179 # We need a database interface 180 dbif = tgis.SQLDatabaseInterfaceConnection() 181 dbif.connect() 182 183 sp = tgis.open_old_stds(strds, "strds", dbif) 184 maps = sp.get_registered_maps_as_objects(where=where, order=order, 185 dbif=dbif) 186 dbif.close() 187 if not maps: 188 gscript.fatal(_("Space time raster dataset <%s> is empty") % sp.get_id()) 189 190 # Setup flags are disabled due to test issues 191 flags = "" 192 #if output_cat_label is True: 193 # flags += "f" 194 #if output_color is True: 195 # flags += "r" 196 #if output_cat is True: 197 # flags += "i" 198 if vcat is True: 199 flags += "v" 200 201 # Configure the r.what module 202 if points: 203 r_what = pymod.Module("r.what", map="dummy", 204 output="dummy", run_=False, 205 separator=separator, points=points, 206 overwrite=overwrite, flags=flags, 207 null_value=null_value, 208 quiet=True) 209 elif coordinates: 210 # Create a list of values 211 coord_list = coordinates.split(",") 212 r_what = pymod.Module("r.what", map="dummy", 213 output="dummy", run_=False, 214 separator=separator, 215 coordinates=coord_list, 216 overwrite=overwrite, flags=flags, 217 null_value=null_value, 218 quiet=True) 219 elif use_stdin: 220 r_what = pymod.Module("r.what", map="dummy", 221 output="dummy", run_=False, 222 separator=separator, 223 stdin_=coordinates_stdin, 224 overwrite=overwrite, flags=flags, 225 null_value=null_value, 226 quiet=True) 227 else: 228 gscript.error(_("Please specify points or coordinates")) 229 230 if len(maps) < nprocs: 231 nprocs = len(maps) 232 233 # The module queue for parallel execution 234 process_queue = pymod.ParallelModuleQueue(int(nprocs)) 235 num_maps = len(maps) 236 237 # 400 Maps is the absolute maximum in r.what 238 # We need to determie the number of maps that can be processed 239 # in parallel 240 241 # First estimate the number of maps per process. We use 400 maps 242 # simultaniously as maximum for a single process 243 244 num_loops = int(num_maps / (400 * nprocs)) 245 remaining_maps = num_maps % (400 * nprocs) 246 247 if num_loops == 0: 248 num_loops = 1 249 remaining_maps = 0 250 251 # Compute the number of maps for each process 252 maps_per_loop = int((num_maps - remaining_maps) / num_loops) 253 maps_per_process = int(maps_per_loop / nprocs) 254 remaining_maps_per_loop = maps_per_loop % nprocs 255 256 # We put the output files in an ordered list 257 output_files = [] 258 output_time_list = [] 259 260 count = 0 261 for loop in range(num_loops): 262 file_name = gscript.tempfile() + "_%i"%(loop) 263 count = process_loop(nprocs, maps, file_name, count, maps_per_process, 264 remaining_maps_per_loop, output_files, 265 output_time_list, r_what, process_queue) 266 267 process_queue.wait() 268 269 gscript.verbose("Number of raster map layers remaining for sampling %i"%(remaining_maps)) 270 if remaining_maps > 0: 271 # Use a single process if less then 100 maps 272 if remaining_maps <= 100: 273 map_names = [] 274 for i in range(remaining_maps): 275 map = maps[count] 276 map_names.append(map.get_id()) 277 count += 1 278 mod = copy.deepcopy(r_what) 279 mod(map=map_names, output=file_name) 280 process_queue.put(mod) 281 else: 282 maps_per_process = int(remaining_maps / nprocs) 283 remaining_maps_per_loop = remaining_maps % nprocs 284 285 file_name = "out_remain" 286 process_loop(nprocs, maps, file_name, count, maps_per_process, 287 remaining_maps_per_loop, output_files, 288 output_time_list, r_what, process_queue) 289 290 # Wait for unfinished processes 291 process_queue.wait() 292 293 # Out the output files in the correct order together 294 if layout == "row": 295 one_point_per_row_output(separator, output_files, output_time_list, 296 output, write_header, site_input, vcat) 297 elif layout == "col": 298 one_point_per_col_output(separator, output_files, output_time_list, 299 output, write_header, site_input, vcat) 300 else: 301 one_point_per_timerow_output(separator, output_files, output_time_list, 302 output, write_header, site_input, vcat) 303 304############################################################################ 305 306def one_point_per_row_output(separator, output_files, output_time_list, 307 output, write_header, site_input, vcat): 308 """Write one point per row 309 output is of type: x,y,start,end,value 310 """ 311 # open the output file for writing 312 out_file = open(output, 'w') if output != "-" else sys.stdout 313 314 if write_header is True: 315 out_str = "" 316 if vcat: 317 out_str += "cat{sep}" 318 if site_input: 319 out_str += "x{sep}y{sep}site{sep}start{sep}end{sep}value\n" 320 else: 321 out_str += "x{sep}y{sep}start{sep}end{sep}value\n" 322 out_file.write(out_str.format(sep=separator)) 323 324 for count in range(len(output_files)): 325 file_name = output_files[count] 326 gscript.verbose(_("Transforming r.what output file %s"%(file_name))) 327 map_list = output_time_list[count] 328 in_file = open(file_name, "r") 329 for line in in_file: 330 line = line.split(separator) 331 if vcat: 332 cat = line[0] 333 x = line[1] 334 y = line[2] 335 values = line[4:] 336 if site_input: 337 site = line[3] 338 values = line[5:] 339 340 else: 341 x = line[0] 342 y = line[1] 343 if site_input: 344 site = line[2] 345 values = line[3:] 346 347 for i in range(len(values)): 348 start, end = map_list[i].get_temporal_extent_as_tuple() 349 if vcat: 350 cat_str = "{ca}{sep}".format(ca=cat, sep=separator) 351 else: 352 cat_str = "" 353 if site_input: 354 coor_string = "%(x)10.10f%(sep)s%(y)10.10f%(sep)s%(site_name)s%(sep)s"\ 355 %({"x":float(x),"y":float(y),"site_name":str(site),"sep":separator}) 356 else: 357 coor_string = "%(x)10.10f%(sep)s%(y)10.10f%(sep)s"\ 358 %({"x":float(x),"y":float(y),"sep":separator}) 359 time_string = "%(start)s%(sep)s%(end)s%(sep)s%(val)s\n"\ 360 %({"start":str(start), "end":str(end), 361 "val":(values[i].strip()),"sep":separator}) 362 363 out_file.write(cat_str + coor_string + time_string) 364 365 in_file.close() 366 367 if out_file is not sys.stdout: 368 out_file.close() 369 370############################################################################ 371 372def one_point_per_col_output(separator, output_files, output_time_list, 373 output, write_header, site_input, vcat): 374 """Write one point per col 375 output is of type: 376 start,end,point_1 value,point_2 value,...,point_n value 377 378 Each row represents a single raster map, hence a single time stamp 379 """ 380 # open the output file for writing 381 out_file = open(output, 'w') if output != "-" else sys.stdout 382 383 first = True 384 for count in range(len(output_files)): 385 file_name = output_files[count] 386 gscript.verbose(_("Transforming r.what output file %s"%(file_name))) 387 map_list = output_time_list[count] 388 in_file = open(file_name, "r") 389 lines = in_file.readlines() 390 391 matrix = [] 392 for line in lines: 393 matrix.append(line.split(separator)) 394 395 num_cols = len(matrix[0]) 396 397 if first is True: 398 if write_header is True: 399 out_str = "start%(sep)send"%({"sep":separator}) 400 401 # Define different separator for coordinates and sites 402 if separator == ',': 403 coor_sep = ';' 404 else: 405 coor_sep = ',' 406 407 for row in matrix: 408 if vcat: 409 cat = row[0] 410 x = row[1] 411 y = row[2] 412 out_str += "{sep}{cat}{csep}{x:10.10f}{csep}" \ 413 "{y:10.10f}".format(cat=cat, x=float(x), 414 y=float(y), 415 sep=separator, 416 csep=coor_sep) 417 if site_input: 418 site = row[3] 419 out_str += "{sep}{site}".format(sep=coor_sep, 420 site=site) 421 else: 422 x = row[0] 423 y = row[1] 424 out_str += "{sep}{x:10.10f}{csep}" \ 425 "{y:10.10f}".format(x=float(x), y=float(y), 426 sep=separator, 427 csep=coor_sep) 428 if site_input: 429 site = row[2] 430 out_str += "{sep}{site}".format(sep=coor_sep, 431 site=site) 432 433 out_file.write(out_str + "\n") 434 435 first = False 436 437 if vcat: 438 ncol = 4 439 else: 440 ncol = 3 441 for col in range(num_cols - ncol): 442 start, end = output_time_list[count][col].get_temporal_extent_as_tuple() 443 time_string = "%(start)s%(sep)s%(end)s"\ 444 %({"start":str(start), "end":str(end), 445 "sep":separator}) 446 out_file.write(time_string) 447 for row in range(len(matrix)): 448 value = matrix[row][col + ncol] 449 out_file.write("%(sep)s%(value)s"\ 450 %({"sep":separator, 451 "value":value.strip()})) 452 out_file.write("\n") 453 454 in_file.close() 455 if out_file is not sys.stdout: 456 out_file.close() 457 458############################################################################ 459 460def one_point_per_timerow_output(separator, output_files, output_time_list, 461 output, write_header, site_input, vcat): 462 """Use the original layout of the r.what output and print instead of 463 the raster names, the time stamps as header 464 465 One point per line for all time stamps: 466 x|y|1991-01-01 00:00:00;1991-01-02 00:00:00|1991-01-02 00:00:00;1991-01-03 00:00:00|1991-01-03 00:00:00;1991-01-04 00:00:00|1991-01-04 00:00:00;1991-01-05 00:00:00 467 3730731.49590371|5642483.51236521|6|8|7|7 468 3581249.04638104|5634411.97526282|5|8|7|7 469 """ 470 out_file = open(output, 'w') if output != "-" else sys.stdout 471 472 matrix = [] 473 header = "" 474 475 first = True 476 for count in range(len(output_files)): 477 file_name = output_files[count] 478 gscript.verbose("Transforming r.what output file %s"%(file_name)) 479 map_list = output_time_list[count] 480 in_file = open(file_name, "r") 481 482 if write_header: 483 if first is True: 484 if vcat: 485 header = "cat{sep}".format(sep=separator) 486 else: 487 header = "" 488 if site_input: 489 header += "x%(sep)sy%(sep)ssite"%({"sep":separator}) 490 else: 491 header += "x%(sep)sy"%({"sep":separator}) 492 for map in map_list: 493 start, end = map.get_temporal_extent_as_tuple() 494 time_string = "%(sep)s%(start)s;%(end)s"\ 495 %({"start":str(start), "end":str(end), 496 "sep":separator}) 497 header += time_string 498 499 lines = in_file.readlines() 500 501 for i in range(len(lines)): 502 cols = lines[i].split(separator) 503 504 if first is True: 505 if vcat and site_input: 506 matrix.append(cols[:4]) 507 elif vcat or site_input: 508 matrix.append(cols[:3]) 509 else: 510 matrix.append(cols[:2]) 511 512 if vcat: 513 matrix[i] = matrix[i] + cols[4:] 514 else: 515 matrix[i] = matrix[i] + cols[3:] 516 517 first = False 518 519 in_file.close() 520 521 if write_header: 522 out_file.write(header + "\n") 523 524 gscript.verbose(_("Writing the output file <%s>"%(output))) 525 for row in matrix: 526 first = True 527 for col in row: 528 value = col.strip() 529 530 if first is False: 531 out_file.write("%s"%(separator)) 532 out_file.write(value) 533 534 first = False 535 536 out_file.write("\n") 537 if out_file is not sys.stdout: 538 out_file.close() 539 540############################################################################ 541 542def process_loop(nprocs, maps, file_name, count, maps_per_process, 543 remaining_maps_per_loop, output_files, 544 output_time_list, r_what, process_queue): 545 """Call r.what in parallel subprocesses""" 546 first = True 547 for process in range(nprocs): 548 num = maps_per_process 549 # Add the remaining maps to the first process 550 if first is True: 551 num += remaining_maps_per_loop 552 first = False 553 554 # Temporary output file 555 final_file_name = file_name + "_%i"%(process) 556 output_files.append(final_file_name) 557 558 map_names = [] 559 map_list = [] 560 for i in range(num): 561 map = maps[count] 562 map_names.append(map.get_id()) 563 map_list.append(map) 564 count += 1 565 566 output_time_list.append(map_list) 567 568 gscript.verbose(_("Process maps %(samp_start)i to %(samp_end)i (of %(total)i)"\ 569 %({"samp_start":count-len(map_names)+1, 570 "samp_end":count, "total":len(maps)}))) 571 mod = copy.deepcopy(r_what) 572 mod(map=map_names, output=final_file_name) 573 #print(mod.get_bash()) 574 process_queue.put(mod) 575 576 return count 577 578############################################################################ 579 580if __name__ == "__main__": 581 options, flags = gscript.parser() 582 main(options, flags) 583