| 1 | #!/usr/bin/env python
|
|---|
| 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 |
|
|---|
| 122 | import sys
|
|---|
| 123 | import copy
|
|---|
| 124 | import grass.script as gscript
|
|---|
| 125 |
|
|---|
| 126 | import pdb
|
|---|
| 127 |
|
|---|
| 128 | ############################################################################
|
|---|
| 129 |
|
|---|
| 130 | def 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 |
|
|---|
| 306 | def 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 |
|
|---|
| 372 | def 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 |
|
|---|
| 460 | def 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 |
|
|---|
| 542 | def 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 |
|
|---|
| 580 | if __name__ == "__main__":
|
|---|
| 581 | options, flags = gscript.parser()
|
|---|
| 582 | main(options, flags)
|
|---|
