source: branches/numpy/anuga/abstract_2d_finite_volumes/util.py @ 6883

Last change on this file since 6883 was 6689, checked in by rwilson, 16 years ago

Back-merge from Numeric trunk to numpy branch.

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