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

Last change on this file since 6171 was 6171, checked in by ole, 15 years ago

Removed use of 'num' as a variable

File size: 98.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
8import sys
9import os
10
11from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep,getcwd
12from os.path import exists, basename, split,join
13from warnings import warn
14from shutil import copy
15
16from anuga.utilities.numerical_tools import ensure_numeric
17from Scientific.IO.NetCDF import NetCDFFile
18   
19from anuga.geospatial_data.geospatial_data import ensure_absolute
20from math import sqrt, atan, degrees
21
22# FIXME (Ole): Temporary short cuts -
23# FIXME (Ole): remove and update scripts where they are used
24from anuga.utilities.system_tools import get_revision_number
25from anuga.utilities.system_tools import store_version_info
26
27from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
28
29import Numeric as num
30
31
32##
33# @brief Read time history of data from NetCDF file, return callable object.
34# @param filename  Name of .sww or .tms file.
35# @param domain Associated domain object.
36# @param quantities Name of quantity to be interpolated or a list of names.
37# @param interpolation_points List of absolute UTM coordinates for points
38#                             (N x 2) or geospatial object or
39#                             points file name at which values are sought.
40# @param time_thinning
41# @param verbose True if this function is to be verbose.
42# @param use_cache True means that caching of intermediate result is attempted.
43# @param boundary_polygon
44# @return A callable object.
45def file_function(filename,
46                  domain=None,
47                  quantities=None,
48                  interpolation_points=None,
49                  time_thinning=1,
50                  verbose=False,
51                  use_cache=False,
52                  boundary_polygon=None):
53    """Read time history of spatial data from NetCDF file and return
54    a callable object.
55
56    Input variables:
57   
58    filename - Name of sww or tms file
59       
60       If the file has extension 'sww' then it is assumed to be spatio-temporal
61       or temporal and the callable object will have the form f(t,x,y) or f(t)
62       depending on whether the file contains spatial data
63
64       If the file has extension 'tms' then it is assumed to be temporal only
65       and the callable object will have the form f(t)
66
67       Either form will return interpolated values based on the input file
68       using the underlying interpolation_function.
69
70    domain - Associated domain object   
71       If domain is specified, model time (domain.starttime)
72       will be checked and possibly modified.
73   
74       All times are assumed to be in UTC
75       
76       All spatial information is assumed to be in absolute UTM coordinates.
77
78    quantities - the name of the quantity to be interpolated or a
79                 list of quantity names. The resulting function will return
80                 a tuple of values - one for each quantity
81                 If quantities are None, the default quantities are
82                 ['stage', 'xmomentum', 'ymomentum']
83                 
84
85    interpolation_points - list of absolute UTM coordinates for points (N x 2)
86    or geospatial object or points file name at which values are sought
87
88    time_thinning -
89
90    verbose -
91
92    use_cache: True means that caching of intermediate result of
93               Interpolation_function is attempted
94
95    boundary_polygon -
96
97   
98    See Interpolation function in anuga.fit_interpolate.interpolation for
99    further documentation
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
106    # Take into account:
107    # - domain's georef
108    # - sww file's georef
109    # - interpolation points as absolute UTM coordinates
110
111    if quantities is None:
112        if verbose:
113            msg = 'Quantities specified in file_function are None,'
114            msg += ' so I will use stage, xmomentum, and ymomentum in that order'
115            print msg
116        quantities = ['stage', 'xmomentum', 'ymomentum']
117
118    # Use domain's startime if available
119    if domain is not None:   
120        domain_starttime = domain.get_starttime()
121    else:
122        domain_starttime = None
123
124    # Build arguments and keyword arguments for use with caching or apply.
125    args = (filename,)
126
127    # FIXME (Ole): Caching this function will not work well
128    # if domain is passed in as instances change hash code.
129    # Instead we pass in those attributes that are needed (and return them
130    # if modified)
131    kwargs = {'quantities': quantities,
132              'interpolation_points': interpolation_points,
133              'domain_starttime': domain_starttime,
134              'time_thinning': time_thinning,                   
135              'verbose': verbose,
136              'boundary_polygon': boundary_polygon}
137
138    # Call underlying engine with or without caching
139    if use_cache is True:
140        try:
141            from caching import cache
142        except:
143            msg = 'Caching was requested, but caching module'+\
144                  'could not be imported'
145            raise msg
146
147        f, starttime = cache(_file_function,
148                             args, kwargs,
149                             dependencies=[filename],
150                             compression=False,                 
151                             verbose=verbose)
152    else:
153        f, starttime = apply(_file_function,
154                             args, kwargs)
155
156    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
157    #structure
158
159    f.starttime = starttime
160    f.filename = filename
161   
162    if domain is not None:
163        #Update domain.startime if it is *earlier* than starttime from file
164        if starttime > domain.starttime:
165            msg = 'WARNING: Start time as specified in domain (%f)'\
166                  %domain.starttime
167            msg += ' is earlier than the starttime of file %s (%f).'\
168                     %(filename, starttime)
169            msg += ' Modifying domain starttime accordingly.'
170           
171            if verbose: print msg
172
173            domain.set_starttime(starttime) #Modifying model time
174
175            if verbose: print 'Domain starttime is now set to %f'\
176               %domain.starttime
177    return f
178
179
180##
181# @brief ??
182# @param filename  Name of .sww or .tms file.
183# @param domain Associated domain object.
184# @param quantities Name of quantity to be interpolated or a list of names.
185# @param interpolation_points List of absolute UTM coordinates for points
186#                             (N x 2) or geospatial object or
187#                             points file name at which values are sought.
188# @param time_thinning
189# @param verbose True if this function is to be verbose.
190# @param use_cache True means that caching of intermediate result is attempted.
191# @param boundary_polygon
192def _file_function(filename,
193                   quantities=None,
194                   interpolation_points=None,
195                   domain_starttime=None,
196                   time_thinning=1,
197                   verbose=False,
198                   boundary_polygon=None):
199    """Internal function
200   
201    See file_function for documentatiton
202    """
203
204    assert type(filename) == type(''),\
205               'First argument to File_function must be a string'
206
207    try:
208        fid = open(filename)
209    except Exception, e:
210        msg = 'File "%s" could not be opened: Error="%s"' % (filename, e)
211        raise msg
212
213    # read first line of file, guess file type
214    line = fid.readline()
215    fid.close()
216
217    if line[:3] == 'CDF':
218        return get_netcdf_file_function(filename,
219                                        quantities,
220                                        interpolation_points,
221                                        domain_starttime,
222                                        time_thinning=time_thinning,
223                                        verbose=verbose,
224                                        boundary_polygon=boundary_polygon)
225    else:
226        # FIXME (Ole): Could add csv file here to address Ted Rigby's
227        # suggestion about reading hydrographs.
228        # This may also deal with the gist of ticket:289
229        raise 'Must be a NetCDF File'
230
231
232##
233# @brief ??
234# @param filename  Name of .sww or .tms file.
235# @param quantity_names Name of quantity to be interpolated or a list of names.
236# @param interpolation_points List of absolute UTM coordinates for points
237#                             (N x 2) or geospatial object or
238#                             points file name at which values are sought.
239# @param domain_starttime Start time from domain object.
240# @param time_thinning ??
241# @param verbose True if this function is to be verbose.
242# @param boundary_polygon ??
243# @return A callable object.
244def get_netcdf_file_function(filename,
245                             quantity_names=None,
246                             interpolation_points=None,
247                             domain_starttime=None,                           
248                             time_thinning=1,                             
249                             verbose=False,
250                             boundary_polygon=None):
251    """Read time history of spatial data from NetCDF sww file and
252    return a callable object f(t,x,y)
253    which will return interpolated values based on the input file.
254
255    Model time (domain_starttime)
256    will be checked, possibly modified and returned
257   
258    All times are assumed to be in UTC
259
260    See Interpolation function for further documetation
261    """
262
263    # FIXME: Check that model origin is the same as file's origin
264    # (both in UTM coordinates)
265    # If not - modify those from file to match domain
266    # (origin should be passed in)
267    # Take this code from e.g. dem2pts in data_manager.py
268    # FIXME: Use geo_reference to read and write xllcorner...
269
270    import time, calendar, types
271    from anuga.config import time_format
272
273    # Open NetCDF file
274    if verbose: print 'Reading', filename
275
276    fid = NetCDFFile(filename, netcdf_mode_r)
277
278    if type(quantity_names) == types.StringType:
279        quantity_names = [quantity_names]       
280
281    if quantity_names is None or len(quantity_names) < 1:
282        msg = 'No quantities are specified in file_function'
283        raise Exception, msg
284 
285    if interpolation_points is not None:
286        interpolation_points = ensure_absolute(interpolation_points)
287        msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1]
288        assert interpolation_points.shape[1] == 2, msg
289
290    # Now assert that requested quantitites (and the independent ones)
291    # are present in file
292    missing = []
293    for quantity in ['time'] + quantity_names:
294        if not fid.variables.has_key(quantity):
295            missing.append(quantity)
296
297    if len(missing) > 0:
298        msg = 'Quantities %s could not be found in file %s'\
299              % (str(missing), filename)
300        fid.close()
301        raise Exception, msg
302
303    # Decide whether this data has a spatial dimension
304    spatial = True
305    for quantity in ['x', 'y']:
306        if not fid.variables.has_key(quantity):
307            spatial = False
308
309    if filename[-3:] == 'tms' and spatial is True:
310        msg = 'Files of type tms must not contain spatial  information'
311        raise msg
312
313    if filename[-3:] == 'sww' and spatial is False:
314        msg = 'Files of type sww must contain spatial information'       
315        raise msg
316
317    if filename[-3:] == 'sts' and spatial is False:
318        #What if mux file only contains one point
319        msg = 'Files of type sts must contain spatial information'       
320        raise msg
321
322    if filename[-3:] == 'sts' and boundary_polygon is None:
323        #What if mux file only contains one point
324        msg = 'Files of type sts require boundary polygon'       
325        raise msg
326
327    # Get first timestep
328    try:
329        starttime = fid.starttime[0]
330    except ValueError:
331        msg = 'Could not read starttime from file %s' %filename
332        raise msg
333
334    # Get variables
335    # if verbose: print 'Get variables'   
336    time = fid.variables['time'][:]   
337
338    # Get time independent stuff
339    if spatial:
340        # Get origin
341        xllcorner = fid.xllcorner[0]
342        yllcorner = fid.yllcorner[0]
343        zone = fid.zone[0]       
344
345        x = fid.variables['x'][:]
346        y = fid.variables['y'][:]
347        if filename[-3:] == 'sww':
348            triangles = fid.variables['volumes'][:]
349
350        x = num.reshape(x, (len(x),1))
351        y = num.reshape(y, (len(y),1))
352        vertex_coordinates = num.concatenate((x,y), axis=1) #m x 2 array
353
354        if boundary_polygon is not None:
355            #removes sts points that do not lie on boundary
356            boundary_polygon=ensure_numeric(boundary_polygon)
357            boundary_polygon[:,0] -= xllcorner
358            boundary_polygon[:,1] -= yllcorner
359            temp=[]
360            boundary_id=[]
361            gauge_id=[]
362            for i in range(len(boundary_polygon)):
363                for j in range(len(x)):
364                    if num.allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
365                        #FIX ME:
366                        #currently gauges lat and long is stored as float and
367                        #then cast to double. This cuases slight repositioning
368                        #of vertex_coordinates.
369                        temp.append(boundary_polygon[i])
370                        gauge_id.append(j)
371                        boundary_id.append(i)
372                        break
373            gauge_neighbour_id=[]
374            for i in range(len(boundary_id)-1):
375                if boundary_id[i]+1==boundary_id[i+1]:
376                    gauge_neighbour_id.append(i+1)
377                else:
378                    gauge_neighbour_id.append(-1)
379            if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 \
380               and boundary_id[0]==0:
381                gauge_neighbour_id.append(0)
382            else:
383                gauge_neighbour_id.append(-1)
384            gauge_neighbour_id=ensure_numeric(gauge_neighbour_id)
385
386            if len(num.compress(gauge_neighbour_id>=0,gauge_neighbour_id)) \
387               != len(temp)-1:
388                msg='incorrect number of segments'
389                raise msg
390            vertex_coordinates=ensure_numeric(temp)
391            if len(vertex_coordinates)==0:
392                msg = 'None of the sts gauges fall on the boundary'
393                raise msg
394        else:
395            gauge_neighbour_id=None
396
397        if interpolation_points is not None:
398            # Adjust for georef
399            interpolation_points[:,0] -= xllcorner
400            interpolation_points[:,1] -= yllcorner       
401    else:
402        gauge_neighbour_id=None
403       
404    if domain_starttime is not None:
405        # If domain_startime is *later* than starttime,
406        # move time back - relative to domain's time
407        if domain_starttime > starttime:
408            time = time - domain_starttime + starttime
409
410        # FIXME Use method in geo to reconcile
411        # if spatial:
412        # assert domain.geo_reference.xllcorner == xllcorner
413        # assert domain.geo_reference.yllcorner == yllcorner
414        # assert domain.geo_reference.zone == zone       
415       
416    if verbose:
417        print 'File_function data obtained from: %s' %filename
418        print '  References:'
419        #print '    Datum ....' #FIXME
420        if spatial:
421            print '    Lower left corner: [%f, %f]'\
422                  %(xllcorner, yllcorner)
423        print '    Start time:   %f' %starttime               
424       
425   
426    # Produce values for desired data points at
427    # each timestep for each quantity
428    quantities = {}
429    for i, name in enumerate(quantity_names):
430        quantities[name] = fid.variables[name][:]
431        if boundary_polygon is not None:
432            #removes sts points that do not lie on boundary
433            quantities[name] = num.take(quantities[name],gauge_id,1)
434           
435    # Close sww, tms or sts netcdf file         
436    fid.close()
437
438    from anuga.fit_interpolate.interpolate import Interpolation_function
439
440    if not spatial:
441        vertex_coordinates = triangles = interpolation_points = None
442    if filename[-3:] == 'sts':#added
443        triangles = None
444        #vertex coordinates is position of urs gauges
445
446    # Return Interpolation_function instance as well as
447    # starttime for use to possible modify that of domain
448    return (Interpolation_function(time,
449                                   quantities,
450                                   quantity_names,
451                                   vertex_coordinates,
452                                   triangles,
453                                   interpolation_points,
454                                   time_thinning=time_thinning,
455                                   verbose=verbose,
456                                   gauge_neighbour_id=gauge_neighbour_id),
457            starttime)
458
459    # NOTE (Ole): Caching Interpolation function is too slow as
460    # the very long parameters need to be hashed.
461
462
463##
464# @brief Replace multiple substrings in a string.
465# @param text The string to operate on.
466# @param dictionary A dict containing replacements, key->value.
467# @return The new string.
468def multiple_replace(text, dictionary):
469    """Multiple replace of words in text
470
471    text:       String to be modified
472    dictionary: Mapping of words that are to be substituted
473
474    Python Cookbook 3.14 page 88 and page 90
475    http://code.activestate.com/recipes/81330/
476    """
477
478    import re
479   
480    #Create a regular expression from all of the dictionary keys
481    #matching only entire words
482    regex = re.compile(r'\b'+ \
483                       r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \
484                       r'\b' )
485
486    #For each match, lookup the corresponding value in the dictionary
487    return regex.sub(lambda match: dictionary[match.group(0)], text)
488
489
490##
491# @brief Apply arbitrary expressions to the values of a dict.
492# @param expression A string expression to apply.
493# @param dictionary The dictionary to apply the expression to.
494def apply_expression_to_dictionary(expression, dictionary):
495    """Apply arbitrary expression to values of dictionary
496
497    Given an expression in terms of the keys, replace key by the
498    corresponding values and evaluate.
499
500    expression: Arbitrary, e.g. arithmetric, expression relating keys
501                from dictionary.
502
503    dictionary: Mapping between symbols in expression and objects that
504                will be evaluated by expression.
505                Values in dictionary must support operators given in
506                expression e.g. by overloading
507
508    Due to a limitation with Numeric, this can not evaluate 0/0
509    In general, the user can fix by adding 1e-30 to the numerator.
510    SciPy core can handle this situation.
511    """
512
513    import types
514    import re
515
516    assert isinstance(expression, basestring)
517    assert type(dictionary) == types.DictType
518
519    #Convert dictionary values to textual representations suitable for eval
520    D = {}
521    for key in dictionary:
522        D[key] = 'dictionary["%s"]' % key
523
524    #Perform substitution of variables   
525    expression = multiple_replace(expression, D)
526
527    #Evaluate and return
528    try:
529        return eval(expression)
530    except NameError, e:
531        msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)
532        raise NameError, msg
533    except ValueError, e:
534        msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)
535        raise ValueError, msg
536
537
538##
539# @brief Format a float into a string.
540# @param value Float value to format.
541# @param format The format to use (%.2f is default).
542# @return The formatted float as a string.
543def get_textual_float(value, format = '%.2f'):
544    """Get textual representation of floating point numbers
545    and accept None as valid entry
546
547    format is a string - default = '%.2f'
548    """
549
550    if value is None:
551        return 'None'
552    else:
553        try:
554            float(value)
555        except:
556            # May this is a vector
557            if len(value) > 1:
558                s = '('
559                for v in value:
560                    s += get_textual_float(v, format) + ', '
561                   
562                s = s[:-2] + ')' # Strip trailing comma and close
563                return s
564            else:
565                raise 'Illegal input to get_textual_float:', value
566        else:
567            return format % float(value)
568
569
570#################################################################################
571# OBSOLETE STUFF
572#################################################################################
573
574# @note TEMP
575def angle(v1, v2):
576    """Temporary Interface to new location"""
577
578    import anuga.utilities.numerical_tools as NT       
579   
580    msg = 'angle has moved from util.py.  '
581    msg += 'Please use "from anuga.utilities.numerical_tools import angle"'
582    warn(msg, DeprecationWarning) 
583
584    return NT.angle(v1, v2)
585
586   
587# @note TEMP
588def anglediff(v0, v1):
589    """Temporary Interface to new location"""
590
591    import anuga.utilities.numerical_tools as NT
592   
593    msg = 'anglediff has moved from util.py.  '
594    msg += 'Please use "from anuga.utilities.numerical_tools import anglediff"'
595    warn(msg, DeprecationWarning) 
596
597    return NT.anglediff(v0, v1)   
598
599   
600# @note TEMP
601def mean(x):
602    """Temporary Interface to new location"""
603
604    import anuga.utilities.numerical_tools as NT   
605   
606    msg = 'mean has moved from util.py.  '
607    msg += 'Please use "from anuga.utilities.numerical_tools import mean"'
608    warn(msg, DeprecationWarning) 
609
610    return NT.mean(x)   
611
612
613# @note TEMP
614def point_on_line(*args, **kwargs):
615    """Temporary Interface to new location"""
616
617    msg = 'point_on_line has moved from util.py.  '
618    msg += 'Please use "from anuga.utilities.polygon import point_on_line"'
619    warn(msg, DeprecationWarning) 
620
621    return utilities.polygon.point_on_line(*args, **kwargs)     
622
623   
624# @note TEMP
625def inside_polygon(*args, **kwargs):
626    """Temporary Interface to new location"""
627
628    print 'inside_polygon has moved from util.py.  ',
629    print 'Please use "from anuga.utilities.polygon import inside_polygon"'
630
631    return utilities.polygon.inside_polygon(*args, **kwargs)   
632
633   
634# @note TEMP
635def outside_polygon(*args, **kwargs):
636    """Temporary Interface to new location"""
637
638    print 'outside_polygon has moved from util.py.  ',
639    print 'Please use "from anuga.utilities.polygon import outside_polygon"'
640
641    return utilities.polygon.outside_polygon(*args, **kwargs)   
642
643
644# @note TEMP
645def separate_points_by_polygon(*args, **kwargs):
646    """Temporary Interface to new location"""
647
648    print 'separate_points_by_polygon has moved from util.py.  ',
649    print 'Please use "from anuga.utilities.polygon import ' \
650          'separate_points_by_polygon"'
651
652    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)   
653
654
655# @note TEMP
656def read_polygon(*args, **kwargs):
657    """Temporary Interface to new location"""
658
659    print 'read_polygon has moved from util.py.  ',
660    print 'Please use "from anuga.utilities.polygon import read_polygon"'
661
662    return utilities.polygon.read_polygon(*args, **kwargs)   
663
664
665# @note TEMP
666def populate_polygon(*args, **kwargs):
667    """Temporary Interface to new location"""
668
669    print 'populate_polygon has moved from util.py.  ',
670    print 'Please use "from anuga.utilities.polygon import populate_polygon"'
671
672    return utilities.polygon.populate_polygon(*args, **kwargs)   
673
674
675#################################################################################
676# End of obsolete stuff ?
677#################################################################################
678
679# @note TEMP
680def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='',
681                         verbose=False):
682    """Temporary Interface to new location"""
683    from anuga.shallow_water.data_manager import start_screen_catcher \
684         as dm_start_screen_catcher
685
686    print 'start_screen_catcher has moved from util.py.  ',
687    print 'Please use "from anuga.shallow_water.data_manager import ' \
688          'start_screen_catcher"'
689   
690    return dm_start_screen_catcher(dir_name, myid='', numprocs='',
691                                   extra_info='', verbose=False)
692
693
694##
695# @brief Read a .sww file and plot the time series.
696# @param swwfiles Dictionary of .sww files.
697# @param gauge_filename Name of gauge file.
698# @param production_dirs ??
699# @param report If True, write figures to report directory.
700# @param reportname Name of generated report (if any).
701# @param plot_quantity List containing quantities to plot.
702# @param generate_fig If True, generate figures as well as CSV files.
703# @param surface If True, then generate solution surface with 3d plot.
704# @param time_min Beginning of user defined time range for plotting purposes.
705# @param time_max End of user defined time range for plotting purposes.
706# @param time_thinning ??
707# @param time_unit ??
708# @param title_on If True, export standard graphics with title.
709# @param use_cache If True, use caching.
710# @param verbose If True, this function is verbose.
711def sww2timeseries(swwfiles,
712                   gauge_filename,
713                   production_dirs,
714                   report=None,
715                   reportname=None,
716                   plot_quantity=None,
717                   generate_fig=False,
718                   surface=None,
719                   time_min=None,
720                   time_max=None,
721                   time_thinning=1,                   
722                   time_unit=None,
723                   title_on=None,
724                   use_cache=False,
725                   verbose=False):
726    """ Read sww file and plot the time series for the
727    prescribed quantities at defined gauge locations and
728    prescribed time range.
729
730    Input variables:
731
732    swwfiles        - dictionary of sww files with label_ids (used in
733                      generating latex output. It will be part of
734                      the directory name of file_loc (typically the timestamp).
735                      Helps to differentiate latex files for different
736                      simulations for a particular scenario. 
737                    - assume that all conserved quantities have been stored
738                    - assume each sww file has been simulated with same timestep
739   
740    gauge_filename  - name of file containing gauge data
741                        - easting, northing, name , elevation?
742                    - OR (this is not yet done)
743                        - structure which can be converted to a Numeric array,
744                          such as a geospatial data object
745                     
746    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
747                                                'boundaries': 'urs boundary'}
748                      this will use the second part as the label and the
749                      first part as the ?
750                      #FIXME: Is it a list or a dictionary
751                      # This is probably obsolete by now
752                     
753    report          - if True, then write figures to report_figures directory in
754                      relevant production directory
755                    - if False, figures are already stored with sww file
756                    - default to False
757
758    reportname      - name for report if wishing to generate report
759   
760    plot_quantity   - list containing quantities to plot, they must
761                      be the name of an existing quantity or one of
762                      the following possibilities
763                    - possibilities:
764                        - stage; 'stage'
765                        - depth; 'depth'
766                        - speed; calculated as absolute momentum
767                         (pointwise) divided by depth; 'speed'
768                        - bearing; calculated as the angle of the momentum
769                          vector (xmomentum, ymomentum) from the North; 'bearing'
770                        - absolute momentum; calculated as
771                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
772                        - x momentum; 'xmomentum'
773                        - y momentum; 'ymomentum'
774                    - default will be ['stage', 'speed', 'bearing']
775
776    generate_fig     - if True, generate figures as well as csv file
777                     - if False, csv files created only
778                     
779    surface          - if True, then generate solution surface with 3d plot
780                       and save to current working directory
781                     - default = False
782   
783    time_min         - beginning of user defined time range for plotting purposes
784                        - default will be first available time found in swwfile
785                       
786    time_max         - end of user defined time range for plotting purposes
787                        - default will be last available time found in swwfile
788                       
789    title_on        - if True, export standard graphics with title
790                    - if False, export standard graphics without title
791
792
793    Output:
794   
795    - time series data stored in .csv for later use if required.
796      Name = gauges_timeseries followed by gauge name
797    - latex file will be generated in same directory as where script is
798      run (usually production scenario directory.
799      Name = latexoutputlabel_id.tex
800
801    Other important information:
802   
803    It is assumed that the used has stored all the conserved quantities
804    and elevation during the scenario run, i.e.
805    ['stage', 'elevation', 'xmomentum', 'ymomentum']
806    If this has not occurred then sww2timeseries will not work.
807
808
809    Usage example
810    texname = sww2timeseries({project.boundary_name + '.sww': ''},
811                             project.polygons_dir + sep + 'boundary_extent.csv',
812                             project.anuga_dir,
813                             report = False,
814                             plot_quantity = ['stage', 'speed', 'bearing'],
815                             time_min = None,
816                             time_max = None,
817                             title_on = True,   
818                             verbose = True)
819   
820    """
821
822    msg = 'NOTE: A new function is available to create csv files from sww '
823    msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
824    msg += ' PLUS another new function to create graphs from csv files called '
825    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
826    print msg
827   
828    k = _sww2timeseries(swwfiles,
829                        gauge_filename,
830                        production_dirs,
831                        report,
832                        reportname,
833                        plot_quantity,
834                        generate_fig,
835                        surface,
836                        time_min,
837                        time_max,
838                        time_thinning,                       
839                        time_unit,
840                        title_on,
841                        use_cache,
842                        verbose)
843    return k
844
845
846##
847# @brief Read a .sww file and plot the time series.
848# @param swwfiles Dictionary of .sww files.
849# @param gauge_filename Name of gauge file.
850# @param production_dirs ??
851# @param report If True, write figures to report directory.
852# @param reportname Name of generated report (if any).
853# @param plot_quantity List containing quantities to plot.
854# @param generate_fig If True, generate figures as well as CSV files.
855# @param surface If True, then generate solution surface with 3d plot.
856# @param time_min Beginning of user defined time range for plotting purposes.
857# @param time_max End of user defined time range for plotting purposes.
858# @param time_thinning ??
859# @param time_unit ??
860# @param title_on If True, export standard graphics with title.
861# @param use_cache If True, use caching.
862# @param verbose If True, this function is verbose.
863def _sww2timeseries(swwfiles,
864                    gauge_filename,
865                    production_dirs,
866                    report = None,
867                    reportname = None,
868                    plot_quantity = None,
869                    generate_fig = False,
870                    surface = None,
871                    time_min = None,
872                    time_max = None,
873                    time_thinning = 1,                   
874                    time_unit = None,
875                    title_on = None,
876                    use_cache = False,
877                    verbose = False):   
878       
879    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
880   
881    try:
882        fid = open(gauge_filename)
883    except Exception, e:
884        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
885        raise msg
886
887    if report is None:
888        report = False
889       
890    if plot_quantity is None:
891        plot_quantity = ['depth', 'speed']
892    else:
893        assert type(plot_quantity) == list, 'plot_quantity must be a list'
894        check_list(plot_quantity)
895
896    if surface is None:
897        surface = False
898
899    if time_unit is None:
900        time_unit = 'hours'
901   
902    if title_on is None:
903        title_on = True
904   
905    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
906
907    gauges, locations, elev = get_gauges_from_file(gauge_filename)
908
909    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
910
911    file_loc = []
912    f_list = []
913    label_id = []
914    leg_label = []
915    themaxT = 0.0
916    theminT = 0.0
917
918    for swwfile in swwfiles.keys():
919        try:
920            fid = open(swwfile)
921        except Exception, e:
922            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
923            raise msg
924
925        print 'swwfile', swwfile
926
927        # Extract parent dir name and use as label
928        path, _ = os.path.split(swwfile)
929        _, label = os.path.split(path)       
930       
931        #print 'label', label
932        leg_label.append(label)
933
934        f = file_function(swwfile,
935                          quantities = sww_quantity,
936                          interpolation_points = gauges,
937                          time_thinning = time_thinning,
938                          verbose = verbose,
939                          use_cache = use_cache)
940
941        # determine which gauges are contained in sww file
942        count = 0
943        gauge_index = []
944        for k, g in enumerate(gauges):
945            if f(0.0, point_id = k)[2] > 1.0e6:
946                count += 1
947                if count == 1: print 'Gauges not contained here:'
948                print locations[k]
949            else:
950                gauge_index.append(k)
951
952        if len(gauge_index) > 0:
953            print 'Gauges contained here: \n',
954        else:
955            print 'No gauges contained here. \n'
956        for i in range(len(gauge_index)):
957             print locations[gauge_index[i]]
958             
959        index = swwfile.rfind(sep)
960        file_loc.append(swwfile[:index+1])
961        label_id.append(swwfiles[swwfile])
962       
963        f_list.append(f)
964        maxT = max(f.get_time())
965        minT = min(f.get_time())
966        if maxT > themaxT: themaxT = maxT
967        if minT > theminT: theminT = minT
968
969    if time_min is None:
970        time_min = theminT # min(T)
971    else:
972        if time_min < theminT: # min(T):
973            msg = 'Minimum time entered not correct - please try again'
974            raise Exception, msg
975
976    if time_max is None:
977        time_max = themaxT # max(T)
978    else:
979        if time_max > themaxT: # max(T):
980            msg = 'Maximum time entered not correct - please try again'
981            raise Exception, msg
982
983    if verbose and len(gauge_index) > 0:
984         print 'Inputs OK - going to generate figures'
985
986    if len(gauge_index) <> 0:
987        texfile, elev_output = \
988            generate_figures(plot_quantity, file_loc, report, reportname,
989                             surface, leg_label, f_list, gauges, locations,
990                             elev, gauge_index, production_dirs, time_min,
991                             time_max, time_unit, title_on, label_id,
992                             generate_fig, verbose)
993    else:
994        texfile = ''
995        elev_output = []
996
997    return texfile, elev_output
998
999
1000##
1001# @brief Read gauge info from a file.
1002# @param filename The name of the file to read.
1003# @return A (gauges, gaugelocation, elev) tuple.
1004def get_gauges_from_file(filename):
1005    """ Read in gauge information from file
1006    """
1007
1008    from os import sep, getcwd, access, F_OK, mkdir
1009
1010    # Get data from the gauge file
1011    fid = open(filename)
1012    lines = fid.readlines()
1013    fid.close()
1014   
1015    gauges = []
1016    gaugelocation = []
1017    elev = []
1018
1019    # Check header information   
1020    line1 = lines[0]
1021    line11 = line1.split(',')
1022
1023    if isinstance(line11[0], str) is True:
1024        # We have found text in the first line
1025        east_index = None
1026        north_index = None
1027        name_index = None
1028        elev_index = None
1029
1030        for i in range(len(line11)):
1031            if line11[i].strip().lower() == 'easting':   east_index = i
1032            if line11[i].strip().lower() == 'northing':  north_index = i
1033            if line11[i].strip().lower() == 'name':      name_index = i
1034            if line11[i].strip().lower() == 'elevation': elev_index = i
1035
1036        if east_index < len(line11) and north_index < len(line11):
1037            pass
1038        else:
1039            msg = 'WARNING: %s does not contain correct header information' \
1040                  % filename
1041            msg += 'The header must be: easting, northing, name, elevation'
1042            raise Exception, msg
1043
1044        if elev_index is None: 
1045            raise Exception
1046   
1047        if name_index is None: 
1048            raise Exception
1049
1050        lines = lines[1:] # Remove header from data
1051    else:
1052        # No header, assume that this is a simple easting, northing file
1053
1054        msg = 'There was no header in file %s and the number of columns is %d' \
1055              % (filename, len(line11))
1056        msg += '- was assuming two columns corresponding to Easting and Northing'
1057        assert len(line11) == 2, msg
1058
1059        east_index = 0
1060        north_index = 1
1061
1062        N = len(lines)
1063        elev = [-9999]*N
1064        gaugelocation = range(N)
1065       
1066    # Read in gauge data
1067    for line in lines:
1068        fields = line.split(',')
1069
1070        gauges.append([float(fields[east_index]), float(fields[north_index])])
1071
1072        if len(fields) > 2:
1073            elev.append(float(fields[elev_index]))
1074            loc = fields[name_index]
1075            gaugelocation.append(loc.strip('\n'))
1076
1077    return gauges, gaugelocation, elev
1078
1079
1080##
1081# @brief Check that input quantities in quantity list are legal.
1082# @param quantity Quantity list to check.
1083# @note Raises an exception of list is not legal.
1084def check_list(quantity):
1085    """ Check that input quantities in quantity list are possible
1086    """
1087    import sys
1088    from sets import Set as set
1089
1090    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
1091                    'ymomentum', 'speed', 'bearing', 'elevation']
1092
1093    # convert all quanitiy names to lowercase
1094    for i,j in enumerate(quantity):
1095        quantity[i] = quantity[i].lower()
1096
1097    # check that all names in 'quantity' appear in 'all_quantity'
1098    p = list(set(quantity).difference(set(all_quantity)))
1099    if len(p) != 0:
1100        msg = 'Quantities %s do not exist - please try again' %p
1101        raise Exception, msg
1102
1103
1104##
1105# @brief Calculate velocity bearing from North.
1106# @param uh ??
1107# @param vh ??
1108# @return The calculated bearing.
1109def calc_bearing(uh, vh):
1110    """ Calculate velocity bearing from North
1111    """
1112    #FIXME (Ole): I reckon we should refactor this one to use
1113    #             the function angle() in utilities/numerical_tools
1114    #
1115    #             It will be a simple matter of
1116    # * converting from radians to degrees
1117    # * moving the reference direction from [1,0] to North
1118    # * changing from counter clockwise to clocwise.
1119       
1120    angle = degrees(atan(vh/(uh+1.e-15)))
1121
1122    if (0 < angle < 90.0):
1123        if vh > 0:
1124            bearing = 90.0 - abs(angle)
1125        if vh < 0:
1126            bearing = 270.0 - abs(angle)
1127   
1128    if (-90 < angle < 0):
1129        if vh < 0:
1130            bearing = 90.0 - (angle)
1131        if vh > 0:
1132            bearing = 270.0 - (angle)
1133    if angle == 0: bearing = 0.0
1134
1135    return bearing
1136
1137
1138##
1139# @brief Generate figures from quantities and gauges for each sww file.
1140# @param plot_quantity  ??
1141# @param file_loc ??
1142# @param report ??
1143# @param reportname ??
1144# @param surface ??
1145# @param leg_label ??
1146# @param f_list ??
1147# @param gauges ??
1148# @param locations ??
1149# @param elev ??
1150# @param gauge_index ??
1151# @param production_dirs ??
1152# @param time_min ??
1153# @param time_max ??
1154# @param time_unit ??
1155# @param title_on ??
1156# @param label_id ??
1157# @param generate_fig ??
1158# @param verbose??
1159# @return (texfile2, elev_output)
1160def generate_figures(plot_quantity, file_loc, report, reportname, surface,
1161                     leg_label, f_list, gauges, locations, elev, gauge_index,
1162                     production_dirs, time_min, time_max, time_unit,
1163                     title_on, label_id, generate_fig, verbose):
1164    """ Generate figures based on required quantities and gauges for
1165    each sww file
1166    """
1167    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
1168
1169    if generate_fig is True:
1170        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
1171             xlabel, ylabel, title, close, subplot
1172   
1173        if surface is True:
1174            import pylab as p1
1175            import mpl3d.mplot3d as p3
1176       
1177    if report == True:   
1178        texdir = getcwd()+sep+'report'+sep
1179        if access(texdir,F_OK) == 0:
1180            mkdir (texdir)
1181        if len(label_id) == 1:
1182            label_id1 = label_id[0].replace(sep,'')
1183            label_id2 = label_id1.replace('_','')
1184            texfile = texdir + reportname + '%s' % label_id2
1185            texfile2 = reportname + '%s' % label_id2
1186            texfilename = texfile + '.tex'
1187            fid = open(texfilename, 'w')
1188
1189            if verbose: print '\n Latex output printed to %s \n' %texfilename
1190        else:
1191            texfile = texdir+reportname
1192            texfile2 = reportname
1193            texfilename = texfile + '.tex' 
1194            fid = open(texfilename, 'w')
1195
1196            if verbose: print '\n Latex output printed to %s \n' %texfilename
1197    else:
1198        texfile = ''
1199        texfile2 = ''
1200
1201    p = len(f_list)
1202    n = []
1203    n0 = 0
1204    for i in range(len(f_list)):
1205        n.append(len(f_list[i].get_time()))
1206        if n[i] > n0: n0 = n[i] 
1207    n0 = int(n0)
1208    m = len(locations)
1209    model_time = num.zeros((n0, m, p), num.Float) 
1210    stages = num.zeros((n0, m, p), num.Float)
1211    elevations = num.zeros((n0, m, p), num.Float) 
1212    momenta = num.zeros((n0, m, p), num.Float)
1213    xmom = num.zeros((n0, m, p), num.Float)
1214    ymom = num.zeros((n0, m, p), num.Float)
1215    speed = num.zeros((n0, m, p), num.Float)
1216    bearings = num.zeros((n0, m, p), num.Float)
1217    due_east = 90.0*num.ones((n0, 1), num.Float)
1218    due_west = 270.0*num.ones((n0, 1), num.Float)
1219    depths = num.zeros((n0, m, p), num.Float)
1220    eastings = num.zeros((n0, m, p), num.Float)
1221    min_stages = []
1222    max_stages = []
1223    min_momentums = []   
1224    max_momentums = []
1225    max_xmomentums = []
1226    max_ymomentums = []
1227    min_xmomentums = []
1228    min_ymomentums = []
1229    max_speeds = []
1230    min_speeds = []   
1231    max_depths = []
1232    model_time_plot3d = num.zeros((n0, m), num.Float)
1233    stages_plot3d = num.zeros((n0, m), num.Float)
1234    eastings_plot3d = num.zeros((n0, m),num.Float)
1235    if time_unit is 'mins': scale = 60.0
1236    if time_unit is 'hours': scale = 3600.0
1237
1238    ##### loop over each swwfile #####
1239    for j, f in enumerate(f_list):
1240        if verbose: print 'swwfile %d of %d' % (j, len(f_list))
1241
1242        starttime = f.starttime
1243        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
1244        fid_compare = open(comparefile, 'w')
1245        file0 = file_loc[j] + 'gauges_t0.csv'
1246        fid_0 = open(file0, 'w')
1247
1248        ##### loop over each gauge #####
1249        for k in gauge_index:
1250            if verbose: print 'Gauge %d of %d' % (k, len(gauges))
1251
1252            g = gauges[k]
1253            min_stage = 10
1254            max_stage = 0
1255            max_momentum = max_xmomentum = max_ymomentum = 0
1256            min_momentum = min_xmomentum = min_ymomentum = 100
1257            max_speed = 0
1258            min_speed = 0           
1259            max_depth = 0           
1260            gaugeloc = str(locations[k])
1261            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
1262                       + gaugeloc + '.csv'
1263            fid_out = open(thisfile, 'w')
1264            s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
1265            fid_out.write(s)
1266
1267            #### generate quantities #######
1268            for i, t in enumerate(f.get_time()):
1269                if time_min <= t <= time_max:
1270                    w = f(t, point_id = k)[0]
1271                    z = f(t, point_id = k)[1]
1272                    uh = f(t, point_id = k)[2]
1273                    vh = f(t, point_id = k)[3]
1274                    depth = w-z     
1275                    m = sqrt(uh*uh + vh*vh)
1276                    if depth < 0.001:
1277                        vel = 0.0
1278                    else:
1279                        vel = m / (depth + 1.e-6/depth) 
1280                    bearing = calc_bearing(uh, vh)                   
1281                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
1282                    stages[i,k,j] = w
1283                    elevations[i,k,j] = z
1284                    xmom[i,k,j] = uh
1285                    ymom[i,k,j] = vh
1286                    momenta[i,k,j] = m
1287                    speed[i,k,j] = vel
1288                    bearings[i,k,j] = bearing
1289                    depths[i,k,j] = depth
1290                    thisgauge = gauges[k]
1291                    eastings[i,k,j] = thisgauge[0]
1292                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
1293                            % (t, w, m, vel, z, uh, vh, bearing)
1294                    fid_out.write(s)
1295                    if t == 0:
1296                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
1297                        fid_0.write(s)
1298                    if t/60.0 <= 13920: tindex = i
1299                    if w > max_stage: max_stage = w
1300                    if w < min_stage: min_stage = w
1301                    if m > max_momentum: max_momentum = m
1302                    if m < min_momentum: min_momentum = m                   
1303                    if uh > max_xmomentum: max_xmomentum = uh
1304                    if vh > max_ymomentum: max_ymomentum = vh
1305                    if uh < min_xmomentum: min_xmomentum = uh
1306                    if vh < min_ymomentum: min_ymomentum = vh
1307                    if vel > max_speed: max_speed = vel
1308                    if vel < min_speed: min_speed = vel                   
1309                    if z > 0 and depth > max_depth: max_depth = depth
1310                   
1311                   
1312            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
1313                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
1314            fid_compare.write(s)
1315            max_stages.append(max_stage)
1316            min_stages.append(min_stage)
1317            max_momentums.append(max_momentum)
1318            max_xmomentums.append(max_xmomentum)
1319            max_ymomentums.append(max_ymomentum)
1320            min_xmomentums.append(min_xmomentum)
1321            min_ymomentums.append(min_ymomentum)
1322            min_momentums.append(min_momentum)           
1323            max_depths.append(max_depth)
1324            max_speeds.append(max_speed)
1325            min_speeds.append(min_speed)           
1326            #### finished generating quantities for each swwfile #####
1327       
1328        model_time_plot3d[:,:] = model_time[:,:,j]
1329        stages_plot3d[:,:] = stages[:,:,j]
1330        eastings_plot3d[:,] = eastings[:,:,j]
1331           
1332        if surface is True:
1333            print 'Printing surface figure'
1334            for i in range(2):
1335                fig = p1.figure(10)
1336                ax = p3.Axes3D(fig)
1337                if len(gauges) > 80:
1338                    ax.plot_surface(model_time[:,:,j],
1339                                    eastings[:,:,j],
1340                                    stages[:,:,j])
1341                else:
1342                    ax.plot3D(num.ravel(eastings[:,:,j]),
1343                              num.ravel(model_time[:,:,j]),
1344                              num.ravel(stages[:,:,j]))
1345                ax.set_xlabel('time')
1346                ax.set_ylabel('x')
1347                ax.set_zlabel('stage')
1348                fig.add_axes(ax)
1349                p1.show()
1350                surfacefig = 'solution_surface%s' % leg_label[j]
1351                p1.savefig(surfacefig)
1352                p1.close()
1353           
1354    #### finished generating quantities for all swwfiles #####
1355
1356    # x profile for given time
1357    if surface is True:
1358        figure(11)
1359        plot(eastings[tindex,:,j], stages[tindex,:,j])
1360        xlabel('x')
1361        ylabel('stage')
1362        profilefig = 'solution_xprofile' 
1363        savefig('profilefig')
1364
1365    elev_output = []
1366    if generate_fig is True:
1367        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
1368                           max(max_depths)*1.1])
1369        stage_axis = axis([starttime/scale, time_max/scale,
1370                           min(min_stages), max(max_stages)*1.1])
1371        vel_axis = axis([starttime/scale, time_max/scale,
1372                         min(min_speeds), max(max_speeds)*1.1])
1373        mom_axis = axis([starttime/scale, time_max/scale,
1374                         min(min_momentums), max(max_momentums)*1.1])
1375        xmom_axis = axis([starttime/scale, time_max/scale,
1376                          min(min_xmomentums), max(max_xmomentums)*1.1])
1377        ymom_axis = axis([starttime/scale, time_max/scale,
1378                          min(min_ymomentums), max(max_ymomentums)*1.1])
1379        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
1380        nn = len(plot_quantity)
1381        no_cols = 2
1382       
1383        if len(label_id) > 1: graphname_report = []
1384        pp = 1
1385        div = 11.
1386        cc = 0
1387        for k in gauge_index:
1388            g = gauges[k]
1389            count1 = 0
1390            if report == True and len(label_id) > 1:
1391                s = '\\begin{figure}[ht] \n' \
1392                    '\\centering \n' \
1393                    '\\begin{tabular}{cc} \n'
1394                fid.write(s)
1395            if len(label_id) > 1: graphname_report = []
1396
1397            #### generate figures for each gauge ####
1398            for j, f in enumerate(f_list):
1399                ion()
1400                hold(True)
1401                count = 0
1402                where1 = 0
1403                where2 = 0
1404                word_quantity = ''
1405                if report == True and len(label_id) == 1:
1406                    s = '\\begin{figure}[hbt] \n' \
1407                        '\\centering \n' \
1408                        '\\begin{tabular}{cc} \n'
1409                    fid.write(s)
1410                   
1411                for which_quantity in plot_quantity:
1412                    count += 1
1413                    where1 += 1
1414                    figure(count, frameon = False)
1415                    if which_quantity == 'depth':
1416                        plot(model_time[0:n[j]-1,k,j],
1417                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
1418                        units = 'm'
1419                        axis(depth_axis)
1420                    if which_quantity == 'stage':
1421                        if elevations[0,k,j] <= 0:
1422                            plot(model_time[0:n[j]-1,k,j],
1423                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
1424                            axis(stage_axis)
1425                        else:
1426                            plot(model_time[0:n[j]-1,k,j],
1427                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
1428                            #axis(depth_axis)                 
1429                        units = 'm'
1430                    if which_quantity == 'momentum':
1431                        plot(model_time[0:n[j]-1,k,j],
1432                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
1433                        axis(mom_axis)
1434                        units = 'm^2 / sec'
1435                    if which_quantity == 'xmomentum':
1436                        plot(model_time[0:n[j]-1,k,j],
1437                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
1438                        axis(xmom_axis)
1439                        units = 'm^2 / sec'
1440                    if which_quantity == 'ymomentum':
1441                        plot(model_time[0:n[j]-1,k,j],
1442                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
1443                        axis(ymom_axis)
1444                        units = 'm^2 / sec'
1445                    if which_quantity == 'speed':
1446                        plot(model_time[0:n[j]-1,k,j],
1447                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
1448                        axis(vel_axis)
1449                        units = 'm / sec'
1450                    if which_quantity == 'bearing':
1451                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
1452                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
1453                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
1454                        units = 'degrees from North'
1455                        #ax = axis([time_min, time_max, 0.0, 360.0])
1456                        legend(('Bearing','West','East'))
1457
1458                    if time_unit is 'mins': xlabel('time (mins)')
1459                    if time_unit is 'hours': xlabel('time (hours)')
1460                    #if which_quantity == 'stage' \
1461                    #   and elevations[0:n[j]-1,k,j] > 0:
1462                    #    ylabel('%s (%s)' %('depth', units))
1463                    #else:
1464                    #    ylabel('%s (%s)' %(which_quantity, units))
1465                        #ylabel('%s (%s)' %('wave height', units))
1466                    ylabel('%s (%s)' %(which_quantity, units))
1467                    if len(label_id) > 1: legend((leg_label),loc='upper right')
1468
1469                    #gaugeloc1 = gaugeloc.replace(' ','')
1470                    #gaugeloc2 = gaugeloc1.replace('_','')
1471                    gaugeloc2 = str(locations[k]).replace(' ','')
1472                    graphname = '%sgauge%s_%s' %(file_loc[j],
1473                                                 gaugeloc2,
1474                                                 which_quantity)
1475
1476                    if report == True and len(label_id) > 1:
1477                        figdir = getcwd()+sep+'report_figures'+sep
1478                        if access(figdir,F_OK) == 0 :
1479                            mkdir (figdir)
1480                        latex_file_loc = figdir.replace(sep,altsep) 
1481                        # storing files in production directory   
1482                        graphname_latex = '%sgauge%s%s' \
1483                                          % (latex_file_loc, gaugeloc2,
1484                                             which_quantity)
1485                        # giving location in latex output file
1486                        graphname_report_input = '%sgauge%s%s' % \
1487                                                 ('..' + altsep + 
1488                                                      'report_figures' + altsep,
1489                                                  gaugeloc2, which_quantity)
1490                        graphname_report.append(graphname_report_input)
1491                       
1492                        # save figures in production directory for report
1493                        savefig(graphname_latex)
1494
1495                    if report == True:
1496                        figdir = getcwd() + sep + 'report_figures' + sep
1497                        if access(figdir,F_OK) == 0:
1498                            mkdir(figdir)
1499                        latex_file_loc = figdir.replace(sep,altsep)   
1500
1501                        if len(label_id) == 1: 
1502                            # storing files in production directory 
1503                            graphname_latex = '%sgauge%s%s%s' % \
1504                                              (latex_file_loc, gaugeloc2,
1505                                               which_quantity, label_id2)
1506                            # giving location in latex output file
1507                            graphname_report = '%sgauge%s%s%s' % \
1508                                               ('..' + altsep +
1509                                                    'report_figures' + altsep,
1510                                                gaugeloc2, which_quantity,
1511                                                label_id2)
1512                            s = '\includegraphics' \
1513                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1514                                (graphname_report, '.png')
1515                            fid.write(s)
1516                            if where1 % 2 == 0:
1517                                s = '\\\\ \n'
1518                                where1 = 0
1519                            else:
1520                                s = '& \n'
1521                            fid.write(s)
1522                            savefig(graphname_latex)
1523                   
1524                    if title_on == True:
1525                        title('%s scenario: %s at %s gauge' % \
1526                              (label_id, which_quantity, gaugeloc2))
1527                        #title('Gauge %s (MOST elevation %.2f, ' \
1528                        #      'ANUGA elevation %.2f)' % \
1529                        #      (gaugeloc2, elevations[10,k,0],
1530                        #       elevations[10,k,1]))
1531
1532                    savefig(graphname) # save figures with sww file
1533
1534                if report == True and len(label_id) == 1:
1535                    for i in range(nn-1):
1536                        if nn > 2:
1537                            if plot_quantity[i] == 'stage' \
1538                               and elevations[0,k,j] > 0:
1539                                word_quantity += 'depth' + ', '
1540                            else:
1541                                word_quantity += plot_quantity[i] + ', '
1542                        else:
1543                            if plot_quantity[i] == 'stage' \
1544                               and elevations[0,k,j] > 0:
1545                                word_quantity += 'depth' + ', '
1546                            else:
1547                                word_quantity += plot_quantity[i]
1548                       
1549                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
1550                        word_quantity += ' and ' + 'depth'
1551                    else:
1552                        word_quantity += ' and ' + plot_quantity[nn-1]
1553                    caption = 'Time series for %s at %s location ' \
1554                              '(elevation %.2fm)' % \
1555                              (word_quantity, locations[k], elev[k])
1556                    if elev[k] == 0.0:
1557                        caption = 'Time series for %s at %s location ' \
1558                                  '(elevation %.2fm)' % \
1559                                  (word_quantity, locations[k],
1560                                   elevations[0,k,j])
1561                        east = gauges[0]
1562                        north = gauges[1]
1563                        elev_output.append([locations[k], east, north,
1564                                            elevations[0,k,j]])
1565                    label = '%sgauge%s' % (label_id2, gaugeloc2)
1566                    s = '\end{tabular} \n' \
1567                        '\\caption{%s} \n' \
1568                        '\label{fig:%s} \n' \
1569                        '\end{figure} \n \n' % (caption, label)
1570                    fid.write(s)
1571                    cc += 1
1572                    if cc % 6 == 0: fid.write('\\clearpage \n')
1573                    savefig(graphname_latex)               
1574                   
1575            if report == True and len(label_id) > 1:
1576                for i in range(nn-1):
1577                    if nn > 2:
1578                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1579                            word_quantity += 'depth' + ','
1580                        else:
1581                            word_quantity += plot_quantity[i] + ', '
1582                    else:
1583                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1584                            word_quantity += 'depth'
1585                        else:
1586                            word_quantity += plot_quantity[i]
1587                    where1 = 0
1588                    count1 += 1
1589                    index = j*len(plot_quantity)
1590                    for which_quantity in plot_quantity:
1591                        where1 += 1
1592                        s = '\includegraphics' \
1593                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1594                            (graphname_report[index], '.png')
1595                        index += 1
1596                        fid.write(s)
1597                        if where1 % 2 == 0:
1598                            s = '\\\\ \n'
1599                            where1 = 0
1600                        else:
1601                            s = '& \n'
1602                        fid.write(s)
1603                word_quantity += ' and ' + plot_quantity[nn-1]           
1604                label = 'gauge%s' %(gaugeloc2) 
1605                caption = 'Time series for %s at %s location ' \
1606                          '(elevation %.2fm)' % \
1607                          (word_quantity, locations[k], elev[k])
1608                if elev[k] == 0.0:
1609                        caption = 'Time series for %s at %s location ' \
1610                                  '(elevation %.2fm)' % \
1611                                  (word_quantity, locations[k],
1612                                   elevations[0,k,j])
1613                        thisgauge = gauges[k]
1614                        east = thisgauge[0]
1615                        north = thisgauge[1]
1616                        elev_output.append([locations[k], east, north,
1617                                            elevations[0,k,j]])
1618                       
1619                s = '\end{tabular} \n' \
1620                    '\\caption{%s} \n' \
1621                    '\label{fig:%s} \n' \
1622                    '\end{figure} \n \n' % (caption, label)
1623                fid.write(s)
1624                if float((k+1)/div - pp) == 0.:
1625                    fid.write('\\clearpage \n')
1626                    pp += 1
1627                #### finished generating figures ###
1628
1629            close('all')
1630       
1631    return texfile2, elev_output
1632
1633
1634# FIXME (DSG): Add unit test, make general, not just 2 files,
1635# but any number of files.
1636##
1637# @brief ??
1638# @param dir_name ??
1639# @param filename1 ??
1640# @param filename2 ??
1641# @return ??
1642# @note TEMP
1643def copy_code_files(dir_name, filename1, filename2):
1644    """Temporary Interface to new location"""
1645
1646    from anuga.shallow_water.data_manager import \
1647                    copy_code_files as dm_copy_code_files
1648    print 'copy_code_files has moved from util.py.  ',
1649    print 'Please use "from anuga.shallow_water.data_manager \
1650                                        import copy_code_files"'
1651   
1652    return dm_copy_code_files(dir_name, filename1, filename2)
1653
1654
1655##
1656# @brief Create a nested sub-directory path.
1657# @param root_directory The base diretory path.
1658# @param directories An iterable of sub-directory names.
1659# @return The final joined directory path.
1660# @note If each sub-directory doesn't exist, it will be created.
1661def add_directories(root_directory, directories):
1662    """
1663    Add the first sub-directory in 'directories' to root_directory.
1664    Then add the second sub-directory to the accumulating path and so on.
1665
1666    Return the path of the final directory.
1667
1668    This is handy for specifying and creating a directory where data will go.
1669    """
1670    dir = root_directory
1671    for new_dir in directories:
1672        dir = os.path.join(dir, new_dir)
1673        if not access(dir,F_OK):
1674            mkdir(dir)
1675    return dir
1676
1677
1678##
1679# @brief
1680# @param filename
1681# @param separator_value
1682# @return
1683# @note TEMP
1684def get_data_from_file(filename, separator_value=','):
1685    """Temporary Interface to new location"""
1686    from anuga.shallow_water.data_manager import \
1687                        get_data_from_file as dm_get_data_from_file
1688    print 'get_data_from_file has moved from util.py'
1689    print 'Please use "from anuga.shallow_water.data_manager \
1690                                     import get_data_from_file"'
1691   
1692    return dm_get_data_from_file(filename,separator_value = ',')
1693
1694
1695##
1696# @brief
1697# @param verbose
1698# @param kwargs
1699# @return
1700# @note TEMP
1701def store_parameters(verbose=False,**kwargs):
1702    """Temporary Interface to new location"""
1703   
1704    from anuga.shallow_water.data_manager \
1705                    import store_parameters as dm_store_parameters
1706    print 'store_parameters has moved from util.py.'
1707    print 'Please use "from anuga.shallow_water.data_manager \
1708                                     import store_parameters"'
1709   
1710    return dm_store_parameters(verbose=False,**kwargs)
1711
1712
1713##
1714# @brief Remove vertices that are not associated with any triangle.
1715# @param verts An iterable (or array) of points.
1716# @param triangles An iterable of 3 element tuples.
1717# @param number_of_full_nodes ??
1718# @return (verts, triangles) where 'verts' has been updated.
1719def remove_lone_verts(verts, triangles, number_of_full_nodes=None):
1720    """Removes vertices that are not associated with any triangles.
1721
1722    verts is a list/array of points.
1723    triangles is a list of 3 element tuples.  Each tuple represents a triangle.
1724    number_of_full_nodes relate to parallelism when a mesh has an
1725        extra layer of ghost points.
1726    """
1727
1728    verts = ensure_numeric(verts)
1729    triangles = ensure_numeric(triangles)
1730   
1731    N = len(verts)
1732   
1733    # initialise the array to easily find the index of the first loner
1734    # ie, if N=3 -> [6,5,4]
1735    loners=num.arange(2*N, N, -1)
1736    for t in triangles:
1737        for vert in t:
1738            loners[vert]= vert # all non-loners will have loners[i]=i
1739
1740    lone_start = 2*N - max(loners) # The index of the first loner
1741
1742    if lone_start-1 == N:
1743        # no loners
1744        pass
1745    elif min(loners[lone_start:N]) > N:
1746        # All the loners are at the end of the vert array
1747        verts = verts[0:lone_start]
1748    else:
1749        # change the loners list so it can be used to modify triangles
1750        # Remove the loners from verts
1751        # Could've used X=compress(less(loners,N),loners)
1752        # verts=num.take(verts,X)  to Remove the loners from verts
1753        # but I think it would use more memory
1754        new_i = lone_start      # point at first loner - 'shuffle down' target
1755        for i in range(lone_start, N):
1756            if loners[i] >= N:  # [i] is a loner, leave alone
1757                pass
1758            else:               # a non-loner, move down
1759                loners[i] = new_i
1760                verts[new_i] = verts[i]
1761                new_i += 1
1762        verts = verts[0:new_i]
1763
1764        # Modify the triangles
1765        #print "loners", loners
1766        #print "triangles before", triangles
1767        triangles = num.choose(triangles,loners)
1768        #print "triangles after", triangles
1769    return verts, triangles
1770
1771
1772##
1773# @brief Compute centroid values from vertex values
1774# @param x Values at vertices of triangular mesh.
1775# @param triangles Nx3 integer array pointing to vertex information.
1776# @return [N] array of centroid values.
1777def get_centroid_values(x, triangles):
1778    """Compute centroid values from vertex values
1779   
1780    x: Values at vertices of triangular mesh
1781    triangles: Nx3 integer array pointing to vertex information
1782    for each of the N triangels. Elements of triangles are
1783    indices into x
1784    """
1785       
1786    xc = num.zeros(triangles.shape[0], num.Float) # Space for centroid info
1787   
1788    for k in range(triangles.shape[0]):
1789        # Indices of vertices
1790        i0 = triangles[k][0]
1791        i1 = triangles[k][1]
1792        i2 = triangles[k][2]       
1793       
1794        xc[k] = (x[i0] + x[i1] + x[i2])/3
1795
1796    return xc
1797
1798
1799# @note TEMP
1800def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]},
1801                                output_dir='',
1802                                base_name='',
1803                                plot_numbers=['3-5'],
1804                                quantities=['speed','stage','momentum'],
1805                                assess_all_csv_files=True,
1806                                extra_plot_name='test'):
1807
1808    msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs '
1809    msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import ' \
1810           'csv2timeseries_graphs"'
1811    raise Exception, msg
1812
1813    return csv2timeseries_graphs(directories_dic,
1814                                 output_dir,
1815                                 base_name,
1816                                 plot_numbers,
1817                                 quantities,
1818                                 extra_plot_name,
1819                                 assess_all_csv_files)
1820
1821
1822##
1823# @brief Plot time series from CSV files.
1824# @param directories_dic
1825# @param output_dir
1826# @param base_name
1827# @param plot_numbers
1828# @param quantities
1829# @param extra_plot_name
1830# @param assess_all_csv_files
1831# @param create_latex
1832# @param verbose
1833# @note Assumes that 'elevation' is in the CSV file(s).
1834def csv2timeseries_graphs(directories_dic={},
1835                          output_dir='',
1836                          base_name=None,
1837                          plot_numbers='',
1838                          quantities=['stage'],
1839                          extra_plot_name='',
1840                          assess_all_csv_files=True,
1841                          create_latex=False,
1842                          verbose=False):
1843                               
1844    """
1845    Read in csv files that have the right header information and
1846    plot time series such as Stage, Speed, etc. Will also plot several
1847    time series on one plot. Filenames must follow this convention,
1848    <base_name><plot_number>.csv eg gauge_timeseries3.csv
1849   
1850    NOTE: relies that 'elevation' is in the csv file!
1851
1852    Each file represents a location and within each file there are
1853    time, quantity columns.
1854   
1855    For example:   
1856    if "directories_dic" defines 4 directories and in each directories
1857    there is a csv files corresponding to the right "plot_numbers",
1858    this will create a plot with 4 lines one for each directory AND
1859    one plot for each "quantities".  ??? FIXME: unclear.
1860   
1861    Usage:
1862        csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1863                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1864                            output_dir='fixed_wave'+sep,
1865                            base_name='gauge_timeseries_',
1866                            plot_numbers='',
1867                            quantities=['stage','speed'],
1868                            extra_plot_name='',
1869                            assess_all_csv_files=True,                           
1870                            create_latex=False,
1871                            verbose=True)
1872            this will create one plot for stage with both 'slide' and
1873            'fixed_wave' lines on it for stage and speed for each csv
1874            file with 'gauge_timeseries_' as the prefix. The graghs
1875            will be in the output directory 'fixed_wave' and the graph
1876            axis will be determined by assessing all the
1877   
1878    ANOTHER EXAMPLE
1879        new_csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1880                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1881                            output_dir='fixed_wave'+sep,
1882                            base_name='gauge_timeseries_',
1883                            plot_numbers=['1-3'],
1884                            quantities=['stage','speed'],
1885                            extra_plot_name='',
1886                            assess_all_csv_files=False,                           
1887                            create_latex=False,
1888                            verbose=True)
1889        This will plot csv files called gauge_timeseries_1.csv and
1890        gauge_timeseries3.csv from both 'slide' and 'fixed_wave' directories
1891        to 'fixed_wave'. There will be 4 plots created two speed and two stage
1892        one for each csv file. There will be two lines on each of these plots.
1893        And the axis will have been determined from only these files, had
1894        assess_all_csv_files = True all csv file with 'gauges_timeseries_' prefix
1895        would of been assessed.
1896   
1897    ANOTHER EXAMPLE   
1898         csv2timeseries_graphs({'J:'+sep+'anuga_validation'+sep:['new',20,-.1],
1899                                   'J:'+sep+'conical_island'+sep:['test',0,0]},
1900                                   output_dir='',
1901                                   plot_numbers=['1','3'],
1902                                   quantities=['stage','depth','bearing'],
1903                                   base_name='gauge_b',
1904                                   assess_all_csv_files=True,
1905                                  verbose=True)   
1906       
1907            This will produce one plot for each quantity (therefore 3) in the
1908            current directory, each plot will have 2 lines on them. The first
1909            plot named 'new' will have the time offseted by 20secs and the stage
1910            height adjusted by -0.1m
1911       
1912    Inputs:
1913        directories_dic: dictionary of directory with values (plot
1914                         legend name for directory), (start time of
1915                         the time series) and the (value to add to
1916                         stage if needed). For example
1917                         {dir1:['Anuga_ons',5000, 0],
1918                          dir2:['b_emoth',5000,1.5],
1919                          dir3:['b_ons',5000,1.5]}
1920                         Having multiple directories defined will plot them on
1921                         one plot, therefore there will be 3 lines on each of
1922                         these plot. If you only want one line per plot call
1923                         csv2timeseries_graph separately for each directory,
1924                         eg only have one directory in the 'directories_dic' in
1925                         each call.
1926                         
1927        output_dir: directory for the plot outputs. Only important to define when
1928                    you have more than one directory in your directories_dic, if
1929                    you have not defined it and you have multiple directories in
1930                    'directories_dic' there will be plots in each directory,
1931                    however only one directory will contain the complete
1932                    plot/graphs.
1933       
1934        base_name: Is used a couple of times.
1935                   1) to find the csv files to be plotted if there is no
1936                      'plot_numbers' then csv files with 'base_name' are plotted
1937                   2) in the title of the plots, the length of base_name is
1938                      removed from the front of the filename to be used in the
1939                      title.
1940                   This could be changed if needed.
1941                   Note is ignored if assess_all_csv_files=True
1942       
1943        plot_numbers: a String list of numbers to plot. For example
1944                      [0-4,10,15-17] will read and attempt to plot
1945                      the follow 0,1,2,3,4,10,15,16,17
1946                      NOTE: if no plot numbers this will create one plot per
1947                            quantity, per gauge
1948
1949        quantities: Will get available quantities from the header in the csv
1950                    file.  Quantities must be one of these.
1951                    NOTE: ALL QUANTITY NAMES MUST BE lower case!
1952                   
1953        extra_plot_name: A string that is appended to the end of the
1954                         output filename.
1955                   
1956        assess_all_csv_files: if true it will read ALL csv file with
1957                             "base_name", regardless of 'plot_numbers'
1958                              and determine a uniform set of axes for
1959                              Stage, Speed and Momentum. IF FALSE it
1960                              will only read the csv file within the
1961                             'plot_numbers'
1962                             
1963        create_latex: NOT IMPLEMENTED YET!! sorry Jane....
1964       
1965    OUTPUTS: saves the plots to
1966              <output_dir><base_name><plot_number><extra_plot_name>.png
1967    """
1968
1969    try: 
1970        import pylab
1971    except ImportError:
1972        msg='csv2timeseries_graphs needs pylab to be installed correctly'
1973        raise msg
1974            #ANUGA don't need pylab to work so the system doesn't
1975            #rely on pylab being installed
1976        return
1977
1978    from os import sep
1979    from anuga.shallow_water.data_manager import \
1980                               get_all_files_with_extension, csv2dict
1981
1982    seconds_in_hour = 3600
1983    seconds_in_minutes = 60
1984   
1985    quantities_label={}
1986#    quantities_label['time'] = 'time (hours)'
1987    quantities_label['time'] = 'time (minutes)'
1988    quantities_label['stage'] = 'wave height (m)'
1989    quantities_label['speed'] = 'speed (m/s)'
1990    quantities_label['momentum'] = 'momentum (m^2/sec)'
1991    quantities_label['depth'] = 'water depth (m)'
1992    quantities_label['xmomentum'] = 'momentum (m^2/sec)'
1993    quantities_label['ymomentum'] = 'momentum (m^2/sec)'
1994    quantities_label['bearing'] = 'degrees (o)'
1995    quantities_label['elevation'] = 'elevation (m)'
1996   
1997    if extra_plot_name != '':
1998        extra_plot_name = '_' + extra_plot_name
1999
2000    new_plot_numbers=[]
2001    #change plot_numbers to list, eg ['0-4','10']
2002    #to ['0','1','2','3','4','10']
2003    for i, num_string in enumerate(plot_numbers):
2004        if '-' in num_string: 
2005            start = int(num_string[:num_string.rfind('-')])
2006            end = int(num_string[num_string.rfind('-') + 1:]) + 1
2007            for x in range(start, end):
2008                new_plot_numbers.append(str(x))
2009        else:
2010            new_plot_numbers.append(num_string)
2011
2012    #finds all the files that fit the specs provided and return a list of them
2013    #so to help find a uniform max and min for the plots...
2014    list_filenames=[]
2015    all_csv_filenames=[]
2016    if verbose: print 'Determining files to access for axes ranges.'
2017   
2018    for i,directory in enumerate(directories_dic.keys()):
2019        all_csv_filenames.append(get_all_files_with_extension(directory,
2020                                                              base_name, '.csv'))
2021
2022        filenames=[]
2023        if plot_numbers == '': 
2024            list_filenames.append(get_all_files_with_extension(directory,
2025                                                               base_name,'.csv'))
2026        else:
2027            for number in new_plot_numbers:
2028                filenames.append(base_name + number)
2029            list_filenames.append(filenames)
2030
2031    #use all the files to get the values for the plot axis
2032    max_start_time= -1000.
2033    min_start_time = 100000 
2034   
2035    if verbose: print 'Determining uniform axes' 
2036
2037    #this entire loop is to determine the min and max range for the
2038    #axes of the plots
2039
2040#    quantities.insert(0,'elevation')
2041    quantities.insert(0,'time')
2042
2043    directory_quantity_value={}
2044#    quantity_value={}
2045    min_quantity_value={}
2046    max_quantity_value={}
2047
2048    for i, directory in enumerate(directories_dic.keys()):
2049        filename_quantity_value = {}
2050        if assess_all_csv_files == False:
2051            which_csv_to_assess = list_filenames[i]
2052        else:
2053            #gets list of filenames for directory "i"
2054            which_csv_to_assess = all_csv_filenames[i]
2055       
2056        for j, filename in enumerate(which_csv_to_assess):
2057            quantity_value = {}
2058
2059            dir_filename = join(directory,filename)
2060            attribute_dic, title_index_dic = csv2dict(dir_filename + '.csv')
2061            directory_start_time = directories_dic[directory][1]
2062            directory_add_tide = directories_dic[directory][2]
2063
2064            if verbose: print 'reading: %s.csv' %dir_filename
2065
2066            #add time to get values
2067            for k, quantity in enumerate(quantities):
2068                quantity_value[quantity] = [float(x) for
2069                                                x in attribute_dic[quantity]]
2070
2071                #add tide to stage if provided
2072                if quantity == 'stage':
2073                     quantity_value[quantity] = num.array(quantity_value[quantity]) \
2074                                                          + directory_add_tide
2075
2076                #condition to find max and mins for all the plots
2077                # populate the list with something when i=0 and j=0 and
2078                # then compare to the other values to determine abs max and min
2079                if i==0 and j==0:
2080                    min_quantity_value[quantity], \
2081                        max_quantity_value[quantity] = \
2082                            get_min_max_values(quantity_value[quantity])
2083
2084                    if quantity != 'time':
2085                        min_quantity_value[quantity] = \
2086                            min_quantity_value[quantity] *1.1
2087                        max_quantity_value[quantity] = \
2088                            max_quantity_value[quantity] *1.1
2089                else:
2090                    min, max = get_min_max_values(quantity_value[quantity])
2091               
2092                    # min and max are multipled by "1+increase_axis" to get axes
2093                    # that are slighty bigger than the max and mins
2094                    # so the plots look good.
2095
2096                    increase_axis = (max-min)*0.05
2097                    if min <= min_quantity_value[quantity]:
2098                        if quantity == 'time': 
2099                            min_quantity_value[quantity] = min
2100                        else:
2101                            if round(min,2) == 0.00:
2102                                min_quantity_value[quantity] = -increase_axis
2103#                                min_quantity_value[quantity] = -2.
2104                                #min_quantity_value[quantity] = \
2105                                #    -max_quantity_value[quantity]*increase_axis
2106                            else:
2107#                                min_quantity_value[quantity] = \
2108#                                    min*(1+increase_axis)
2109                                min_quantity_value[quantity]=min-increase_axis
2110                   
2111                    if max > max_quantity_value[quantity]: 
2112                        if quantity == 'time': 
2113                            max_quantity_value[quantity] = max
2114                        else:
2115                            max_quantity_value[quantity] = max + increase_axis
2116#                            max_quantity_value[quantity]=max*(1+increase_axis)
2117
2118            #set the time... ???
2119            if min_start_time > directory_start_time: 
2120                min_start_time = directory_start_time
2121            if max_start_time < directory_start_time: 
2122                max_start_time = directory_start_time
2123           
2124            filename_quantity_value[filename]=quantity_value
2125           
2126        directory_quantity_value[directory]=filename_quantity_value
2127   
2128    #final step to unifrom axis for the graphs
2129    quantities_axis={}
2130   
2131    for i, quantity in enumerate(quantities):
2132        quantities_axis[quantity] = (float(min_start_time) \
2133                                         / float(seconds_in_minutes),
2134                                     (float(max_quantity_value['time']) \
2135                                          + float(max_start_time)) \
2136                                              / float(seconds_in_minutes),
2137                                     min_quantity_value[quantity],
2138                                     max_quantity_value[quantity])
2139
2140        if verbose and (quantity != 'time' and quantity != 'elevation'): 
2141            print 'axis for quantity %s are x:(%s to %s)%s and y:(%s to %s)%s' \
2142                  % (quantity, 
2143                     quantities_axis[quantity][0],
2144                     quantities_axis[quantity][1],
2145                     quantities_label['time'],
2146                     quantities_axis[quantity][2],
2147                     quantities_axis[quantity][3],
2148                     quantities_label[quantity])
2149
2150    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
2151
2152    if verbose: print 'Now start to plot \n'
2153   
2154    i_max = len(directories_dic.keys())
2155    legend_list_dic = {}
2156    legend_list = []
2157    for i, directory in enumerate(directories_dic.keys()):
2158        if verbose: print 'Plotting in %s %s' % (directory, new_plot_numbers)
2159
2160        # FIXME THIS SORT IS VERY IMPORTANT
2161        # Without it the assigned plot numbers may not work correctly
2162        # there must be a better way
2163        list_filenames[i].sort()
2164        for j, filename in enumerate(list_filenames[i]):
2165            if verbose: print 'Starting %s' % filename 
2166
2167            directory_name = directories_dic[directory][0]
2168            directory_start_time = directories_dic[directory][1]
2169            directory_add_tide = directories_dic[directory][2]
2170           
2171            # create an if about the start time and tide height if don't exist
2172            attribute_dic, title_index_dic = csv2dict(directory + sep
2173                                                      + filename + '.csv')
2174            #get data from dict in to list
2175            #do maths to list by changing to array
2176            t = (num.array(directory_quantity_value[directory][filename]['time'])
2177                     + directory_start_time) / seconds_in_minutes
2178
2179            #finds the maximum elevation, used only as a test
2180            # and as info in the graphs
2181            max_ele=-100000
2182            min_ele=100000
2183            elevation = [float(x) for x in attribute_dic["elevation"]]
2184           
2185            min_ele, max_ele = get_min_max_values(elevation)
2186           
2187            if min_ele != max_ele:
2188                print "Note! Elevation changes in %s" %dir_filename
2189
2190            # creates a dictionary with keys that is the filename and attributes
2191            # are a list of lists containing 'directory_name' and 'elevation'.
2192            # This is used to make the contents for the legends in the graphs,
2193            # this is the name of the model and the elevation.  All in this
2194            # great one liner from DG. If the key 'filename' doesn't exist it
2195            # creates the entry if the entry exist it appends to the key.
2196
2197            legend_list_dic.setdefault(filename,[]) \
2198                .append([directory_name, round(max_ele, 3)])
2199
2200            # creates a LIST for the legend on the last iteration of the
2201            # directories which is when "legend_list_dic" has been fully
2202            # populated. Creates a list of strings which is used in the legend
2203            # only runs on the last iteration for all the gauges(csv) files
2204            # empties the list before creating it
2205
2206            if i == i_max - 1:
2207                legend_list = []
2208   
2209                for name_and_elevation in legend_list_dic[filename]:
2210                    legend_list.append('%s (elevation = %sm)'\
2211                                       % (name_and_elevation[0],
2212                                          name_and_elevation[1]))
2213           
2214            #skip time and elevation so it is not plotted!
2215            for k, quantity in enumerate(quantities):
2216                if quantity != 'time' and quantity != 'elevation':
2217                    pylab.figure(int(k*100+j))
2218                    pylab.ylabel(quantities_label[quantity])
2219                    pylab.plot(t,
2220                               directory_quantity_value[directory]\
2221                                                       [filename][quantity],
2222                               c = cstr[i], linewidth=1)
2223                    pylab.xlabel(quantities_label['time'])
2224                    pylab.axis(quantities_axis[quantity])
2225                    pylab.legend(legend_list,loc='upper right')
2226                   
2227                    pylab.title('%s at %s gauge'
2228                                % (quantity, filename[len(base_name):]))
2229
2230                    if output_dir == '':
2231                        figname = '%s%s%s_%s%s.png' \
2232                                  % (directory, sep, filename, quantity,
2233                                     extra_plot_name)
2234                    else:
2235                        figname = '%s%s%s_%s%s.png' \
2236                                  % (output_dir, sep, filename, quantity,
2237                                     extra_plot_name)
2238
2239                    if verbose: print 'saving figure here %s' %figname
2240
2241                    pylab.savefig(figname)
2242           
2243    if verbose: print 'Closing all plots'
2244
2245    pylab.close('all')
2246    del pylab
2247
2248    if verbose: print 'Finished closing plots'
2249
2250##
2251# @brief Return min and max of an iterable.
2252# @param list The iterable to return min & max of.
2253# @return (min, max) of 'list'.
2254def get_min_max_values(list=None):
2255    """
2256    Returns the min and max of the list it was provided.
2257    """
2258
2259    if list == None: print 'List must be provided'
2260       
2261    return min(list), max(list)
2262
2263
2264##
2265# @brief Get runup around a point in a CSV file.
2266# @param gauge_filename gauge file name.
2267# @param sww_filename SWW file name.
2268# @param runup_filename Name of file to report into.
2269# @param size ??
2270# @param verbose ??
2271def get_runup_data_for_locations_from_file(gauge_filename,
2272                                           sww_filename,
2273                                           runup_filename,
2274                                           size=10,
2275                                           verbose=False):
2276    """this will read a csv file with the header x,y. Then look in a square
2277    'size'x2 around this position for the 'max_inundaiton_height' in the
2278    'sww_filename' and report the findings in the 'runup_filename'.
2279   
2280    WARNING: NO TESTS!
2281    """
2282
2283    from anuga.shallow_water.data_manager import get_all_directories_with_name,\
2284                                                 get_maximum_inundation_data,\
2285                                                 csv2dict
2286                                                 
2287    file = open(runup_filename, "w")
2288    file.write("easting,northing,runup \n ")
2289    file.close()
2290   
2291    #read gauge csv file to dictionary
2292    attribute_dic, title_index_dic = csv2dict(gauge_filename)
2293    northing = [float(x) for x in attribute_dic["y"]]
2294    easting = [float(x) for x in attribute_dic["x"]]
2295
2296    print 'Reading %s' %sww_filename
2297
2298    runup_locations=[]
2299    for i, x in enumerate(northing):
2300        poly = [[int(easting[i]+size),int(northing[i]+size)],
2301                [int(easting[i]+size),int(northing[i]-size)],
2302                [int(easting[i]-size),int(northing[i]-size)],
2303                [int(easting[i]-size),int(northing[i]+size)]]
2304       
2305        run_up, x_y = get_maximum_inundation_data(filename=sww_filename,
2306                                                  polygon=poly,
2307                                                  verbose=False) 
2308
2309        #if no runup will return 0 instead of NONE
2310        if run_up==None: run_up=0
2311        if x_y==None: x_y=[0,0]
2312       
2313        if verbose:
2314            print 'maximum inundation runup near %s is %s meters' %(x_y,run_up)
2315       
2316        #writes to file
2317        file = open(runup_filename, "a")
2318        temp = '%s,%s,%s \n' % (x_y[0], x_y[1], run_up)
2319        file.write(temp)
2320        file.close()
2321
2322
2323##
2324# @brief ??
2325# @param sww_file ??
2326# @param gauge_file ??
2327# @param out_name ??
2328# @param quantities ??
2329# @param verbose ??
2330# @param use_cache ??
2331def sww2csv_gauges(sww_file,
2332                   gauge_file,
2333                   out_name='gauge_',
2334                   quantities=['stage', 'depth', 'elevation',
2335                               'xmomentum', 'ymomentum'],
2336                   verbose=False,
2337                   use_cache = True):
2338    """
2339   
2340    Inputs:
2341       
2342        NOTE: if using csv2timeseries_graphs after creating csv file,
2343        it is essential to export quantities 'depth' and 'elevation'.
2344        'depth' is good to analyse gauges on land and elevation is used
2345        automatically by csv2timeseries_graphs in the legend.
2346       
2347        sww_file: path to any sww file
2348       
2349        gauge_file: Assumes that it follows this format
2350            name, easting, northing, elevation
2351            point1, 100.3, 50.2, 10.0
2352            point2, 10.3, 70.3, 78.0
2353       
2354        NOTE: order of column can change but names eg 'easting', 'elevation'
2355        must be the same! ALL lowercaps!
2356
2357        out_name: prefix for output file name (default is 'gauge_')
2358       
2359    Outputs:
2360        one file for each gauge/point location in the points file. They
2361        will be named with this format in the same directory as the 'sww_file'
2362            <out_name><name>.csv
2363        eg gauge_point1.csv if <out_name> not supplied
2364           myfile_2_point1.csv if <out_name> ='myfile_2_'
2365           
2366           
2367        They will all have a header
2368   
2369    Usage: sww2csv_gauges(sww_file='test1.sww',
2370                          quantities = ['stage', 'elevation','depth','bearing'],
2371                          gauge_file='gauge.txt')   
2372   
2373    Interpolate the quantities at a given set of locations, given
2374    an sww file.
2375    The results are written to a csv file.
2376
2377    In the future let points be a points file.
2378    And the user choose the quantities.
2379
2380    This is currently quite specific.
2381    If it needs to be more general, change things.
2382
2383    This is really returning speed, not velocity.
2384    """
2385   
2386    from csv import reader,writer
2387    from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
2388    import string
2389    from anuga.shallow_water.data_manager import get_all_swwfiles
2390
2391#    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
2392    #print "points",points
2393
2394    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
2395    assert type(out_name) == type(''), 'Output filename prefix must be a string'
2396   
2397    try:
2398        point_reader = reader(file(gauge_file))
2399    except Exception, e:
2400        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e)
2401        raise msg
2402
2403    if verbose: print '\n Gauges obtained from: %s \n' %gauge_file
2404   
2405    point_reader = reader(file(gauge_file))
2406    points = []
2407    point_name = []
2408   
2409    #read point info from file
2410    for i,row in enumerate(point_reader):
2411        #read header and determine the column numbers to read correcty.
2412        if i==0:
2413            for j,value in enumerate(row):
2414                if value.strip()=='easting':easting=j
2415                if value.strip()=='northing':northing=j
2416                if value.strip()=='name':name=j
2417                if value.strip()=='elevation':elevation=j
2418        else:
2419            points.append([float(row[easting]),float(row[northing])])
2420            point_name.append(row[name])
2421       
2422    #convert to array for file_function
2423    points_array = num.array(points,num.Float)
2424       
2425    points_array = ensure_absolute(points_array)
2426
2427    dir_name, base = os.path.split(sww_file)   
2428
2429    #need to get current directory so when path and file
2430    #are "joined" below the directory is correct
2431    if dir_name == '':
2432        dir_name =getcwd()
2433       
2434    if access(sww_file,R_OK):
2435        if verbose: print 'File %s exists' %(sww_file)
2436    else:
2437        msg = 'File "%s" could not be opened: no read permission' % sww_file
2438        raise msg
2439
2440    sww_files = get_all_swwfiles(look_in_dir=dir_name,
2441                                 base_name=base,
2442                                 verbose=verbose)
2443   
2444    #to make all the quantities lower case for file_function
2445    quantities = [quantity.lower() for quantity in quantities]
2446
2447    # what is quantities are needed from sww file to calculate output quantities
2448    # also
2449
2450    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
2451   
2452    for sww_file in sww_files:
2453        sww_file = join(dir_name, sww_file+'.sww')
2454        callable_sww = file_function(sww_file,
2455                                     quantities=core_quantities,
2456                                     interpolation_points=points_array,
2457                                     verbose=verbose,
2458                                     use_cache=use_cache)
2459    gauge_file = out_name
2460
2461    heading = [quantity for quantity in quantities]
2462    heading.insert(0,'time')
2463
2464    #create a list of csv writers for all the points and write header
2465    points_writer = []
2466    for i,point in enumerate(points):
2467        points_writer.append(writer(file(dir_name + sep + gauge_file
2468                                         + point_name[i] + '.csv', "wb")))
2469        points_writer[i].writerow(heading)
2470   
2471    if verbose: print 'Writing csv files'
2472
2473    for time in callable_sww.get_time():
2474        for point_i, point in enumerate(points_array):
2475            #add domain starttime to relative time.
2476            points_list = [time + callable_sww.starttime]
2477            point_quantities = callable_sww(time,point_i)
2478           
2479            for quantity in quantities:
2480                if quantity == NAN:
2481                    print 'quantity does not exist in' % callable_sww.get_name
2482                else:
2483                    if quantity == 'stage':
2484                        points_list.append(point_quantities[0])
2485                       
2486                    if quantity == 'elevation':
2487                        points_list.append(point_quantities[1])
2488                       
2489                    if quantity == 'xmomentum':
2490                        points_list.append(point_quantities[2])
2491                       
2492                    if quantity == 'ymomentum':
2493                        points_list.append(point_quantities[3])
2494                       
2495                    if quantity == 'depth':
2496                        points_list.append(point_quantities[0] 
2497                                           - point_quantities[1])
2498
2499                    if quantity == 'momentum':
2500                        momentum = sqrt(point_quantities[2]**2 
2501                                        + point_quantities[3]**2)
2502                        points_list.append(momentum)
2503                       
2504                    if quantity == 'speed':
2505                        #if depth is less than 0.001 then speed = 0.0
2506                        if point_quantities[0] - point_quantities[1] < 0.001:
2507                            vel = 0.0
2508                        else:
2509                            if point_quantities[2] < 1.0e6:
2510                                momentum = sqrt(point_quantities[2]**2
2511                                                + point_quantities[3]**2)
2512    #                            vel = momentum/depth             
2513                                vel = momentum / (point_quantities[0] 
2514                                                  - point_quantities[1])
2515    #                            vel = momentum/(depth + 1.e-6/depth)
2516                            else:
2517                                momentum = 0
2518                                vel = 0
2519                           
2520                        points_list.append(vel)
2521                       
2522                    if quantity == 'bearing':
2523                        points_list.append(calc_bearing(point_quantities[2],
2524                                                        point_quantities[3]))
2525
2526            points_writer[point_i].writerow(points_list)
2527       
2528
2529##
2530# @brief Get a wave height at a certain depth given wave height at another depth.
2531# @param d1 The first depth.
2532# @param d2 The second depth.
2533# @param h1 Wave ampitude at d1
2534# @param verbose True if this function is to be verbose.
2535# @return The wave height at d2.
2536def greens_law(d1, d2, h1, verbose=False):
2537    """Green's Law
2538
2539    Green's Law allows an approximation of wave amplitude at
2540    a given depth based on the fourh root of the ratio of two depths
2541    and the amplitude at another given depth.
2542
2543    Note, wave amplitude is equal to stage.
2544   
2545    Inputs:
2546
2547    d1, d2 - the two depths
2548    h1     - the wave amplitude at d1
2549    h2     - the derived amplitude at d2
2550
2551    h2 = h1 (d1/d2)^(1/4), where d2 cannot equal 0.
2552   
2553    """
2554
2555    d1 = ensure_numeric(d1)
2556    d2 = ensure_numeric(d2)
2557    h1 = ensure_numeric(h1)
2558
2559    if d1 <= 0.0:
2560        msg = 'the first depth, d1 (%f), must be strictly positive' % (d1)
2561        raise Exception(msg)
2562
2563    if d2 <= 0.0:
2564        msg = 'the second depth, d2 (%f), must be strictly positive' % (d2)
2565        raise Exception(msg)
2566   
2567    if h1 <= 0.0:
2568        msg = 'the wave amplitude, h1 (%f), must be strictly positive' % (h1)
2569        raise Exception(msg)
2570   
2571    h2 = h1*(d1/d2)**0.25
2572
2573    assert h2 > 0
2574   
2575    return h2
2576       
2577
2578##
2579# @brief Get the square-root of a value.
2580# @param s The value to get the square-root of.
2581# @return The square-root of 's'.
2582def square_root(s):
2583    return sqrt(s)
2584
2585
Note: See TracBrowser for help on using the repository browser.