"""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 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 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. #If we use the suggested Point_set class for interpolation points #here it would be easier #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 from utilities.numerical_tools import ensure_numeric #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_numeric(interpolation_points, Float) #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) ''' 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(swwfile, gauge_filename, label_id, report = None, plot_quantity = 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: swwfile - name of sww file - assume that all conserved quantities have been stored gauge_filename - name of file containing gauge data - name, easting, northing - OR (this is not yet done) - structure which can be converted to a Numeric array, such as a geospatial data object label_id - 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. 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 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' - velocity; calculated as absolute momentum (pointwise) divided by depth; 'velocity' - 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', 'velocity', 'bearing'] 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(swwfile, gauge_filename, label_id, report, plot_quantity, time_min, time_max, title_on, verbose) return k def _sww2timeseries(swwfile, gauge_filename, label_id, report = None, plot_quantity = None, time_min = None, time_max = None, title_on = None, verbose = False): assert type(swwfile) == type(''),\ 'The sww filename must be a string' try: fid = open(swwfile) except Exception, e: msg = 'File "%s" could not be opened: Error="%s"'\ %(swwfile, e) raise msg index = swwfile.rfind(sep) file_loc = swwfile[:index+1] 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 assert type(plot_quantity) == list,\ 'plot_quantity must be a list' if plot_quantity is None: plot_quantity = ['depth', 'velocity', 'bearing'] else: check_list(plot_quantity) if title_on is None: title_on = True assert type(label_id) == type(''),\ 'label_id to sww2timeseries must be a string' if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename gauges, locations = get_gauges_from_file(gauge_filename) sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum'] f = file_function(swwfile, quantities = sww_quantity, interpolation_points = gauges, verbose = True, use_cache = True) T = f.get_time() if max(f.quantities['xmomentum']) > 1.e10: msg = 'Not all conserved quantities available from sww file. \n sww2timeseries requires all conserved quantities stored in sww file' raise Exception, msg if time_min is None: time_min = min(T) else: if time_min < min(T): msg = 'Minimum time entered not correct - please try again' raise Exception, msg if time_max is None: time_max = max(T) else: if time_max > 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, f, gauges, locations, time_min, time_max, title_on, label_id, verbose) def get_gauges_from_file(filename): fid = open(filename) lines = fid.readlines() fid.close() gauges = [] gaugelocation = [] 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 for line in lines[1:]: fields = line.split(',') gauges.append([float(fields[east_index]), float(fields[north_index])]) loc = fields[name_index] gaugelocation.append(loc.strip('\n')) return gauges, gaugelocation def check_list(quantity): all_quantity = ['stage', 'depth', 'momentum', 'xmomentum', 'ymomentum', 'velocity', '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, f, gauges, locations, time_min, time_max, title_on, label_id, verbose): from math import sqrt, atan, degrees from Numeric import ones 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 filename = file_loc.split(sep) if report == True: label_id1 = label_id.replace(sep,'') label_id2 = label_id1.replace('_','') texfile = 'latexoutput%s' %(label_id2) texfilename = texfile + '.tex' if verbose: print '\n Latex output printed to %s \n' %texfilename fid = open(texfilename, 'w') s = '\\begin{center} \n \\begin{tabular}{|l|l|l|}\hline \n \\bf{Gauge Name} & \\bf{Easting} & \\bf{Northing} \\\\ \hline \n' fid.write(s) gauges1 = gauges for name, gauges1 in zip(locations, gauges1): east = gauges1[0] north = gauges1[1] s = '%s & %.2f & %.2f \\\\ \hline \n' %(name.replace('_',' '), east, north) fid.write(s) s = '\\end{tabular} \n \\end{center} \n \n' fid.write(s) else: texfile = '' c = 0 ##### loop over each gauge ##### for k, g in enumerate(gauges): if verbose: print 'Gauge %d of %d' %(k, len(gauges)) model_time = [] stages = [] elevations = [] momenta = [] xmom = [] ymom = [] velocity = [] bearings = [] depths = [] max_depth = 0 max_momentum = 0 max_velocity = 0 #### 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] gaugeloc = locations[k] depth = w-z m = sqrt(uh*uh + vh*vh) vel = m / (depth + 1.e-30) bearing = calc_bearing(uh, vh) model_time.append(t) stages.append(w) elevations.append(z) xmom.append(uh) ymom.append(vh) momenta.append(m) velocity.append(vel) bearings.append(bearing) depths.append(depth) if depth > max_depth: max_depth = w-z if m > max_momentum: max_momentum = m if vel > max_velocity: max_velocity = vel #### finished generating quantities ##### #### generate figures ### ion() hold(False) for i in range(len(plot_quantity)): which_quantity = plot_quantity[i] figure(i+1, frameon = False) if which_quantity == 'depth': plot(model_time, depths, '-') units = 'm' if which_quantity == 'stage': plot(model_time, stages, '-') units = 'm' if which_quantity == 'momentum': plot(model_time, momenta, '-') units = 'm^2 / sec' if which_quantity == 'xmomentum': plot(model_time, xmom, '-') units = 'm^2 / sec' if which_quantity == 'ymomentum': plot(model_time, ymom, '-') units = 'm^2 / sec' if which_quantity == 'velocity': plot(model_time, velocity, '-') units = 'm / sec' if which_quantity == 'bearing': due_east = 90.0*ones([len(model_time)]) due_west = 270.0*ones([len(model_time)]) plot(model_time, bearings, '-', model_time, due_west, '-.', model_time, due_east, '-.') units = 'degrees from North' ax = axis([time_min, time_max, 0, 360]) legend(('Bearing','West','East')) xlabel('time (secs)') ylabel('%s (%s)' %(which_quantity, units)) gaugeloc1 = gaugeloc.replace(' ','') gaugeloc2 = gaugeloc1.replace('_','') graphname = '%sgauge%s_%s' %(file_loc, gaugeloc2, which_quantity) if report == True: 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%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # giving location in latex output file graphname_report = '%sgauge%s%s%s' %('report_figures'+altsep, gaugeloc2, which_quantity, label_id2) label = '%s%sgauge%s' %(label_id2, which_quantity, gaugeloc2) caption = 'Time series for %s at %s gauge location' %(which_quantity, gaugeloc.replace('_',' ')) s = '\\begin{figure}[hbt] \n \\centerline{\includegraphics[width=100mm, height=75mm]{%s%s}} \n' %(graphname_report, '.png') fid.write(s) s = '\\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label) fid.write(s) c += 1 if c % 10 == 0: fid.write('\\clearpage \n') savefig(graphname_latex) # save figures in production directory for report generation if title_on == True: title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc)) savefig(graphname) # save figures with sww file thisfile = file_loc+sep+'gauges_time_series'+'_'+gaugeloc+'.csv' fid_out = open(thisfile, 'w') s = 'Time, Depth, Momentum, Velocity \n' fid_out.write(s) for i_t, i_d, i_m, i_vel in zip(model_time, depths, momenta, velocity): s = '%.2f, %.2f, %.2f, %.2f\n' %(i_t, i_d, i_m, i_vel) fid_out.write(s) #### finished generating figures ### close('all') return texfile