"""This module contains various auxiliary function used by pyvolution. It is also a clearing house for functions that may later earn a module of their own. """ import utilities.polygon from warnings import warn from geospatial_data.geospatial_data import ensure_absolute def file_function(filename, domain = None, quantities = None, interpolation_points = None, verbose = False, use_cache = False): """Read time history of spatial data from NetCDF file and return a callable object. Input variables: filename - Name of sww or tms file If the file has extension 'sww' then it is assumed to be spatio-temporal or temporal and the callable object will have the form f(t,x,y) or f(t) depending on whether the file contains spatial data If the file has extension 'tms' then it is assumed to be temporal only and the callable object will have the form f(t) Either form will return interpolated values based on the input file using the underlying interpolation_function. domain - Associated domain object If domain is specified, model time (domain.starttime) will be checked and possibly modified. All times are assumed to be in UTC All spatial information is assumed to be in absolute UTM coordinates. quantities - the name of the quantity to be interpolated or a list of quantity names. The resulting function will return a tuple of values - one for each quantity interpolation_points - list of absolute UTM coordinates for points (N x 2) or geospatial object or points file name at which values are sought use_cache: True means that caching of intermediate result of Interpolation_function is attempted See Interpolation function for further documentation """ if use_cache is True: try: from caching import cache except: msg = 'Caching was requested, but caching module'+\ 'could not be imported' raise msg f = cache(_file_function, filename, {'domain': domain, 'quantities': quantities, 'interpolation_points': interpolation_points, 'verbose': verbose}, dependencies = [filename], compression = False, verbose = verbose) #FIXME (Ole): Pass cache arguments, such as compression, in some sort of #structure else: f = _file_function(filename, domain, quantities, interpolation_points, verbose) return f def _file_function(filename, domain = None, quantities = None, interpolation_points = None, verbose = False): """Internal function See file_function for documentatiton """ #FIXME (OLE): Should check origin of domain against that of file #In fact, this is where origin should be converted to that of domain #Also, check that file covers domain fully. #Take into account: #- domain's georef #- sww file's georef #- interpolation points as absolute UTM coordinates assert type(filename) == type(''),\ 'First argument to File_function must be a string' try: fid = open(filename) except Exception, e: msg = 'File "%s" could not be opened: Error="%s"'\ %(filename, e) raise msg line = fid.readline() fid.close() if quantities is None: if domain is not None: quantities = domain.conserved_quantities if line[:3] == 'CDF': return get_netcdf_file_function(filename, domain, quantities, interpolation_points, verbose = verbose) else: raise 'Must be a NetCDF File' def get_netcdf_file_function(filename, domain=None, quantity_names=None, interpolation_points=None, verbose = False): """Read time history of spatial data from NetCDF sww file and return a callable object f(t,x,y) which will return interpolated values based on the input file. If domain is specified, model time (domain.starttime) will be checked and possibly modified All times are assumed to be in UTC See Interpolation function for further documetation """ #FIXME: Check that model origin is the same as file's origin #(both in UTM coordinates) #If not - modify those from file to match domain #Take this code from e.g. dem2pts in data_manager.py #FIXME: Use geo_reference to read and write xllcorner... #FIXME: Maybe move caching out to this level rather than at the #Interpolation_function level (below) import time, calendar, types from config import time_format from Scientific.IO.NetCDF import NetCDFFile from Numeric import array, zeros, Float, alltrue, concatenate, reshape #Open NetCDF file if verbose: print 'Reading', filename fid = NetCDFFile(filename, 'r') if type(quantity_names) == types.StringType: quantity_names = [quantity_names] if quantity_names is None or len(quantity_names) < 1: #If no quantities are specified get quantities from file #x, y, time are assumed as the independent variables so #they are excluded as potentiol quantities quantity_names = [] for name in fid.variables.keys(): if name not in ['x', 'y', 'time']: quantity_names.append(name) if len(quantity_names) < 1: msg = 'ERROR: At least one independent value must be specified' raise msg if interpolation_points is not None: interpolation_points = ensure_absolute(interpolation_points) msg = 'Points must by N x 2. I got %d' %interpolation_points.shape[1] assert interpolation_points.shape[1] == 2, msg #Now assert that requested quantitites (and the independent ones) #are present in file missing = [] for quantity in ['time'] + quantity_names: if not fid.variables.has_key(quantity): missing.append(quantity) if len(missing) > 0: msg = 'Quantities %s could not be found in file %s'\ %(str(missing), filename) fid.close() raise msg #Decide whether this data has a spatial dimension spatial = True for quantity in ['x', 'y']: if not fid.variables.has_key(quantity): spatial = False if filename[-3:] == 'tms' and spatial is True: msg = 'Files of type tms must not contain spatial information' raise msg if filename[-3:] == 'sww' and spatial is False: msg = 'Files of type sww must contain spatial information' raise msg #Get first timestep try: starttime = fid.starttime[0] except ValueError: msg = 'Could not read starttime from file %s' %filename raise msg #Get variables if verbose: print 'Get variables' time = fid.variables['time'][:] if spatial: #Get origin xllcorner = fid.xllcorner[0] yllcorner = fid.yllcorner[0] zone = fid.zone[0] x = fid.variables['x'][:] y = fid.variables['y'][:] triangles = fid.variables['volumes'][:] x = reshape(x, (len(x),1)) y = reshape(y, (len(y),1)) vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array if interpolation_points is not None: #Adjust for georef interpolation_points[:,0] -= xllcorner interpolation_points[:,1] -= yllcorner if domain is not None: #Update domain.startime if it is *earlier* than starttime if starttime > domain.starttime: msg = 'WARNING: Start time as specified in domain (%f)'\ %domain.starttime msg += ' is earlier than the starttime of file %s (%f).'\ %(filename, starttime) msg += ' Modifying domain starttime accordingly.' if verbose: print msg domain.starttime = starttime #Modifying model time if verbose: print 'Domain starttime is now set to %f'\ %domain.starttime #If domain.startime is *later* than starttime, #move time back - relative to domain's time if domain.starttime > starttime: time = time - domain.starttime + starttime #FIXME Use method in geo to reconcile #if spatial: #assert domain.geo_reference.xllcorner == xllcorner #assert domain.geo_reference.yllcorner == yllcorner #assert domain.geo_reference.zone == zone if verbose: print 'File_function data obtained from: %s' %filename print ' References:' #print ' Datum ....' #FIXME if spatial: print ' Lower left corner: [%f, %f]'\ %(xllcorner, yllcorner) print ' Start time: %f' %starttime #Produce values for desired data points at #each timestep quantities = {} for i, name in enumerate(quantity_names): quantities[name] = fid.variables[name][:] fid.close() #from least_squares import Interpolation_function from fit_interpolate.interpolate import Interpolation_function if not spatial: vertex_coordinates = triangles = interpolation_points = None return Interpolation_function(time, quantities, quantity_names, vertex_coordinates, triangles, interpolation_points, verbose = verbose) def multiple_replace(text, dictionary): """Multiple replace of words in text text: String to be modified dictionary: Mapping of words that are to be substituted Python Cookbook 3.14 page 88 and page 90 """ import re #Create a regular expression from all of the dictionary keys #matching only entire words regex = re.compile(r'\b'+ \ r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \ r'\b' ) #For each match, lookup the corresponding value in the dictionary return regex.sub(lambda match: dictionary[match.group(0)], text) def apply_expression_to_dictionary(expression, dictionary):#dictionary): """Apply arbitrary expression to values of dictionary Given an expression in terms of the keys, replace key by the corresponding values and evaluate. expression: Arbitrary, e.g. arithmetric, expression relating keys from dictionary. dictionary: Mapping between symbols in expression and objects that will be evaluated by expression. Values in dictionary must support operators given in expression e.g. by overloading due to a limitation with Numeric, this can not evaluate 0/0 In general, the user can fix by adding 1e-30 to the numerator. SciPy core can handle this situation. """ import types import re assert isinstance(expression, basestring) assert type(dictionary) == types.DictType #Convert dictionary values to textual representations suitable for eval D = {} for key in dictionary: D[key] = 'dictionary["%s"]' %key #Perform substitution of variables expression = multiple_replace(expression, D) #Evaluate and return try: return eval(expression) except NameError, e: msg = 'Expression "%s" could not be evaluated: %s' %(expression, e) raise NameError, msg except ValueError, e: msg = 'Expression "%s" could not be evaluated: %s' %(expression, e) raise ValueError, msg #################################### ####OBSOLETE STUFF def angle(v1, v2): """Temporary Interface to new location""" import utilities.numerical_tools as NT msg = 'angle has moved from util.py. ' msg += 'Please use "from utilities.numerical_tools import angle"' warn(msg, DeprecationWarning) return NT.angle(v1, v2) def anglediff(v0, v1): """Temporary Interface to new location""" import utilities.numerical_tools as NT msg = 'anglediff has moved from util.py. ' msg += 'Please use "from utilities.numerical_tools import anglediff"' warn(msg, DeprecationWarning) return NT.anglediff(v0, v1) def mean(x): """Temporary Interface to new location""" import utilities.numerical_tools as NT msg = 'mean has moved from util.py. ' msg += 'Please use "from utilities.numerical_tools import mean"' warn(msg, DeprecationWarning) return NT.mean(x) def point_on_line(*args, **kwargs): """Temporary Interface to new location""" msg = 'point_on_line has moved from util.py. ' msg += 'Please use "from utilities.polygon import point_on_line"' warn(msg, DeprecationWarning) return utilities.polygon.point_on_line(*args, **kwargs) def inside_polygon(*args, **kwargs): """Temporary Interface to new location""" print 'inside_polygon has moved from util.py. ', print 'Please use "from utilities.polygon import inside_polygon"' return utilities.polygon.inside_polygon(*args, **kwargs) def outside_polygon(*args, **kwargs): """Temporary Interface to new location""" print 'outside_polygon has moved from util.py. ', print 'Please use "from utilities.polygon import outside_polygon"' return utilities.polygon.outside_polygon(*args, **kwargs) def separate_points_by_polygon(*args, **kwargs): """Temporary Interface to new location""" print 'separate_points_by_polygon has moved from util.py. ', print 'Please use "from utilities.polygon import separate_points_by_polygon"' return utilities.polygon.separate_points_by_polygon(*args, **kwargs) def read_polygon(*args, **kwargs): """Temporary Interface to new location""" print 'read_polygon has moved from util.py. ', print 'Please use "from utilities.polygon import read_polygon"' return utilities.polygon.read_polygon(*args, **kwargs) def populate_polygon(*args, **kwargs): """Temporary Interface to new location""" print 'populate_polygon has moved from util.py. ', print 'Please use "from utilities.polygon import populate_polygon"' return utilities.polygon.populate_polygon(*args, **kwargs) ##################### end of obsolete stuff ? ############ ''' this simply catches the screen output and files it to a file ''' from os import sep class Screen_Catcher: def __init__(self, filename): # self.data = '' self.filename = filename def write(self, stuff): fid = open(self.filename, 'a') # fid = open(self.filename, 'w') # self.data = self.data + stuff fid.write(stuff) fid.close() def get_version_info(): """gets the version number of the SVN """ import os, sys # Create dummy info info = 'Revision: Version info could not be obtained.' info += 'A command line version of svn and access to the ' info += 'repository is necessary and the output must ' info += 'contain a line starting with "Revision:"' try: fid = os.popen('svn info') except: msg = 'svn is not recognised' warn(msg, UserWarning) else: lines = fid.readlines() fid.close() for line in lines: if line.startswith('Revision:'): info = line break return info #if sys.platform == 'win32': # msg = 'does not work in Windows .... yet' # raise OSError, msg #else: # fid = os.popen('svn -info') # info = fid.readlines() # fid.close() # return info def sww2timeseries(swwfiles, gauge_filename, production_dirs, report = None, reportname = None, plot_quantity = None, surface = None, time_min = None, time_max = None, title_on = None, verbose = 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 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'] 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. """ k = _sww2timeseries(swwfiles, gauge_filename, production_dirs, report, reportname, plot_quantity, surface, time_min, time_max, title_on, verbose) return k def _sww2timeseries(swwfiles, gauge_filename, production_dirs, report = None, reportname = None, plot_quantity = None, surface = None, time_min = None, time_max = None, title_on = None, verbose = False): 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', 'bearing'] else: assert type(plot_quantity) == list,\ 'plot_quantity must be a list' check_list(plot_quantity) if surface is None: surface = False if title_on is None: title_on = True if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename gauges, locations, elev = get_gauges_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 f = file_function(swwfile, quantities = sww_quantity, interpolation_points = gauges, verbose = True, use_cache = True) index = swwfile.rfind(sep) file_loc.append(swwfile[:index+1]) label_id.append(swwfiles[swwfile]) leg_label.append(production_dirs[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: print 'Inputs OK - going to generate figures' return generate_figures(plot_quantity, file_loc, report, reportname, surface, leg_label, f_list, gauges, locations, elev, production_dirs, time_min, time_max, title_on, label_id, verbose) #Fixme - Use geospatial to read this file - it's an xya file def get_gauges_from_file(filename): from os import sep, getcwd, access, F_OK, mkdir fid = open(filename) lines = fid.readlines() fid.close() gauges = [] gaugelocation = [] elev = [] line1 = lines[0] line11 = line1.split(',') for i in range(len(line11)): if line11[i].strip('\n').strip(' ') == 'Easting': east_index = i if line11[i].strip('\n').strip(' ') == 'Northing': north_index = i if line11[i].strip('\n').strip(' ') == 'Name': name_index = i if line11[i].strip('\n').strip(' ') == 'Elevation': elev_index = i for line in lines[1:]: fields = line.split(',') gauges.append([float(fields[east_index]), float(fields[north_index])]) elev.append(float(fields[elev_index])) loc = fields[name_index] gaugelocation.append(loc.strip('\n')) return gauges, gaugelocation, elev def check_list(quantity): all_quantity = ['stage', 'depth', 'momentum', 'xmomentum', 'ymomentum', 'speed', 'bearing', 'elevation'] p = list(set(quantity).difference(set(all_quantity))) if len(p) <> 0: msg = 'Quantities %s do not exist - please try again' %p raise Exception, msg return def calc_bearing(uh, vh): from math import atan, degrees angle = degrees(atan(vh/(uh+1.e-15))) if (0 < angle < 90.0): if vh > 0: bearing = 90.0 - abs(angle) if vh < 0: bearing = 270.0 - abs(angle) if (-90 < angle < 0): if vh < 0: bearing = 90.0 - abs(angle) if vh > 0: bearing = 270.0 - abs(angle) if angle == 0: bearing = 0.0 return bearing def generate_figures(plot_quantity, file_loc, report, reportname, surface, leg_label, f_list, gauges, locations, elev, production_dirs, time_min, time_max, title_on, label_id, verbose): from math import sqrt, atan, degrees from Numeric import ones, allclose, zeros, Float, ravel from os import sep, altsep, getcwd, mkdir, access, F_OK, environ from pylab import ion, hold, plot, axis, figure, legend, savefig, \ xlabel, ylabel, title, close, subplot 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' if verbose: print '\n Latex output printed to %s \n' %texfilename fid = open(texfilename, 'w') else: texfile = texdir+reportname texfile2 = reportname texfilename = texfile + '.tex' if verbose: print '\n Latex output printed to %s \n' %texfilename fid = open(texfilename, 'w') else: texfile = '' 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 = zeros((n0,m,p), Float) stages = zeros((n0,m,p), Float) elevations = zeros((n0,m,p), Float) momenta = zeros((n0,m,p), Float) xmom = zeros((n0,m,p), Float) ymom = zeros((n0,m,p), Float) speed = zeros((n0,m,p), Float) bearings = zeros((n0,m,p), Float) depths = zeros((n0,m,p), Float) eastings = zeros((n0,m,p), Float) min_stages = [] max_stages = [] max_momentums = [] max_speeds = [] c = 0 model_time_plot3d = zeros((n0,m), Float) stages_plot3d = zeros((n0,m), Float) eastings_plot3d = zeros((n0,m),Float) ##### loop over each swwfile ##### for j, f in enumerate(f_list): if verbose: print 'swwfile %d of %d' %(j, len(f_list)) comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv' fid_compare = open(comparefile, 'w') ##### loop over each gauge ##### for k, g in enumerate(gauges): if verbose: print 'Gauge %d of %d' %(k, len(gauges)) min_stage = 10 max_stage = 0 max_momentum = 0 max_speed = 0 gaugeloc = locations[k] thisfile = file_loc[j]+sep+'gauges_time_series'+'_'+gaugeloc+'.csv' fid_out = open(thisfile, 'w') s = 'Time, Stage, Momentum, Speed \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-30) bearing = calc_bearing(uh, vh) model_time[i,k,j] = 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\n' %(t, w, m, vel) fid_out.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 vel > max_speed: max_speed = vel s = '%.2f, %.2f, %.2f, %s\n' %(max_stage, min_stage, z, leg_label[j]) fid_compare.write(s) max_stages.append(max_stage) min_stages.append(min_stage) max_momentums.append(max_momentum) max_speeds.append(max_speed) #### finished generating quantities for each swwfile ##### model_time_plot3d[:,:] = model_time[:,:,j] stages_plot3d[:,:] = stages[:,:,j] eastings_plot3d[:,] = eastings[:,:,j] if surface == True: print '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.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j]) #ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),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 == True: figure(11) plot(eastings[tindex,:,j],stages[tindex,:,j]) xlabel('x') ylabel('stage') profilefig = 'solution_xprofile' savefig('profilefig') stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1]) vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1]) mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1]) cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k'] nn = len(plot_quantity) no_cols = 2 elev_output = [] if len(label_id) > 1: graphname_report = [] for k, g in enumerate(gauges): count1 = 0 if report == True and len(label_id) > 1: s = '\\begin{figure}[hbt] \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' if which_quantity == 'stage': plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j]) axis(stage_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(mom_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(mom_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': due_east = 90.0*ones(shape(model_time[0:n[j]-1,k,j],Float)) due_west = 270.0*ones(shape(model_time[0:n[j]-1,k,j],Float)) plot(model_time[0:n[j]-1,k,j], bearings, '-', model_time[0:n[j]-1,k,j], due_west, '-.', model_time[0:n[j]-1,k,j], due_east, '-.') units = 'degrees from North' ax = axis([time_min, time_max, 0.0, 360.0]) legend(('Bearing','West','East')) xlabel('time (mins)') ylabel('%s (%s)' %(which_quantity, units)) if len(label_id) > 1: legend((leg_label),loc='upper right') gaugeloc1 = gaugeloc.replace(' ','') #gaugeloc2 = gaugeloc1.replace('_','') gaugeloc2 = 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) graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file graphname_report.append(graphname_report_input) savefig(graphname_latex) # save figures in production directory for report generation 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: graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file 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, gaugeloc)) savefig(graphname) # save figures with sww file if report == True and len(label_id) == 1: for i in range(nn-1): if nn > 2: word_quantity += plot_quantity[i] + ', ' else: word_quantity += plot_quantity[i] word_quantity += ' and ' + plot_quantity[nn-1] caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' ')) 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) c += 1 if c % 25 == 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: word_quantity += plot_quantity[i] + ', ' 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) c += 1 if c % 25 == 0: fid.write('\\clearpage \n') #### finished generating figures ### close('all') return texfile2, elev_output