[7672] | 1 | """Gauge functions |
---|
| 2 | |
---|
[7679] | 3 | High-level functions for converting gauge and sww files into timeseries plots. |
---|
[7672] | 4 | |
---|
| 5 | |
---|
| 6 | Copyright 2010 |
---|
| 7 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou, James Hudson |
---|
| 8 | Geoscience Australia |
---|
| 9 | """ |
---|
| 10 | |
---|
| 11 | import numpy as num |
---|
| 12 | |
---|
| 13 | from anuga.geospatial_data.geospatial_data import ensure_absolute |
---|
[7691] | 14 | from util import check_list, calc_bearing |
---|
[7679] | 15 | from file_function import file_function |
---|
[7672] | 16 | |
---|
| 17 | import os |
---|
| 18 | |
---|
| 19 | from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep, getcwd |
---|
| 20 | from os.path import exists, split, join |
---|
| 21 | import anuga.utilities.log as log |
---|
| 22 | |
---|
| 23 | from math import sqrt |
---|
| 24 | |
---|
[7693] | 25 | def _quantities2csv(quantities, point_quantities, centroids, point_i): |
---|
| 26 | points_list = [] |
---|
| 27 | |
---|
| 28 | for quantity in quantities: |
---|
| 29 | #core quantities that are exported from the interpolator |
---|
| 30 | if quantity == 'stage': |
---|
| 31 | points_list.append(point_quantities[0]) |
---|
| 32 | |
---|
| 33 | if quantity == 'elevation': |
---|
| 34 | points_list.append(point_quantities[1]) |
---|
| 35 | |
---|
| 36 | if quantity == 'xmomentum': |
---|
| 37 | points_list.append(point_quantities[2]) |
---|
| 38 | |
---|
| 39 | if quantity == 'ymomentum': |
---|
| 40 | points_list.append(point_quantities[3]) |
---|
| 41 | |
---|
| 42 | #derived quantities that are calculated from the core ones |
---|
| 43 | if quantity == 'depth': |
---|
| 44 | points_list.append(point_quantities[0] |
---|
| 45 | - point_quantities[1]) |
---|
| 46 | |
---|
| 47 | if quantity == 'momentum': |
---|
| 48 | momentum = sqrt(point_quantities[2]**2 |
---|
| 49 | + point_quantities[3]**2) |
---|
| 50 | points_list.append(momentum) |
---|
| 51 | |
---|
| 52 | if quantity == 'speed': |
---|
| 53 | #if depth is less than 0.001 then speed = 0.0 |
---|
| 54 | if point_quantities[0] - point_quantities[1] < 0.001: |
---|
| 55 | vel = 0.0 |
---|
| 56 | else: |
---|
| 57 | if point_quantities[2] < 1.0e6: |
---|
| 58 | momentum = sqrt(point_quantities[2]**2 |
---|
| 59 | + point_quantities[3]**2) |
---|
| 60 | vel = momentum / (point_quantities[0] |
---|
| 61 | - point_quantities[1]) |
---|
| 62 | else: |
---|
| 63 | momentum = 0 |
---|
| 64 | vel = 0 |
---|
| 65 | |
---|
| 66 | points_list.append(vel) |
---|
| 67 | |
---|
| 68 | if quantity == 'bearing': |
---|
| 69 | points_list.append(calc_bearing(point_quantities[2], |
---|
| 70 | point_quantities[3])) |
---|
| 71 | if quantity == 'xcentroid': |
---|
| 72 | points_list.append(centroids[point_i][0]) |
---|
| 73 | |
---|
| 74 | if quantity == 'ycentroid': |
---|
| 75 | points_list.append(centroids[point_i][1]) |
---|
| 76 | |
---|
| 77 | return points_list |
---|
| 78 | |
---|
| 79 | |
---|
[7672] | 80 | ## |
---|
| 81 | # @brief Take gauge readings, given a gauge file and a sww mesh |
---|
| 82 | # |
---|
| 83 | # Use this function to take a timeseries sample, given a list of gauge points |
---|
| 84 | # @param sww_file sww file to use as input |
---|
| 85 | # @param gauge_file gauge file as input, containing list of gauge points to sample |
---|
| 86 | # @param out_name output file prefix |
---|
| 87 | # @param quantities which quantities in the sww file we want to export |
---|
| 88 | # @param verbose show extra logging information for debug purposes |
---|
| 89 | # @param use_cache cache requests if possible, for speed |
---|
| 90 | # @param output_centroids Set to true to output the values at the centroid of the mesh triangle |
---|
[7679] | 91 | def sww2csv_gauges(sww_file, |
---|
[7672] | 92 | gauge_file, |
---|
| 93 | out_name='gauge_', |
---|
| 94 | quantities=['stage', 'depth', 'elevation', |
---|
| 95 | 'xmomentum', 'ymomentum'], |
---|
| 96 | verbose=False, |
---|
| 97 | use_cache=True, |
---|
| 98 | output_centroids=False): |
---|
| 99 | """ |
---|
| 100 | |
---|
| 101 | Inputs: |
---|
| 102 | NOTE: if using csv2timeseries_graphs after creating csv file, |
---|
| 103 | it is essential to export quantities 'depth' and 'elevation'. |
---|
| 104 | 'depth' is good to analyse gauges on land and elevation is used |
---|
| 105 | automatically by csv2timeseries_graphs in the legend. |
---|
| 106 | |
---|
| 107 | sww_file: path to any sww file |
---|
| 108 | |
---|
| 109 | gauge_file: Assumes that it follows this format |
---|
| 110 | name, easting, northing, elevation |
---|
| 111 | point1, 100.3, 50.2, 10.0 |
---|
| 112 | point2, 10.3, 70.3, 78.0 |
---|
| 113 | |
---|
| 114 | NOTE: order of column can change but names eg 'easting', 'elevation' |
---|
| 115 | must be the same! ALL lowercaps! |
---|
| 116 | |
---|
| 117 | out_name: prefix for output file name (default is 'gauge_') |
---|
| 118 | |
---|
| 119 | Outputs: |
---|
| 120 | one file for each gauge/point location in the points file. They |
---|
| 121 | will be named with this format in the same directory as the 'sww_file' |
---|
| 122 | <out_name><name>.csv |
---|
| 123 | eg gauge_point1.csv if <out_name> not supplied |
---|
| 124 | myfile_2_point1.csv if <out_name> ='myfile_2_' |
---|
| 125 | |
---|
| 126 | They will all have a header |
---|
| 127 | |
---|
| 128 | Usage: sww2csv_gauges(sww_file='test1.sww', |
---|
| 129 | quantities = ['stage', 'elevation','depth','bearing'], |
---|
| 130 | gauge_file='gauge.txt') |
---|
| 131 | |
---|
| 132 | Interpolate the quantities at a given set of locations, given |
---|
| 133 | an sww file. |
---|
| 134 | The results are written to a csv file. |
---|
| 135 | |
---|
| 136 | In the future let points be a points file. |
---|
| 137 | And the user choose the quantities. |
---|
| 138 | |
---|
| 139 | This is currently quite specific. |
---|
| 140 | If it needs to be more general, change things. |
---|
| 141 | |
---|
| 142 | This is really returning speed, not velocity. |
---|
| 143 | """ |
---|
| 144 | |
---|
| 145 | from csv import reader,writer |
---|
| 146 | from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN |
---|
| 147 | import string |
---|
[7780] | 148 | from anuga.utilities.file_utils import get_all_swwfiles |
---|
[7685] | 149 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
[7672] | 150 | |
---|
| 151 | assert type(gauge_file) == type(''), 'Gauge filename must be a string' |
---|
| 152 | assert type(out_name) == type(''), 'Output filename prefix must be a string' |
---|
| 153 | |
---|
| 154 | try: |
---|
| 155 | point_reader = reader(file(gauge_file)) |
---|
| 156 | except Exception, e: |
---|
| 157 | msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e) |
---|
| 158 | raise msg |
---|
| 159 | |
---|
| 160 | if verbose: log.critical('Gauges obtained from: %s' % gauge_file) |
---|
| 161 | |
---|
| 162 | point_reader = reader(file(gauge_file)) |
---|
| 163 | points = [] |
---|
| 164 | point_name = [] |
---|
| 165 | |
---|
| 166 | # read point info from file |
---|
| 167 | for i,row in enumerate(point_reader): |
---|
| 168 | # read header and determine the column numbers to read correctly. |
---|
| 169 | if i==0: |
---|
| 170 | for j,value in enumerate(row): |
---|
| 171 | if value.strip()=='easting':easting=j |
---|
| 172 | if value.strip()=='northing':northing=j |
---|
| 173 | if value.strip()=='name':name=j |
---|
| 174 | if value.strip()=='elevation':elevation=j |
---|
| 175 | else: |
---|
| 176 | points.append([float(row[easting]),float(row[northing])]) |
---|
| 177 | point_name.append(row[name]) |
---|
| 178 | |
---|
| 179 | #convert to array for file_function |
---|
| 180 | points_array = num.array(points,num.float) |
---|
| 181 | |
---|
| 182 | points_array = ensure_absolute(points_array) |
---|
| 183 | |
---|
| 184 | dir_name, base = os.path.split(sww_file) |
---|
| 185 | |
---|
| 186 | #need to get current directory so when path and file |
---|
| 187 | #are "joined" below the directory is correct |
---|
| 188 | if dir_name == '': |
---|
| 189 | dir_name =getcwd() |
---|
| 190 | |
---|
| 191 | if access(sww_file,R_OK): |
---|
| 192 | if verbose: log.critical('File %s exists' % sww_file) |
---|
| 193 | else: |
---|
| 194 | msg = 'File "%s" could not be opened: no read permission' % sww_file |
---|
| 195 | raise msg |
---|
| 196 | |
---|
| 197 | sww_files = get_all_swwfiles(look_in_dir=dir_name, |
---|
| 198 | base_name=base, |
---|
| 199 | verbose=verbose) |
---|
| 200 | |
---|
| 201 | # fudge to get SWW files in 'correct' order, oldest on the left |
---|
| 202 | sww_files.sort() |
---|
| 203 | |
---|
| 204 | if verbose: |
---|
| 205 | log.critical('sww files=%s' % sww_files) |
---|
| 206 | |
---|
| 207 | #to make all the quantities lower case for file_function |
---|
| 208 | quantities = [quantity.lower() for quantity in quantities] |
---|
| 209 | |
---|
| 210 | # what is quantities are needed from sww file to calculate output quantities |
---|
| 211 | # also |
---|
| 212 | |
---|
| 213 | core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
| 214 | |
---|
| 215 | gauge_file = out_name |
---|
| 216 | |
---|
| 217 | heading = [quantity for quantity in quantities] |
---|
| 218 | heading.insert(0,'time') |
---|
| 219 | heading.insert(1,'hours') |
---|
| 220 | |
---|
| 221 | if verbose: log.critical('Writing csv files') |
---|
| 222 | |
---|
| 223 | quake_offset_time = None |
---|
| 224 | |
---|
| 225 | for sww_file in sww_files: |
---|
| 226 | sww_file = join(dir_name, sww_file+'.sww') |
---|
| 227 | callable_sww = file_function(sww_file, |
---|
| 228 | quantities=core_quantities, |
---|
| 229 | interpolation_points=points_array, |
---|
| 230 | verbose=verbose, |
---|
| 231 | use_cache=use_cache, |
---|
[7673] | 232 | output_centroids = output_centroids) |
---|
[7672] | 233 | |
---|
| 234 | if quake_offset_time is None: |
---|
| 235 | quake_offset_time = callable_sww.starttime |
---|
| 236 | |
---|
[7693] | 237 | for point_i, point in enumerate(points_array): |
---|
[7694] | 238 | is_opened = False |
---|
[7672] | 239 | for time in callable_sww.get_time(): |
---|
[7693] | 240 | #add domain starttime to relative time. |
---|
| 241 | quake_time = time + quake_offset_time |
---|
| 242 | point_quantities = callable_sww(time, point_i) # __call__ is overridden |
---|
[7675] | 243 | |
---|
[7693] | 244 | if point_quantities[0] != NAN: |
---|
| 245 | if is_opened == False: |
---|
| 246 | points_writer = writer(file(dir_name + sep + gauge_file |
---|
| 247 | + point_name[point_i] + '.csv', "wb")) |
---|
| 248 | points_writer.writerow(heading) |
---|
| 249 | is_opened = True |
---|
| 250 | points_list = [quake_time, quake_time/3600.] + _quantities2csv(quantities, point_quantities, callable_sww.centroids, point_i) |
---|
| 251 | points_writer.writerow(points_list) |
---|
| 252 | else: |
---|
[7694] | 253 | if verbose: |
---|
| 254 | msg = 'gauge' + point_name[point_i] + 'falls off the mesh in file ' + sww_file + '.' |
---|
| 255 | log.warning(msg) |
---|
[7672] | 256 | ## |
---|
| 257 | # @brief Read a .sww file and plot the time series. |
---|
| 258 | # @param swwfiles Dictionary of .sww files. |
---|
| 259 | # @param gauge_filename Name of gauge file. |
---|
| 260 | # @param production_dirs ?? |
---|
| 261 | # @param report If True, write figures to report directory. |
---|
| 262 | # @param reportname Name of generated report (if any). |
---|
| 263 | # @param plot_quantity List containing quantities to plot. |
---|
| 264 | # @param generate_fig If True, generate figures as well as CSV files. |
---|
| 265 | # @param surface If True, then generate solution surface with 3d plot. |
---|
| 266 | # @param time_min Beginning of user defined time range for plotting purposes. |
---|
| 267 | # @param time_max End of user defined time range for plotting purposes. |
---|
| 268 | # @param time_thinning ?? |
---|
| 269 | # @param time_unit ?? |
---|
| 270 | # @param title_on If True, export standard graphics with title. |
---|
| 271 | # @param use_cache If True, use caching. |
---|
| 272 | # @param verbose If True, this function is verbose. |
---|
[7679] | 273 | def sww2timeseries(swwfiles, |
---|
[7672] | 274 | gauge_filename, |
---|
| 275 | production_dirs, |
---|
| 276 | report=None, |
---|
| 277 | reportname=None, |
---|
| 278 | plot_quantity=None, |
---|
| 279 | generate_fig=False, |
---|
| 280 | surface=None, |
---|
| 281 | time_min=None, |
---|
| 282 | time_max=None, |
---|
| 283 | time_thinning=1, |
---|
| 284 | time_unit=None, |
---|
| 285 | title_on=None, |
---|
| 286 | use_cache=False, |
---|
[7673] | 287 | verbose=False, |
---|
[7685] | 288 | output_centroids=False): |
---|
[7672] | 289 | """ Read sww file and plot the time series for the |
---|
| 290 | prescribed quantities at defined gauge locations and |
---|
| 291 | prescribed time range. |
---|
| 292 | |
---|
| 293 | Input variables: |
---|
| 294 | |
---|
| 295 | swwfiles - dictionary of sww files with label_ids (used in |
---|
| 296 | generating latex output. It will be part of |
---|
| 297 | the directory name of file_loc (typically the timestamp). |
---|
| 298 | Helps to differentiate latex files for different |
---|
| 299 | simulations for a particular scenario. |
---|
| 300 | - assume that all conserved quantities have been stored |
---|
| 301 | - assume each sww file has been simulated with same timestep |
---|
| 302 | |
---|
| 303 | gauge_filename - name of file containing gauge data |
---|
| 304 | - easting, northing, name , elevation? |
---|
| 305 | - OR (this is not yet done) |
---|
| 306 | - structure which can be converted to a numeric array, |
---|
| 307 | such as a geospatial data object |
---|
| 308 | |
---|
| 309 | production_dirs - A list of list, example {20061101_121212: '1 in 10000', |
---|
| 310 | 'boundaries': 'urs boundary'} |
---|
| 311 | this will use the second part as the label and the |
---|
| 312 | first part as the ? |
---|
| 313 | #FIXME: Is it a list or a dictionary |
---|
| 314 | # This is probably obsolete by now |
---|
| 315 | |
---|
| 316 | report - if True, then write figures to report_figures directory in |
---|
| 317 | relevant production directory |
---|
| 318 | - if False, figures are already stored with sww file |
---|
| 319 | - default to False |
---|
| 320 | |
---|
| 321 | reportname - name for report if wishing to generate report |
---|
| 322 | |
---|
| 323 | plot_quantity - list containing quantities to plot, they must |
---|
| 324 | be the name of an existing quantity or one of |
---|
| 325 | the following possibilities |
---|
| 326 | - possibilities: |
---|
| 327 | - stage; 'stage' |
---|
| 328 | - depth; 'depth' |
---|
| 329 | - speed; calculated as absolute momentum |
---|
| 330 | (pointwise) divided by depth; 'speed' |
---|
| 331 | - bearing; calculated as the angle of the momentum |
---|
| 332 | vector (xmomentum, ymomentum) from the North; 'bearing' |
---|
| 333 | - absolute momentum; calculated as |
---|
| 334 | sqrt(xmomentum^2 + ymomentum^2); 'momentum' |
---|
| 335 | - x momentum; 'xmomentum' |
---|
| 336 | - y momentum; 'ymomentum' |
---|
| 337 | - default will be ['stage', 'speed', 'bearing'] |
---|
| 338 | |
---|
| 339 | generate_fig - if True, generate figures as well as csv file |
---|
| 340 | - if False, csv files created only |
---|
| 341 | |
---|
| 342 | surface - if True, then generate solution surface with 3d plot |
---|
| 343 | and save to current working directory |
---|
| 344 | - default = False |
---|
| 345 | |
---|
| 346 | time_min - beginning of user defined time range for plotting purposes |
---|
| 347 | - default will be first available time found in swwfile |
---|
| 348 | |
---|
| 349 | time_max - end of user defined time range for plotting purposes |
---|
| 350 | - default will be last available time found in swwfile |
---|
| 351 | |
---|
| 352 | title_on - if True, export standard graphics with title |
---|
| 353 | - if False, export standard graphics without title |
---|
| 354 | |
---|
| 355 | |
---|
| 356 | Output: |
---|
| 357 | |
---|
| 358 | - time series data stored in .csv for later use if required. |
---|
| 359 | Name = gauges_timeseries followed by gauge name |
---|
| 360 | - latex file will be generated in same directory as where script is |
---|
| 361 | run (usually production scenario directory. |
---|
| 362 | Name = latexoutputlabel_id.tex |
---|
| 363 | |
---|
| 364 | Other important information: |
---|
| 365 | |
---|
| 366 | It is assumed that the used has stored all the conserved quantities |
---|
| 367 | and elevation during the scenario run, i.e. |
---|
| 368 | ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
| 369 | If this has not occurred then sww2timeseries will not work. |
---|
| 370 | |
---|
| 371 | |
---|
| 372 | Usage example |
---|
| 373 | texname = sww2timeseries({project.boundary_name + '.sww': ''}, |
---|
| 374 | project.polygons_dir + sep + 'boundary_extent.csv', |
---|
| 375 | project.anuga_dir, |
---|
| 376 | report = False, |
---|
| 377 | plot_quantity = ['stage', 'speed', 'bearing'], |
---|
| 378 | time_min = None, |
---|
| 379 | time_max = None, |
---|
| 380 | title_on = True, |
---|
| 381 | verbose = True) |
---|
| 382 | |
---|
| 383 | """ |
---|
| 384 | |
---|
| 385 | msg = 'NOTE: A new function is available to create csv files from sww ' |
---|
| 386 | msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util' |
---|
| 387 | msg += ' PLUS another new function to create graphs from csv files called ' |
---|
| 388 | msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util' |
---|
| 389 | log.critical(msg) |
---|
| 390 | |
---|
| 391 | k = _sww2timeseries(swwfiles, |
---|
| 392 | gauge_filename, |
---|
| 393 | production_dirs, |
---|
| 394 | report, |
---|
| 395 | reportname, |
---|
| 396 | plot_quantity, |
---|
| 397 | generate_fig, |
---|
| 398 | surface, |
---|
| 399 | time_min, |
---|
| 400 | time_max, |
---|
| 401 | time_thinning, |
---|
| 402 | time_unit, |
---|
| 403 | title_on, |
---|
| 404 | use_cache, |
---|
[7673] | 405 | verbose, |
---|
[7685] | 406 | output_centroids = output_centroids) |
---|
[7672] | 407 | return k |
---|
| 408 | |
---|
| 409 | |
---|
| 410 | ## |
---|
| 411 | # @brief Read a .sww file and plot the time series. |
---|
| 412 | # @param swwfiles Dictionary of .sww files. |
---|
| 413 | # @param gauge_filename Name of gauge file. |
---|
| 414 | # @param production_dirs ?? |
---|
| 415 | # @param report If True, write figures to report directory. |
---|
| 416 | # @param reportname Name of generated report (if any). |
---|
| 417 | # @param plot_quantity List containing quantities to plot. |
---|
| 418 | # @param generate_fig If True, generate figures as well as CSV files. |
---|
| 419 | # @param surface If True, then generate solution surface with 3d plot. |
---|
| 420 | # @param time_min Beginning of user defined time range for plotting purposes. |
---|
| 421 | # @param time_max End of user defined time range for plotting purposes. |
---|
| 422 | # @param time_thinning ?? |
---|
| 423 | # @param time_unit ?? |
---|
| 424 | # @param title_on If True, export standard graphics with title. |
---|
| 425 | # @param use_cache If True, use caching. |
---|
| 426 | # @param verbose If True, this function is verbose. |
---|
| 427 | def _sww2timeseries(swwfiles, |
---|
| 428 | gauge_filename, |
---|
| 429 | production_dirs, |
---|
| 430 | report = None, |
---|
| 431 | reportname = None, |
---|
| 432 | plot_quantity = None, |
---|
| 433 | generate_fig = False, |
---|
| 434 | surface = None, |
---|
| 435 | time_min = None, |
---|
| 436 | time_max = None, |
---|
| 437 | time_thinning = 1, |
---|
| 438 | time_unit = None, |
---|
| 439 | title_on = None, |
---|
| 440 | use_cache = False, |
---|
[7673] | 441 | verbose = False, |
---|
[7685] | 442 | output_centroids = False): |
---|
[7672] | 443 | |
---|
| 444 | # FIXME(Ole): Shouldn't print statements here be governed by verbose? |
---|
| 445 | assert type(gauge_filename) == type(''), 'Gauge filename must be a string' |
---|
| 446 | |
---|
| 447 | try: |
---|
| 448 | fid = open(gauge_filename) |
---|
| 449 | except Exception, e: |
---|
| 450 | msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e) |
---|
| 451 | raise msg |
---|
| 452 | |
---|
| 453 | if report is None: |
---|
| 454 | report = False |
---|
| 455 | |
---|
| 456 | if plot_quantity is None: |
---|
| 457 | plot_quantity = ['depth', 'speed'] |
---|
| 458 | else: |
---|
| 459 | assert type(plot_quantity) == list, 'plot_quantity must be a list' |
---|
| 460 | check_list(plot_quantity) |
---|
| 461 | |
---|
| 462 | if surface is None: |
---|
| 463 | surface = False |
---|
| 464 | |
---|
| 465 | if time_unit is None: |
---|
| 466 | time_unit = 'hours' |
---|
| 467 | |
---|
| 468 | if title_on is None: |
---|
| 469 | title_on = True |
---|
| 470 | |
---|
| 471 | if verbose: log.critical('Gauges obtained from: %s' % gauge_filename) |
---|
| 472 | |
---|
[7679] | 473 | gauges, locations, elev = gauge_get_from_file(gauge_filename) |
---|
[7672] | 474 | |
---|
| 475 | sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
| 476 | |
---|
| 477 | file_loc = [] |
---|
| 478 | f_list = [] |
---|
| 479 | label_id = [] |
---|
| 480 | leg_label = [] |
---|
| 481 | themaxT = 0.0 |
---|
| 482 | theminT = 0.0 |
---|
| 483 | |
---|
| 484 | for swwfile in swwfiles.keys(): |
---|
| 485 | try: |
---|
| 486 | fid = open(swwfile) |
---|
| 487 | except Exception, e: |
---|
| 488 | msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e) |
---|
| 489 | raise msg |
---|
| 490 | |
---|
| 491 | if verbose: |
---|
| 492 | log.critical('swwfile = %s' % swwfile) |
---|
| 493 | |
---|
| 494 | # Extract parent dir name and use as label |
---|
| 495 | path, _ = os.path.split(swwfile) |
---|
| 496 | _, label = os.path.split(path) |
---|
| 497 | |
---|
| 498 | leg_label.append(label) |
---|
| 499 | |
---|
| 500 | f = file_function(swwfile, |
---|
| 501 | quantities = sww_quantity, |
---|
| 502 | interpolation_points = gauges, |
---|
| 503 | time_thinning = time_thinning, |
---|
| 504 | verbose = verbose, |
---|
[7673] | 505 | use_cache = use_cache, |
---|
[7685] | 506 | output_centroids = output_centroids) |
---|
[7672] | 507 | |
---|
| 508 | # determine which gauges are contained in sww file |
---|
| 509 | count = 0 |
---|
| 510 | gauge_index = [] |
---|
| 511 | for k, g in enumerate(gauges): |
---|
| 512 | if f(0.0, point_id = k)[2] > 1.0e6: |
---|
| 513 | count += 1 |
---|
| 514 | if count == 1: log.critical('Gauges not contained here:') |
---|
| 515 | log.critical(locations[k]) |
---|
| 516 | else: |
---|
| 517 | gauge_index.append(k) |
---|
| 518 | |
---|
| 519 | if len(gauge_index) > 0: |
---|
| 520 | log.critical('Gauges contained here:') |
---|
| 521 | else: |
---|
| 522 | log.critical('No gauges contained here.') |
---|
| 523 | for i in range(len(gauge_index)): |
---|
| 524 | log.critical(locations[gauge_index[i]]) |
---|
| 525 | |
---|
| 526 | index = swwfile.rfind(sep) |
---|
| 527 | file_loc.append(swwfile[:index+1]) |
---|
| 528 | label_id.append(swwfiles[swwfile]) |
---|
| 529 | |
---|
| 530 | f_list.append(f) |
---|
| 531 | maxT = max(f.get_time()) |
---|
| 532 | minT = min(f.get_time()) |
---|
| 533 | if maxT > themaxT: themaxT = maxT |
---|
| 534 | if minT > theminT: theminT = minT |
---|
| 535 | |
---|
| 536 | if time_min is None: |
---|
| 537 | time_min = theminT # min(T) |
---|
| 538 | else: |
---|
| 539 | if time_min < theminT: # min(T): |
---|
| 540 | msg = 'Minimum time entered not correct - please try again' |
---|
| 541 | raise Exception, msg |
---|
| 542 | |
---|
| 543 | if time_max is None: |
---|
| 544 | time_max = themaxT # max(T) |
---|
| 545 | else: |
---|
| 546 | if time_max > themaxT: # max(T): |
---|
| 547 | msg = 'Maximum time entered not correct - please try again' |
---|
| 548 | raise Exception, msg |
---|
| 549 | |
---|
| 550 | if verbose and len(gauge_index) > 0: |
---|
| 551 | log.critical('Inputs OK - going to generate figures') |
---|
| 552 | |
---|
| 553 | if len(gauge_index) <> 0: |
---|
| 554 | texfile, elev_output = \ |
---|
[7690] | 555 | _generate_figures(plot_quantity, file_loc, report, reportname, |
---|
[7672] | 556 | surface, leg_label, f_list, gauges, locations, |
---|
| 557 | elev, gauge_index, production_dirs, time_min, |
---|
| 558 | time_max, time_unit, title_on, label_id, |
---|
| 559 | generate_fig, verbose) |
---|
| 560 | else: |
---|
| 561 | texfile = '' |
---|
| 562 | elev_output = [] |
---|
| 563 | |
---|
| 564 | return texfile, elev_output |
---|
| 565 | |
---|
| 566 | |
---|
| 567 | ## |
---|
| 568 | # @brief Read gauge info from a file. |
---|
| 569 | # @param filename The name of the file to read. |
---|
| 570 | # @return A (gauges, gaugelocation, elev) tuple. |
---|
[7679] | 571 | def gauge_get_from_file(filename): |
---|
[7672] | 572 | """ Read in gauge information from file |
---|
| 573 | """ |
---|
| 574 | |
---|
| 575 | from os import sep, getcwd, access, F_OK, mkdir |
---|
| 576 | |
---|
| 577 | # Get data from the gauge file |
---|
| 578 | fid = open(filename) |
---|
| 579 | lines = fid.readlines() |
---|
| 580 | fid.close() |
---|
| 581 | |
---|
| 582 | gauges = [] |
---|
| 583 | gaugelocation = [] |
---|
| 584 | elev = [] |
---|
| 585 | |
---|
| 586 | # Check header information |
---|
| 587 | line1 = lines[0] |
---|
| 588 | line11 = line1.split(',') |
---|
| 589 | |
---|
| 590 | if isinstance(line11[0], str) is True: |
---|
| 591 | # We have found text in the first line |
---|
| 592 | east_index = None |
---|
| 593 | north_index = None |
---|
| 594 | name_index = None |
---|
| 595 | elev_index = None |
---|
| 596 | |
---|
| 597 | for i in range(len(line11)): |
---|
| 598 | if line11[i].strip().lower() == 'easting': east_index = i |
---|
| 599 | if line11[i].strip().lower() == 'northing': north_index = i |
---|
| 600 | if line11[i].strip().lower() == 'name': name_index = i |
---|
| 601 | if line11[i].strip().lower() == 'elevation': elev_index = i |
---|
| 602 | |
---|
| 603 | if east_index < len(line11) and north_index < len(line11): |
---|
| 604 | pass |
---|
| 605 | else: |
---|
| 606 | msg = 'WARNING: %s does not contain correct header information' \ |
---|
| 607 | % filename |
---|
| 608 | msg += 'The header must be: easting, northing, name, elevation' |
---|
| 609 | raise Exception, msg |
---|
| 610 | |
---|
| 611 | if elev_index is None: |
---|
| 612 | raise Exception |
---|
| 613 | |
---|
| 614 | if name_index is None: |
---|
| 615 | raise Exception |
---|
| 616 | |
---|
| 617 | lines = lines[1:] # Remove header from data |
---|
| 618 | else: |
---|
| 619 | # No header, assume that this is a simple easting, northing file |
---|
| 620 | |
---|
| 621 | msg = 'There was no header in file %s and the number of columns is %d' \ |
---|
| 622 | % (filename, len(line11)) |
---|
| 623 | msg += '- was assuming two columns corresponding to Easting and Northing' |
---|
| 624 | assert len(line11) == 2, msg |
---|
| 625 | |
---|
| 626 | east_index = 0 |
---|
| 627 | north_index = 1 |
---|
| 628 | |
---|
| 629 | N = len(lines) |
---|
| 630 | elev = [-9999]*N |
---|
| 631 | gaugelocation = range(N) |
---|
| 632 | |
---|
| 633 | # Read in gauge data |
---|
| 634 | for line in lines: |
---|
| 635 | fields = line.split(',') |
---|
| 636 | |
---|
| 637 | gauges.append([float(fields[east_index]), float(fields[north_index])]) |
---|
| 638 | |
---|
| 639 | if len(fields) > 2: |
---|
| 640 | elev.append(float(fields[elev_index])) |
---|
| 641 | loc = fields[name_index] |
---|
| 642 | gaugelocation.append(loc.strip('\n')) |
---|
| 643 | |
---|
| 644 | return gauges, gaugelocation, elev |
---|
[7690] | 645 | |
---|
| 646 | |
---|
| 647 | ## |
---|
| 648 | # @brief Generate figures from quantities and gauges for each sww file. |
---|
| 649 | # This is a legacy function, used only by sww2timeseries |
---|
| 650 | # @param plot_quantity ?? |
---|
| 651 | # @param file_loc ?? |
---|
| 652 | # @param report ?? |
---|
| 653 | # @param reportname ?? |
---|
| 654 | # @param surface ?? |
---|
| 655 | # @param leg_label ?? |
---|
| 656 | # @param f_list ?? |
---|
| 657 | # @param gauges ?? |
---|
| 658 | # @param locations ?? |
---|
| 659 | # @param elev ?? |
---|
| 660 | # @param gauge_index ?? |
---|
| 661 | # @param production_dirs ?? |
---|
| 662 | # @param time_min ?? |
---|
| 663 | # @param time_max ?? |
---|
| 664 | # @param time_unit ?? |
---|
| 665 | # @param title_on ?? |
---|
| 666 | # @param label_id ?? |
---|
| 667 | # @param generate_fig ?? |
---|
| 668 | # @param verbose?? |
---|
| 669 | # @return (texfile2, elev_output) |
---|
| 670 | def _generate_figures(plot_quantity, file_loc, report, reportname, surface, |
---|
| 671 | leg_label, f_list, gauges, locations, elev, gauge_index, |
---|
| 672 | production_dirs, time_min, time_max, time_unit, |
---|
| 673 | title_on, label_id, generate_fig, verbose): |
---|
| 674 | """ Generate figures based on required quantities and gauges for |
---|
| 675 | each sww file |
---|
| 676 | """ |
---|
| 677 | from os import sep, altsep, getcwd, mkdir, access, F_OK, environ |
---|
| 678 | |
---|
| 679 | if generate_fig is True: |
---|
| 680 | from pylab import ion, hold, plot, axis, figure, legend, savefig, \ |
---|
| 681 | xlabel, ylabel, title, close, subplot |
---|
| 682 | |
---|
| 683 | if surface is True: |
---|
| 684 | import pylab as p1 |
---|
| 685 | import mpl3d.mplot3d as p3 |
---|
| 686 | |
---|
| 687 | if report == True: |
---|
| 688 | texdir = getcwd()+sep+'report'+sep |
---|
| 689 | if access(texdir,F_OK) == 0: |
---|
| 690 | mkdir (texdir) |
---|
| 691 | if len(label_id) == 1: |
---|
| 692 | label_id1 = label_id[0].replace(sep,'') |
---|
| 693 | label_id2 = label_id1.replace('_','') |
---|
| 694 | texfile = texdir + reportname + '%s' % label_id2 |
---|
| 695 | texfile2 = reportname + '%s' % label_id2 |
---|
| 696 | texfilename = texfile + '.tex' |
---|
| 697 | fid = open(texfilename, 'w') |
---|
| 698 | |
---|
| 699 | if verbose: log.critical('Latex output printed to %s' % texfilename) |
---|
| 700 | else: |
---|
| 701 | texfile = texdir+reportname |
---|
| 702 | texfile2 = reportname |
---|
| 703 | texfilename = texfile + '.tex' |
---|
| 704 | fid = open(texfilename, 'w') |
---|
| 705 | |
---|
| 706 | if verbose: log.critical('Latex output printed to %s' % texfilename) |
---|
| 707 | else: |
---|
| 708 | texfile = '' |
---|
| 709 | texfile2 = '' |
---|
| 710 | |
---|
| 711 | p = len(f_list) |
---|
| 712 | n = [] |
---|
| 713 | n0 = 0 |
---|
| 714 | for i in range(len(f_list)): |
---|
| 715 | n.append(len(f_list[i].get_time())) |
---|
| 716 | if n[i] > n0: n0 = n[i] |
---|
| 717 | n0 = int(n0) |
---|
| 718 | m = len(locations) |
---|
| 719 | model_time = num.zeros((n0, m, p), num.float) |
---|
| 720 | stages = num.zeros((n0, m, p), num.float) |
---|
| 721 | elevations = num.zeros((n0, m, p), num.float) |
---|
| 722 | momenta = num.zeros((n0, m, p), num.float) |
---|
| 723 | xmom = num.zeros((n0, m, p), num.float) |
---|
| 724 | ymom = num.zeros((n0, m, p), num.float) |
---|
| 725 | speed = num.zeros((n0, m, p), num.float) |
---|
| 726 | bearings = num.zeros((n0, m, p), num.float) |
---|
| 727 | due_east = 90.0*num.ones((n0, 1), num.float) |
---|
| 728 | due_west = 270.0*num.ones((n0, 1), num.float) |
---|
| 729 | depths = num.zeros((n0, m, p), num.float) |
---|
| 730 | eastings = num.zeros((n0, m, p), num.float) |
---|
| 731 | min_stages = [] |
---|
| 732 | max_stages = [] |
---|
| 733 | min_momentums = [] |
---|
| 734 | max_momentums = [] |
---|
| 735 | max_xmomentums = [] |
---|
| 736 | max_ymomentums = [] |
---|
| 737 | min_xmomentums = [] |
---|
| 738 | min_ymomentums = [] |
---|
| 739 | max_speeds = [] |
---|
| 740 | min_speeds = [] |
---|
| 741 | max_depths = [] |
---|
| 742 | model_time_plot3d = num.zeros((n0, m), num.float) |
---|
| 743 | stages_plot3d = num.zeros((n0, m), num.float) |
---|
| 744 | eastings_plot3d = num.zeros((n0, m),num.float) |
---|
| 745 | if time_unit is 'mins': scale = 60.0 |
---|
| 746 | if time_unit is 'hours': scale = 3600.0 |
---|
| 747 | |
---|
| 748 | ##### loop over each swwfile ##### |
---|
| 749 | for j, f in enumerate(f_list): |
---|
| 750 | if verbose: log.critical('swwfile %d of %d' % (j, len(f_list))) |
---|
| 751 | |
---|
| 752 | starttime = f.starttime |
---|
| 753 | comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv' |
---|
| 754 | fid_compare = open(comparefile, 'w') |
---|
| 755 | file0 = file_loc[j] + 'gauges_t0.csv' |
---|
| 756 | fid_0 = open(file0, 'w') |
---|
| 757 | |
---|
| 758 | ##### loop over each gauge ##### |
---|
| 759 | for k in gauge_index: |
---|
| 760 | if verbose: log.critical('Gauge %d of %d' % (k, len(gauges))) |
---|
| 761 | |
---|
| 762 | g = gauges[k] |
---|
| 763 | min_stage = 10 |
---|
| 764 | max_stage = 0 |
---|
| 765 | max_momentum = max_xmomentum = max_ymomentum = 0 |
---|
| 766 | min_momentum = min_xmomentum = min_ymomentum = 100 |
---|
| 767 | max_speed = 0 |
---|
| 768 | min_speed = 0 |
---|
| 769 | max_depth = 0 |
---|
| 770 | gaugeloc = str(locations[k]) |
---|
| 771 | thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \ |
---|
| 772 | + gaugeloc + '.csv' |
---|
| 773 | if j == 0: |
---|
| 774 | fid_out = open(thisfile, 'w') |
---|
| 775 | s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n' |
---|
| 776 | fid_out.write(s) |
---|
| 777 | |
---|
| 778 | #### generate quantities ####### |
---|
| 779 | for i, t in enumerate(f.get_time()): |
---|
| 780 | if time_min <= t <= time_max: |
---|
| 781 | w = f(t, point_id = k)[0] |
---|
| 782 | z = f(t, point_id = k)[1] |
---|
| 783 | uh = f(t, point_id = k)[2] |
---|
| 784 | vh = f(t, point_id = k)[3] |
---|
| 785 | depth = w-z |
---|
| 786 | m = sqrt(uh*uh + vh*vh) |
---|
| 787 | if depth < 0.001: |
---|
| 788 | vel = 0.0 |
---|
| 789 | else: |
---|
| 790 | vel = m / (depth + 1.e-6/depth) |
---|
| 791 | bearing = calc_bearing(uh, vh) |
---|
| 792 | model_time[i,k,j] = (t + starttime)/scale #t/60.0 |
---|
| 793 | stages[i,k,j] = w |
---|
| 794 | elevations[i,k,j] = z |
---|
| 795 | xmom[i,k,j] = uh |
---|
| 796 | ymom[i,k,j] = vh |
---|
| 797 | momenta[i,k,j] = m |
---|
| 798 | speed[i,k,j] = vel |
---|
| 799 | bearings[i,k,j] = bearing |
---|
| 800 | depths[i,k,j] = depth |
---|
| 801 | thisgauge = gauges[k] |
---|
| 802 | eastings[i,k,j] = thisgauge[0] |
---|
| 803 | s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \ |
---|
| 804 | % (t, w, m, vel, z, uh, vh, bearing) |
---|
| 805 | fid_out.write(s) |
---|
| 806 | if t == 0: |
---|
| 807 | s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w) |
---|
| 808 | fid_0.write(s) |
---|
| 809 | if t/60.0 <= 13920: tindex = i |
---|
| 810 | if w > max_stage: max_stage = w |
---|
| 811 | if w < min_stage: min_stage = w |
---|
| 812 | if m > max_momentum: max_momentum = m |
---|
| 813 | if m < min_momentum: min_momentum = m |
---|
| 814 | if uh > max_xmomentum: max_xmomentum = uh |
---|
| 815 | if vh > max_ymomentum: max_ymomentum = vh |
---|
| 816 | if uh < min_xmomentum: min_xmomentum = uh |
---|
| 817 | if vh < min_ymomentum: min_ymomentum = vh |
---|
| 818 | if vel > max_speed: max_speed = vel |
---|
| 819 | if vel < min_speed: min_speed = vel |
---|
| 820 | if z > 0 and depth > max_depth: max_depth = depth |
---|
| 821 | |
---|
| 822 | |
---|
| 823 | s = '%.2f, %.2f, %.2f, %.2f, %s\n' \ |
---|
| 824 | % (max_stage, min_stage, z, thisgauge[0], leg_label[j]) |
---|
| 825 | fid_compare.write(s) |
---|
| 826 | max_stages.append(max_stage) |
---|
| 827 | min_stages.append(min_stage) |
---|
| 828 | max_momentums.append(max_momentum) |
---|
| 829 | max_xmomentums.append(max_xmomentum) |
---|
| 830 | max_ymomentums.append(max_ymomentum) |
---|
| 831 | min_xmomentums.append(min_xmomentum) |
---|
| 832 | min_ymomentums.append(min_ymomentum) |
---|
| 833 | min_momentums.append(min_momentum) |
---|
| 834 | max_depths.append(max_depth) |
---|
| 835 | max_speeds.append(max_speed) |
---|
| 836 | min_speeds.append(min_speed) |
---|
| 837 | #### finished generating quantities for each swwfile ##### |
---|
| 838 | |
---|
| 839 | model_time_plot3d[:,:] = model_time[:,:,j] |
---|
| 840 | stages_plot3d[:,:] = stages[:,:,j] |
---|
| 841 | eastings_plot3d[:,] = eastings[:,:,j] |
---|
| 842 | |
---|
| 843 | if surface is True: |
---|
| 844 | log.critical('Printing surface figure') |
---|
| 845 | for i in range(2): |
---|
| 846 | fig = p1.figure(10) |
---|
| 847 | ax = p3.Axes3D(fig) |
---|
| 848 | if len(gauges) > 80: |
---|
| 849 | ax.plot_surface(model_time[:,:,j], |
---|
| 850 | eastings[:,:,j], |
---|
| 851 | stages[:,:,j]) |
---|
| 852 | else: |
---|
| 853 | ax.plot3D(num.ravel(eastings[:,:,j]), |
---|
| 854 | num.ravel(model_time[:,:,j]), |
---|
| 855 | num.ravel(stages[:,:,j])) |
---|
| 856 | ax.set_xlabel('time') |
---|
| 857 | ax.set_ylabel('x') |
---|
| 858 | ax.set_zlabel('stage') |
---|
| 859 | fig.add_axes(ax) |
---|
| 860 | p1.show() |
---|
| 861 | surfacefig = 'solution_surface%s' % leg_label[j] |
---|
| 862 | p1.savefig(surfacefig) |
---|
| 863 | p1.close() |
---|
| 864 | |
---|
| 865 | #### finished generating quantities for all swwfiles ##### |
---|
| 866 | |
---|
| 867 | # x profile for given time |
---|
| 868 | if surface is True: |
---|
| 869 | figure(11) |
---|
| 870 | plot(eastings[tindex,:,j], stages[tindex,:,j]) |
---|
| 871 | xlabel('x') |
---|
| 872 | ylabel('stage') |
---|
| 873 | profilefig = 'solution_xprofile' |
---|
| 874 | savefig('profilefig') |
---|
| 875 | |
---|
| 876 | elev_output = [] |
---|
| 877 | if generate_fig is True: |
---|
| 878 | depth_axis = axis([starttime/scale, time_max/scale, -0.1, |
---|
| 879 | max(max_depths)*1.1]) |
---|
| 880 | stage_axis = axis([starttime/scale, time_max/scale, |
---|
| 881 | min(min_stages), max(max_stages)*1.1]) |
---|
| 882 | vel_axis = axis([starttime/scale, time_max/scale, |
---|
| 883 | min(min_speeds), max(max_speeds)*1.1]) |
---|
| 884 | mom_axis = axis([starttime/scale, time_max/scale, |
---|
| 885 | min(min_momentums), max(max_momentums)*1.1]) |
---|
| 886 | xmom_axis = axis([starttime/scale, time_max/scale, |
---|
| 887 | min(min_xmomentums), max(max_xmomentums)*1.1]) |
---|
| 888 | ymom_axis = axis([starttime/scale, time_max/scale, |
---|
| 889 | min(min_ymomentums), max(max_ymomentums)*1.1]) |
---|
| 890 | cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k'] |
---|
| 891 | nn = len(plot_quantity) |
---|
| 892 | no_cols = 2 |
---|
| 893 | |
---|
| 894 | if len(label_id) > 1: graphname_report = [] |
---|
| 895 | pp = 1 |
---|
| 896 | div = 11. |
---|
| 897 | cc = 0 |
---|
| 898 | for k in gauge_index: |
---|
| 899 | g = gauges[k] |
---|
| 900 | count1 = 0 |
---|
| 901 | if report == True and len(label_id) > 1: |
---|
| 902 | s = '\\begin{figure}[ht] \n' \ |
---|
| 903 | '\\centering \n' \ |
---|
| 904 | '\\begin{tabular}{cc} \n' |
---|
| 905 | fid.write(s) |
---|
| 906 | if len(label_id) > 1: graphname_report = [] |
---|
| 907 | |
---|
| 908 | #### generate figures for each gauge #### |
---|
| 909 | for j, f in enumerate(f_list): |
---|
| 910 | ion() |
---|
| 911 | hold(True) |
---|
| 912 | count = 0 |
---|
| 913 | where1 = 0 |
---|
| 914 | where2 = 0 |
---|
| 915 | word_quantity = '' |
---|
| 916 | if report == True and len(label_id) == 1: |
---|
| 917 | s = '\\begin{figure}[hbt] \n' \ |
---|
| 918 | '\\centering \n' \ |
---|
| 919 | '\\begin{tabular}{cc} \n' |
---|
| 920 | fid.write(s) |
---|
| 921 | |
---|
| 922 | for which_quantity in plot_quantity: |
---|
| 923 | count += 1 |
---|
| 924 | where1 += 1 |
---|
| 925 | figure(count, frameon = False) |
---|
| 926 | if which_quantity == 'depth': |
---|
| 927 | plot(model_time[0:n[j]-1,k,j], |
---|
| 928 | depths[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 929 | units = 'm' |
---|
| 930 | axis(depth_axis) |
---|
| 931 | if which_quantity == 'stage': |
---|
| 932 | if elevations[0,k,j] <= 0: |
---|
| 933 | plot(model_time[0:n[j]-1,k,j], |
---|
| 934 | stages[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 935 | axis(stage_axis) |
---|
| 936 | else: |
---|
| 937 | plot(model_time[0:n[j]-1,k,j], |
---|
| 938 | depths[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 939 | #axis(depth_axis) |
---|
| 940 | units = 'm' |
---|
| 941 | if which_quantity == 'momentum': |
---|
| 942 | plot(model_time[0:n[j]-1,k,j], |
---|
| 943 | momenta[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 944 | axis(mom_axis) |
---|
| 945 | units = 'm^2 / sec' |
---|
| 946 | if which_quantity == 'xmomentum': |
---|
| 947 | plot(model_time[0:n[j]-1,k,j], |
---|
| 948 | xmom[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 949 | axis(xmom_axis) |
---|
| 950 | units = 'm^2 / sec' |
---|
| 951 | if which_quantity == 'ymomentum': |
---|
| 952 | plot(model_time[0:n[j]-1,k,j], |
---|
| 953 | ymom[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 954 | axis(ymom_axis) |
---|
| 955 | units = 'm^2 / sec' |
---|
| 956 | if which_quantity == 'speed': |
---|
| 957 | plot(model_time[0:n[j]-1,k,j], |
---|
| 958 | speed[0:n[j]-1,k,j], '-', c = cstr[j]) |
---|
| 959 | axis(vel_axis) |
---|
| 960 | units = 'm / sec' |
---|
| 961 | if which_quantity == 'bearing': |
---|
| 962 | plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-', |
---|
| 963 | model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', |
---|
| 964 | model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.') |
---|
| 965 | units = 'degrees from North' |
---|
| 966 | #ax = axis([time_min, time_max, 0.0, 360.0]) |
---|
| 967 | legend(('Bearing','West','East')) |
---|
| 968 | |
---|
| 969 | if time_unit is 'mins': xlabel('time (mins)') |
---|
| 970 | if time_unit is 'hours': xlabel('time (hours)') |
---|
| 971 | #if which_quantity == 'stage' \ |
---|
| 972 | # and elevations[0:n[j]-1,k,j] > 0: |
---|
| 973 | # ylabel('%s (%s)' %('depth', units)) |
---|
| 974 | #else: |
---|
| 975 | # ylabel('%s (%s)' %(which_quantity, units)) |
---|
| 976 | #ylabel('%s (%s)' %('wave height', units)) |
---|
| 977 | ylabel('%s (%s)' %(which_quantity, units)) |
---|
| 978 | if len(label_id) > 1: legend((leg_label),loc='upper right') |
---|
| 979 | |
---|
| 980 | #gaugeloc1 = gaugeloc.replace(' ','') |
---|
| 981 | #gaugeloc2 = gaugeloc1.replace('_','') |
---|
| 982 | gaugeloc2 = str(locations[k]).replace(' ','') |
---|
| 983 | graphname = '%sgauge%s_%s' %(file_loc[j], |
---|
| 984 | gaugeloc2, |
---|
| 985 | which_quantity) |
---|
| 986 | |
---|
| 987 | if report == True and len(label_id) > 1: |
---|
| 988 | figdir = getcwd()+sep+'report_figures'+sep |
---|
| 989 | if access(figdir,F_OK) == 0 : |
---|
| 990 | mkdir (figdir) |
---|
| 991 | latex_file_loc = figdir.replace(sep,altsep) |
---|
| 992 | # storing files in production directory |
---|
| 993 | graphname_latex = '%sgauge%s%s' \ |
---|
| 994 | % (latex_file_loc, gaugeloc2, |
---|
| 995 | which_quantity) |
---|
| 996 | # giving location in latex output file |
---|
| 997 | graphname_report_input = '%sgauge%s%s' % \ |
---|
| 998 | ('..' + altsep + |
---|
| 999 | 'report_figures' + altsep, |
---|
| 1000 | gaugeloc2, which_quantity) |
---|
| 1001 | graphname_report.append(graphname_report_input) |
---|
| 1002 | |
---|
| 1003 | # save figures in production directory for report |
---|
| 1004 | savefig(graphname_latex) |
---|
| 1005 | |
---|
| 1006 | if report == True: |
---|
| 1007 | figdir = getcwd() + sep + 'report_figures' + sep |
---|
| 1008 | if access(figdir,F_OK) == 0: |
---|
| 1009 | mkdir(figdir) |
---|
| 1010 | latex_file_loc = figdir.replace(sep,altsep) |
---|
| 1011 | |
---|
| 1012 | if len(label_id) == 1: |
---|
| 1013 | # storing files in production directory |
---|
| 1014 | graphname_latex = '%sgauge%s%s%s' % \ |
---|
| 1015 | (latex_file_loc, gaugeloc2, |
---|
| 1016 | which_quantity, label_id2) |
---|
| 1017 | # giving location in latex output file |
---|
| 1018 | graphname_report = '%sgauge%s%s%s' % \ |
---|
| 1019 | ('..' + altsep + |
---|
| 1020 | 'report_figures' + altsep, |
---|
| 1021 | gaugeloc2, which_quantity, |
---|
| 1022 | label_id2) |
---|
| 1023 | s = '\includegraphics' \ |
---|
| 1024 | '[width=0.49\linewidth, height=50mm]{%s%s}' % \ |
---|
| 1025 | (graphname_report, '.png') |
---|
| 1026 | fid.write(s) |
---|
| 1027 | if where1 % 2 == 0: |
---|
| 1028 | s = '\\\\ \n' |
---|
| 1029 | where1 = 0 |
---|
| 1030 | else: |
---|
| 1031 | s = '& \n' |
---|
| 1032 | fid.write(s) |
---|
| 1033 | savefig(graphname_latex) |
---|
| 1034 | |
---|
| 1035 | if title_on == True: |
---|
| 1036 | title('%s scenario: %s at %s gauge' % \ |
---|
| 1037 | (label_id, which_quantity, gaugeloc2)) |
---|
| 1038 | #title('Gauge %s (MOST elevation %.2f, ' \ |
---|
| 1039 | # 'ANUGA elevation %.2f)' % \ |
---|
| 1040 | # (gaugeloc2, elevations[10,k,0], |
---|
| 1041 | # elevations[10,k,1])) |
---|
| 1042 | |
---|
| 1043 | savefig(graphname) # save figures with sww file |
---|
| 1044 | |
---|
| 1045 | if report == True and len(label_id) == 1: |
---|
| 1046 | for i in range(nn-1): |
---|
| 1047 | if nn > 2: |
---|
| 1048 | if plot_quantity[i] == 'stage' \ |
---|
| 1049 | and elevations[0,k,j] > 0: |
---|
| 1050 | word_quantity += 'depth' + ', ' |
---|
| 1051 | else: |
---|
| 1052 | word_quantity += plot_quantity[i] + ', ' |
---|
| 1053 | else: |
---|
| 1054 | if plot_quantity[i] == 'stage' \ |
---|
| 1055 | and elevations[0,k,j] > 0: |
---|
| 1056 | word_quantity += 'depth' + ', ' |
---|
| 1057 | else: |
---|
| 1058 | word_quantity += plot_quantity[i] |
---|
| 1059 | |
---|
| 1060 | if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0: |
---|
| 1061 | word_quantity += ' and ' + 'depth' |
---|
| 1062 | else: |
---|
| 1063 | word_quantity += ' and ' + plot_quantity[nn-1] |
---|
| 1064 | caption = 'Time series for %s at %s location ' \ |
---|
| 1065 | '(elevation %.2fm)' % \ |
---|
| 1066 | (word_quantity, locations[k], elev[k]) |
---|
| 1067 | if elev[k] == 0.0: |
---|
| 1068 | caption = 'Time series for %s at %s location ' \ |
---|
| 1069 | '(elevation %.2fm)' % \ |
---|
| 1070 | (word_quantity, locations[k], |
---|
| 1071 | elevations[0,k,j]) |
---|
| 1072 | east = gauges[0] |
---|
| 1073 | north = gauges[1] |
---|
| 1074 | elev_output.append([locations[k], east, north, |
---|
| 1075 | elevations[0,k,j]]) |
---|
| 1076 | label = '%sgauge%s' % (label_id2, gaugeloc2) |
---|
| 1077 | s = '\end{tabular} \n' \ |
---|
| 1078 | '\\caption{%s} \n' \ |
---|
| 1079 | '\label{fig:%s} \n' \ |
---|
| 1080 | '\end{figure} \n \n' % (caption, label) |
---|
| 1081 | fid.write(s) |
---|
| 1082 | cc += 1 |
---|
| 1083 | if cc % 6 == 0: fid.write('\\clearpage \n') |
---|
| 1084 | savefig(graphname_latex) |
---|
| 1085 | |
---|
| 1086 | if report == True and len(label_id) > 1: |
---|
| 1087 | for i in range(nn-1): |
---|
| 1088 | if nn > 2: |
---|
| 1089 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
| 1090 | word_quantity += 'depth' + ',' |
---|
| 1091 | else: |
---|
| 1092 | word_quantity += plot_quantity[i] + ', ' |
---|
| 1093 | else: |
---|
| 1094 | if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: |
---|
| 1095 | word_quantity += 'depth' |
---|
| 1096 | else: |
---|
| 1097 | word_quantity += plot_quantity[i] |
---|
| 1098 | where1 = 0 |
---|
| 1099 | count1 += 1 |
---|
| 1100 | index = j*len(plot_quantity) |
---|
| 1101 | for which_quantity in plot_quantity: |
---|
| 1102 | where1 += 1 |
---|
| 1103 | s = '\includegraphics' \ |
---|
| 1104 | '[width=0.49\linewidth, height=50mm]{%s%s}' % \ |
---|
| 1105 | (graphname_report[index], '.png') |
---|
| 1106 | index += 1 |
---|
| 1107 | fid.write(s) |
---|
| 1108 | if where1 % 2 == 0: |
---|
| 1109 | s = '\\\\ \n' |
---|
| 1110 | where1 = 0 |
---|
| 1111 | else: |
---|
| 1112 | s = '& \n' |
---|
| 1113 | fid.write(s) |
---|
| 1114 | word_quantity += ' and ' + plot_quantity[nn-1] |
---|
| 1115 | label = 'gauge%s' %(gaugeloc2) |
---|
| 1116 | caption = 'Time series for %s at %s location ' \ |
---|
| 1117 | '(elevation %.2fm)' % \ |
---|
| 1118 | (word_quantity, locations[k], elev[k]) |
---|
| 1119 | if elev[k] == 0.0: |
---|
| 1120 | caption = 'Time series for %s at %s location ' \ |
---|
| 1121 | '(elevation %.2fm)' % \ |
---|
| 1122 | (word_quantity, locations[k], |
---|
| 1123 | elevations[0,k,j]) |
---|
| 1124 | thisgauge = gauges[k] |
---|
| 1125 | east = thisgauge[0] |
---|
| 1126 | north = thisgauge[1] |
---|
| 1127 | elev_output.append([locations[k], east, north, |
---|
| 1128 | elevations[0,k,j]]) |
---|
| 1129 | |
---|
| 1130 | s = '\end{tabular} \n' \ |
---|
| 1131 | '\\caption{%s} \n' \ |
---|
| 1132 | '\label{fig:%s} \n' \ |
---|
| 1133 | '\end{figure} \n \n' % (caption, label) |
---|
| 1134 | fid.write(s) |
---|
| 1135 | if float((k+1)/div - pp) == 0.: |
---|
| 1136 | fid.write('\\clearpage \n') |
---|
| 1137 | pp += 1 |
---|
| 1138 | #### finished generating figures ### |
---|
| 1139 | |
---|
| 1140 | close('all') |
---|
| 1141 | |
---|
| 1142 | return texfile2, elev_output |
---|