"""Gauge functions High-level functions for converting gauge and sww files into timeseries plots. Copyright 2010 Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou, James Hudson Geoscience Australia """ import numpy as num from anuga.geospatial_data.geospatial_data import ensure_absolute from util import check_list, calc_bearing from file_function import file_function import os from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep, getcwd from os.path import exists, split, join import anuga.utilities.log as log from math import sqrt def _quantities2csv(quantities, point_quantities, centroids, point_i): points_list = [] for quantity in quantities: #core quantities that are exported from the interpolator if quantity == 'stage': points_list.append(point_quantities[0]) if quantity == 'elevation': points_list.append(point_quantities[1]) if quantity == 'xmomentum': points_list.append(point_quantities[2]) if quantity == 'ymomentum': points_list.append(point_quantities[3]) #derived quantities that are calculated from the core ones if quantity == 'depth': points_list.append(point_quantities[0] - point_quantities[1]) if quantity == 'momentum': momentum = sqrt(point_quantities[2]**2 + point_quantities[3]**2) points_list.append(momentum) if quantity == 'speed': #if depth is less than 0.001 then speed = 0.0 if point_quantities[0] - point_quantities[1] < 0.001: vel = 0.0 else: if point_quantities[2] < 1.0e6: momentum = sqrt(point_quantities[2]**2 + point_quantities[3]**2) vel = momentum / (point_quantities[0] - point_quantities[1]) else: momentum = 0 vel = 0 points_list.append(vel) if quantity == 'bearing': points_list.append(calc_bearing(point_quantities[2], point_quantities[3])) if quantity == 'xcentroid': points_list.append(centroids[point_i][0]) if quantity == 'ycentroid': points_list.append(centroids[point_i][1]) return points_list ## # @brief Take gauge readings, given a gauge file and a sww mesh # # Use this function to take a timeseries sample, given a list of gauge points # @param sww_file sww file to use as input # @param gauge_file gauge file as input, containing list of gauge points to sample # @param out_name output file prefix # @param quantities which quantities in the sww file we want to export # @param verbose show extra logging information for debug purposes # @param use_cache cache requests if possible, for speed # @param output_centroids Set to true to output the values at the centroid of the mesh triangle def sww2csv_gauges(sww_file, gauge_file, out_name='gauge_', quantities=['stage', 'depth', 'elevation', 'xmomentum', 'ymomentum'], verbose=False, use_cache=True, output_centroids=False): """ Inputs: NOTE: if using csv2timeseries_graphs after creating csv file, it is essential to export quantities 'depth' and 'elevation'. 'depth' is good to analyse gauges on land and elevation is used automatically by csv2timeseries_graphs in the legend. sww_file: path to any sww file gauge_file: Assumes that it follows this format name, easting, northing, elevation point1, 100.3, 50.2, 10.0 point2, 10.3, 70.3, 78.0 NOTE: order of column can change but names eg 'easting', 'elevation' must be the same! ALL lowercaps! out_name: prefix for output file name (default is 'gauge_') Outputs: one file for each gauge/point location in the points file. They will be named with this format in the same directory as the 'sww_file' .csv eg gauge_point1.csv if not supplied myfile_2_point1.csv if ='myfile_2_' They will all have a header Usage: sww2csv_gauges(sww_file='test1.sww', quantities = ['stage', 'elevation','depth','bearing'], gauge_file='gauge.txt') Interpolate the quantities at a given set of locations, given an sww file. The results are written to a csv file. In the future let points be a points file. And the user choose the quantities. This is currently quite specific. If it needs to be more general, change things. This is really returning speed, not velocity. """ from csv import reader,writer from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN import string from anuga.utilities.file_utils import get_all_swwfiles from anuga.abstract_2d_finite_volumes.util import file_function assert type(gauge_file) == type(''), 'Gauge filename must be a string' assert type(out_name) == type(''), 'Output filename prefix must be a string' try: point_reader = reader(file(gauge_file)) except Exception, e: msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e) raise msg if verbose: log.critical('Gauges obtained from: %s' % gauge_file) point_reader = reader(file(gauge_file)) points = [] point_name = [] # read point info from file for i,row in enumerate(point_reader): # read header and determine the column numbers to read correctly. if i==0: for j,value in enumerate(row): if value.strip()=='easting':easting=j if value.strip()=='northing':northing=j if value.strip()=='name':name=j if value.strip()=='elevation':elevation=j else: points.append([float(row[easting]),float(row[northing])]) point_name.append(row[name]) #convert to array for file_function points_array = num.array(points,num.float) points_array = ensure_absolute(points_array) dir_name, base = os.path.split(sww_file) #need to get current directory so when path and file #are "joined" below the directory is correct if dir_name == '': dir_name =getcwd() if access(sww_file,R_OK): if verbose: log.critical('File %s exists' % sww_file) else: msg = 'File "%s" could not be opened: no read permission' % sww_file raise msg sww_files = get_all_swwfiles(look_in_dir=dir_name, base_name=base, verbose=verbose) # fudge to get SWW files in 'correct' order, oldest on the left sww_files.sort() if verbose: log.critical('sww files=%s' % sww_files) #to make all the quantities lower case for file_function quantities = [quantity.lower() for quantity in quantities] # what is quantities are needed from sww file to calculate output quantities # also core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] gauge_file = out_name heading = [quantity for quantity in quantities] heading.insert(0,'time') heading.insert(1,'hours') if verbose: log.critical('Writing csv files') quake_offset_time = None for sww_file in sww_files: sww_file = join(dir_name, sww_file+'.sww') callable_sww = file_function(sww_file, quantities=core_quantities, interpolation_points=points_array, verbose=verbose, use_cache=use_cache, output_centroids = output_centroids) if quake_offset_time is None: quake_offset_time = callable_sww.starttime for point_i, point in enumerate(points_array): is_opened = False for time in callable_sww.get_time(): #add domain starttime to relative time. quake_time = time + quake_offset_time point_quantities = callable_sww(time, point_i) # __call__ is overridden if point_quantities[0] != NAN: if is_opened == False: points_writer = writer(file(dir_name + sep + gauge_file + point_name[point_i] + '.csv', "wb")) points_writer.writerow(heading) is_opened = True points_list = [quake_time, quake_time/3600.] + _quantities2csv(quantities, point_quantities, callable_sww.centroids, point_i) points_writer.writerow(points_list) else: if verbose: msg = 'gauge' + point_name[point_i] + 'falls off the mesh in file ' + sww_file + '.' log.warning(msg) ## # @brief Read a .sww file and plot the time series. # @param swwfiles Dictionary of .sww files. # @param gauge_filename Name of gauge file. # @param production_dirs ?? # @param report If True, write figures to report directory. # @param reportname Name of generated report (if any). # @param plot_quantity List containing quantities to plot. # @param generate_fig If True, generate figures as well as CSV files. # @param surface If True, then generate solution surface with 3d plot. # @param time_min Beginning of user defined time range for plotting purposes. # @param time_max End of user defined time range for plotting purposes. # @param time_thinning ?? # @param time_unit ?? # @param title_on If True, export standard graphics with title. # @param use_cache If True, use caching. # @param verbose If True, this function is verbose. def sww2timeseries(swwfiles, gauge_filename, production_dirs, report=None, reportname=None, plot_quantity=None, generate_fig=False, surface=None, time_min=None, time_max=None, time_thinning=1, time_unit=None, title_on=None, use_cache=False, verbose=False, output_centroids=False): """ Read sww file and plot the time series for the prescribed quantities at defined gauge locations and prescribed time range. Input variables: swwfiles - dictionary of sww files with label_ids (used in generating latex output. It will be part of the directory name of file_loc (typically the timestamp). Helps to differentiate latex files for different simulations for a particular scenario. - assume that all conserved quantities have been stored - assume each sww file has been simulated with same timestep gauge_filename - name of file containing gauge data - easting, northing, name , elevation? - OR (this is not yet done) - structure which can be converted to a numeric array, such as a geospatial data object production_dirs - A list of list, example {20061101_121212: '1 in 10000', 'boundaries': 'urs boundary'} this will use the second part as the label and the first part as the ? #FIXME: Is it a list or a dictionary # This is probably obsolete by now report - if True, then write figures to report_figures directory in relevant production directory - if False, figures are already stored with sww file - default to False reportname - name for report if wishing to generate report plot_quantity - list containing quantities to plot, they must be the name of an existing quantity or one of the following possibilities - possibilities: - stage; 'stage' - depth; 'depth' - speed; calculated as absolute momentum (pointwise) divided by depth; 'speed' - bearing; calculated as the angle of the momentum vector (xmomentum, ymomentum) from the North; 'bearing' - absolute momentum; calculated as sqrt(xmomentum^2 + ymomentum^2); 'momentum' - x momentum; 'xmomentum' - y momentum; 'ymomentum' - default will be ['stage', 'speed', 'bearing'] generate_fig - if True, generate figures as well as csv file - if False, csv files created only surface - if True, then generate solution surface with 3d plot and save to current working directory - default = False time_min - beginning of user defined time range for plotting purposes - default will be first available time found in swwfile time_max - end of user defined time range for plotting purposes - default will be last available time found in swwfile title_on - if True, export standard graphics with title - if False, export standard graphics without title Output: - time series data stored in .csv for later use if required. Name = gauges_timeseries followed by gauge name - latex file will be generated in same directory as where script is run (usually production scenario directory. Name = latexoutputlabel_id.tex Other important information: It is assumed that the used has stored all the conserved quantities and elevation during the scenario run, i.e. ['stage', 'elevation', 'xmomentum', 'ymomentum'] If this has not occurred then sww2timeseries will not work. Usage example texname = sww2timeseries({project.boundary_name + '.sww': ''}, project.polygons_dir + sep + 'boundary_extent.csv', project.anuga_dir, report = False, plot_quantity = ['stage', 'speed', 'bearing'], time_min = None, time_max = None, title_on = True, verbose = True) """ msg = 'NOTE: A new function is available to create csv files from sww ' msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util' msg += ' PLUS another new function to create graphs from csv files called ' msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util' log.critical(msg) k = _sww2timeseries(swwfiles, gauge_filename, production_dirs, report, reportname, plot_quantity, generate_fig, surface, time_min, time_max, time_thinning, time_unit, title_on, use_cache, verbose, output_centroids = output_centroids) return k ## # @brief Read a .sww file and plot the time series. # @param swwfiles Dictionary of .sww files. # @param gauge_filename Name of gauge file. # @param production_dirs ?? # @param report If True, write figures to report directory. # @param reportname Name of generated report (if any). # @param plot_quantity List containing quantities to plot. # @param generate_fig If True, generate figures as well as CSV files. # @param surface If True, then generate solution surface with 3d plot. # @param time_min Beginning of user defined time range for plotting purposes. # @param time_max End of user defined time range for plotting purposes. # @param time_thinning ?? # @param time_unit ?? # @param title_on If True, export standard graphics with title. # @param use_cache If True, use caching. # @param verbose If True, this function is verbose. def _sww2timeseries(swwfiles, gauge_filename, production_dirs, report = None, reportname = None, plot_quantity = None, generate_fig = False, surface = None, time_min = None, time_max = None, time_thinning = 1, time_unit = None, title_on = None, use_cache = False, verbose = False, output_centroids = False): # FIXME(Ole): Shouldn't print statements here be governed by verbose? assert type(gauge_filename) == type(''), 'Gauge filename must be a string' try: fid = open(gauge_filename) except Exception, e: msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e) raise msg if report is None: report = False if plot_quantity is None: plot_quantity = ['depth', 'speed'] else: assert type(plot_quantity) == list, 'plot_quantity must be a list' check_list(plot_quantity) if surface is None: surface = False if time_unit is None: time_unit = 'hours' if title_on is None: title_on = True if verbose: log.critical('Gauges obtained from: %s' % gauge_filename) gauges, locations, elev = gauge_get_from_file(gauge_filename) sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] file_loc = [] f_list = [] label_id = [] leg_label = [] themaxT = 0.0 theminT = 0.0 for swwfile in swwfiles.keys(): try: fid = open(swwfile) except Exception, e: msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e) raise msg if verbose: log.critical('swwfile = %s' % swwfile) # Extract parent dir name and use as label path, _ = os.path.split(swwfile) _, label = os.path.split(path) leg_label.append(label) f = file_function(swwfile, quantities = sww_quantity, interpolation_points = gauges, time_thinning = time_thinning, verbose = verbose, use_cache = use_cache, output_centroids = output_centroids) # determine which gauges are contained in sww file count = 0 gauge_index = [] for k, g in enumerate(gauges): if f(0.0, point_id = k)[2] > 1.0e6: count += 1 if count == 1: log.critical('Gauges not contained here:') log.critical(locations[k]) else: gauge_index.append(k) if len(gauge_index) > 0: log.critical('Gauges contained here:') else: log.critical('No gauges contained here.') for i in range(len(gauge_index)): log.critical(locations[gauge_index[i]]) index = swwfile.rfind(sep) file_loc.append(swwfile[:index+1]) label_id.append(swwfiles[swwfile]) f_list.append(f) maxT = max(f.get_time()) minT = min(f.get_time()) if maxT > themaxT: themaxT = maxT if minT > theminT: theminT = minT if time_min is None: time_min = theminT # min(T) else: if time_min < theminT: # min(T): msg = 'Minimum time entered not correct - please try again' raise Exception, msg if time_max is None: time_max = themaxT # max(T) else: if time_max > themaxT: # max(T): msg = 'Maximum time entered not correct - please try again' raise Exception, msg if verbose and len(gauge_index) > 0: log.critical('Inputs OK - going to generate figures') if len(gauge_index) <> 0: texfile, elev_output = \ _generate_figures(plot_quantity, file_loc, report, reportname, surface, leg_label, f_list, gauges, locations, elev, gauge_index, production_dirs, time_min, time_max, time_unit, title_on, label_id, generate_fig, verbose) else: texfile = '' elev_output = [] return texfile, elev_output ## # @brief Read gauge info from a file. # @param filename The name of the file to read. # @return A (gauges, gaugelocation, elev) tuple. def gauge_get_from_file(filename): """ Read in gauge information from file """ from os import sep, getcwd, access, F_OK, mkdir # Get data from the gauge file fid = open(filename) lines = fid.readlines() fid.close() gauges = [] gaugelocation = [] elev = [] # Check header information line1 = lines[0] line11 = line1.split(',') if isinstance(line11[0], str) is True: # We have found text in the first line east_index = None north_index = None name_index = None elev_index = None for i in range(len(line11)): if line11[i].strip().lower() == 'easting': east_index = i if line11[i].strip().lower() == 'northing': north_index = i if line11[i].strip().lower() == 'name': name_index = i if line11[i].strip().lower() == 'elevation': elev_index = i if east_index < len(line11) and north_index < len(line11): pass else: msg = 'WARNING: %s does not contain correct header information' \ % filename msg += 'The header must be: easting, northing, name, elevation' raise Exception, msg if elev_index is None: raise Exception if name_index is None: raise Exception lines = lines[1:] # Remove header from data else: # No header, assume that this is a simple easting, northing file msg = 'There was no header in file %s and the number of columns is %d' \ % (filename, len(line11)) msg += '- was assuming two columns corresponding to Easting and Northing' assert len(line11) == 2, msg east_index = 0 north_index = 1 N = len(lines) elev = [-9999]*N gaugelocation = range(N) # Read in gauge data for line in lines: fields = line.split(',') gauges.append([float(fields[east_index]), float(fields[north_index])]) if len(fields) > 2: elev.append(float(fields[elev_index])) loc = fields[name_index] gaugelocation.append(loc.strip('\n')) return gauges, gaugelocation, elev ## # @brief Generate figures from quantities and gauges for each sww file. # This is a legacy function, used only by sww2timeseries # @param plot_quantity ?? # @param file_loc ?? # @param report ?? # @param reportname ?? # @param surface ?? # @param leg_label ?? # @param f_list ?? # @param gauges ?? # @param locations ?? # @param elev ?? # @param gauge_index ?? # @param production_dirs ?? # @param time_min ?? # @param time_max ?? # @param time_unit ?? # @param title_on ?? # @param label_id ?? # @param generate_fig ?? # @param verbose?? # @return (texfile2, elev_output) def _generate_figures(plot_quantity, file_loc, report, reportname, surface, leg_label, f_list, gauges, locations, elev, gauge_index, production_dirs, time_min, time_max, time_unit, title_on, label_id, generate_fig, verbose): """ Generate figures based on required quantities and gauges for each sww file """ from os import sep, altsep, getcwd, mkdir, access, F_OK, environ if generate_fig is True: from pylab import ion, hold, plot, axis, figure, legend, savefig, \ xlabel, ylabel, title, close, subplot if surface is True: import pylab as p1 import mpl3d.mplot3d as p3 if report == True: texdir = getcwd()+sep+'report'+sep if access(texdir,F_OK) == 0: mkdir (texdir) if len(label_id) == 1: label_id1 = label_id[0].replace(sep,'') label_id2 = label_id1.replace('_','') texfile = texdir + reportname + '%s' % label_id2 texfile2 = reportname + '%s' % label_id2 texfilename = texfile + '.tex' fid = open(texfilename, 'w') if verbose: log.critical('Latex output printed to %s' % texfilename) else: texfile = texdir+reportname texfile2 = reportname texfilename = texfile + '.tex' fid = open(texfilename, 'w') if verbose: log.critical('Latex output printed to %s' % texfilename) else: texfile = '' texfile2 = '' p = len(f_list) n = [] n0 = 0 for i in range(len(f_list)): n.append(len(f_list[i].get_time())) if n[i] > n0: n0 = n[i] n0 = int(n0) m = len(locations) model_time = num.zeros((n0, m, p), num.float) stages = num.zeros((n0, m, p), num.float) elevations = num.zeros((n0, m, p), num.float) momenta = num.zeros((n0, m, p), num.float) xmom = num.zeros((n0, m, p), num.float) ymom = num.zeros((n0, m, p), num.float) speed = num.zeros((n0, m, p), num.float) bearings = num.zeros((n0, m, p), num.float) due_east = 90.0*num.ones((n0, 1), num.float) due_west = 270.0*num.ones((n0, 1), num.float) depths = num.zeros((n0, m, p), num.float) eastings = num.zeros((n0, m, p), num.float) min_stages = [] max_stages = [] min_momentums = [] max_momentums = [] max_xmomentums = [] max_ymomentums = [] min_xmomentums = [] min_ymomentums = [] max_speeds = [] min_speeds = [] max_depths = [] model_time_plot3d = num.zeros((n0, m), num.float) stages_plot3d = num.zeros((n0, m), num.float) eastings_plot3d = num.zeros((n0, m),num.float) if time_unit is 'mins': scale = 60.0 if time_unit is 'hours': scale = 3600.0 ##### loop over each swwfile ##### for j, f in enumerate(f_list): if verbose: log.critical('swwfile %d of %d' % (j, len(f_list))) starttime = f.starttime comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv' fid_compare = open(comparefile, 'w') file0 = file_loc[j] + 'gauges_t0.csv' fid_0 = open(file0, 'w') ##### loop over each gauge ##### for k in gauge_index: if verbose: log.critical('Gauge %d of %d' % (k, len(gauges))) g = gauges[k] min_stage = 10 max_stage = 0 max_momentum = max_xmomentum = max_ymomentum = 0 min_momentum = min_xmomentum = min_ymomentum = 100 max_speed = 0 min_speed = 0 max_depth = 0 gaugeloc = str(locations[k]) thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \ + gaugeloc + '.csv' if j == 0: fid_out = open(thisfile, 'w') s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n' fid_out.write(s) #### generate quantities ####### for i, t in enumerate(f.get_time()): if time_min <= t <= time_max: w = f(t, point_id = k)[0] z = f(t, point_id = k)[1] uh = f(t, point_id = k)[2] vh = f(t, point_id = k)[3] depth = w-z m = sqrt(uh*uh + vh*vh) if depth < 0.001: vel = 0.0 else: vel = m / (depth + 1.e-6/depth) bearing = calc_bearing(uh, vh) model_time[i,k,j] = (t + starttime)/scale #t/60.0 stages[i,k,j] = w elevations[i,k,j] = z xmom[i,k,j] = uh ymom[i,k,j] = vh momenta[i,k,j] = m speed[i,k,j] = vel bearings[i,k,j] = bearing depths[i,k,j] = depth thisgauge = gauges[k] eastings[i,k,j] = thisgauge[0] s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \ % (t, w, m, vel, z, uh, vh, bearing) fid_out.write(s) if t == 0: s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w) fid_0.write(s) if t/60.0 <= 13920: tindex = i if w > max_stage: max_stage = w if w < min_stage: min_stage = w if m > max_momentum: max_momentum = m if m < min_momentum: min_momentum = m if uh > max_xmomentum: max_xmomentum = uh if vh > max_ymomentum: max_ymomentum = vh if uh < min_xmomentum: min_xmomentum = uh if vh < min_ymomentum: min_ymomentum = vh if vel > max_speed: max_speed = vel if vel < min_speed: min_speed = vel if z > 0 and depth > max_depth: max_depth = depth s = '%.2f, %.2f, %.2f, %.2f, %s\n' \ % (max_stage, min_stage, z, thisgauge[0], leg_label[j]) fid_compare.write(s) max_stages.append(max_stage) min_stages.append(min_stage) max_momentums.append(max_momentum) max_xmomentums.append(max_xmomentum) max_ymomentums.append(max_ymomentum) min_xmomentums.append(min_xmomentum) min_ymomentums.append(min_ymomentum) min_momentums.append(min_momentum) max_depths.append(max_depth) max_speeds.append(max_speed) min_speeds.append(min_speed) #### finished generating quantities for each swwfile ##### model_time_plot3d[:,:] = model_time[:,:,j] stages_plot3d[:,:] = stages[:,:,j] eastings_plot3d[:,] = eastings[:,:,j] if surface is True: log.critical('Printing surface figure') for i in range(2): fig = p1.figure(10) ax = p3.Axes3D(fig) if len(gauges) > 80: ax.plot_surface(model_time[:,:,j], eastings[:,:,j], stages[:,:,j]) else: ax.plot3D(num.ravel(eastings[:,:,j]), num.ravel(model_time[:,:,j]), num.ravel(stages[:,:,j])) ax.set_xlabel('time') ax.set_ylabel('x') ax.set_zlabel('stage') fig.add_axes(ax) p1.show() surfacefig = 'solution_surface%s' % leg_label[j] p1.savefig(surfacefig) p1.close() #### finished generating quantities for all swwfiles ##### # x profile for given time if surface is True: figure(11) plot(eastings[tindex,:,j], stages[tindex,:,j]) xlabel('x') ylabel('stage') profilefig = 'solution_xprofile' savefig('profilefig') elev_output = [] if generate_fig is True: depth_axis = axis([starttime/scale, time_max/scale, -0.1, max(max_depths)*1.1]) stage_axis = axis([starttime/scale, time_max/scale, min(min_stages), max(max_stages)*1.1]) vel_axis = axis([starttime/scale, time_max/scale, min(min_speeds), max(max_speeds)*1.1]) mom_axis = axis([starttime/scale, time_max/scale, min(min_momentums), max(max_momentums)*1.1]) xmom_axis = axis([starttime/scale, time_max/scale, min(min_xmomentums), max(max_xmomentums)*1.1]) ymom_axis = axis([starttime/scale, time_max/scale, min(min_ymomentums), max(max_ymomentums)*1.1]) cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k'] nn = len(plot_quantity) no_cols = 2 if len(label_id) > 1: graphname_report = [] pp = 1 div = 11. cc = 0 for k in gauge_index: g = gauges[k] count1 = 0 if report == True and len(label_id) > 1: s = '\\begin{figure}[ht] \n' \ '\\centering \n' \ '\\begin{tabular}{cc} \n' fid.write(s) if len(label_id) > 1: graphname_report = [] #### generate figures for each gauge #### for j, f in enumerate(f_list): ion() hold(True) count = 0 where1 = 0 where2 = 0 word_quantity = '' if report == True and len(label_id) == 1: s = '\\begin{figure}[hbt] \n' \ '\\centering \n' \ '\\begin{tabular}{cc} \n' fid.write(s) for which_quantity in plot_quantity: count += 1 where1 += 1 figure(count, frameon = False) if which_quantity == 'depth': plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j]) units = 'm' axis(depth_axis) if which_quantity == 'stage': if elevations[0,k,j] <= 0: plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j]) axis(stage_axis) else: plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j]) #axis(depth_axis) units = 'm' if which_quantity == 'momentum': plot(model_time[0:n[j]-1,k,j], momenta[0:n[j]-1,k,j], '-', c = cstr[j]) axis(mom_axis) units = 'm^2 / sec' if which_quantity == 'xmomentum': plot(model_time[0:n[j]-1,k,j], xmom[0:n[j]-1,k,j], '-', c = cstr[j]) axis(xmom_axis) units = 'm^2 / sec' if which_quantity == 'ymomentum': plot(model_time[0:n[j]-1,k,j], ymom[0:n[j]-1,k,j], '-', c = cstr[j]) axis(ymom_axis) units = 'm^2 / sec' if which_quantity == 'speed': plot(model_time[0:n[j]-1,k,j], speed[0:n[j]-1,k,j], '-', c = cstr[j]) axis(vel_axis) units = 'm / sec' if which_quantity == 'bearing': plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-', model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.') units = 'degrees from North' #ax = axis([time_min, time_max, 0.0, 360.0]) legend(('Bearing','West','East')) if time_unit is 'mins': xlabel('time (mins)') if time_unit is 'hours': xlabel('time (hours)') #if which_quantity == 'stage' \ # and elevations[0:n[j]-1,k,j] > 0: # ylabel('%s (%s)' %('depth', units)) #else: # ylabel('%s (%s)' %(which_quantity, units)) #ylabel('%s (%s)' %('wave height', units)) ylabel('%s (%s)' %(which_quantity, units)) if len(label_id) > 1: legend((leg_label),loc='upper right') #gaugeloc1 = gaugeloc.replace(' ','') #gaugeloc2 = gaugeloc1.replace('_','') gaugeloc2 = str(locations[k]).replace(' ','') graphname = '%sgauge%s_%s' %(file_loc[j], gaugeloc2, which_quantity) if report == True and len(label_id) > 1: figdir = getcwd()+sep+'report_figures'+sep if access(figdir,F_OK) == 0 : mkdir (figdir) latex_file_loc = figdir.replace(sep,altsep) # storing files in production directory graphname_latex = '%sgauge%s%s' \ % (latex_file_loc, gaugeloc2, which_quantity) # giving location in latex output file graphname_report_input = '%sgauge%s%s' % \ ('..' + altsep + 'report_figures' + altsep, gaugeloc2, which_quantity) graphname_report.append(graphname_report_input) # save figures in production directory for report savefig(graphname_latex) if report == True: figdir = getcwd() + sep + 'report_figures' + sep if access(figdir,F_OK) == 0: mkdir(figdir) latex_file_loc = figdir.replace(sep,altsep) if len(label_id) == 1: # storing files in production directory graphname_latex = '%sgauge%s%s%s' % \ (latex_file_loc, gaugeloc2, which_quantity, label_id2) # giving location in latex output file graphname_report = '%sgauge%s%s%s' % \ ('..' + altsep + 'report_figures' + altsep, gaugeloc2, which_quantity, label_id2) s = '\includegraphics' \ '[width=0.49\linewidth, height=50mm]{%s%s}' % \ (graphname_report, '.png') fid.write(s) if where1 % 2 == 0: s = '\\\\ \n' where1 = 0 else: s = '& \n' fid.write(s) savefig(graphname_latex) if title_on == True: title('%s scenario: %s at %s gauge' % \ (label_id, which_quantity, gaugeloc2)) #title('Gauge %s (MOST elevation %.2f, ' \ # 'ANUGA elevation %.2f)' % \ # (gaugeloc2, elevations[10,k,0], # elevations[10,k,1])) savefig(graphname) # save figures with sww file if report == True and len(label_id) == 1: for i in range(nn-1): if nn > 2: if plot_quantity[i] == 'stage' \ and elevations[0,k,j] > 0: word_quantity += 'depth' + ', ' else: word_quantity += plot_quantity[i] + ', ' else: if plot_quantity[i] == 'stage' \ and elevations[0,k,j] > 0: word_quantity += 'depth' + ', ' else: word_quantity += plot_quantity[i] if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0: word_quantity += ' and ' + 'depth' else: word_quantity += ' and ' + plot_quantity[nn-1] caption = 'Time series for %s at %s location ' \ '(elevation %.2fm)' % \ (word_quantity, locations[k], elev[k]) if elev[k] == 0.0: caption = 'Time series for %s at %s location ' \ '(elevation %.2fm)' % \ (word_quantity, locations[k], elevations[0,k,j]) east = gauges[0] north = gauges[1] elev_output.append([locations[k], east, north, elevations[0,k,j]]) label = '%sgauge%s' % (label_id2, gaugeloc2) s = '\end{tabular} \n' \ '\\caption{%s} \n' \ '\label{fig:%s} \n' \ '\end{figure} \n \n' % (caption, label) fid.write(s) cc += 1 if cc % 6 == 0: fid.write('\\clearpage \n') savefig(graphname_latex) if report == True and len(label_id) > 1: for i in range(nn-1): if nn > 2: if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: word_quantity += 'depth' + ',' else: word_quantity += plot_quantity[i] + ', ' else: if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0: word_quantity += 'depth' else: word_quantity += plot_quantity[i] where1 = 0 count1 += 1 index = j*len(plot_quantity) for which_quantity in plot_quantity: where1 += 1 s = '\includegraphics' \ '[width=0.49\linewidth, height=50mm]{%s%s}' % \ (graphname_report[index], '.png') index += 1 fid.write(s) if where1 % 2 == 0: s = '\\\\ \n' where1 = 0 else: s = '& \n' fid.write(s) word_quantity += ' and ' + plot_quantity[nn-1] label = 'gauge%s' %(gaugeloc2) caption = 'Time series for %s at %s location ' \ '(elevation %.2fm)' % \ (word_quantity, locations[k], elev[k]) if elev[k] == 0.0: caption = 'Time series for %s at %s location ' \ '(elevation %.2fm)' % \ (word_quantity, locations[k], elevations[0,k,j]) thisgauge = gauges[k] east = thisgauge[0] north = thisgauge[1] elev_output.append([locations[k], east, north, elevations[0,k,j]]) s = '\end{tabular} \n' \ '\\caption{%s} \n' \ '\label{fig:%s} \n' \ '\end{figure} \n \n' % (caption, label) fid.write(s) if float((k+1)/div - pp) == 0.: fid.write('\\clearpage \n') pp += 1 #### finished generating figures ### close('all') return texfile2, elev_output