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

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

(1) updates to Dampier script based on Perth script (2) minor updates to Onslow report

File size: 45.2 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(dir_name, 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    dir_name = dir_name
510    if access(dir_name,F_OK) == 0:
511        print 'Make directory %s' %dir_name
512        mkdir (dir_name,0777)
513    screen_output_name = dir_name + "screen_output_%d_%d.txt" %(myid,numprocs)
514    screen_error_name = dir_name + "screen_error_%d_%d.txt" %(myid,numprocs)
515   
516    #used to catch screen output to file
517    sys.stdout = Screen_Catcher(screen_output_name)
518    sys.stderr = Screen_Catcher(screen_error_name)
519
520class Screen_Catcher:
521    """this simply catches the screen output and stores it to file defined by
522    start_screen_catcher (above)
523    """
524   
525    def __init__(self, filename):
526        self.filename = filename
527 
528        if exists(self.filename)is True:
529            remove(self.filename)
530            print'Old existing file "%s" has been deleted' %(self.filename)
531
532    def write(self, stuff):
533        fid = open(self.filename, 'a')
534        fid.write(stuff)
535
536def get_version_info():
537    """gets the version number of the SVN
538    NOTE: doesn't work on all systems eg GA's cyclone (64 bit linux cluster)
539    """
540   
541    import os, sys
542
543    # Create dummy info
544    info = 'Revision: Version info could not be obtained.'
545    info += 'A command line version of svn and access to the '
546    info += 'repository is necessary and the output must '
547    info += 'contain a line starting with "Revision:"'
548
549    try:
550        fid = os.popen('svn info')
551    except:
552        msg = 'svn is not recognised'
553        warn(msg, UserWarning)
554    else:   
555        lines = fid.readlines()
556        fid.close()
557        for line in lines:
558            if line.startswith('Revision:'):
559                info = line
560                break
561       
562    return info
563   
564def sww2timeseries(swwfiles,
565                   gauge_filename,
566                   production_dirs,
567                   report = None,
568                   reportname = None,
569                   plot_quantity = None,
570                   surface = None,
571                   time_min = None,
572                   time_max = None,
573                   title_on = None,
574                   verbose = False):
575   
576    """ Read sww file and plot the time series for the
577    prescribed quantities at defined gauge locations and
578    prescribed time range.
579
580    Input variables:
581
582    swwfiles        - dictionary of sww files with label_ids (used in
583                      generating latex output. It will be part of
584                      the directory name of file_loc (typically the timestamp).
585                      Helps to differentiate latex files for different simulations
586                      for a particular scenario. 
587                    - assume that all conserved quantities have been stored
588                    - assume each sww file has been simulated with same timestep
589   
590    gauge_filename  - name of file containing gauge data
591                        - easting, northing, name , elevation?
592                    - OR (this is not yet done)
593                        - structure which can be converted to a Numeric array,
594                          such as a geospatial data object
595                     
596    report          - if True, then write figures to report_figures directory in
597                      relevant production directory
598                    - if False, figures are already stored with sww file
599                    - default to False
600
601    reportname      - name for report if wishing to generate report
602   
603    plot_quantity   - list containing quantities to plot, they must
604                      be the name of an existing quantity or one of
605                      the following possibilities
606                    - possibilities:
607                        - stage; 'stage'
608                        - depth; 'depth'
609                        - speed; calculated as absolute momentum
610                         (pointwise) divided by depth; 'speed'
611                        - bearing; calculated as the angle of the momentum
612                          vector (xmomentum, ymomentum) from the North; 'bearing'
613                        - absolute momentum; calculated as
614                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
615                        - x momentum; 'xmomentum'
616                        - y momentum; 'ymomentum'
617                    - default will be ['stage', 'speed', 'bearing']
618
619    surface          - if True, then generate solution surface with 3d plot
620                       and save to current working directory
621                     - default = False
622   
623    time_min         - beginning of user defined time range for plotting purposes
624                        - default will be first available time found in swwfile
625                       
626    time_max         - end of user defined time range for plotting purposes
627                        - default will be last available time found in swwfile
628                       
629    title_on        - if True, export standard graphics with title
630                    - if False, export standard graphics without title
631
632
633                     
634    Output:
635   
636    - time series data stored in .csv for later use if required.
637      Name = gauges_timeseries followed by gauge name
638    - latex file will be generated in same directory as where script is
639      run (usually production scenario directory.
640      Name = latexoutputlabel_id.tex
641
642    Other important information:
643   
644    It is assumed that the used has stored all the conserved quantities
645    and elevation during the scenario run, i.e.
646    ['stage', 'elevation', 'xmomentum', 'ymomentum']
647    If this has not occurred then sww2timeseries will not work.
648                     
649    """
650
651   
652    k = _sww2timeseries(swwfiles,
653                        gauge_filename,
654                        production_dirs,
655                        report,
656                        reportname,
657                        plot_quantity,
658                        surface,
659                        time_min,
660                        time_max,
661                        title_on,
662                        verbose)
663
664    return k
665
666def _sww2timeseries(swwfiles,
667                    gauge_filename,
668                    production_dirs,
669                    report = None,
670                    reportname = None,
671                    plot_quantity = None,
672                    surface = None,
673                    time_min = None,
674                    time_max = None,
675                    title_on = None,
676                    verbose = False):   
677       
678    assert type(gauge_filename) == type(''),\
679           'Gauge filename must be a string'
680   
681    try:
682        fid = open(gauge_filename)
683    except Exception, e:
684        msg = 'File "%s" could not be opened: Error="%s"'\
685                  %(gauge_filename, e)
686        raise msg
687
688    if report is None:
689        report = False
690       
691    if plot_quantity is None:
692        plot_quantity = ['depth', 'speed', 'bearing']
693    else:
694        assert type(plot_quantity) == list,\
695               'plot_quantity must be a list'
696        check_list(plot_quantity)
697
698    if surface is None:
699        surface = False
700       
701    if title_on is None:
702        title_on = True
703   
704    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
705    gauges, locations, elev = get_gauges_from_file(gauge_filename)
706
707    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
708
709    file_loc = []
710    f_list = []
711    label_id = []
712    leg_label = []
713    themaxT = 0.0
714    theminT = 0.0
715   
716    for swwfile in swwfiles.keys():
717
718        try:
719            fid = open(swwfile)
720        except Exception, e:
721            msg = 'File "%s" could not be opened: Error="%s"'\
722                  %(swwfile, e)
723            raise msg
724
725        f = file_function(swwfile,
726                          quantities = sww_quantity,
727                          interpolation_points = gauges,
728                          verbose = True,
729                          use_cache = True)
730
731        # determine which gauges are contained in sww file
732        count = 0
733        gauge_index = []
734        print 'swwfile', swwfile
735        for k, g in enumerate(gauges):
736            if f(0.0, point_id = k)[2] > 1.0e6:
737                count += 1
738                if count == 1: print 'Gauges not contained here:'
739                print locations[k]
740            else:
741                gauge_index.append(k)
742
743        if len(gauge_index) > 0:
744            print 'Gauges contained here: \n',
745        else:
746            print 'No gauges contained here. \n'
747        for i in range(len(gauge_index)):
748             print locations[gauge_index[i]]
749             
750        index = swwfile.rfind(sep)
751        file_loc.append(swwfile[:index+1])
752        label_id.append(swwfiles[swwfile])
753        leg_label.append(production_dirs[swwfiles[swwfile]])
754       
755        f_list.append(f)
756        maxT = max(f.get_time())
757        minT = min(f.get_time())
758        if maxT > themaxT: themaxT = maxT
759        if minT > theminT: theminT = minT
760
761    if time_min is None:
762        time_min = theminT # min(T)
763    else:
764        if time_min < theminT: # min(T):
765            msg = 'Minimum time entered not correct - please try again'
766            raise Exception, msg
767
768    if time_max is None:
769        time_max = themaxT # max(T)
770    else:
771        if time_max > themaxT: # max(T):
772            msg = 'Maximum time entered not correct - please try again'
773            raise Exception, msg
774
775    if verbose and len(gauge_index) > 0: print 'Inputs OK - going to generate figures'
776
777    if len(gauge_index) <> 0:
778        texfile, elev_output = generate_figures(plot_quantity, file_loc, report, reportname, surface,
779                                                leg_label, f_list, gauges, locations, elev, gauge_index,
780                                                production_dirs, time_min, time_max, title_on, label_id, verbose)
781    else:
782        texfile = ''
783        elev_output = []
784
785    return texfile, elev_output
786                         
787#Fixme - Use geospatial to read this file - it's an xya file
788#Need to include other information into this filename, so xya + Name - required for report
789def get_gauges_from_file(filename):
790    """ Read in gauge information from file
791    """
792    from os import sep, getcwd, access, F_OK, mkdir
793    fid = open(filename)
794    lines = fid.readlines()
795    fid.close()
796   
797    gauges = []
798    gaugelocation = []
799    elev = []
800    line1 = lines[0]
801    line11 = line1.split(',')
802    east_index = len(line11)+1
803    north_index = len(line11)+1
804    name_index = len(line11)+1
805    elev_index = len(line11)+1
806    for i in range(len(line11)):
807        if line11[i].strip('\n').strip(' ').lower() == 'easting': east_index = i
808        if line11[i].strip('\n').strip(' ').lower() == 'northing': north_index = i
809        if line11[i].strip('\n').strip(' ').lower() == 'name': name_index = i
810        if line11[i].strip('\n').strip(' ').lower() == 'elevation': elev_index = i
811
812    for line in lines[1:]:
813        fields = line.split(',')
814        if east_index < len(line11) and north_index < len(line11):
815            gauges.append([float(fields[east_index]), float(fields[north_index])])
816        else:
817            msg = 'WARNING: %s does not contain location information' %(filename)
818            raise Exception, msg
819        if elev_index < len(line11): elev.append(float(fields[elev_index]))
820        if name_index < len(line11):
821            loc = fields[name_index]
822            gaugelocation.append(loc.strip('\n'))
823
824    return gauges, gaugelocation, elev
825
826def check_list(quantity):
827    """ Check that input quantities in quantity list are possible
828    """
829    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
830                    'ymomentum', 'speed', 'bearing', 'elevation']
831
832       
833                   
834    import sys
835    #if not sys.version.startswith('2.4'):
836        # Backwards compatibility
837    #   from sets import Set as set
838
839    from sets import Set as set
840           
841    for i,j in enumerate(quantity):
842        quantity[i] = quantity[i].lower()
843    p = list(set(quantity).difference(set(all_quantity)))
844    if len(p) <> 0:
845        msg = 'Quantities %s do not exist - please try again' %p
846        raise Exception, msg
847       
848    return 
849
850def calc_bearing(uh, vh):
851    """ Calculate velocity bearing from North
852    """
853    from math import atan, degrees
854   
855    angle = degrees(atan(vh/(uh+1.e-15)))
856    if (0 < angle < 90.0):
857        if vh > 0:
858            bearing = 90.0 - abs(angle)
859        if vh < 0:
860            bearing = 270.0 - abs(angle)
861    if (-90 < angle < 0):
862        if vh < 0:
863            bearing = 90.0 - abs(angle)
864        if vh > 0:
865            bearing = 270.0 - abs(angle)
866    if angle == 0: bearing = 0.0
867           
868    return bearing
869
870def generate_figures(plot_quantity, file_loc, report, reportname, surface,
871                     leg_label, f_list, gauges, locations, elev, gauge_index,
872                     production_dirs, time_min, time_max, title_on, label_id,
873                     verbose):
874    """ Generate figures based on required quantities and gauges for each sww file
875    """
876    from math import sqrt, atan, degrees
877    from Numeric import ones, allclose, zeros, Float, ravel
878    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
879    from pylab import ion, hold, plot, axis, figure, legend, savefig, \
880         xlabel, ylabel, title, close, subplot
881
882    import pylab as p1
883    if surface is True:
884        import mpl3d.mplot3d as p3
885       
886    if report == True:   
887        texdir = getcwd()+sep+'report'+sep
888        if access(texdir,F_OK) == 0:
889            mkdir (texdir)
890        if len(label_id) == 1:
891            label_id1 = label_id[0].replace(sep,'')
892            label_id2 = label_id1.replace('_','')
893            texfile = texdir+reportname+'%s' %(label_id2)
894            texfile2 = reportname+'%s' %(label_id2)
895            texfilename = texfile + '.tex'
896            if verbose: print '\n Latex output printed to %s \n' %texfilename
897            fid = open(texfilename, 'w')
898        else:
899            texfile = texdir+reportname
900            texfile2 = reportname
901            texfilename = texfile + '.tex' 
902            if verbose: print '\n Latex output printed to %s \n' %texfilename
903            fid = open(texfilename, 'w')
904    else:
905        texfile = ''
906        texfile2 = ''
907
908    p = len(f_list)
909    n = []
910    n0 = 0
911    for i in range(len(f_list)):
912        n.append(len(f_list[i].get_time()))
913        if n[i] > n0: n0 = n[i] 
914    n0 = int(n0)
915    m = len(locations)
916    model_time = zeros((n0,m,p), Float) 
917    stages = zeros((n0,m,p), Float)
918    elevations = zeros((n0,m,p), Float) 
919    momenta = zeros((n0,m,p), Float)
920    xmom = zeros((n0,m,p), Float)
921    ymom = zeros((n0,m,p), Float)
922    speed = zeros((n0,m,p), Float)
923    bearings = zeros((n0,m,p), Float)
924    depths = zeros((n0,m,p), Float)
925    eastings = zeros((n0,m,p), Float)
926    min_stages = []
927    max_stages = []
928    max_momentums = []
929    max_speeds = []
930    max_depths = []
931    model_time_plot3d = zeros((n0,m), Float)
932    stages_plot3d = zeros((n0,m), Float)
933    eastings_plot3d = zeros((n0,m),Float)
934    ##### loop over each swwfile #####
935    for j, f in enumerate(f_list):
936        if verbose: print 'swwfile %d of %d' %(j, len(f_list))
937        comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv'
938        fid_compare = open(comparefile, 'w')
939        ##### loop over each gauge #####
940        for k in gauge_index:
941            g = gauges[k]
942            if verbose: print 'Gauge %d of %d' %(k, len(gauges))
943            min_stage = 10
944            max_stage = 0
945            max_momentum = 0
946            max_speed = 0
947            max_depth = 0
948            gaugeloc = locations[k]
949            thisfile = file_loc[j]+sep+'gauges_time_series'+'_'+gaugeloc+'.csv'
950            fid_out = open(thisfile, 'w')
951            s = 'Time, Stage, Momentum, Speed, Elevation \n'
952            fid_out.write(s)
953           
954            #### generate quantities #######
955            for i, t in enumerate(f.get_time()):
956                if time_min <= t <= time_max:
957                    w = f(t, point_id = k)[0]
958                    z = f(t, point_id = k)[1]
959                    uh = f(t, point_id = k)[2]
960                    vh = f(t, point_id = k)[3]
961                    depth = w-z     
962                    m = sqrt(uh*uh + vh*vh)   
963                    if depth < 0.001:
964                        vel = 0.0
965                    else:
966                        vel = m / (depth + 1.e-30) 
967                    bearing = calc_bearing(uh, vh)
968                    model_time[i,k,j] = t/60.0
969                    stages[i,k,j] = w
970                    elevations[i,k,j] = z
971                    xmom[i,k,j] = uh
972                    ymom[i,k,j] = vh
973                    momenta[i,k,j] = m
974                    speed[i,k,j] = vel
975                    bearings[i,k,j] = bearing
976                    depths[i,k,j] = depth
977                    thisgauge = gauges[k]
978                    eastings[i,k,j] = thisgauge[0]
979                    s = '%.2f, %.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel, z)
980                    fid_out.write(s)
981                    if t/60.0 <= 13920: tindex = i
982                    if w > max_stage: max_stage = w
983                    if w < min_stage: min_stage = w
984                    if m > max_momentum: max_momentum = m
985                    if vel > max_speed: max_speed = vel
986                    if z > 0 and depth > max_depth: max_depth = depth
987                   
988                   
989            s = '%.2f, %.2f, %.2f, %.2f, %s\n' %(max_stage, min_stage, z, thisgauge[0], leg_label[j])
990            fid_compare.write(s)
991            max_stages.append(max_stage)
992            min_stages.append(min_stage)
993            max_momentums.append(max_momentum)
994            max_speeds.append(max_speed)
995            max_depths.append(max_depth)
996            #### finished generating quantities for each swwfile #####
997       
998        model_time_plot3d[:,:] = model_time[:,:,j]
999        stages_plot3d[:,:] = stages[:,:,j]
1000        eastings_plot3d[:,] = eastings[:,:,j]
1001           
1002        if surface ==  True:
1003            print 'Printing surface figure'
1004            for i in range(2):
1005                fig = p1.figure(10)
1006                ax = p3.Axes3D(fig)
1007                if len(gauges) > 80:
1008                    ax.plot_surface(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
1009                else:
1010                    #ax.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
1011                    ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j]))
1012                ax.set_xlabel('time')
1013                ax.set_ylabel('x')
1014                ax.set_zlabel('stage')
1015                fig.add_axes(ax)
1016                p1.show()
1017                surfacefig = 'solution_surface%s' %leg_label[j]
1018                p1.savefig(surfacefig)
1019                p1.close()
1020           
1021    #### finished generating quantities for all swwfiles #####
1022
1023    # x profile for given time
1024    if surface == True:
1025        figure(11)
1026        plot(eastings[tindex,:,j],stages[tindex,:,j])
1027        xlabel('x')
1028        ylabel('stage')
1029        profilefig = 'solution_xprofile' 
1030        savefig('profilefig')
1031
1032    depth_axis = axis([time_min/60.0, time_max/60.0, -0.1, max(max_depths)*1.1])
1033    stage_axis = axis([time_min/60.0, time_max/60.0, min(min_stages), max(max_stages)*1.1])
1034    vel_axis = axis([time_min/60.0, time_max/60.0, min(max_speeds), max(max_speeds)*1.1])
1035    mom_axis = axis([time_min/60.0, time_max/60.0, min(max_momentums), max(max_momentums)*1.1]) 
1036   
1037    cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
1038    nn = len(plot_quantity)
1039    no_cols = 2
1040    elev_output = []
1041    if len(label_id) > 1: graphname_report = []
1042    pp = 1
1043    div = 11.
1044    for k in gauge_index:
1045        g = gauges[k]
1046        count1 = 0
1047        if report == True and len(label_id) > 1:
1048            s = '\\begin{figure}[ht] \n \\centering \n \\begin{tabular}{cc} \n'
1049            fid.write(s)
1050        if len(label_id) > 1: graphname_report = []
1051        #### generate figures for each gauge ####
1052        for j, f in enumerate(f_list):
1053            ion()
1054            hold(True)
1055            count = 0
1056            where1 = 0
1057            where2 = 0
1058            word_quantity = ''
1059            if report == True and len(label_id) == 1:
1060                s = '\\begin{figure}[hbt] \n \\centering \n \\begin{tabular}{cc} \n'
1061                fid.write(s)
1062               
1063            for which_quantity in plot_quantity:
1064                count += 1
1065                where1 += 1
1066                figure(count, frameon = False)
1067                if which_quantity == 'depth':
1068                    plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j])
1069                    units = 'm'
1070                    axis(depth_axis)
1071                if which_quantity == 'stage':
1072                    if elevations[0,k,j] < 0:
1073                        plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j])
1074                        axis(stage_axis)
1075                    else:
1076                        plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j])
1077                        axis(depth_axis)                 
1078                    units = 'm'
1079                if which_quantity == 'momentum':
1080                    plot(model_time[0:n[j]-1,k,j], momenta[0:n[j]-1,k,j], '-', c = cstr[j])
1081                    axis(mom_axis)
1082                    units = 'm^2 / sec'
1083                if which_quantity == 'xmomentum':
1084                    plot(model_time[0:n[j]-1,k,j], xmom[0:n[j]-1,k,j], '-', c = cstr[j])
1085                    axis(mom_axis)
1086                    units = 'm^2 / sec'
1087                if which_quantity == 'ymomentum':
1088                    plot(model_time[0:n[j]-1,k,j], ymom[0:n[j]-1,k,j], '-', c = cstr[j])
1089                    axis(mom_axis)
1090                    units = 'm^2 / sec'
1091                if which_quantity == 'speed':
1092                    plot(model_time[0:n[j]-1,k,j], speed[0:n[j]-1,k,j], '-', c = cstr[j])
1093                    axis(vel_axis)
1094                    units = 'm / sec'
1095                if which_quantity == 'bearing':
1096                    due_east = 90.0*ones(shape(model_time[0:n[j]-1,k,j],Float))
1097                    due_west = 270.0*ones(shape(model_time[0:n[j]-1,k,j],Float))
1098                    plot(model_time[0:n[j]-1,k,j], bearings, '-', 
1099                         model_time[0:n[j]-1,k,j], due_west, '-.', 
1100                         model_time[0:n[j]-1,k,j], due_east, '-.')
1101                    units = 'degrees from North'
1102                    ax = axis([time_min, time_max, 0.0, 360.0])
1103                    legend(('Bearing','West','East'))
1104
1105                xlabel('time (mins)')
1106                if which_quantity == 'stage' and elevations[0:n[j]-1,k,j] > 0:
1107                    ylabel('%s (%s)' %('depth', units))
1108                else:
1109                    ylabel('%s (%s)' %(which_quantity, units))
1110                if len(label_id) > 1: legend((leg_label),loc='upper right')
1111
1112                gaugeloc1 = gaugeloc.replace(' ','')
1113                #gaugeloc2 = gaugeloc1.replace('_','')
1114                gaugeloc2 = locations[k].replace(' ','')
1115                graphname = '%sgauge%s_%s' %(file_loc[j], gaugeloc2, which_quantity)
1116
1117                if report == True and len(label_id) > 1:
1118                    figdir = getcwd()+sep+'report_figures'+sep
1119                    if access(figdir,F_OK) == 0 :
1120                        mkdir (figdir)
1121                    latex_file_loc = figdir.replace(sep,altsep) 
1122                    graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory   
1123                    graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file
1124                    graphname_report.append(graphname_report_input)
1125                   
1126                    savefig(graphname_latex) # save figures in production directory for report generation
1127
1128                if report == True:
1129
1130                    figdir = getcwd()+sep+'report_figures'+sep
1131                    if access(figdir,F_OK) == 0 :
1132                        mkdir (figdir)
1133                    latex_file_loc = figdir.replace(sep,altsep)   
1134
1135                    if len(label_id) == 1: 
1136                        graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory 
1137                        graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file
1138                        s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png')
1139                        fid.write(s)
1140                        if where1 % 2 == 0:
1141                            s = '\\\\ \n'
1142                            where1 = 0
1143                        else:
1144                            s = '& \n'
1145                        fid.write(s)
1146                        savefig(graphname_latex)
1147               
1148                if title_on == True:
1149                    title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc2))
1150
1151                savefig(graphname) # save figures with sww file
1152
1153            if report == True and len(label_id) == 1:
1154                for i in range(nn-1):
1155                    if nn > 2:
1156                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1157                            word_quantity += 'depth' + ', '
1158                        else:
1159                            word_quantity += plot_quantity[i] + ', '
1160                    else:
1161                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1162                            word_quantity += 'depth' + ', '
1163                        else:
1164                            word_quantity += plot_quantity[i]
1165                   
1166                if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
1167                    word_quantity += ' and ' + 'depth'
1168                else:
1169                    word_quantity += ' and ' + plot_quantity[nn-1]
1170                caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' '))
1171                if elev[k] == 0.0:
1172                    caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j])
1173                    east = gauges[0]
1174                    north = gauges[1]
1175                    elev_output.append([locations[k],east,north,elevations[0,k,j]])
1176                label = '%sgauge%s' %(label_id2, gaugeloc2)
1177                s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
1178                fid.write(s)
1179                cc += 1
1180                if cc % 6 == 0: fid.write('\\clearpage \n')
1181                savefig(graphname_latex)               
1182               
1183        if report == True and len(label_id) > 1:
1184            for i in range(nn-1):
1185                if nn > 2:
1186                    if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1187                        word_quantity += 'depth' + ','
1188                    else:
1189                        word_quantity += plot_quantity[i] + ', '
1190                else:
1191                    if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1192                        word_quantity += 'depth'
1193                    else:
1194                        word_quantity += plot_quantity[i]
1195                where1 = 0
1196                count1 += 1
1197                index = j*len(plot_quantity)
1198                for which_quantity in plot_quantity:
1199                    where1 += 1
1200                    s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report[index], '.png')
1201                    index += 1
1202                    fid.write(s)
1203                    if where1 % 2 == 0:
1204                        s = '\\\\ \n'
1205                        where1 = 0
1206                    else:
1207                        s = '& \n'
1208                    fid.write(s)
1209            word_quantity += ' and ' + plot_quantity[nn-1]           
1210            label = 'gauge%s' %(gaugeloc2) 
1211            caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k])
1212            if elev[k] == 0.0:
1213                    caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j])
1214                    thisgauge = gauges[k]
1215                    east = thisgauge[0]
1216                    north = thisgauge[1]
1217                    elev_output.append([locations[k],east,north,elevations[0,k,j]])
1218                   
1219            s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
1220            fid.write(s)
1221            if float((k+1)/div - pp) == 0.:
1222                fid.write('\\clearpage \n')
1223                pp += 1
1224           
1225            #### finished generating figures ###
1226
1227        close('all')
1228   
1229    return texfile2, elev_output
1230
1231# FIXME (DSG): Add unit test, make general, not just 2 files,
1232# but any number of files.
1233def copy_code_files(dir_name, filename1, filename2):
1234    """Copies "filename1" and "filename2" to "dir_name". Very useful for
1235    information management """
1236
1237    if access(dir_name,F_OK) == 0:
1238        print 'Make directory %s' %dir_name
1239        mkdir (dir_name,0777)
1240    copy(filename1, dir_name + sep + basename(filename1))
1241    copy(filename2, dir_name + sep + basename(filename2))
1242#    copy (__file__, project.output_run_time_dir + basename(__file__))
1243    print 'Files %s and %s copied' %(filename1, filename2)
1244
1245
1246def add_directories(root_directory, directories):
1247    """
1248    Add the first directory in directories to root_directory.
1249    Then add the second
1250    direcotory to the first directory and so on.
1251
1252    Return the path of the final directory.
1253
1254    This is handy for specifying and creating a directory where data will go.
1255    """
1256    dir = root_directory
1257    for new_dir in directories:
1258        dir = os.path.join(dir, new_dir)
1259        if not access(dir,F_OK):
1260            mkdir (dir)
1261    return dir
1262
1263
1264
Note: See TracBrowser for help on using the repository browser.