Changeset 2795

Timestamp:

May 3, 2006, 4:12:40 PM (18 years ago)

Author:

sexton

Message:

(1) inclusion of sww2timeseries which generates time series figures from sww file and which can output latex file if requested. (2) Script using sww2timeseries for sydney scenario

Files:

• inundation/pyvolution/util.py (modified) (1 diff)
• production/sydney_2006/run_timeseries.py (added)

inundation/pyvolution/util.py

@@ -502 +502 @@
         fid.write(self.data)
         fid.close()
+
+
+from pylab import *
+    
+def sww2timeseries(swwfile,
+                   gauge_filename,
+                   file_loc,
+                   label_id,
+                   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
+                          
+    file_loc        - file location of sww file
+                    - gauge graphs and data written to this directory
+    
+    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
+                    - if False, export graphics as eps and generate
+                      latex output
+                    - latex output to be written to same directory as
+                      where this function is being used from and named
+                      according to scenario name
+                      
+    label_id        - used in generating latex output. It will generally
+                      be part of the directory name of file_loc (the
+                      time stamp, if using that). Helps to differentiate
+                      between runs within a particular scenario.
+                      
+    Output:
+    
+    - time series data stored in .csv for later use if required.
+      Name = gauges_timeseries followed by gauge name 
+    - if title_on = False then output latex file will be generated
+      in same directory as where script is run (usually production
+      scenario direcotry. Name = latex_output_label_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,
+                        file_loc,
+                        label_id,
+                        plot_quantity,
+                        time_min,
+                        time_max,
+                        title_on,
+                        verbose)
+
+    return k
+
+def _sww2timeseries(swwfile,
+                    gauge_filename,
+                    file_loc,
+                    label_id,
+                    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
+
+    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
+        
+    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)
+
+    #from pyvolution.util import file_function
+
+    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+
+    f = file_function(swwfile,
+                      quantities = sww_quantity,
+                      interpolation_points = gauges,
+                      verbose = True,
+                      use_cache = True)
+
+    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(f.T)
+    else:
+        if time_min < min(f.T):
+            msg = 'Minimum time entered not correct - please try again'
+            raise Exception, msg
+
+    if time_max is None:
+        time_max = max(f.T)
+    else:
+        if time_max > max(f.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, 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, f, gauges,
+                     locations, time_min, time_max, title_on, label_id, verbose):
+
+    from math import sqrt, atan, degrees
+    from os import sep, altsep
+
+    filename = file_loc.split(sep)
+        
+    if title_on == False:
+        ext = '.eps'
+        texfilename = 'latex_output_%s.tex' %(label_id.replace(sep,'_'))
+        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, east, north)
+            fid.write(s)
+
+        s = '\\end{tabular} \n \\end{center} \n \n'
+        fid.write(s)
+    else:
+        texfilename = ''
+        ext = '' 
+    
+    ##### 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.T): 
+
+            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, 'g-')
+                units = 'm'
+            if which_quantity == 'stage':
+                plot(model_time, stages, 'g-')
+                units = 'm'
+            if which_quantity == 'momentum':
+                plot(model_time, momenta, 'g-')
+                units = 'm^2 / sec'
+            if which_quantity == 'xmomentum':
+                plot(model_time, xmom, 'g-')
+                units = 'm^2 / sec'
+            if which_quantity == 'ymomentum':
+                plot(model_time, ymom, 'g-')
+                units = 'm^2 / sec'
+            if which_quantity == 'velocity':
+                plot(model_time, velocity, 'g-')
+                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, 'r-', model_time, due_west, '-.g',
+                     model_time, due_east, '-.b')
+                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(' ','_')
+            graphname = '%splot_test_%s_%s%s' %(file_loc, which_quantity, gaugeloc1, ext)
+            latex_file_loc = file_loc.replace(sep,altsep)
+            graphname_latex = '%splot_test_%s_%s%s' %(latex_file_loc, which_quantity, gaugeloc1, ext)
+
+            if title_on == True:
+                title('%s scenario: %s at %s gauge'
+                      %(label_id, which_quantity, gaugeloc))
+            else:
+                label = '%s_%s_gauge_%s' %(label_id.replace(sep,'_'), which_quantity, gaugeloc1) 
+                caption = 'Time series for %s' %(which_quantity)
+                s = '\\begin{figure}[hbt] \n \\centerline{\includegraphics[width=75mm, height=75mm]{%s}} \n' %(graphname_latex)
+                fid.write(s)
+                s = '\\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
+                fid.write(s)
+
+            savefig(graphname)
+        
+        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 texfilename