source: inundation/pyvolution/util.py @ 3190

Last change on this file since 3190 was 3190, checked in by sexton, 18 years ago

MOST and ANUGA comparisons and updates

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