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

Last change on this file since 3900 was 3900, checked in by ole, 16 years ago

Added time_thinning to Interpolation_function and improved diagnostics.

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