source: anuga_core/source/anuga/abstract_2d_finite_volumes/util.py @ 4076

Last change on this file since 4076 was 4076, checked in by sexton, 17 years ago

add unit test for part of sww2timeseries functionality

File size: 43.5 KB
Line 
1"""This module contains various auxiliary function used by pyvolution.
2
3It is also a clearing house for functions that may later earn a module
4of their own.
5"""
6
7import anuga.utilities.polygon
8import sys
9import os
10
11from os import remove, mkdir, access, F_OK, sep
12from os.path import exists, basename
13from warnings import warn
14from shutil import copy
15
16from anuga.geospatial_data.geospatial_data import ensure_absolute
17
18def file_function(filename,
19                  domain=None,
20                  quantities=None,
21                  interpolation_points=None,
22                  time_thinning=1,
23                  verbose=False,
24                  use_cache=False):
25    """Read time history of spatial data from NetCDF file and return
26    a callable object.
27
28    Input variables:
29   
30    filename - Name of sww or tms file
31       
32       If the file has extension 'sww' then it is assumed to be spatio-temporal
33       or temporal and the callable object will have the form f(t,x,y) or f(t)
34       depending on whether the file contains spatial data
35
36       If the file has extension 'tms' then it is assumed to be temporal only
37       and the callable object will have the form f(t)
38
39       Either form will return interpolated values based on the input file
40       using the underlying interpolation_function.
41
42    domain - Associated domain object   
43       If domain is specified, model time (domain.starttime)
44       will be checked and possibly modified.
45   
46       All times are assumed to be in UTC
47       
48       All spatial information is assumed to be in absolute UTM coordinates.
49
50    quantities - the name of the quantity to be interpolated or a
51                 list of quantity names. The resulting function will return
52                 a tuple of values - one for each quantity 
53
54    interpolation_points - list of absolute UTM coordinates for points (N x 2)
55    or geospatial object or points file name at which values are sought
56   
57    use_cache: True means that caching of intermediate result of
58               Interpolation_function is attempted
59
60   
61    See Interpolation function for further documentation
62    """
63
64
65    # Build arguments and keyword arguments for use with caching or apply.
66    args = (filename,)
67   
68    kwargs = {'domain': domain,
69              'quantities': quantities,
70              'interpolation_points': interpolation_points,
71              'time_thinning': time_thinning,                   
72              'verbose': verbose}
73
74
75    # Call underlying engine with or without caching
76    if use_cache is True:
77        try:
78            from caching import cache
79        except:
80            msg = 'Caching was requested, but caching module'+\
81                  'could not be imported'
82            raise msg
83
84        f = cache(_file_function,
85                  args, kwargs,
86                  dependencies=[filename],
87                  compression=False,                 
88                  verbose=verbose)
89
90    else:
91        f = apply(_file_function,
92                  args, kwargs)
93
94
95    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
96    #structure
97       
98
99    return f
100
101
102
103def _file_function(filename,
104                   domain=None,
105                   quantities=None,
106                   interpolation_points=None,
107                   time_thinning=1,                                               
108                   verbose=False):
109    """Internal function
110   
111    See file_function for documentatiton
112    """
113   
114
115    #FIXME (OLE): Should check origin of domain against that of file
116    #In fact, this is where origin should be converted to that of domain
117    #Also, check that file covers domain fully.
118
119    #Take into account:
120    #- domain's georef
121    #- sww file's georef
122    #- interpolation points as absolute UTM coordinates
123
124
125    assert type(filename) == type(''),\
126               'First argument to File_function must be a string'
127
128    try:
129        fid = open(filename)
130    except Exception, e:
131        msg = 'File "%s" could not be opened: Error="%s"'\
132                  %(filename, e)
133        raise msg
134
135    line = fid.readline()
136    fid.close()
137
138    if quantities is None: 
139        if domain is not None:
140            quantities = domain.conserved_quantities
141
142
143
144    if line[:3] == 'CDF':
145        return get_netcdf_file_function(filename, domain, quantities,
146                                        interpolation_points,
147                                        time_thinning=time_thinning,
148                                        verbose=verbose)
149    else:
150        raise 'Must be a NetCDF File'
151
152
153
154def get_netcdf_file_function(filename,
155                             domain=None,
156                             quantity_names=None,
157                             interpolation_points=None,
158                             time_thinning=1,                             
159                             verbose=False):
160    """Read time history of spatial data from NetCDF sww file and
161    return a callable object f(t,x,y)
162    which will return interpolated values based on the input file.
163
164    If domain is specified, model time (domain.starttime)
165    will be checked and possibly modified
166   
167    All times are assumed to be in UTC
168
169    See Interpolation function for further documetation
170   
171    """
172   
173   
174    #FIXME: Check that model origin is the same as file's origin
175    #(both in UTM coordinates)
176    #If not - modify those from file to match domain
177    #Take this code from e.g. dem2pts in data_manager.py
178    #FIXME: Use geo_reference to read and write xllcorner...
179       
180
181    #FIXME: Maybe move caching out to this level rather than at the
182    #Interpolation_function level (below)
183
184    import time, calendar, types
185    from anuga.config import time_format
186    from Scientific.IO.NetCDF import NetCDFFile
187    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
188
189    #Open NetCDF file
190    if verbose: print 'Reading', filename
191    fid = NetCDFFile(filename, 'r')
192
193    if type(quantity_names) == types.StringType:
194        quantity_names = [quantity_names]       
195   
196    if quantity_names is None or len(quantity_names) < 1:
197        #If no quantities are specified get quantities from file
198        #x, y, time are assumed as the independent variables so
199        #they are excluded as potentiol quantities
200        quantity_names = []
201        for name in fid.variables.keys():
202            if name not in ['x', 'y', 'time']:
203                quantity_names.append(name)
204
205    if len(quantity_names) < 1:               
206        msg = 'ERROR: At least one independent value must be specified'
207        raise msg
208
209
210    if interpolation_points is not None:
211        interpolation_points = ensure_absolute(interpolation_points)
212        msg = 'Points must by N x 2. I got %d' %interpolation_points.shape[1]
213        assert interpolation_points.shape[1] == 2, msg
214
215
216    #Now assert that requested quantitites (and the independent ones)
217    #are present in file
218    missing = []
219    for quantity in ['time'] + quantity_names:
220        if not fid.variables.has_key(quantity):
221            missing.append(quantity)
222
223    if len(missing) > 0:
224        msg = 'Quantities %s could not be found in file %s'\
225              %(str(missing), filename)
226        fid.close()
227        raise Exception, msg
228
229    #Decide whether this data has a spatial dimension
230    spatial = True
231    for quantity in ['x', 'y']:
232        if not fid.variables.has_key(quantity):
233            spatial = False
234
235    if filename[-3:] == 'tms' and spatial is True:
236        msg = 'Files of type tms must not contain spatial information'
237        raise msg
238
239    if filename[-3:] == 'sww' and spatial is False:
240        msg = 'Files of type sww must contain spatial information'       
241        raise msg       
242
243    #Get first timestep
244    try:
245        starttime = fid.starttime[0]
246    except ValueError:
247        msg = 'Could not read starttime from file %s' %filename
248        raise msg
249
250    #Get variables
251    #if verbose: print 'Get variables'   
252    time = fid.variables['time'][:]   
253
254    # Get time independent stuff
255    if spatial:
256        #Get origin
257        xllcorner = fid.xllcorner[0]
258        yllcorner = fid.yllcorner[0]
259        zone = fid.zone[0]       
260
261        x = fid.variables['x'][:]
262        y = fid.variables['y'][:]
263        triangles = fid.variables['volumes'][:]
264
265        x = reshape(x, (len(x),1))
266        y = reshape(y, (len(y),1))
267        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
268
269        if interpolation_points is not None:
270            #Adjust for georef
271            interpolation_points[:,0] -= xllcorner
272            interpolation_points[:,1] -= yllcorner       
273       
274
275
276
277    if domain is not None:
278        #Update domain.startime if it is *earlier* than starttime
279        if starttime > domain.starttime:
280            msg = 'WARNING: Start time as specified in domain (%f)'\
281                  %domain.starttime
282            msg += ' is earlier than the starttime of file %s (%f).'\
283                     %(filename, starttime)
284            msg += ' Modifying domain starttime accordingly.'
285           
286            if verbose: print msg
287            domain.starttime = starttime #Modifying model time
288            if verbose: print 'Domain starttime is now set to %f'\
289               %domain.starttime
290
291
292        #If domain.startime is *later* than starttime,
293        #move time back - relative to domain's time
294        if domain.starttime > starttime:
295            time = time - domain.starttime + starttime
296
297        #FIXME Use method in geo to reconcile
298        #if spatial:
299        #assert domain.geo_reference.xllcorner == xllcorner
300        #assert domain.geo_reference.yllcorner == yllcorner
301        #assert domain.geo_reference.zone == zone       
302       
303    if verbose:
304        print 'File_function data obtained from: %s' %filename
305        print '  References:'
306        #print '    Datum ....' #FIXME
307        if spatial:
308            print '    Lower left corner: [%f, %f]'\
309                  %(xllcorner, yllcorner)
310        print '    Start time:   %f' %starttime               
311       
312   
313    # Produce values for desired data points at
314    # each timestep for each quantity
315    quantities = {}
316    for i, name in enumerate(quantity_names):
317        quantities[name] = fid.variables[name][:]
318    fid.close()
319   
320
321    #from least_squares import Interpolation_function
322    from anuga.fit_interpolate.interpolate import Interpolation_function
323
324    if not spatial:
325        vertex_coordinates = triangles = interpolation_points = None         
326
327
328    return Interpolation_function(time,
329                                  quantities,
330                                  quantity_names,
331                                  vertex_coordinates,
332                                  triangles,
333                                  interpolation_points,
334                                  time_thinning=time_thinning,
335                                  verbose=verbose)
336
337
338
339
340
341def multiple_replace(text, dictionary):
342    """Multiple replace of words in text
343
344    text:       String to be modified
345    dictionary: Mapping of words that are to be substituted
346
347    Python Cookbook 3.14 page 88 and page 90
348    """
349
350    import re
351   
352    #Create a regular expression from all of the dictionary keys
353    #matching only entire words
354    regex = re.compile(r'\b'+ \
355                       r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \
356                       r'\b' )
357
358    #For each match, lookup the corresponding value in the dictionary
359    return regex.sub(lambda match: dictionary[match.group(0)], text)
360
361
362
363
364def apply_expression_to_dictionary(expression, dictionary):#dictionary):
365    """Apply arbitrary expression to values of dictionary
366
367    Given an expression in terms of the keys, replace key by the
368    corresponding values and evaluate.
369   
370
371    expression: Arbitrary, e.g. arithmetric, expression relating keys
372                from dictionary.
373
374    dictionary: Mapping between symbols in expression and objects that
375                will be evaluated by expression.
376                Values in dictionary must support operators given in
377                expression e.g. by overloading
378
379    due to a limitation with Numeric, this can not evaluate 0/0
380    In general, the user can fix by adding 1e-30 to the numerator.
381    SciPy core can handle this situation.
382    """
383
384    import types
385    import re
386
387    assert isinstance(expression, basestring)
388    assert type(dictionary) == types.DictType
389
390    #Convert dictionary values to textual representations suitable for eval
391    D = {}
392    for key in dictionary:
393        D[key] = 'dictionary["%s"]' %key
394
395    #Perform substitution of variables   
396    expression = multiple_replace(expression, D)
397
398    #Evaluate and return
399    try:
400        return eval(expression)
401    except NameError, e:
402        msg = 'Expression "%s" could not be evaluated: %s' %(expression, e)
403        raise NameError, msg
404    except ValueError, e:
405        msg = 'Expression "%s" could not be evaluated: %s' %(expression, e)
406        raise ValueError, msg
407   
408
409
410
411####################################
412####OBSOLETE STUFF
413
414
415def angle(v1, v2):
416    """Temporary Interface to new location"""
417
418    import anuga.utilities.numerical_tools as NT       
419   
420    msg = 'angle has moved from util.py.  '
421    msg += 'Please use "from anuga.utilities.numerical_tools import angle"'
422    warn(msg, DeprecationWarning) 
423
424    return NT.angle(v1, v2)
425   
426def anglediff(v0, v1):
427    """Temporary Interface to new location"""
428
429    import anuga.utilities.numerical_tools as NT
430   
431    msg = 'anglediff has moved from util.py.  '
432    msg += 'Please use "from anuga.utilities.numerical_tools import anglediff"'
433    warn(msg, DeprecationWarning) 
434
435    return NT.anglediff(v0, v1)   
436
437   
438def mean(x):
439    """Temporary Interface to new location"""
440
441    import anuga.utilities.numerical_tools as NT   
442   
443    msg = 'mean has moved from util.py.  '
444    msg += 'Please use "from anuga.utilities.numerical_tools import mean"'
445    warn(msg, DeprecationWarning) 
446
447    return NT.mean(x)   
448
449def point_on_line(*args, **kwargs):
450    """Temporary Interface to new location"""
451
452    msg = 'point_on_line has moved from util.py.  '
453    msg += 'Please use "from anuga.utilities.polygon import point_on_line"'
454    warn(msg, DeprecationWarning) 
455
456    return utilities.polygon.point_on_line(*args, **kwargs)     
457   
458def inside_polygon(*args, **kwargs):
459    """Temporary Interface to new location"""
460
461    print 'inside_polygon has moved from util.py.  ',
462    print 'Please use "from anuga.utilities.polygon import inside_polygon"'
463
464    return utilities.polygon.inside_polygon(*args, **kwargs)   
465   
466def outside_polygon(*args, **kwargs):
467    """Temporary Interface to new location"""
468
469    print 'outside_polygon has moved from util.py.  ',
470    print 'Please use "from anuga.utilities.polygon import outside_polygon"'
471
472    return utilities.polygon.outside_polygon(*args, **kwargs)   
473
474
475def separate_points_by_polygon(*args, **kwargs):
476    """Temporary Interface to new location"""
477
478    print 'separate_points_by_polygon has moved from util.py.  ',
479    print 'Please use "from anuga.utilities.polygon import separate_points_by_polygon"'
480
481    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)   
482
483
484
485def read_polygon(*args, **kwargs):
486    """Temporary Interface to new location"""
487
488    print 'read_polygon has moved from util.py.  ',
489    print 'Please use "from anuga.utilities.polygon import read_polygon"'
490
491    return utilities.polygon.read_polygon(*args, **kwargs)   
492
493
494def populate_polygon(*args, **kwargs):
495    """Temporary Interface to new location"""
496
497    print 'populate_polygon has moved from util.py.  ',
498    print 'Please use "from anuga.utilities.polygon import populate_polygon"'
499
500    return utilities.polygon.populate_polygon(*args, **kwargs)   
501
502##################### end of obsolete stuff ? ############
503
504def start_screen_catcher(dirname, myid=0, numprocs=1):
505    """Used to store screen output and errors to file, if run on multiple
506    processes eachprocessor will have its own output and error file.
507    """
508
509    dirname = dirname
510    screen_output_name = dirname + "screen_output_%d_%d.txt" %(myid,numprocs)
511    screen_error_name = dirname + "screen_error_%d_%d.txt" %(myid,numprocs)
512
513    #used to catch screen output to file
514    sys.stdout = Screen_Catcher(screen_output_name)
515    sys.stderr = Screen_Catcher(screen_error_name)
516
517class Screen_Catcher:
518    """this simply catches the screen output and stores it to file defined by
519    start_screen_catcher (above)
520    """
521   
522    def __init__(self, filename):
523        self.filename = filename
524        if exists(self.filename)is True:
525            remove(self.filename)
526            print'Old existing file "%s" has been deleted' %(self.filename)
527
528    def write(self, stuff):
529        fid = open(self.filename, 'a')
530        fid.write(stuff)
531
532def get_version_info():
533    """gets the version number of the SVN
534    NOTE: doesn't work on all systems eg GA's cyclone (64 bit linux cluster)
535    """
536   
537    import os, sys
538
539    # Create dummy info
540    info = 'Revision: Version info could not be obtained.'
541    info += 'A command line version of svn and access to the '
542    info += 'repository is necessary and the output must '
543    info += 'contain a line starting with "Revision:"'
544
545    try:
546        fid = os.popen('svn info')
547    except:
548        msg = 'svn is not recognised'
549        warn(msg, UserWarning)
550    else:   
551        lines = fid.readlines()
552        fid.close()
553        for line in lines:
554            if line.startswith('Revision:'):
555                info = line
556                break
557       
558    return info
559   
560def sww2timeseries(swwfiles,
561                   gauge_filename,
562                   production_dirs,
563                   report = None,
564                   reportname = None,
565                   plot_quantity = None,
566                   surface = None,
567                   time_min = None,
568                   time_max = None,
569                   title_on = None,
570                   verbose = False):
571   
572    """ Read sww file and plot the time series for the
573    prescribed quantities at defined gauge locations and
574    prescribed time range.
575
576    Input variables:
577
578    swwfiles        - dictionary of sww files with label_ids (used in
579                      generating latex output. It will be part of
580                      the directory name of file_loc (typically the timestamp).
581                      Helps to differentiate latex files for different simulations
582                      for a particular scenario. 
583                    - assume that all conserved quantities have been stored
584                    - assume each sww file has been simulated with same timestep
585   
586    gauge_filename  - name of file containing gauge data
587                        - easting, northing, name , elevation?
588                    - OR (this is not yet done)
589                        - structure which can be converted to a Numeric array,
590                          such as a geospatial data object
591                     
592    report          - if True, then write figures to report_figures directory in
593                      relevant production directory
594                    - if False, figures are already stored with sww file
595                    - default to False
596
597    reportname      - name for report if wishing to generate report
598   
599    plot_quantity   - list containing quantities to plot, they must
600                      be the name of an existing quantity or one of
601                      the following possibilities
602                    - possibilities:
603                        - stage; 'stage'
604                        - depth; 'depth'
605                        - speed; calculated as absolute momentum
606                         (pointwise) divided by depth; 'speed'
607                        - bearing; calculated as the angle of the momentum
608                          vector (xmomentum, ymomentum) from the North; 'bearing'
609                        - absolute momentum; calculated as
610                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
611                        - x momentum; 'xmomentum'
612                        - y momentum; 'ymomentum'
613                    - default will be ['stage', 'speed', 'bearing']
614
615    surface          - if True, then generate solution surface with 3d plot
616                       and save to current working directory
617                     - default = False
618   
619    time_min         - beginning of user defined time range for plotting purposes
620                        - default will be first available time found in swwfile
621                       
622    time_max         - end of user defined time range for plotting purposes
623                        - default will be last available time found in swwfile
624                       
625    title_on        - if True, export standard graphics with title
626                    - if False, export standard graphics without title
627
628
629                     
630    Output:
631   
632    - time series data stored in .csv for later use if required.
633      Name = gauges_timeseries followed by gauge name
634    - latex file will be generated in same directory as where script is
635      run (usually production scenario directory.
636      Name = latexoutputlabel_id.tex
637
638    Other important information:
639   
640    It is assumed that the used has stored all the conserved quantities
641    and elevation during the scenario run, i.e.
642    ['stage', 'elevation', 'xmomentum', 'ymomentum']
643    If this has not occurred then sww2timeseries will not work.
644                     
645    """
646
647   
648    k = _sww2timeseries(swwfiles,
649                        gauge_filename,
650                        production_dirs,
651                        report,
652                        reportname,
653                        plot_quantity,
654                        surface,
655                        time_min,
656                        time_max,
657                        title_on,
658                        verbose)
659
660    return k
661
662def _sww2timeseries(swwfiles,
663                    gauge_filename,
664                    production_dirs,
665                    report = None,
666                    reportname = None,
667                    plot_quantity = None,
668                    surface = None,
669                    time_min = None,
670                    time_max = None,
671                    title_on = None,
672                    verbose = False):   
673       
674    assert type(gauge_filename) == type(''),\
675           'Gauge filename must be a string'
676   
677    try:
678        fid = open(gauge_filename)
679    except Exception, e:
680        msg = 'File "%s" could not be opened: Error="%s"'\
681                  %(gauge_filename, e)
682        raise msg
683
684    if report is None:
685        report = False
686       
687    if plot_quantity is None:
688        plot_quantity = ['depth', 'speed', 'bearing']
689    else:
690        assert type(plot_quantity) == list,\
691               'plot_quantity must be a list'
692        check_list(plot_quantity)
693
694    if surface is None:
695        surface = False
696       
697    if title_on is None:
698        title_on = True
699   
700    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
701    gauges, locations, elev = get_gauges_from_file(gauge_filename)
702
703    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
704
705    file_loc = []
706    f_list = []
707    label_id = []
708    leg_label = []
709    themaxT = 0.0
710    theminT = 0.0
711   
712    for swwfile in swwfiles.keys():
713
714        try:
715            fid = open(swwfile)
716        except Exception, e:
717            msg = 'File "%s" could not be opened: Error="%s"'\
718                  %(swwfile, e)
719            raise msg
720
721        f = file_function(swwfile,
722                          quantities = sww_quantity,
723                          interpolation_points = gauges,
724                          verbose = True,
725                          use_cache = True)
726
727        # determine which gauges are contained in sww file
728        count = 0
729        gauge_index = []
730        print 'swwfile', swwfile
731        for k, g in enumerate(gauges):
732            if f(0.0, point_id = k)[2] > 1.0e6:
733                count += 1
734                if count == 1: print 'Gauges not contained here:'
735                print locations[k]
736            else:
737                gauge_index.append(k)
738
739        if len(gauge_index) > 0:
740            print 'Gauges contained here: \n',
741        else:
742            print 'No gauges contained here. \n'
743        for i in range(len(gauge_index)):
744             print locations[gauge_index[i]]
745             
746        index = swwfile.rfind(sep)
747        file_loc.append(swwfile[:index+1])
748        label_id.append(swwfiles[swwfile])
749        leg_label.append(production_dirs[swwfiles[swwfile]])
750       
751        f_list.append(f)
752        maxT = max(f.get_time())
753        minT = min(f.get_time())
754        if maxT > themaxT: themaxT = maxT
755        if minT > theminT: theminT = minT
756
757    if time_min is None:
758        time_min = theminT # min(T)
759    else:
760        if time_min < theminT: # min(T):
761            msg = 'Minimum time entered not correct - please try again'
762            raise Exception, msg
763
764    if time_max is None:
765        time_max = themaxT # max(T)
766    else:
767        if time_max > themaxT: # max(T):
768            msg = 'Maximum time entered not correct - please try again'
769            raise Exception, msg
770
771    if verbose and len(gauge_index) > 0: print 'Inputs OK - going to generate figures'
772
773    if len(gauge_index) <> 0:
774        texfile, elev_output = generate_figures(plot_quantity, file_loc, report, reportname, surface,
775                                                leg_label, f_list, gauges, locations, elev, gauge_index,
776                                                production_dirs, time_min, time_max, title_on, label_id, verbose)
777    else:
778        texfile = ''
779        elev_output = []
780
781    return texfile, elev_output
782                         
783#Fixme - Use geospatial to read this file - it's an xya file
784#Need to include other information into this filename, so xya + Name - required for report
785def get_gauges_from_file(filename):
786    """ Read in gauge information from file
787    """
788    from os import sep, getcwd, access, F_OK, mkdir
789    fid = open(filename)
790    lines = fid.readlines()
791    fid.close()
792   
793    gauges = []
794    gaugelocation = []
795    elev = []
796    line1 = lines[0]
797    line11 = line1.split(',')
798    east_index = len(line11)+1
799    north_index = len(line11)+1
800    name_index = len(line11)+1
801    elev_index = len(line11)+1
802    for i in range(len(line11)):
803        if line11[i].strip('\n').strip(' ').lower() == 'easting': east_index = i
804        if line11[i].strip('\n').strip(' ').lower() == 'northing': north_index = i
805        if line11[i].strip('\n').strip(' ').lower() == 'name': name_index = i
806        if line11[i].strip('\n').strip(' ').lower() == 'elevation': elev_index = i
807
808    for line in lines[1:]:
809        fields = line.split(',')
810        if east_index < len(line11) and north_index < len(line11):
811            gauges.append([float(fields[east_index]), float(fields[north_index])])
812        else:
813            msg = 'WARNING: %s does not contain location information' %(filename)
814            raise Exception, msg
815        if elev_index < len(line11): elev.append(float(fields[elev_index]))
816        if name_index < len(line11):
817            loc = fields[name_index]
818            gaugelocation.append(loc.strip('\n'))
819
820    return gauges, gaugelocation, elev
821
822def check_list(quantity):
823    """ Check that input quantities in quantity list are possible
824    """
825    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
826                    'ymomentum', 'speed', 'bearing', 'elevation']
827
828    import sys
829    if sys.platform == 'win32':
830        for i,j in enumerate(quantity):
831            quantity[i] = quantity[i].lower()
832        p = list(set(quantity).difference(set(all_quantity)))
833        if len(p) <> 0:
834            msg = 'Quantities %s do not exist - please try again' %p
835            raise Exception, msg
836       
837    return 
838
839def calc_bearing(uh, vh):
840    """ Calculate velocity bearing from North
841    """
842    from math import atan, degrees
843   
844    angle = degrees(atan(vh/(uh+1.e-15)))
845    if (0 < angle < 90.0):
846        if vh > 0:
847            bearing = 90.0 - abs(angle)
848        if vh < 0:
849            bearing = 270.0 - abs(angle)
850    if (-90 < angle < 0):
851        if vh < 0:
852            bearing = 90.0 - abs(angle)
853        if vh > 0:
854            bearing = 270.0 - abs(angle)
855    if angle == 0: bearing = 0.0
856           
857    return bearing
858
859def generate_figures(plot_quantity, file_loc, report, reportname, surface,
860                     leg_label, f_list, gauges, locations, elev, gauge_index,
861                     production_dirs, time_min, time_max, title_on, label_id,
862                     verbose):
863    """ Generate figures based on required quantities and gauges for each sww file
864    """
865    from math import sqrt, atan, degrees
866    from Numeric import ones, allclose, zeros, Float, ravel
867    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
868    from pylab import ion, hold, plot, axis, figure, legend, savefig, \
869         xlabel, ylabel, title, close, subplot
870
871    import pylab as p1
872    if surface is True:
873        import mpl3d.mplot3d as p3
874       
875    if report == True:   
876        texdir = getcwd()+sep+'report'+sep
877        if access(texdir,F_OK) == 0:
878            mkdir (texdir)
879        if len(label_id) == 1:
880            label_id1 = label_id[0].replace(sep,'')
881            label_id2 = label_id1.replace('_','')
882            texfile = texdir+reportname+'%s' %(label_id2)
883            texfile2 = reportname+'%s' %(label_id2)
884            texfilename = texfile + '.tex'
885            if verbose: print '\n Latex output printed to %s \n' %texfilename
886            fid = open(texfilename, 'w')
887        else:
888            texfile = texdir+reportname
889            texfile2 = reportname
890            texfilename = texfile + '.tex' 
891            if verbose: print '\n Latex output printed to %s \n' %texfilename
892            fid = open(texfilename, 'w')
893    else:
894        texfile = ''
895        texfile2 = ''
896
897    p = len(f_list)
898    n = []
899    n0 = 0
900    for i in range(len(f_list)):
901        n.append(len(f_list[i].get_time()))
902        if n[i] > n0: n0 = n[i] 
903    n0 = int(n0)
904    m = len(locations)
905    model_time = zeros((n0,m,p), Float) 
906    stages = zeros((n0,m,p), Float)
907    elevations = zeros((n0,m,p), Float) 
908    momenta = zeros((n0,m,p), Float)
909    xmom = zeros((n0,m,p), Float)
910    ymom = zeros((n0,m,p), Float)
911    speed = zeros((n0,m,p), Float)
912    bearings = zeros((n0,m,p), Float)
913    depths = zeros((n0,m,p), Float)
914    eastings = zeros((n0,m,p), Float)
915    min_stages = []
916    max_stages = []
917    max_momentums = []
918    max_speeds = []
919    c = 0
920    model_time_plot3d = zeros((n0,m), Float)
921    stages_plot3d = zeros((n0,m), Float)
922    eastings_plot3d = zeros((n0,m),Float)
923    ##### loop over each swwfile #####
924    for j, f in enumerate(f_list):
925        if verbose: print 'swwfile %d of %d' %(j, len(f_list))
926        comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv'
927        fid_compare = open(comparefile, 'w')
928        ##### loop over each gauge #####
929        for k in gauge_index:
930            g = gauges[k]
931            if verbose: print 'Gauge %d of %d' %(k, len(gauges))
932            min_stage = 10
933            max_stage = 0
934            max_momentum = 0
935            max_speed = 0   
936            gaugeloc = locations[k]
937            thisfile = file_loc[j]+sep+'gauges_time_series'+'_'+gaugeloc+'.csv'
938            fid_out = open(thisfile, 'w')
939            s = 'Time, Stage, Momentum, Speed, Elevation \n'
940            fid_out.write(s)
941           
942            #### generate quantities #######
943            for i, t in enumerate(f.get_time()):
944                if time_min <= t <= time_max:
945                    w = f(t, point_id = k)[0]
946                    z = f(t, point_id = k)[1]
947                    uh = f(t, point_id = k)[2]
948                    vh = f(t, point_id = k)[3]
949                    depth = w-z     
950                    m = sqrt(uh*uh + vh*vh)   
951                    if depth < 0.001:
952                        vel = 0.0
953                    else:
954                        vel = m / (depth + 1.e-30) 
955                    bearing = calc_bearing(uh, vh)
956                    model_time[i,k,j] = t/60.0
957                    stages[i,k,j] = w
958                    elevations[i,k,j] = z
959                    xmom[i,k,j] = uh
960                    ymom[i,k,j] = vh
961                    momenta[i,k,j] = m
962                    speed[i,k,j] = vel
963                    bearings[i,k,j] = bearing
964                    depths[i,k,j] = depth
965                    thisgauge = gauges[k]
966                    eastings[i,k,j] = thisgauge[0]
967                    s = '%.2f, %.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel, z)
968                    fid_out.write(s)
969                    if t/60.0 <= 13920: tindex = i
970                    if w > max_stage: max_stage = w
971                    if w < min_stage: min_stage = w
972                    if m > max_momentum: max_momentum = m
973                    if vel > max_speed: max_speed = vel
974                   
975                   
976            s = '%.2f, %.2f, %.2f, %.2f, %s\n' %(max_stage, min_stage, z, thisgauge[0], leg_label[j])
977            fid_compare.write(s)
978            max_stages.append(max_stage)
979            min_stages.append(min_stage)
980            max_momentums.append(max_momentum)
981            max_speeds.append(max_speed)     
982            #### finished generating quantities for each swwfile #####
983       
984        model_time_plot3d[:,:] = model_time[:,:,j]
985        stages_plot3d[:,:] = stages[:,:,j]
986        eastings_plot3d[:,] = eastings[:,:,j]
987           
988        if surface ==  True:
989            print 'Printing surface figure'
990            for i in range(2):
991                fig = p1.figure(10)
992                ax = p3.Axes3D(fig)
993                if len(gauges) > 80:
994                    ax.plot_surface(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
995                else:
996                    #ax.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
997                    ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j]))
998                ax.set_xlabel('time')
999                ax.set_ylabel('x')
1000                ax.set_zlabel('stage')
1001                fig.add_axes(ax)
1002                p1.show()
1003                surfacefig = 'solution_surface%s' %leg_label[j]
1004                p1.savefig(surfacefig)
1005                p1.close()
1006           
1007    #### finished generating quantities for all swwfiles #####
1008
1009    # x profile for given time
1010    if surface == True:
1011        figure(11)
1012        plot(eastings[tindex,:,j],stages[tindex,:,j])
1013        xlabel('x')
1014        ylabel('stage')
1015        profilefig = 'solution_xprofile' 
1016        savefig('profilefig')
1017               
1018    #stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
1019    #stage_axis = axis([time_min/60.0, time_max/60.0, -3.0, 3.0])   
1020    #vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
1021    #mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1]) 
1022   
1023    cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
1024    nn = len(plot_quantity)
1025    no_cols = 2
1026    elev_output = []
1027    if len(label_id) > 1: graphname_report = []
1028    for k in gauge_index:
1029        g = gauges[k]
1030        count1 = 0
1031        if report == True and len(label_id) > 1:
1032            s = '\\begin{figure}[hbt] \n \\centering \n \\begin{tabular}{cc} \n'
1033            fid.write(s)
1034        if len(label_id) > 1: graphname_report = []
1035        #### generate figures for each gauge ####
1036        for j, f in enumerate(f_list):
1037            ion()
1038            hold(True)
1039            count = 0
1040            where1 = 0
1041            where2 = 0
1042            word_quantity = ''
1043            if report == True and len(label_id) == 1:
1044                s = '\\begin{figure}[hbt] \n \\centering \n \\begin{tabular}{cc} \n'
1045                fid.write(s)
1046               
1047            for which_quantity in plot_quantity:
1048                count += 1
1049                where1 += 1
1050                figure(count, frameon = False)
1051                if which_quantity == 'depth':
1052                    plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j])
1053                    units = 'm'
1054                if which_quantity == 'stage':
1055                    plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j])
1056                    #axis(stage_axis)
1057                    units = 'm'
1058                if which_quantity == 'momentum':
1059                    plot(model_time[0:n[j]-1,k,j], momenta[0:n[j]-1,k,j], '-', c = cstr[j])
1060                    #axis(mom_axis)
1061                    units = 'm^2 / sec'
1062                if which_quantity == 'xmomentum':
1063                    plot(model_time[0:n[j]-1,k,j], xmom[0:n[j]-1,k,j], '-', c = cstr[j])
1064                    #axis(mom_axis)
1065                    units = 'm^2 / sec'
1066                if which_quantity == 'ymomentum':
1067                    plot(model_time[0:n[j]-1,k,j], ymom[0:n[j]-1,k,j], '-', c = cstr[j])
1068                    #axis(mom_axis)
1069                    units = 'm^2 / sec'
1070                if which_quantity == 'speed':
1071                    plot(model_time[0:n[j]-1,k,j], speed[0:n[j]-1,k,j], '-', c = cstr[j])
1072                    #axis(vel_axis)
1073                    units = 'm / sec'
1074                if which_quantity == 'bearing':
1075                    due_east = 90.0*ones(shape(model_time[0:n[j]-1,k,j],Float))
1076                    due_west = 270.0*ones(shape(model_time[0:n[j]-1,k,j],Float))
1077                    plot(model_time[0:n[j]-1,k,j], bearings, '-', 
1078                         model_time[0:n[j]-1,k,j], due_west, '-.', 
1079                         model_time[0:n[j]-1,k,j], due_east, '-.')
1080                    units = 'degrees from North'
1081                    ax = axis([time_min, time_max, 0.0, 360.0])
1082                    legend(('Bearing','West','East'))
1083
1084                xlabel('time (mins)')
1085                ylabel('%s (%s)' %(which_quantity, units))
1086                if len(label_id) > 1: legend((leg_label),loc='upper right')
1087
1088                gaugeloc1 = gaugeloc.replace(' ','')
1089                #gaugeloc2 = gaugeloc1.replace('_','')
1090                gaugeloc2 = locations[k].replace(' ','')
1091                graphname = '%sgauge%s_%s' %(file_loc[j], gaugeloc2, which_quantity)
1092
1093                if report == True and len(label_id) > 1:
1094                    figdir = getcwd()+sep+'report_figures'+sep
1095                    if access(figdir,F_OK) == 0 :
1096                        mkdir (figdir)
1097                    latex_file_loc = figdir.replace(sep,altsep) 
1098                    graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory   
1099                    graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file
1100                    graphname_report.append(graphname_report_input)
1101                   
1102                    savefig(graphname_latex) # save figures in production directory for report generation
1103
1104                if report == True:
1105
1106                    figdir = getcwd()+sep+'report_figures'+sep
1107                    if access(figdir,F_OK) == 0 :
1108                        mkdir (figdir)
1109                    latex_file_loc = figdir.replace(sep,altsep)   
1110
1111                    if len(label_id) == 1: 
1112                        graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory 
1113                        graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file
1114                        s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png')
1115                        fid.write(s)
1116                        if where1 % 2 == 0:
1117                            s = '\\\\ \n'
1118                            where1 = 0
1119                        else:
1120                            s = '& \n'
1121                        fid.write(s)
1122                        savefig(graphname_latex)
1123               
1124                if title_on == True:
1125                    title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc2))
1126
1127                savefig(graphname) # save figures with sww file
1128
1129            if report == True and len(label_id) == 1:
1130                for i in range(nn-1):
1131                    if nn > 2:
1132                        word_quantity += plot_quantity[i] + ', '
1133                    else:
1134                        word_quantity += plot_quantity[i]
1135                   
1136                word_quantity += ' and ' + plot_quantity[nn-1]
1137                caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' '))
1138                if elev[k] == 0.0:
1139                    caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j])
1140                    east = gauges[0]
1141                    north = gauges[1]
1142                    elev_output.append([locations[k],east,north,elevations[0,k,j]])
1143                label = '%sgauge%s' %(label_id2, gaugeloc2)
1144                s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
1145                fid.write(s)
1146                c += 1
1147                if c % 6 == 0: fid.write('\\clearpage \n')
1148                savefig(graphname_latex)               
1149               
1150        if report == True and len(label_id) > 1:
1151            for i in range(nn-1):
1152                if nn > 2:
1153                    word_quantity += plot_quantity[i] + ', '
1154                else:
1155                    word_quantity += plot_quantity[i]
1156                where1 = 0
1157                count1 += 1
1158                index = j*len(plot_quantity)
1159                for which_quantity in plot_quantity:
1160                    where1 += 1
1161                    s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report[index], '.png')
1162                    index += 1
1163                    fid.write(s)
1164                    if where1 % 2 == 0:
1165                        s = '\\\\ \n'
1166                        where1 = 0
1167                    else:
1168                        s = '& \n'
1169                    fid.write(s)
1170            word_quantity += ' and ' + plot_quantity[nn-1]           
1171            label = 'gauge%s' %(gaugeloc2) 
1172            caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k])
1173            if elev[k] == 0.0:
1174                    caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j])
1175                    thisgauge = gauges[k]
1176                    east = thisgauge[0]
1177                    north = thisgauge[1]
1178                    elev_output.append([locations[k],east,north,elevations[0,k,j]])
1179                   
1180            s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
1181            fid.write(s)
1182            c += 1
1183            if c % 6 == 0: fid.write('\\clearpage \n')         
1184           
1185            #### finished generating figures ###
1186
1187        close('all')
1188   
1189    return texfile2, elev_output
1190
1191# FIXME (DSG): Add unit test, make general, not just 2 files,
1192# but any number of files.
1193def copy_code_files(dir_name, filename1, filename2):
1194    """Copies "filename1" and "filename2" to "dir_name". Very useful for
1195    information management """
1196
1197    if access(dir_name,F_OK) == 0:
1198        print 'Make directory %s' %dir_name
1199        mkdir (dir_name,0777)
1200    copy(filename1, dir_name + sep + basename(filename1))
1201    copy(filename2, dir_name + sep + basename(filename2))
1202#    copy (__file__, project.output_run_time_dir + basename(__file__))
1203    print 'Files %s and %s copied' %(filename1, filename2)
1204
1205
1206def add_directories(root_directory, directories):
1207    """
1208    Add the first directory in directories to root_directory.
1209    Then add the second
1210    direcotory to the first directory and so on.
1211
1212    Return the path of the final directory.
1213
1214    This is handy for specifying and creating a directory where data will go.
1215    """
1216    dir = root_directory
1217    for new_dir in directories:
1218        dir = os.path.join(dir, new_dir)
1219        if not access(dir,F_OK):
1220            mkdir (dir)
1221    return dir
1222
1223
1224
Note: See TracBrowser for help on using the repository browser.