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

Last change on this file since 6070 was 6070, checked in by rwilson, 15 years ago

Fixed problem in remove_lone_verts() and added an extra test case.

File size: 98.7 KB
Line 
1#12345678901234567890123456789012345678901234567890123456789012345678901234567890
2"""This module contains various auxiliary function used by pyvolution.
3
4It is also a clearing house for functions that may later earn a module
5of their own.
6"""
7
8import anuga.utilities.polygon
9import sys
10import os
11
12from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep,getcwd
13from os.path import exists, basename, split,join
14from warnings import warn
15from shutil import copy
16
17from anuga.utilities.numerical_tools import ensure_numeric
18from Scientific.IO.NetCDF import NetCDFFile
19from Numeric import arange, choose, zeros, Float, array, allclose, take, compress
20   
21from anuga.geospatial_data.geospatial_data import ensure_absolute
22from math import sqrt, atan, degrees
23
24# FIXME (Ole): Temporary short cuts -
25# FIXME (Ole): remove and update scripts where they are used
26from anuga.utilities.system_tools import get_revision_number
27from anuga.utilities.system_tools import store_version_info
28
29
30##
31# @brief Read time history of data from NetCDF file, return callable object.
32# @param filename  Name of .sww or .tms file.
33# @param domain Associated domain object.
34# @param quantities Name of quantity to be interpolated or a list of names.
35# @param interpolation_points List of absolute UTM coordinates for points
36#                             (N x 2) or geospatial object or
37#                             points file name at which values are sought.
38# @param time_thinning
39# @param verbose True if this function is to be verbose.
40# @param use_cache True means that caching of intermediate result is attempted.
41# @param boundary_polygon
42# @return A callable object.
43def file_function(filename,
44                  domain=None,
45                  quantities=None,
46                  interpolation_points=None,
47                  time_thinning=1,
48                  verbose=False,
49                  use_cache=False,
50                  boundary_polygon=None):
51    """Read time history of spatial data from NetCDF file and return
52    a callable object.
53
54    Input variables:
55   
56    filename - Name of sww or tms file
57       
58       If the file has extension 'sww' then it is assumed to be spatio-temporal
59       or temporal and the callable object will have the form f(t,x,y) or f(t)
60       depending on whether the file contains spatial data
61
62       If the file has extension 'tms' then it is assumed to be temporal only
63       and the callable object will have the form f(t)
64
65       Either form will return interpolated values based on the input file
66       using the underlying interpolation_function.
67
68    domain - Associated domain object   
69       If domain is specified, model time (domain.starttime)
70       will be checked and possibly modified.
71   
72       All times are assumed to be in UTC
73       
74       All spatial information is assumed to be in absolute UTM coordinates.
75
76    quantities - the name of the quantity to be interpolated or a
77                 list of quantity names. The resulting function will return
78                 a tuple of values - one for each quantity
79                 If quantities are None, the default quantities are
80                 ['stage', 'xmomentum', 'ymomentum']
81                 
82
83    interpolation_points - list of absolute UTM coordinates for points (N x 2)
84    or geospatial object or points file name at which values are sought
85
86    time_thinning -
87
88    verbose -
89
90    use_cache: True means that caching of intermediate result of
91               Interpolation_function is attempted
92
93    boundary_polygon -
94
95   
96    See Interpolation function in anuga.fit_interpolate.interpolation for
97    further documentation
98    """
99
100    # FIXME (OLE): Should check origin of domain against that of file
101    # In fact, this is where origin should be converted to that of domain
102    # Also, check that file covers domain fully.
103
104    # Take into account:
105    # - domain's georef
106    # - sww file's georef
107    # - interpolation points as absolute UTM coordinates
108
109    if quantities is None:
110        if verbose:
111            msg = 'Quantities specified in file_function are None,'
112            msg += ' so I will use stage, xmomentum, and ymomentum in that order'
113            print msg
114        quantities = ['stage', 'xmomentum', 'ymomentum']
115
116    # Use domain's startime if available
117    if domain is not None:   
118        domain_starttime = domain.get_starttime()
119    else:
120        domain_starttime = None
121
122    # Build arguments and keyword arguments for use with caching or apply.
123    args = (filename,)
124
125    # FIXME (Ole): Caching this function will not work well
126    # if domain is passed in as instances change hash code.
127    # Instead we pass in those attributes that are needed (and return them
128    # if modified)
129    kwargs = {'quantities': quantities,
130              'interpolation_points': interpolation_points,
131              'domain_starttime': domain_starttime,
132              'time_thinning': time_thinning,                   
133              'verbose': verbose,
134              'boundary_polygon': boundary_polygon}
135
136    # Call underlying engine with or without caching
137    if use_cache is True:
138        try:
139            from caching import cache
140        except:
141            msg = 'Caching was requested, but caching module'+\
142                  'could not be imported'
143            raise msg
144
145        f, starttime = cache(_file_function,
146                             args, kwargs,
147                             dependencies=[filename],
148                             compression=False,                 
149                             verbose=verbose)
150    else:
151        f, starttime = apply(_file_function,
152                             args, kwargs)
153
154    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
155    #structure
156
157    f.starttime = starttime
158    f.filename = filename
159   
160    if domain is not None:
161        #Update domain.startime if it is *earlier* than starttime from file
162        if starttime > domain.starttime:
163            msg = 'WARNING: Start time as specified in domain (%f)'\
164                  %domain.starttime
165            msg += ' is earlier than the starttime of file %s (%f).'\
166                     %(filename, starttime)
167            msg += ' Modifying domain starttime accordingly.'
168           
169            if verbose: print msg
170
171            domain.set_starttime(starttime) #Modifying model time
172
173            if verbose: print 'Domain starttime is now set to %f'\
174               %domain.starttime
175    return f
176
177
178##
179# @brief ??
180# @param filename  Name of .sww or .tms file.
181# @param domain Associated domain object.
182# @param quantities Name of quantity to be interpolated or a list of names.
183# @param interpolation_points List of absolute UTM coordinates for points
184#                             (N x 2) or geospatial object or
185#                             points file name at which values are sought.
186# @param time_thinning
187# @param verbose True if this function is to be verbose.
188# @param use_cache True means that caching of intermediate result is attempted.
189# @param boundary_polygon
190def _file_function(filename,
191                   quantities=None,
192                   interpolation_points=None,
193                   domain_starttime=None,
194                   time_thinning=1,
195                   verbose=False,
196                   boundary_polygon=None):
197    """Internal function
198   
199    See file_function for documentatiton
200    """
201
202    assert type(filename) == type(''),\
203               'First argument to File_function must be a string'
204
205    try:
206        fid = open(filename)
207    except Exception, e:
208        msg = 'File "%s" could not be opened: Error="%s"' % (filename, e)
209        raise msg
210
211    # read first line of file, guess file type
212    line = fid.readline()
213    fid.close()
214
215    if line[:3] == 'CDF':
216        return get_netcdf_file_function(filename,
217                                        quantities,
218                                        interpolation_points,
219                                        domain_starttime,
220                                        time_thinning=time_thinning,
221                                        verbose=verbose,
222                                        boundary_polygon=boundary_polygon)
223    else:
224        # FIXME (Ole): Could add csv file here to address Ted Rigby's
225        # suggestion about reading hydrographs.
226        # This may also deal with the gist of ticket:289
227        raise 'Must be a NetCDF File'
228
229
230##
231# @brief ??
232# @param filename  Name of .sww or .tms file.
233# @param quantity_names Name of quantity to be interpolated or a list of names.
234# @param interpolation_points List of absolute UTM coordinates for points
235#                             (N x 2) or geospatial object or
236#                             points file name at which values are sought.
237# @param domain_starttime Start time from domain object.
238# @param time_thinning ??
239# @param verbose True if this function is to be verbose.
240# @param boundary_polygon ??
241# @return A callable object.
242def get_netcdf_file_function(filename,
243                             quantity_names=None,
244                             interpolation_points=None,
245                             domain_starttime=None,                           
246                             time_thinning=1,                             
247                             verbose=False,
248                             boundary_polygon=None):
249    """Read time history of spatial data from NetCDF sww file and
250    return a callable object f(t,x,y)
251    which will return interpolated values based on the input file.
252
253    Model time (domain_starttime)
254    will be checked, possibly modified and returned
255   
256    All times are assumed to be in UTC
257
258    See Interpolation function for further documetation
259    """
260
261    # FIXME: Check that model origin is the same as file's origin
262    # (both in UTM coordinates)
263    # If not - modify those from file to match domain
264    # (origin should be passed in)
265    # Take this code from e.g. dem2pts in data_manager.py
266    # FIXME: Use geo_reference to read and write xllcorner...
267
268    import time, calendar, types
269    from anuga.config import time_format
270    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
271
272    # Open NetCDF file
273    if verbose: print 'Reading', filename
274
275    fid = NetCDFFile(filename, 'r')
276
277    if type(quantity_names) == types.StringType:
278        quantity_names = [quantity_names]       
279
280    if quantity_names is None or len(quantity_names) < 1:
281        msg = 'No quantities are specified in file_function'
282        raise Exception, msg
283 
284    if interpolation_points is not None:
285        interpolation_points = ensure_absolute(interpolation_points)
286        msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1]
287        assert interpolation_points.shape[1] == 2, msg
288
289    # Now assert that requested quantitites (and the independent ones)
290    # are present in file
291    missing = []
292    for quantity in ['time'] + quantity_names:
293        if not fid.variables.has_key(quantity):
294            missing.append(quantity)
295
296    if len(missing) > 0:
297        msg = 'Quantities %s could not be found in file %s'\
298              % (str(missing), filename)
299        fid.close()
300        raise Exception, msg
301
302    # Decide whether this data has a spatial dimension
303    spatial = True
304    for quantity in ['x', 'y']:
305        if not fid.variables.has_key(quantity):
306            spatial = False
307
308    if filename[-3:] == 'tms' and spatial is True:
309        msg = 'Files of type tms must not contain spatial  information'
310        raise msg
311
312    if filename[-3:] == 'sww' and spatial is False:
313        msg = 'Files of type sww must contain spatial information'       
314        raise msg
315
316    if filename[-3:] == 'sts' and spatial is False:
317        #What if mux file only contains one point
318        msg = 'Files of type sts must contain spatial information'       
319        raise msg
320
321    if filename[-3:] == 'sts' and boundary_polygon is None:
322        #What if mux file only contains one point
323        msg = 'Files of type sts require boundary polygon'       
324        raise msg
325
326    # Get first timestep
327    try:
328        starttime = fid.starttime[0]
329    except ValueError:
330        msg = 'Could not read starttime from file %s' %filename
331        raise msg
332
333    # Get variables
334    # if verbose: print 'Get variables'   
335    time = fid.variables['time'][:]   
336
337    # Get time independent stuff
338    if spatial:
339        # Get origin
340        xllcorner = fid.xllcorner[0]
341        yllcorner = fid.yllcorner[0]
342        zone = fid.zone[0]       
343
344        x = fid.variables['x'][:]
345        y = fid.variables['y'][:]
346        if filename[-3:] == 'sww':
347            triangles = fid.variables['volumes'][:]
348
349        x = reshape(x, (len(x),1))
350        y = reshape(y, (len(y),1))
351        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
352
353        if boundary_polygon is not None:
354            #removes sts points that do not lie on boundary
355            boundary_polygon=ensure_numeric(boundary_polygon)
356            boundary_polygon[:,0] -= xllcorner
357            boundary_polygon[:,1] -= yllcorner
358            temp=[]
359            boundary_id=[]
360            gauge_id=[]
361            for i in range(len(boundary_polygon)):
362                for j in range(len(x)):
363                    if allclose(vertex_coordinates[j],boundary_polygon[i],1e-4):
364                        #FIX ME:
365                        #currently gauges lat and long is stored as float and
366                        #then cast to double. This cuases slight repositioning
367                        #of vertex_coordinates.
368                        temp.append(boundary_polygon[i])
369                        gauge_id.append(j)
370                        boundary_id.append(i)
371                        break
372            gauge_neighbour_id=[]
373            for i in range(len(boundary_id)-1):
374                if boundary_id[i]+1==boundary_id[i+1]:
375                    gauge_neighbour_id.append(i+1)
376                else:
377                    gauge_neighbour_id.append(-1)
378            if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 \
379               and boundary_id[0]==0:
380                gauge_neighbour_id.append(0)
381            else:
382                gauge_neighbour_id.append(-1)
383            gauge_neighbour_id=ensure_numeric(gauge_neighbour_id)
384
385            if len(compress(gauge_neighbour_id>=0,gauge_neighbour_id)) \
386               != len(temp)-1:
387                msg='incorrect number of segments'
388                raise msg
389            vertex_coordinates=ensure_numeric(temp)
390            if len(vertex_coordinates)==0:
391                msg = 'None of the sts gauges fall on the boundary'
392                raise msg
393        else:
394            gauge_neighbour_id=None
395
396        if interpolation_points is not None:
397            # Adjust for georef
398            interpolation_points[:,0] -= xllcorner
399            interpolation_points[:,1] -= yllcorner       
400    else:
401        gauge_neighbour_id=None
402       
403    if domain_starttime is not None:
404        # If domain_startime is *later* than starttime,
405        # move time back - relative to domain's time
406        if domain_starttime > starttime:
407            time = time - domain_starttime + starttime
408
409        # FIXME Use method in geo to reconcile
410        # if spatial:
411        # assert domain.geo_reference.xllcorner == xllcorner
412        # assert domain.geo_reference.yllcorner == yllcorner
413        # assert domain.geo_reference.zone == zone       
414       
415    if verbose:
416        print 'File_function data obtained from: %s' %filename
417        print '  References:'
418        #print '    Datum ....' #FIXME
419        if spatial:
420            print '    Lower left corner: [%f, %f]'\
421                  %(xllcorner, yllcorner)
422        print '    Start time:   %f' %starttime               
423       
424   
425    # Produce values for desired data points at
426    # each timestep for each quantity
427    quantities = {}
428    for i, name in enumerate(quantity_names):
429        quantities[name] = fid.variables[name][:]
430        if boundary_polygon is not None:
431            #removes sts points that do not lie on boundary
432            quantities[name] = take(quantities[name],gauge_id,1)
433           
434    # Close sww, tms or sts netcdf file         
435    fid.close()
436
437    from anuga.fit_interpolate.interpolate import Interpolation_function
438
439    if not spatial:
440        vertex_coordinates = triangles = interpolation_points = None
441    if filename[-3:] == 'sts':#added
442        triangles = None
443        #vertex coordinates is position of urs gauges
444
445    # Return Interpolation_function instance as well as
446    # starttime for use to possible modify that of domain
447    return (Interpolation_function(time,
448                                   quantities,
449                                   quantity_names,
450                                   vertex_coordinates,
451                                   triangles,
452                                   interpolation_points,
453                                   time_thinning=time_thinning,
454                                   verbose=verbose,
455                                   gauge_neighbour_id=gauge_neighbour_id),
456            starttime)
457
458    # NOTE (Ole): Caching Interpolation function is too slow as
459    # the very long parameters need to be hashed.
460
461
462##
463# @brief Replace multiple substrings in a string.
464# @param text The string to operate on.
465# @param dictionary A dict containing replacements, key->value.
466# @return The new string.
467def multiple_replace(text, dictionary):
468    """Multiple replace of words in text
469
470    text:       String to be modified
471    dictionary: Mapping of words that are to be substituted
472
473    Python Cookbook 3.14 page 88 and page 90
474    http://code.activestate.com/recipes/81330/
475    """
476
477    import re
478   
479    #Create a regular expression from all of the dictionary keys
480    #matching only entire words
481    regex = re.compile(r'\b'+ \
482                       r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \
483                       r'\b' )
484
485    #For each match, lookup the corresponding value in the dictionary
486    return regex.sub(lambda match: dictionary[match.group(0)], text)
487
488
489##
490# @brief Apply arbitrary expressions to the values of a dict.
491# @param expression A string expression to apply.
492# @param dictionary The dictionary to apply the expression to.
493def apply_expression_to_dictionary(expression, dictionary):
494    """Apply arbitrary expression to values of dictionary
495
496    Given an expression in terms of the keys, replace key by the
497    corresponding values and evaluate.
498
499    expression: Arbitrary, e.g. arithmetric, expression relating keys
500                from dictionary.
501
502    dictionary: Mapping between symbols in expression and objects that
503                will be evaluated by expression.
504                Values in dictionary must support operators given in
505                expression e.g. by overloading
506
507    Due to a limitation with Numeric, this can not evaluate 0/0
508    In general, the user can fix by adding 1e-30 to the numerator.
509    SciPy core can handle this situation.
510    """
511
512    import types
513    import re
514
515    assert isinstance(expression, basestring)
516    assert type(dictionary) == types.DictType
517
518    #Convert dictionary values to textual representations suitable for eval
519    D = {}
520    for key in dictionary:
521        D[key] = 'dictionary["%s"]' % key
522
523    #Perform substitution of variables   
524    expression = multiple_replace(expression, D)
525
526    #Evaluate and return
527    try:
528        return eval(expression)
529    except NameError, e:
530        msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)
531        raise NameError, msg
532    except ValueError, e:
533        msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)
534        raise ValueError, msg
535
536
537##
538# @brief Format a float into a string.
539# @param value Float value to format.
540# @param format The format to use (%.2f is default).
541# @return The formatted float as a string.
542def get_textual_float(value, format = '%.2f'):
543    """Get textual representation of floating point numbers
544    and accept None as valid entry
545
546    format is a string - default = '%.2f'
547    """
548
549    if value is None:
550        return 'None'
551    else:
552        try:
553            float(value)
554        except:
555            # May this is a vector
556            if len(value) > 1:
557                s = '('
558                for v in value:
559                    s += get_textual_float(v, format) + ', '
560                   
561                s = s[:-2] + ')' # Strip trailing comma and close
562                return s
563            else:
564                raise 'Illegal input to get_textual_float:', value
565        else:
566            return format % float(value)
567
568
569#################################################################################
570# OBSOLETE STUFF
571#################################################################################
572
573# @note TEMP
574def angle(v1, v2):
575    """Temporary Interface to new location"""
576
577    import anuga.utilities.numerical_tools as NT       
578   
579    msg = 'angle has moved from util.py.  '
580    msg += 'Please use "from anuga.utilities.numerical_tools import angle"'
581    warn(msg, DeprecationWarning) 
582
583    return NT.angle(v1, v2)
584
585   
586# @note TEMP
587def anglediff(v0, v1):
588    """Temporary Interface to new location"""
589
590    import anuga.utilities.numerical_tools as NT
591   
592    msg = 'anglediff has moved from util.py.  '
593    msg += 'Please use "from anuga.utilities.numerical_tools import anglediff"'
594    warn(msg, DeprecationWarning) 
595
596    return NT.anglediff(v0, v1)   
597
598   
599# @note TEMP
600def mean(x):
601    """Temporary Interface to new location"""
602
603    import anuga.utilities.numerical_tools as NT   
604   
605    msg = 'mean has moved from util.py.  '
606    msg += 'Please use "from anuga.utilities.numerical_tools import mean"'
607    warn(msg, DeprecationWarning) 
608
609    return NT.mean(x)   
610
611
612# @note TEMP
613def point_on_line(*args, **kwargs):
614    """Temporary Interface to new location"""
615
616    msg = 'point_on_line has moved from util.py.  '
617    msg += 'Please use "from anuga.utilities.polygon import point_on_line"'
618    warn(msg, DeprecationWarning) 
619
620    return utilities.polygon.point_on_line(*args, **kwargs)     
621
622   
623# @note TEMP
624def inside_polygon(*args, **kwargs):
625    """Temporary Interface to new location"""
626
627    print 'inside_polygon has moved from util.py.  ',
628    print 'Please use "from anuga.utilities.polygon import inside_polygon"'
629
630    return utilities.polygon.inside_polygon(*args, **kwargs)   
631
632   
633# @note TEMP
634def outside_polygon(*args, **kwargs):
635    """Temporary Interface to new location"""
636
637    print 'outside_polygon has moved from util.py.  ',
638    print 'Please use "from anuga.utilities.polygon import outside_polygon"'
639
640    return utilities.polygon.outside_polygon(*args, **kwargs)   
641
642
643# @note TEMP
644def separate_points_by_polygon(*args, **kwargs):
645    """Temporary Interface to new location"""
646
647    print 'separate_points_by_polygon has moved from util.py.  ',
648    print 'Please use "from anuga.utilities.polygon import ' \
649          'separate_points_by_polygon"'
650
651    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)   
652
653
654# @note TEMP
655def read_polygon(*args, **kwargs):
656    """Temporary Interface to new location"""
657
658    print 'read_polygon has moved from util.py.  ',
659    print 'Please use "from anuga.utilities.polygon import read_polygon"'
660
661    return utilities.polygon.read_polygon(*args, **kwargs)   
662
663
664# @note TEMP
665def populate_polygon(*args, **kwargs):
666    """Temporary Interface to new location"""
667
668    print 'populate_polygon has moved from util.py.  ',
669    print 'Please use "from anuga.utilities.polygon import populate_polygon"'
670
671    return utilities.polygon.populate_polygon(*args, **kwargs)   
672
673
674#################################################################################
675# End of obsolete stuff ?
676#################################################################################
677
678# @note TEMP
679def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='',
680                         verbose=False):
681    """Temporary Interface to new location"""
682    from anuga.shallow_water.data_manager import start_screen_catcher \
683         as dm_start_screen_catcher
684
685    print 'start_screen_catcher has moved from util.py.  ',
686    print 'Please use "from anuga.shallow_water.data_manager import ' \
687          'start_screen_catcher"'
688   
689    return dm_start_screen_catcher(dir_name, myid='', numprocs='',
690                                   extra_info='', verbose=False)
691
692
693##
694# @brief Read a .sww file and plot the time series.
695# @param swwfiles Dictionary of .sww files.
696# @param gauge_filename Name of gauge file.
697# @param production_dirs ??
698# @param report If True, write figures to report directory.
699# @param reportname Name of generated report (if any).
700# @param plot_quantity List containing quantities to plot.
701# @param generate_fig If True, generate figures as well as CSV files.
702# @param surface If True, then generate solution surface with 3d plot.
703# @param time_min Beginning of user defined time range for plotting purposes.
704# @param time_max End of user defined time range for plotting purposes.
705# @param time_thinning ??
706# @param time_unit ??
707# @param title_on If True, export standard graphics with title.
708# @param use_cache If True, use caching.
709# @param verbose If True, this function is verbose.
710def sww2timeseries(swwfiles,
711                   gauge_filename,
712                   production_dirs,
713                   report=None,
714                   reportname=None,
715                   plot_quantity=None,
716                   generate_fig=False,
717                   surface=None,
718                   time_min=None,
719                   time_max=None,
720                   time_thinning=1,                   
721                   time_unit=None,
722                   title_on=None,
723                   use_cache=False,
724                   verbose=False):
725    """ Read sww file and plot the time series for the
726    prescribed quantities at defined gauge locations and
727    prescribed time range.
728
729    Input variables:
730
731    swwfiles        - dictionary of sww files with label_ids (used in
732                      generating latex output. It will be part of
733                      the directory name of file_loc (typically the timestamp).
734                      Helps to differentiate latex files for different
735                      simulations for a particular scenario. 
736                    - assume that all conserved quantities have been stored
737                    - assume each sww file has been simulated with same timestep
738   
739    gauge_filename  - name of file containing gauge data
740                        - easting, northing, name , elevation?
741                    - OR (this is not yet done)
742                        - structure which can be converted to a Numeric array,
743                          such as a geospatial data object
744                     
745    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
746                                                'boundaries': 'urs boundary'}
747                      this will use the second part as the label and the
748                      first part as the ?
749                      #FIXME: Is it a list or a dictionary
750                      # This is probably obsolete by now
751                     
752    report          - if True, then write figures to report_figures directory in
753                      relevant production directory
754                    - if False, figures are already stored with sww file
755                    - default to False
756
757    reportname      - name for report if wishing to generate report
758   
759    plot_quantity   - list containing quantities to plot, they must
760                      be the name of an existing quantity or one of
761                      the following possibilities
762                    - possibilities:
763                        - stage; 'stage'
764                        - depth; 'depth'
765                        - speed; calculated as absolute momentum
766                         (pointwise) divided by depth; 'speed'
767                        - bearing; calculated as the angle of the momentum
768                          vector (xmomentum, ymomentum) from the North; 'bearing'
769                        - absolute momentum; calculated as
770                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
771                        - x momentum; 'xmomentum'
772                        - y momentum; 'ymomentum'
773                    - default will be ['stage', 'speed', 'bearing']
774
775    generate_fig     - if True, generate figures as well as csv file
776                     - if False, csv files created only
777                     
778    surface          - if True, then generate solution surface with 3d plot
779                       and save to current working directory
780                     - default = False
781   
782    time_min         - beginning of user defined time range for plotting purposes
783                        - default will be first available time found in swwfile
784                       
785    time_max         - end of user defined time range for plotting purposes
786                        - default will be last available time found in swwfile
787                       
788    title_on        - if True, export standard graphics with title
789                    - if False, export standard graphics without title
790
791
792    Output:
793   
794    - time series data stored in .csv for later use if required.
795      Name = gauges_timeseries followed by gauge name
796    - latex file will be generated in same directory as where script is
797      run (usually production scenario directory.
798      Name = latexoutputlabel_id.tex
799
800    Other important information:
801   
802    It is assumed that the used has stored all the conserved quantities
803    and elevation during the scenario run, i.e.
804    ['stage', 'elevation', 'xmomentum', 'ymomentum']
805    If this has not occurred then sww2timeseries will not work.
806
807
808    Usage example
809    texname = sww2timeseries({project.boundary_name + '.sww': ''},
810                             project.polygons_dir + sep + 'boundary_extent.csv',
811                             project.anuga_dir,
812                             report = False,
813                             plot_quantity = ['stage', 'speed', 'bearing'],
814                             time_min = None,
815                             time_max = None,
816                             title_on = True,   
817                             verbose = True)
818   
819    """
820
821    msg = 'NOTE: A new function is available to create csv files from sww '
822    msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
823    msg += ' PLUS another new function to create graphs from csv files called '
824    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
825    print msg
826   
827    k = _sww2timeseries(swwfiles,
828                        gauge_filename,
829                        production_dirs,
830                        report,
831                        reportname,
832                        plot_quantity,
833                        generate_fig,
834                        surface,
835                        time_min,
836                        time_max,
837                        time_thinning,                       
838                        time_unit,
839                        title_on,
840                        use_cache,
841                        verbose)
842    return k
843
844
845##
846# @brief Read a .sww file and plot the time series.
847# @param swwfiles Dictionary of .sww files.
848# @param gauge_filename Name of gauge file.
849# @param production_dirs ??
850# @param report If True, write figures to report directory.
851# @param reportname Name of generated report (if any).
852# @param plot_quantity List containing quantities to plot.
853# @param generate_fig If True, generate figures as well as CSV files.
854# @param surface If True, then generate solution surface with 3d plot.
855# @param time_min Beginning of user defined time range for plotting purposes.
856# @param time_max End of user defined time range for plotting purposes.
857# @param time_thinning ??
858# @param time_unit ??
859# @param title_on If True, export standard graphics with title.
860# @param use_cache If True, use caching.
861# @param verbose If True, this function is verbose.
862def _sww2timeseries(swwfiles,
863                    gauge_filename,
864                    production_dirs,
865                    report = None,
866                    reportname = None,
867                    plot_quantity = None,
868                    generate_fig = False,
869                    surface = None,
870                    time_min = None,
871                    time_max = None,
872                    time_thinning = 1,                   
873                    time_unit = None,
874                    title_on = None,
875                    use_cache = False,
876                    verbose = False):   
877       
878    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
879   
880    try:
881        fid = open(gauge_filename)
882    except Exception, e:
883        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
884        raise msg
885
886    if report is None:
887        report = False
888       
889    if plot_quantity is None:
890        plot_quantity = ['depth', 'speed']
891    else:
892        assert type(plot_quantity) == list, 'plot_quantity must be a list'
893        check_list(plot_quantity)
894
895    if surface is None:
896        surface = False
897
898    if time_unit is None:
899        time_unit = 'hours'
900   
901    if title_on is None:
902        title_on = True
903   
904    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
905
906    gauges, locations, elev = get_gauges_from_file(gauge_filename)
907
908    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
909
910    file_loc = []
911    f_list = []
912    label_id = []
913    leg_label = []
914    themaxT = 0.0
915    theminT = 0.0
916
917    for swwfile in swwfiles.keys():
918        try:
919            fid = open(swwfile)
920        except Exception, e:
921            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
922            raise msg
923
924        print 'swwfile', swwfile
925
926        # Extract parent dir name and use as label
927        path, _ = os.path.split(swwfile)
928        _, label = os.path.split(path)       
929       
930        #print 'label', label
931        leg_label.append(label)
932
933        f = file_function(swwfile,
934                          quantities = sww_quantity,
935                          interpolation_points = gauges,
936                          time_thinning = time_thinning,
937                          verbose = verbose,
938                          use_cache = use_cache)
939
940        # determine which gauges are contained in sww file
941        count = 0
942        gauge_index = []
943        for k, g in enumerate(gauges):
944            if f(0.0, point_id = k)[2] > 1.0e6:
945                count += 1
946                if count == 1: print 'Gauges not contained here:'
947                print locations[k]
948            else:
949                gauge_index.append(k)
950
951        if len(gauge_index) > 0:
952            print 'Gauges contained here: \n',
953        else:
954            print 'No gauges contained here. \n'
955        for i in range(len(gauge_index)):
956             print locations[gauge_index[i]]
957             
958        index = swwfile.rfind(sep)
959        file_loc.append(swwfile[:index+1])
960        label_id.append(swwfiles[swwfile])
961       
962        f_list.append(f)
963        maxT = max(f.get_time())
964        minT = min(f.get_time())
965        if maxT > themaxT: themaxT = maxT
966        if minT > theminT: theminT = minT
967
968    if time_min is None:
969        time_min = theminT # min(T)
970    else:
971        if time_min < theminT: # min(T):
972            msg = 'Minimum time entered not correct - please try again'
973            raise Exception, msg
974
975    if time_max is None:
976        time_max = themaxT # max(T)
977    else:
978        if time_max > themaxT: # max(T):
979            msg = 'Maximum time entered not correct - please try again'
980            raise Exception, msg
981
982    if verbose and len(gauge_index) > 0:
983         print 'Inputs OK - going to generate figures'
984
985    if len(gauge_index) <> 0:
986        texfile, elev_output = \
987            generate_figures(plot_quantity, file_loc, report, reportname,
988                             surface, leg_label, f_list, gauges, locations,
989                             elev, gauge_index, production_dirs, time_min,
990                             time_max, time_unit, title_on, label_id,
991                             generate_fig, verbose)
992    else:
993        texfile = ''
994        elev_output = []
995
996    return texfile, elev_output
997
998
999##
1000# @brief Read gauge info from a file.
1001# @param filename The name of the file to read.
1002# @return A (gauges, gaugelocation, elev) tuple.
1003def get_gauges_from_file(filename):
1004    """ Read in gauge information from file
1005    """
1006
1007    from os import sep, getcwd, access, F_OK, mkdir
1008
1009    # Get data from the gauge file
1010    fid = open(filename)
1011    lines = fid.readlines()
1012    fid.close()
1013   
1014    gauges = []
1015    gaugelocation = []
1016    elev = []
1017
1018    # Check header information   
1019    line1 = lines[0]
1020    line11 = line1.split(',')
1021
1022    if isinstance(line11[0], str) is True:
1023        # We have found text in the first line
1024        east_index = None
1025        north_index = None
1026        name_index = None
1027        elev_index = None
1028
1029        for i in range(len(line11)):
1030            if line11[i].strip().lower() == 'easting':   east_index = i
1031            if line11[i].strip().lower() == 'northing':  north_index = i
1032            if line11[i].strip().lower() == 'name':      name_index = i
1033            if line11[i].strip().lower() == 'elevation': elev_index = i
1034
1035        if east_index < len(line11) and north_index < len(line11):
1036            pass
1037        else:
1038            msg = 'WARNING: %s does not contain correct header information' \
1039                  % filename
1040            msg += 'The header must be: easting, northing, name, elevation'
1041            raise Exception, msg
1042
1043        if elev_index is None: 
1044            raise Exception
1045   
1046        if name_index is None: 
1047            raise Exception
1048
1049        lines = lines[1:] # Remove header from data
1050    else:
1051        # No header, assume that this is a simple easting, northing file
1052
1053        msg = 'There was no header in file %s and the number of columns is %d' \
1054              % (filename, len(line11))
1055        msg += '- was assuming two columns corresponding to Easting and Northing'
1056        assert len(line11) == 2, msg
1057
1058        east_index = 0
1059        north_index = 1
1060
1061        N = len(lines)
1062        elev = [-9999]*N
1063        gaugelocation = range(N)
1064       
1065    # Read in gauge data
1066    for line in lines:
1067        fields = line.split(',')
1068
1069        gauges.append([float(fields[east_index]), float(fields[north_index])])
1070
1071        if len(fields) > 2:
1072            elev.append(float(fields[elev_index]))
1073            loc = fields[name_index]
1074            gaugelocation.append(loc.strip('\n'))
1075
1076    return gauges, gaugelocation, elev
1077
1078
1079##
1080# @brief Check that input quantities in quantity list are legal.
1081# @param quantity Quantity list to check.
1082# @note Raises an exception of list is not legal.
1083def check_list(quantity):
1084    """ Check that input quantities in quantity list are possible
1085    """
1086    import sys
1087    from sets import Set as set
1088
1089    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
1090                    'ymomentum', 'speed', 'bearing', 'elevation']
1091
1092    # convert all quanitiy names to lowercase
1093    for i,j in enumerate(quantity):
1094        quantity[i] = quantity[i].lower()
1095
1096    # check that all names in 'quantity' appear in 'all_quantity'
1097    p = list(set(quantity).difference(set(all_quantity)))
1098    if len(p) != 0:
1099        msg = 'Quantities %s do not exist - please try again' %p
1100        raise Exception, msg
1101
1102
1103##
1104# @brief Calculate velocity bearing from North.
1105# @param uh ??
1106# @param vh ??
1107# @return The calculated bearing.
1108def calc_bearing(uh, vh):
1109    """ Calculate velocity bearing from North
1110    """
1111    #FIXME (Ole): I reckon we should refactor this one to use
1112    #             the function angle() in utilities/numerical_tools
1113    #
1114    #             It will be a simple matter of
1115    # * converting from radians to degrees
1116    # * moving the reference direction from [1,0] to North
1117    # * changing from counter clockwise to clocwise.
1118       
1119    angle = degrees(atan(vh/(uh+1.e-15)))
1120
1121    if (0 < angle < 90.0):
1122        if vh > 0:
1123            bearing = 90.0 - abs(angle)
1124        if vh < 0:
1125            bearing = 270.0 - abs(angle)
1126   
1127    if (-90 < angle < 0):
1128        if vh < 0:
1129            bearing = 90.0 - (angle)
1130        if vh > 0:
1131            bearing = 270.0 - (angle)
1132    if angle == 0: bearing = 0.0
1133
1134    return bearing
1135
1136
1137##
1138# @brief Generate figures from quantities and gauges for each sww file.
1139# @param plot_quantity  ??
1140# @param file_loc ??
1141# @param report ??
1142# @param reportname ??
1143# @param surface ??
1144# @param leg_label ??
1145# @param f_list ??
1146# @param gauges ??
1147# @param locations ??
1148# @param elev ??
1149# @param gauge_index ??
1150# @param production_dirs ??
1151# @param time_min ??
1152# @param time_max ??
1153# @param time_unit ??
1154# @param title_on ??
1155# @param label_id ??
1156# @param generate_fig ??
1157# @param verbose??
1158# @return (texfile2, elev_output)
1159def generate_figures(plot_quantity, file_loc, report, reportname, surface,
1160                     leg_label, f_list, gauges, locations, elev, gauge_index,
1161                     production_dirs, time_min, time_max, time_unit,
1162                     title_on, label_id, generate_fig, verbose):
1163    """ Generate figures based on required quantities and gauges for
1164    each sww file
1165    """
1166    from Numeric import ones, allclose, zeros, Float, ravel
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 = zeros((n0, m, p), Float) 
1210    stages = zeros((n0, m, p), Float)
1211    elevations = zeros((n0, m, p), Float) 
1212    momenta = zeros((n0, m, p), Float)
1213    xmom = zeros((n0, m, p), Float)
1214    ymom = zeros((n0, m, p), Float)
1215    speed = zeros((n0, m, p), Float)
1216    bearings = zeros((n0, m, p), Float)
1217    due_east = 90.0*ones((n0, 1), Float)
1218    due_west = 270.0*ones((n0, 1), Float)
1219    depths = zeros((n0, m, p), Float)
1220    eastings = zeros((n0, m, p), 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 = zeros((n0, m), Float)
1233    stages_plot3d = zeros((n0, m), Float)
1234    eastings_plot3d = zeros((n0, m),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(ravel(eastings[:,:,j]),
1343                              ravel(model_time[:,:,j]),
1344                              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=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=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 - first '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 = choose(triangles,loners)
1768        #print "triangles after", triangles
1769    return verts, triangles
1770
1771 
1772def get_centroid_values(x, triangles):
1773    """Compute centroid values from vertex values
1774   
1775    x: Values at vertices of triangular mesh
1776    triangles: Nx3 integer array pointing to vertex information
1777    for each of the N triangels. Elements of triangles are
1778    indices into x
1779    """
1780
1781       
1782    xc = zeros(triangles.shape[0], Float) # Space for centroid info
1783   
1784    for k in range(triangles.shape[0]):
1785        # Indices of vertices
1786        i0 = triangles[k][0]
1787        i1 = triangles[k][1]
1788        i2 = triangles[k][2]       
1789       
1790        xc[k] = (x[i0] + x[i1] + x[i2])/3
1791
1792
1793    return xc
1794
1795# @note TEMP
1796def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]},
1797                                output_dir='',
1798                                base_name='',
1799                                plot_numbers=['3-5'],
1800                                quantities=['speed','stage','momentum'],
1801                                assess_all_csv_files=True,
1802                                extra_plot_name='test' ):
1803
1804    msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs ',
1805    msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import csv2timeseries_graphs"'
1806
1807    raise Exception, msg
1808    return csv2timeseries_graphs(directories_dic,
1809                                 output_dir,
1810                                 base_name,
1811                                 plot_numbers,
1812                                 quantities,
1813                                 extra_plot_name,
1814                                 assess_all_csv_files
1815                                 )
1816   
1817def csv2timeseries_graphs(directories_dic={},
1818                            output_dir='',
1819                            base_name=None,
1820                            plot_numbers='',
1821                            quantities=['stage'],
1822                            extra_plot_name='',
1823                            assess_all_csv_files=True,                           
1824                            create_latex=False,
1825                            verbose=False):
1826                               
1827    """
1828    Read in csv files that have the right header information and
1829    plot time series such as Stage, Speed, etc. Will also plot several
1830    time series on one plot. Filenames must follow this convention,
1831    <base_name><plot_number>.csv eg gauge_timeseries3.csv
1832   
1833    NOTE: relies that 'elevation' is in the csv file!
1834
1835    Each file represents a location and within each file there are
1836    time, quantity columns.
1837   
1838    For example:   
1839    if "directories_dic" defines 4 directories and in each directories
1840    there is a csv files corresponding to the right "plot_numbers",
1841    this will create a plot with 4 lines one for each directory AND
1842    one plot for each "quantities".  ??? FIXME: unclear.
1843   
1844    Usage:
1845        csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1846                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1847                            output_dir='fixed_wave'+sep,
1848                            base_name='gauge_timeseries_',
1849                            plot_numbers='',
1850                            quantities=['stage','speed'],
1851                            extra_plot_name='',
1852                            assess_all_csv_files=True,                           
1853                            create_latex=False,
1854                            verbose=True)
1855            this will create one plot for stage with both 'slide' and
1856            'fixed_wave' lines on it for stage and speed for each csv
1857            file with 'gauge_timeseries_' as the prefix. The graghs
1858            will be in the output directory 'fixed_wave' and the graph
1859            axis will be determined by assessing all the
1860   
1861    ANOTHER EXAMPLE
1862        new_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=['1-3'],
1867                            quantities=['stage','speed'],
1868                            extra_plot_name='',
1869                            assess_all_csv_files=False,                           
1870                            create_latex=False,
1871                            verbose=True)
1872        This will plot csv files called gauge_timeseries_1.csv and
1873        gauge_timeseries3.csv from both 'slide' and 'fixed_wave' directories
1874        to 'fixed_wave'. There will be 4 plots created two speed and two stage
1875        one for each csv file. There will be two lines on each of these plots.
1876        And the axis will have been determined from only these files, had
1877        assess_all_csv_files = True all csv file with 'gauges_timeseries_' prefix
1878        would of been assessed.
1879   
1880    ANOTHER EXAMPLE   
1881         csv2timeseries_graphs({'J:'+sep+'anuga_validation'+sep:['new',20,-.1],
1882                                   'J:'+sep+'conical_island'+sep:['test',0,0]},
1883                                   output_dir='',
1884                                   plot_numbers=['1','3'],
1885                                   quantities=['stage','depth','bearing'],
1886                                   base_name='gauge_b',
1887                                   assess_all_csv_files=True,
1888                                  verbose=True)   
1889       
1890            This will produce one plot for each quantity (therefore 3) in the current directory,
1891            each plot will have 2 lines on them. The first plot named 'new' will have the time
1892            offseted by 20secs and the stage height adjusted by -0.1m
1893       
1894    Inputs:
1895        directories_dic: dictionary of directory with values (plot
1896                         legend name for directory), (start time of
1897                         the time series) and the (value to add to
1898                         stage if needed). For example
1899                        {dir1:['Anuga_ons',5000, 0],
1900                         dir2:['b_emoth',5000,1.5],
1901                         dir3:['b_ons',5000,1.5]}
1902                         Having multiple directories defined will plot them on
1903                         one plot, therefore there will be 3 lines on each of these
1904                         plot. If you only want one line per plot call csv2timeseries_graph
1905                         separately for each directory, eg only have one directory in the
1906                         'directories_dic' in each call.
1907                         
1908        output_dir: directory for the plot outputs. Only important to define
1909                    when you have more than one directory in your directories_dic,
1910                    if you have not defined it and you have multiple directories in
1911                    'directories_dic' there will be plots in each directory however only
1912                    one directory will contain the complete plot/graphs.
1913       
1914        base_name: Is used a couple of times. 1) to find the csv files to be plotted
1915                   if there is no 'plot_numbers' then csv files with 'base_name' are
1916                   plotted and 2) in the title of the plots, the lenght of base_name is
1917                   removed from the front of the filename to be used in the title.
1918                   This could be changed if needed.
1919                   Note is ignored if assess_all_csv_files=True
1920       
1921        plot_numbers: a String list of numbers to plot. For example
1922                      [0-4,10,15-17] will read and attempt to plot
1923                       the follow 0,1,2,3,4,10,15,16,17
1924                       NOTE: if no plot numbers this will create
1925                       one plot per quantity, per gauge
1926        quantities: Will get available quantities from the header in the csv file.
1927                    quantities must be one of these.
1928                    NOTE: ALL QUANTITY NAMES MUST BE lower case!
1929                   
1930        extra_plot_name: A string that is appended to the end of the
1931                         output filename.
1932                   
1933        assess_all_csv_files: if true it will read ALL csv file with
1934                             "base_name", regardless of 'plot_numbers'
1935                              and determine a uniform set of axes for
1936                              Stage, Speed and Momentum. IF FALSE it
1937                              will only read the csv file within the
1938                             'plot_numbers'
1939                             
1940        create_latex: NOT IMPLEMENTED YET!! sorry Jane....
1941       
1942    OUTPUTS: saves the plots to
1943              <output_dir><base_name><plot_number><extra_plot_name>.png
1944             
1945     
1946    """
1947    try: 
1948        import pylab
1949    except ImportError:
1950        msg='csv2timeseries_graphs needs pylab to be installed correctly'
1951        raise msg
1952            #ANUGA don't need pylab to work so the system doesn't
1953            #rely on pylab being installed
1954        return
1955
1956    from os import sep
1957    from anuga.shallow_water.data_manager import \
1958                               get_all_files_with_extension, csv2dict
1959
1960    seconds_in_hour = 3600
1961    seconds_in_minutes = 60
1962   
1963    quantities_label={}
1964#    quantities_label['time'] = 'time (hours)'
1965    quantities_label['time'] = 'time (minutes)'
1966    quantities_label['stage'] = 'wave height (m)'
1967    quantities_label['speed'] = 'speed (m/s)'
1968    quantities_label['momentum'] = 'momentum (m^2/sec)'
1969    quantities_label['depth'] = 'water depth (m)'
1970    quantities_label['xmomentum'] = 'momentum (m^2/sec)'
1971    quantities_label['ymomentum'] = 'momentum (m^2/sec)'
1972    quantities_label['bearing'] = 'degrees (o)'
1973    quantities_label['elevation'] = 'elevation (m)'
1974
1975   
1976    if extra_plot_name != '':
1977        extra_plot_name='_'+extra_plot_name
1978
1979    new_plot_numbers=[]
1980    #change plot_numbers to list, eg ['0-4','10']
1981    #to ['0','1','2','3','4','10']
1982    for i, num_string in enumerate(plot_numbers):
1983        if '-' in num_string: 
1984            start = int(num_string[:num_string.rfind('-')])
1985            end = int(num_string[num_string.rfind('-')+1:])+1
1986            for x in range(start, end):
1987                new_plot_numbers.append(str(x))
1988        else:
1989            new_plot_numbers.append(num_string)
1990
1991    #finds all the files that fit the specs provided and return a list of them
1992    #so to help find a uniform max and min for the plots...
1993    list_filenames=[]
1994    all_csv_filenames=[]
1995    if verbose: print 'Determining files to access for axes ranges \n'
1996   
1997    #print directories_dic.keys(), base_name
1998 
1999    for i,directory in enumerate(directories_dic.keys()):
2000        all_csv_filenames.append(get_all_files_with_extension(directory,
2001                                  base_name,'.csv'))
2002
2003        filenames=[]
2004        if plot_numbers == '': 
2005            list_filenames.append(get_all_files_with_extension(directory,
2006                                  base_name,'.csv'))
2007        else:
2008            for number in new_plot_numbers:
2009#                print 'number!!!', base_name, number
2010                filenames.append(base_name+number)
2011#                print filenames
2012            list_filenames.append(filenames)
2013
2014
2015
2016    #print "list_filenames", list_filenames
2017
2018    #use all the files to get the values for the plot axis
2019    max_start_time= -1000.
2020    min_start_time = 100000 
2021   
2022   
2023    if verbose: print 'Determining uniform axes \n' 
2024    #this entire loop is to determine the min and max range for the
2025    #axes of the plots
2026
2027#    quantities.insert(0,'elevation')
2028    quantities.insert(0,'time')
2029
2030    directory_quantity_value={}
2031#    quantity_value={}
2032    min_quantity_value={}
2033    max_quantity_value={}
2034
2035   
2036    for i, directory in enumerate(directories_dic.keys()):
2037        filename_quantity_value={}
2038        if assess_all_csv_files==False:
2039            which_csv_to_assess = list_filenames[i]
2040        else:
2041            #gets list of filenames for directory "i"
2042            which_csv_to_assess = all_csv_filenames[i]
2043#        print'IN DIR', list_filenames[i]
2044       
2045
2046
2047       
2048        for j, filename in enumerate(which_csv_to_assess):
2049            quantity_value={}
2050
2051            dir_filename=join(directory,filename)
2052            attribute_dic, title_index_dic = csv2dict(dir_filename+
2053                                                       '.csv')
2054            directory_start_time = directories_dic[directory][1]
2055            directory_add_tide = directories_dic[directory][2]
2056
2057            if verbose: print 'reading: %s.csv' %dir_filename
2058#            print 'keys',attribute_dic.keys()
2059            #add time to get values
2060            for k, quantity in enumerate(quantities):
2061                quantity_value[quantity] = [float(x) for x in attribute_dic[quantity]]
2062
2063                #add tide to stage if provided
2064                if quantity == 'stage':
2065                     quantity_value[quantity]=array(quantity_value[quantity])+directory_add_tide
2066
2067                #condition to find max and mins for all the plots
2068                # populate the list with something when i=0 and j=0 and
2069                # then compare to the other values to determine abs max and min
2070                if i==0 and j==0:
2071
2072                    min_quantity_value[quantity], \
2073                    max_quantity_value[quantity] = get_min_max_values(quantity_value[quantity])
2074
2075                    if quantity != 'time':
2076                        min_quantity_value[quantity] = min_quantity_value[quantity] *1.1
2077                        max_quantity_value[quantity] = max_quantity_value[quantity] *1.1
2078
2079                   
2080#                    print '1 min,max',i,j,k,quantity, min_quantity_value[quantity],max_quantity_value[quantity],directory, filename
2081                else:
2082#                    print 'min,max',i,j,k,quantity, min_quantity_value[quantity],max_quantity_value[quantity],directory, filename
2083                    min, max = get_min_max_values(quantity_value[quantity])
2084#                    print "MIN",min, max
2085               
2086                    #min and max are multipled by "1+increase_axis" to get axes that are slighty bigger
2087                    # than the max and mins so the plots look good.
2088
2089                    increase_axis = (max-min)*0.05
2090#                    print quantity, "MIN MAX", max, min
2091                    if min<=min_quantity_value[quantity]:
2092                        if quantity == 'time': 
2093                            min_quantity_value[quantity]=min
2094                        else:
2095                            if round(min,2) == 0.00:
2096                                min_quantity_value[quantity]=-increase_axis
2097#                                min_quantity_value[quantity]=-2.
2098                                #min_quantity_value[quantity]= -max_quantity_value[quantity]*increase_axis
2099                            else:
2100#                                min_quantity_value[quantity]=min*(1+increase_axis)
2101                                min_quantity_value[quantity]=min-increase_axis
2102#                        print quantity, min_quantity_value[quantity]
2103                   
2104                    if max>max_quantity_value[quantity]: 
2105                        if quantity == 'time': 
2106                            max_quantity_value[quantity]=max
2107                        else:
2108                            max_quantity_value[quantity]=max+increase_axis
2109#                            max_quantity_value[quantity]=max*(1+increase_axis)
2110#                        print quantity, max_quantity_value[quantity],increase_axis
2111               
2112#                print 'min,maj',quantity, min_quantity_value[quantity],max_quantity_value[quantity]
2113
2114           
2115           
2116           
2117            #set the time... ???
2118            if min_start_time > directory_start_time: 
2119                min_start_time = directory_start_time
2120            if max_start_time < directory_start_time: 
2121                max_start_time = directory_start_time
2122            #print 'start_time' , max_start_time, min_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)/float(seconds_in_minutes),
2133                                        (float(max_quantity_value['time'])+float(max_start_time))\
2134                                              /float(seconds_in_minutes),
2135                                         min_quantity_value[quantity], 
2136                                         max_quantity_value[quantity])
2137        if verbose and (quantity != 'time' and quantity != 'elevation'): 
2138            print 'axis for quantity %s are x:(%s to %s)%s and y:(%s to %s)%s' %(quantity, 
2139                               quantities_axis[quantity][0],
2140                               quantities_axis[quantity][1],
2141                               quantities_label['time'],
2142                               quantities_axis[quantity][2],
2143                               quantities_axis[quantity][3],
2144                               quantities_label[quantity])
2145   
2146        #print  quantities_axis[quantity]
2147
2148    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
2149
2150    if verbose: print 'Now start to plot \n'
2151   
2152    i_max = len(directories_dic.keys())
2153    legend_list_dic={}
2154    legend_list =[]
2155    for i, directory in enumerate(directories_dic.keys()):
2156       
2157        if verbose: print'Plotting in %s %s' %(directory, new_plot_numbers)
2158#        print 'LIST',list_filenames
2159        #FIXME THIS SORT IS VERY IMPORTANT, without it the assigned plot numbers may not work correctly
2160        #there must be a better way
2161        list_filenames[i].sort()
2162        for j, filename in enumerate(list_filenames[i]):
2163#            print'IN plot', list_filenames[i]
2164           
2165            if verbose: print'Starting %s' %filename 
2166            directory_name = directories_dic[directory][0]
2167            directory_start_time = directories_dic[directory][1]
2168            directory_add_tide = directories_dic[directory][2]
2169           
2170            #create an if about the start time and tide hieght
2171            #if they don't exist
2172            #print 'i %s,j %s, number %s, file %s' %(i,j,number,file)
2173            attribute_dic, title_index_dic = csv2dict(directory+sep+filename+'.csv')
2174            #get data from dict in to list
2175            #do maths to list by changing to array
2176            t=(array(directory_quantity_value[directory][filename]['time'])+directory_start_time)/seconds_in_minutes
2177
2178            #finds the maximum elevation, used only as a test
2179            # and as info in the graphs
2180            max_ele=-100000
2181            min_ele=100000
2182            elevation = [float(x) for x in attribute_dic["elevation"]]
2183           
2184            min_ele, max_ele = get_min_max_values(elevation)
2185           
2186            if min_ele != max_ele:
2187                print "Note! Elevation changes in %s" %dir_filename
2188
2189            # creates a dictionary with keys that is the filename and attributes are a list of
2190            # lists containing 'directory_name' and 'elevation'. This is used to make the contents
2191            # for the legends in the graphs, this is the name of the model and the elevation.
2192            # All in this great one liner from DG. If the key 'filename' doesn't exist it creates the
2193            # entry if the entry exist it appends to the key.
2194           
2195            legend_list_dic.setdefault(filename,[]).append([directory_name,round(max_ele,3)])
2196
2197            # creates a LIST for the legend on the last iteration of the directories
2198            # which is when "legend_list_dic" has been fully populated. Creates a list of strings
2199            # which is used in the legend
2200            # only runs on the last iteration for all the gauges(csv) files
2201            # empties the list before creating it
2202            if i==i_max-1:
2203                legend_list=[]
2204   
2205                #print 'DIC',legend_list_dic
2206                for name_and_elevation in legend_list_dic[filename]:
2207                    legend_list.append('%s (elevation = %sm)'%(name_and_elevation[0],name_and_elevation[1]))
2208           
2209            #print 'filename',filename, quantities
2210            #skip time and elevation so it is not plotted!
2211            for k, quantity in enumerate(quantities):
2212                if quantity != 'time' and quantity != 'elevation':
2213                   
2214                    num=int(k*100+j)
2215                    pylab.figure(num)
2216                    pylab.ylabel(quantities_label[quantity])
2217                    pylab.plot(t, directory_quantity_value[directory][filename][quantity], c = cstr[i], linewidth=1)
2218                    pylab.xlabel(quantities_label['time'])
2219                    pylab.axis(quantities_axis[quantity])
2220                    pylab.legend(legend_list,loc='upper right')
2221                   
2222                    pylab.title('%s at %s gauge' %(quantity,filename[len(base_name):]))
2223                    if output_dir == '':
2224                        figname = '%s%s%s_%s%s.png' %(directory,sep,
2225                                            filename,
2226                                            quantity,
2227                                            extra_plot_name)
2228                    else:
2229                        figname = '%s%s%s_%s%s.png' %(output_dir,sep,
2230                                            filename,
2231                                            quantity,
2232                                            extra_plot_name)
2233                    if verbose: print 'saving figure here %s' %figname
2234                    pylab.savefig(figname)
2235           
2236    if verbose: print 'Closing all plots'
2237    pylab.close('all')
2238    del pylab
2239    if verbose: print 'Finished closing plots'
2240
2241def get_min_max_values(list=None):
2242    """
2243    Returns the min and max of the list it was provided.
2244    """
2245    if list == None: print 'List must be provided'
2246       
2247    return min(list), max(list)
2248
2249
2250
2251def get_runup_data_for_locations_from_file(gauge_filename,
2252                                           sww_filename,
2253                                           runup_filename,
2254                                           size=10,
2255                                           verbose=False):
2256
2257    """this will read a csv file with the header x,y. Then look in a square
2258    'size'x2 around this position for the 'max_inundaiton_height' in the
2259    'sww_filename' and report the findings in the 'runup_filename
2260   
2261    WARNING: NO TESTS!
2262    """
2263
2264    from anuga.shallow_water.data_manager import get_all_directories_with_name,\
2265                                                 get_maximum_inundation_data,\
2266                                                 csv2dict
2267                                                 
2268    file = open(runup_filename,"w")
2269    file.write("easting,northing,runup \n ")
2270    file.close()
2271   
2272    #read gauge csv file to dictionary
2273    attribute_dic, title_index_dic = csv2dict(gauge_filename)
2274    northing = [float(x) for x in attribute_dic["y"]]
2275    easting = [float(x) for x in attribute_dic["x"]]
2276
2277    print 'Reading %s' %sww_filename
2278
2279    runup_locations=[]
2280    for i, x in enumerate(northing):
2281#        print 'easting,northing',i,easting[i],northing[i]
2282        poly = [[int(easting[i]+size),int(northing[i]+size)],
2283                [int(easting[i]+size),int(northing[i]-size)],
2284                [int(easting[i]-size),int(northing[i]-size)],
2285                [int(easting[i]-size),int(northing[i]+size)]]
2286       
2287        run_up, x_y = get_maximum_inundation_data(filename=sww_filename,
2288                                              polygon=poly,
2289                                              verbose=False) 
2290        #if no runup will return 0 instead of NONE
2291        if run_up==None: run_up=0
2292        if x_y==None: x_y=[0,0]
2293       
2294        if verbose:
2295            print 'maximum inundation runup near %s is %s meters' %(x_y,run_up)
2296       
2297        #writes to file
2298        file = open(runup_filename,"a")
2299        temp = '%s,%s,%s \n' %(x_y[0], x_y[1], run_up)
2300        file.write(temp)
2301        file.close()
2302       
2303
2304def sww2csv_gauges(sww_file,
2305                   gauge_file,
2306                   out_name='gauge_',
2307                   quantities = ['stage', 'depth', 'elevation',
2308                                 'xmomentum', 'ymomentum'],
2309                   verbose=False,
2310                   use_cache = True):
2311    """
2312   
2313    Inputs:
2314       
2315        NOTE: if using csv2timeseries_graphs after creating csv file,
2316        it is essential to export quantities 'depth' and 'elevation'.
2317        'depth' is good to analyse gauges on land and elevation is used
2318        automatically by csv2timeseries_graphs in the legend.
2319       
2320        sww_file: path to any sww file
2321       
2322        gauge_file: Assumes that it follows this format
2323            name, easting, northing, elevation
2324            point1, 100.3, 50.2, 10.0
2325            point2, 10.3, 70.3, 78.0
2326       
2327        NOTE: order of column can change but names eg 'easting', 'elevation'
2328        must be the same! ALL lowercaps!
2329
2330        out_name: prefix for output file name (default is 'gauge_')
2331       
2332    Outputs:
2333        one file for each gauge/point location in the points file. They
2334        will be named with this format in the same directory as the 'sww_file'
2335            <out_name><name>.csv
2336        eg gauge_point1.csv if <out_name> not supplied
2337           myfile_2_point1.csv if <out_name> ='myfile_2_'
2338           
2339           
2340        They will all have a header
2341   
2342    Usage: sww2csv_gauges(sww_file='test1.sww',
2343                          quantities = ['stage', 'elevation','depth','bearing'],
2344                          gauge_file='gauge.txt')   
2345   
2346    Interpolate the quantities at a given set of locations, given
2347    an sww file.
2348    The results are written to a csv file.
2349
2350    In the future let points be a points file.
2351    And the user choose the quantities.
2352
2353    This is currently quite specific.
2354    If it needs to be more general, change things.
2355
2356    This is really returning speed, not velocity.
2357     
2358   
2359    """
2360   
2361    from csv import reader,writer
2362    from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
2363    from Numeric import array, resize, shape, Float, zeros, take, argsort, argmin
2364    import string
2365    from anuga.shallow_water.data_manager import get_all_swwfiles
2366#    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
2367    #print "points",points
2368
2369    assert type(gauge_file) == type(''),\
2370           'Gauge filename must be a string'
2371           
2372    assert type(out_name) == type(''),\
2373           'Output filename prefix must be a string'
2374   
2375    try:
2376    #    fid = open(gauge_file)
2377        point_reader = reader(file(gauge_file))
2378
2379    except Exception, e:
2380        msg = 'File "%s" could not be opened: Error="%s"'\
2381                  %(gauge_file, e)
2382        raise msg
2383    if verbose: print '\n Gauges obtained from: %s \n' %gauge_file
2384   
2385   
2386    point_reader = reader(file(gauge_file))
2387    points = []
2388    point_name = []
2389   
2390    #read point info from file
2391    for i,row in enumerate(point_reader):
2392#        print 'i',i,'row',row
2393        #read header and determine the column numbers to read correcty.
2394        if i==0:
2395            for j,value in enumerate(row):
2396#                print 'j',j,value, row
2397                if value.strip()=='easting':easting=j
2398                if value.strip()=='northing':northing=j
2399                if value.strip()=='name':name=j
2400                if value.strip()=='elevation':elevation=j
2401        else:
2402#            print i,'easting',easting,'northing',northing, row[easting]
2403            points.append([float(row[easting]),float(row[northing])])
2404            point_name.append(row[name])
2405       
2406    #convert to array for file_function
2407    points_array = array(points,Float)
2408       
2409    points_array = ensure_absolute(points_array)
2410
2411    dir_name, base = os.path.split(sww_file)   
2412    #print 'dirname',dir_name, base
2413    #need to get current directory so when path and file
2414    #are "joined" below the directory is correct
2415    if dir_name == '':
2416        dir_name =getcwd()
2417       
2418    if access(sww_file,R_OK):
2419        if verbose: print 'File %s exists' %(sww_file)
2420    else:
2421        msg = 'File "%s" could not be opened: no read permission'\
2422               %(sww_file)
2423        raise msg
2424
2425    sww_files = get_all_swwfiles(look_in_dir=dir_name,
2426                                 base_name=base,
2427                                 verbose=verbose)
2428   
2429    #to make all the quantities lower case for file_function
2430    quantities = [quantity.lower() for quantity in quantities]
2431
2432    # what is quantities are needed from sww file to calculate output quantities
2433    # also
2434
2435    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
2436   
2437    for sww_file in sww_files:
2438   
2439#        print 'sww_file',sww_file
2440        sww_file = join(dir_name, sww_file+'.sww')
2441#        print 'sww_file',sww_file, core_quantities
2442       
2443        callable_sww = file_function(sww_file,
2444                                 quantities=core_quantities,
2445                                 interpolation_points=points_array,
2446                                 verbose=verbose,
2447                                 use_cache=use_cache)
2448       
2449    gauge_file = out_name
2450
2451    heading = [quantity for quantity in quantities]
2452    heading.insert(0,'time')
2453
2454#    print 'start time', callable_sww.starttime, heading, quantities
2455
2456    #create a list of csv writers for all the points and write header
2457    points_writer = []
2458    for i,point in enumerate(points):
2459        #print 'gauge file:',dir_name+sep+'gauge_'+point_name[i]+'.csv'
2460        points_writer.append(writer(file(dir_name+sep+gauge_file+point_name[i]+'.csv', "wb")))
2461        points_writer[i].writerow(heading)
2462   
2463    if verbose: print 'Writing csv files'
2464
2465    for time in callable_sww.get_time():
2466
2467        for point_i, point in enumerate(points_array):
2468            #add domain starttime to relative time.
2469            points_list = [time+callable_sww.starttime]
2470#            print'time',time,'point_i',point_i,point, points_array
2471            point_quantities = callable_sww(time,point_i)
2472#            print "quantities", point_quantities
2473           
2474            for quantity in quantities:
2475                if quantity==NAN:
2476                    print 'quantity does not exist in' %callable_sww.get_name
2477                else:
2478                    if quantity == 'stage':
2479                        points_list.append(point_quantities[0])
2480                       
2481                    if quantity == 'elevation':
2482                        points_list.append(point_quantities[1])
2483                       
2484                    if quantity == 'xmomentum':
2485                        points_list.append(point_quantities[2])
2486                       
2487                    if quantity == 'ymomentum':
2488                        points_list.append(point_quantities[3])
2489                       
2490                    if quantity == 'depth':
2491                        points_list.append(point_quantities[0] - point_quantities[1])
2492                       
2493                       
2494                    if quantity == 'momentum':
2495                        momentum = sqrt(point_quantities[2]**2 +\
2496                                        point_quantities[3]**2)
2497                        points_list.append(momentum)
2498                       
2499                    if quantity == 'speed':
2500                        #if depth is less than 0.001 then speed = 0.0
2501                        if point_quantities[0] - point_quantities[1] < 0.001:
2502                            vel = 0.0
2503                        else:
2504                            if point_quantities[2] < 1.0e6:
2505                                momentum = sqrt(point_quantities[2]**2 +\
2506                                                point_quantities[3]**2)
2507    #                            vel = momentum/depth             
2508                                vel = momentum/(point_quantities[0] - point_quantities[1])
2509    #                            vel = momentum/(depth + 1.e-6/depth)
2510                            else:
2511                                momentum = 0
2512                                vel = 0
2513                           
2514                        points_list.append(vel)
2515                       
2516                    if quantity == 'bearing':
2517                        points_list.append(calc_bearing(point_quantities[2],
2518                                                        point_quantities[3]))
2519                       
2520                #print 'list',points_list
2521               
2522            points_writer[point_i].writerow(points_list)
2523       
2524
2525def greens_law(d1,d2,h1,verbose=False):
2526    """
2527
2528    Green's Law allows an approximation of wave amplitude at
2529    a given depth based on the fourh root of the ratio of two depths
2530    and the amplitude at another given depth.
2531
2532    Note, wave amplitude is equal to stage.
2533   
2534    Inputs:
2535
2536    d1, d2 - the two depths
2537    h1     - the wave amplitude at d1
2538    h2     - the derived amplitude at d2
2539
2540    h2 = h1 (d1/d2)^(1/4), where d2 cannot equal 0.
2541   
2542    """
2543
2544    d1 = ensure_numeric(d1)
2545    d2 = ensure_numeric(d2)
2546    h1 = ensure_numeric(h1)
2547
2548    if d1 <= 0.0:
2549        msg = 'the first depth, d1 (%f), must be strictly positive' %(d1)
2550        raise Exception(msg)
2551
2552    if d2 <= 0.0:
2553        msg = 'the second depth, d2 (%f), must be strictly positive' %(d2)
2554        raise Exception(msg)
2555   
2556    if h1 <= 0.0:
2557        msg = 'the wave amplitude, h1 (%f), must be strictly positive' %(h1)
2558        raise Exception(msg)
2559   
2560    h2 = h1*(d1/d2)**0.25
2561
2562    assert h2 > 0
2563   
2564    return h2
2565       
2566
2567def square_root(s):
2568    return sqrt(s)
Note: See TracBrowser for help on using the repository browser.