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

Last change on this file since 3850 was 3850, checked in by ole, 17 years ago

Introduced a covariance measure to verify Okushiri timeseries plus
some incidental work

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