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

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

Initial commit of numpy changes. Still a long way to go.

File size: 100.4 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, 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    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
928   
929    try:
930        fid = open(gauge_filename)
931    except Exception, e:
932        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
933        raise msg
934
935    if report is None:
936        report = False
937       
938    if plot_quantity is None:
939        plot_quantity = ['depth', 'speed']
940    else:
941        assert type(plot_quantity) == list, 'plot_quantity must be a list'
942        check_list(plot_quantity)
943
944    if surface is None:
945        surface = False
946
947    if time_unit is None:
948        time_unit = 'hours'
949   
950    if title_on is None:
951        title_on = True
952   
953    if verbose: print '\n Gauges obtained from: %s \n' %gauge_filename
954
955    gauges, locations, elev = get_gauges_from_file(gauge_filename)
956
957    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
958
959    file_loc = []
960    f_list = []
961    label_id = []
962    leg_label = []
963    themaxT = 0.0
964    theminT = 0.0
965
966    for swwfile in swwfiles.keys():
967        try:
968            fid = open(swwfile)
969        except Exception, e:
970            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
971            raise msg
972
973        print 'swwfile', swwfile
974
975        # Extract parent dir name and use as label
976        path, _ = os.path.split(swwfile)
977        _, label = os.path.split(path)       
978       
979        #print 'label', label
980        leg_label.append(label)
981
982        f = file_function(swwfile,
983                          quantities = sww_quantity,
984                          interpolation_points = gauges,
985                          time_thinning = time_thinning,
986                          verbose = verbose,
987                          use_cache = use_cache)
988
989        # determine which gauges are contained in sww file
990        count = 0
991        gauge_index = []
992        for k, g in enumerate(gauges):
993            if f(0.0, point_id = k)[2] > 1.0e6:
994                count += 1
995                if count == 1: print 'Gauges not contained here:'
996                print locations[k]
997            else:
998                gauge_index.append(k)
999
1000        if len(gauge_index) > 0:
1001            print 'Gauges contained here: \n',
1002        else:
1003            print 'No gauges contained here. \n'
1004        for i in range(len(gauge_index)):
1005             print locations[gauge_index[i]]
1006             
1007        index = swwfile.rfind(sep)
1008        file_loc.append(swwfile[:index+1])
1009        label_id.append(swwfiles[swwfile])
1010       
1011        f_list.append(f)
1012        maxT = max(f.get_time())
1013        minT = min(f.get_time())
1014        if maxT > themaxT: themaxT = maxT
1015        if minT > theminT: theminT = minT
1016
1017    if time_min is None:
1018        time_min = theminT # min(T)
1019    else:
1020        if time_min < theminT: # min(T):
1021            msg = 'Minimum time entered not correct - please try again'
1022            raise Exception, msg
1023
1024    if time_max is None:
1025        time_max = themaxT # max(T)
1026    else:
1027        if time_max > themaxT: # max(T):
1028            msg = 'Maximum time entered not correct - please try again'
1029            raise Exception, msg
1030
1031    if verbose and len(gauge_index) > 0:
1032         print 'Inputs OK - going to generate figures'
1033
1034    if len(gauge_index) <> 0:
1035        texfile, elev_output = \
1036            generate_figures(plot_quantity, file_loc, report, reportname,
1037                             surface, leg_label, f_list, gauges, locations,
1038                             elev, gauge_index, production_dirs, time_min,
1039                             time_max, time_unit, title_on, label_id,
1040                             generate_fig, verbose)
1041    else:
1042        texfile = ''
1043        elev_output = []
1044
1045    return texfile, elev_output
1046
1047
1048##
1049# @brief Read gauge info from a file.
1050# @param filename The name of the file to read.
1051# @return A (gauges, gaugelocation, elev) tuple.
1052def get_gauges_from_file(filename):
1053    """ Read in gauge information from file
1054    """
1055
1056    from os import sep, getcwd, access, F_OK, mkdir
1057
1058    # Get data from the gauge file
1059    fid = open(filename)
1060    lines = fid.readlines()
1061    fid.close()
1062   
1063    gauges = []
1064    gaugelocation = []
1065    elev = []
1066
1067    # Check header information   
1068    line1 = lines[0]
1069    line11 = line1.split(',')
1070
1071    if isinstance(line11[0], str) is True:
1072        # We have found text in the first line
1073        east_index = None
1074        north_index = None
1075        name_index = None
1076        elev_index = None
1077
1078        for i in range(len(line11)):
1079            if line11[i].strip().lower() == 'easting':   east_index = i
1080            if line11[i].strip().lower() == 'northing':  north_index = i
1081            if line11[i].strip().lower() == 'name':      name_index = i
1082            if line11[i].strip().lower() == 'elevation': elev_index = i
1083
1084        if east_index < len(line11) and north_index < len(line11):
1085            pass
1086        else:
1087            msg = 'WARNING: %s does not contain correct header information' \
1088                  % filename
1089            msg += 'The header must be: easting, northing, name, elevation'
1090            raise Exception, msg
1091
1092        if elev_index is None: 
1093            raise Exception
1094   
1095        if name_index is None: 
1096            raise Exception
1097
1098        lines = lines[1:] # Remove header from data
1099    else:
1100        # No header, assume that this is a simple easting, northing file
1101
1102        msg = 'There was no header in file %s and the number of columns is %d' \
1103              % (filename, len(line11))
1104        msg += '- was assuming two columns corresponding to Easting and Northing'
1105        assert len(line11) == 2, msg
1106
1107        east_index = 0
1108        north_index = 1
1109
1110        N = len(lines)
1111        elev = [-9999]*N
1112        gaugelocation = range(N)
1113       
1114    # Read in gauge data
1115    for line in lines:
1116        fields = line.split(',')
1117
1118        gauges.append([float(fields[east_index]), float(fields[north_index])])
1119
1120        if len(fields) > 2:
1121            elev.append(float(fields[elev_index]))
1122            loc = fields[name_index]
1123            gaugelocation.append(loc.strip('\n'))
1124
1125    return gauges, gaugelocation, elev
1126
1127
1128##
1129# @brief Check that input quantities in quantity list are legal.
1130# @param quantity Quantity list to check.
1131# @note Raises an exception of list is not legal.
1132def check_list(quantity):
1133    """ Check that input quantities in quantity list are possible
1134    """
1135    import sys
1136    from sets import Set as set
1137
1138    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
1139                    'ymomentum', 'speed', 'bearing', 'elevation']
1140
1141    # convert all quanitiy names to lowercase
1142    for i,j in enumerate(quantity):
1143        quantity[i] = quantity[i].lower()
1144
1145    # check that all names in 'quantity' appear in 'all_quantity'
1146    p = list(set(quantity).difference(set(all_quantity)))
1147    if len(p) != 0:
1148        msg = 'Quantities %s do not exist - please try again' %p
1149        raise Exception, msg
1150
1151
1152##
1153# @brief Calculate velocity bearing from North.
1154# @param uh ??
1155# @param vh ??
1156# @return The calculated bearing.
1157def calc_bearing(uh, vh):
1158    """ Calculate velocity bearing from North
1159    """
1160    #FIXME (Ole): I reckon we should refactor this one to use
1161    #             the function angle() in utilities/numerical_tools
1162    #
1163    #             It will be a simple matter of
1164    # * converting from radians to degrees
1165    # * moving the reference direction from [1,0] to North
1166    # * changing from counter clockwise to clocwise.
1167       
1168    angle = degrees(atan(vh/(uh+1.e-15)))
1169
1170    if (0 < angle < 90.0):
1171        if vh > 0:
1172            bearing = 90.0 - abs(angle)
1173        if vh < 0:
1174            bearing = 270.0 - abs(angle)
1175   
1176    if (-90 < angle < 0):
1177        if vh < 0:
1178            bearing = 90.0 - (angle)
1179        if vh > 0:
1180            bearing = 270.0 - (angle)
1181    if angle == 0: bearing = 0.0
1182
1183    return bearing
1184
1185
1186##
1187# @brief Generate figures from quantities and gauges for each sww file.
1188# @param plot_quantity  ??
1189# @param file_loc ??
1190# @param report ??
1191# @param reportname ??
1192# @param surface ??
1193# @param leg_label ??
1194# @param f_list ??
1195# @param gauges ??
1196# @param locations ??
1197# @param elev ??
1198# @param gauge_index ??
1199# @param production_dirs ??
1200# @param time_min ??
1201# @param time_max ??
1202# @param time_unit ??
1203# @param title_on ??
1204# @param label_id ??
1205# @param generate_fig ??
1206# @param verbose??
1207# @return (texfile2, elev_output)
1208def generate_figures(plot_quantity, file_loc, report, reportname, surface,
1209                     leg_label, f_list, gauges, locations, elev, gauge_index,
1210                     production_dirs, time_min, time_max, time_unit,
1211                     title_on, label_id, generate_fig, verbose):
1212    """ Generate figures based on required quantities and gauges for
1213    each sww file
1214    """
1215    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
1216
1217    if generate_fig is True:
1218        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
1219             xlabel, ylabel, title, close, subplot
1220   
1221        if surface is True:
1222            import pylab as p1
1223            import mpl3d.mplot3d as p3
1224       
1225    if report == True:   
1226        texdir = getcwd()+sep+'report'+sep
1227        if access(texdir,F_OK) == 0:
1228            mkdir (texdir)
1229        if len(label_id) == 1:
1230            label_id1 = label_id[0].replace(sep,'')
1231            label_id2 = label_id1.replace('_','')
1232            texfile = texdir + reportname + '%s' % label_id2
1233            texfile2 = reportname + '%s' % label_id2
1234            texfilename = texfile + '.tex'
1235            fid = open(texfilename, 'w')
1236
1237            if verbose: print '\n Latex output printed to %s \n' %texfilename
1238        else:
1239            texfile = texdir+reportname
1240            texfile2 = reportname
1241            texfilename = texfile + '.tex' 
1242            fid = open(texfilename, 'w')
1243
1244            if verbose: print '\n Latex output printed to %s \n' %texfilename
1245    else:
1246        texfile = ''
1247        texfile2 = ''
1248
1249    p = len(f_list)
1250    n = []
1251    n0 = 0
1252    for i in range(len(f_list)):
1253        n.append(len(f_list[i].get_time()))
1254        if n[i] > n0: n0 = n[i] 
1255    n0 = int(n0)
1256    m = len(locations)
1257    model_time = num.zeros((n0, m, p), num.float) 
1258    stages = num.zeros((n0, m, p), num.float)
1259    elevations = num.zeros((n0, m, p), num.float) 
1260    momenta = num.zeros((n0, m, p), num.float)
1261    xmom = num.zeros((n0, m, p), num.float)
1262    ymom = num.zeros((n0, m, p), num.float)
1263    speed = num.zeros((n0, m, p), num.float)
1264    bearings = num.zeros((n0, m, p), num.float)
1265    due_east = 90.0*num.ones((n0, 1), num.float)
1266    due_west = 270.0*num.ones((n0, 1), num.float)
1267    depths = num.zeros((n0, m, p), num.float)
1268    eastings = num.zeros((n0, m, p), num.float)
1269    min_stages = []
1270    max_stages = []
1271    min_momentums = []   
1272    max_momentums = []
1273    max_xmomentums = []
1274    max_ymomentums = []
1275    min_xmomentums = []
1276    min_ymomentums = []
1277    max_speeds = []
1278    min_speeds = []   
1279    max_depths = []
1280    model_time_plot3d = num.zeros((n0, m), num.float)
1281    stages_plot3d = num.zeros((n0, m), num.float)
1282    eastings_plot3d = num.zeros((n0, m),num.float)
1283    if time_unit is 'mins': scale = 60.0
1284    if time_unit is 'hours': scale = 3600.0
1285
1286    ##### loop over each swwfile #####
1287    for j, f in enumerate(f_list):
1288        if verbose: print 'swwfile %d of %d' % (j, len(f_list))
1289
1290        starttime = f.starttime
1291        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
1292        fid_compare = open(comparefile, 'w')
1293        file0 = file_loc[j] + 'gauges_t0.csv'
1294        fid_0 = open(file0, 'w')
1295
1296        ##### loop over each gauge #####
1297        for k in gauge_index:
1298            if verbose: print 'Gauge %d of %d' % (k, len(gauges))
1299
1300            g = gauges[k]
1301            min_stage = 10
1302            max_stage = 0
1303            max_momentum = max_xmomentum = max_ymomentum = 0
1304            min_momentum = min_xmomentum = min_ymomentum = 100
1305            max_speed = 0
1306            min_speed = 0           
1307            max_depth = 0           
1308            gaugeloc = str(locations[k])
1309            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
1310                       + gaugeloc + '.csv'
1311            fid_out = open(thisfile, 'w')
1312            s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
1313            fid_out.write(s)
1314
1315            #### generate quantities #######
1316            for i, t in enumerate(f.get_time()):
1317                if time_min <= t <= time_max:
1318                    w = f(t, point_id = k)[0]
1319                    z = f(t, point_id = k)[1]
1320                    uh = f(t, point_id = k)[2]
1321                    vh = f(t, point_id = k)[3]
1322                    depth = w-z     
1323                    m = sqrt(uh*uh + vh*vh)
1324                    if depth < 0.001:
1325                        vel = 0.0
1326                    else:
1327                        vel = m / (depth + 1.e-6/depth) 
1328                    bearing = calc_bearing(uh, vh)                   
1329                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
1330                    stages[i,k,j] = w
1331                    elevations[i,k,j] = z
1332                    xmom[i,k,j] = uh
1333                    ymom[i,k,j] = vh
1334                    momenta[i,k,j] = m
1335                    speed[i,k,j] = vel
1336                    bearings[i,k,j] = bearing
1337                    depths[i,k,j] = depth
1338                    thisgauge = gauges[k]
1339                    eastings[i,k,j] = thisgauge[0]
1340                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
1341                            % (t, w, m, vel, z, uh, vh, bearing)
1342                    fid_out.write(s)
1343                    if t == 0:
1344                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
1345                        fid_0.write(s)
1346                    if t/60.0 <= 13920: tindex = i
1347                    if w > max_stage: max_stage = w
1348                    if w < min_stage: min_stage = w
1349                    if m > max_momentum: max_momentum = m
1350                    if m < min_momentum: min_momentum = m                   
1351                    if uh > max_xmomentum: max_xmomentum = uh
1352                    if vh > max_ymomentum: max_ymomentum = vh
1353                    if uh < min_xmomentum: min_xmomentum = uh
1354                    if vh < min_ymomentum: min_ymomentum = vh
1355                    if vel > max_speed: max_speed = vel
1356                    if vel < min_speed: min_speed = vel                   
1357                    if z > 0 and depth > max_depth: max_depth = depth
1358                   
1359                   
1360            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
1361                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
1362            fid_compare.write(s)
1363            max_stages.append(max_stage)
1364            min_stages.append(min_stage)
1365            max_momentums.append(max_momentum)
1366            max_xmomentums.append(max_xmomentum)
1367            max_ymomentums.append(max_ymomentum)
1368            min_xmomentums.append(min_xmomentum)
1369            min_ymomentums.append(min_ymomentum)
1370            min_momentums.append(min_momentum)           
1371            max_depths.append(max_depth)
1372            max_speeds.append(max_speed)
1373            min_speeds.append(min_speed)           
1374            #### finished generating quantities for each swwfile #####
1375       
1376        model_time_plot3d[:,:] = model_time[:,:,j]
1377        stages_plot3d[:,:] = stages[:,:,j]
1378        eastings_plot3d[:,] = eastings[:,:,j]
1379           
1380        if surface is True:
1381            print 'Printing surface figure'
1382            for i in range(2):
1383                fig = p1.figure(10)
1384                ax = p3.Axes3D(fig)
1385                if len(gauges) > 80:
1386                    ax.plot_surface(model_time[:,:,j],
1387                                    eastings[:,:,j],
1388                                    stages[:,:,j])
1389                else:
1390                    ax.plot3D(num.ravel(eastings[:,:,j]),
1391                              num.ravel(model_time[:,:,j]),
1392                              num.ravel(stages[:,:,j]))
1393                ax.set_xlabel('time')
1394                ax.set_ylabel('x')
1395                ax.set_zlabel('stage')
1396                fig.add_axes(ax)
1397                p1.show()
1398                surfacefig = 'solution_surface%s' % leg_label[j]
1399                p1.savefig(surfacefig)
1400                p1.close()
1401           
1402    #### finished generating quantities for all swwfiles #####
1403
1404    # x profile for given time
1405    if surface is True:
1406        figure(11)
1407        plot(eastings[tindex,:,j], stages[tindex,:,j])
1408        xlabel('x')
1409        ylabel('stage')
1410        profilefig = 'solution_xprofile' 
1411        savefig('profilefig')
1412
1413    elev_output = []
1414    if generate_fig is True:
1415        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
1416                           max(max_depths)*1.1])
1417        stage_axis = axis([starttime/scale, time_max/scale,
1418                           min(min_stages), max(max_stages)*1.1])
1419        vel_axis = axis([starttime/scale, time_max/scale,
1420                         min(min_speeds), max(max_speeds)*1.1])
1421        mom_axis = axis([starttime/scale, time_max/scale,
1422                         min(min_momentums), max(max_momentums)*1.1])
1423        xmom_axis = axis([starttime/scale, time_max/scale,
1424                          min(min_xmomentums), max(max_xmomentums)*1.1])
1425        ymom_axis = axis([starttime/scale, time_max/scale,
1426                          min(min_ymomentums), max(max_ymomentums)*1.1])
1427        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
1428        nn = len(plot_quantity)
1429        no_cols = 2
1430       
1431        if len(label_id) > 1: graphname_report = []
1432        pp = 1
1433        div = 11.
1434        cc = 0
1435        for k in gauge_index:
1436            g = gauges[k]
1437            count1 = 0
1438            if report == True and len(label_id) > 1:
1439                s = '\\begin{figure}[ht] \n' \
1440                    '\\centering \n' \
1441                    '\\begin{tabular}{cc} \n'
1442                fid.write(s)
1443            if len(label_id) > 1: graphname_report = []
1444
1445            #### generate figures for each gauge ####
1446            for j, f in enumerate(f_list):
1447                ion()
1448                hold(True)
1449                count = 0
1450                where1 = 0
1451                where2 = 0
1452                word_quantity = ''
1453                if report == True and len(label_id) == 1:
1454                    s = '\\begin{figure}[hbt] \n' \
1455                        '\\centering \n' \
1456                        '\\begin{tabular}{cc} \n'
1457                    fid.write(s)
1458                   
1459                for which_quantity in plot_quantity:
1460                    count += 1
1461                    where1 += 1
1462                    figure(count, frameon = False)
1463                    if which_quantity == 'depth':
1464                        plot(model_time[0:n[j]-1,k,j],
1465                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
1466                        units = 'm'
1467                        axis(depth_axis)
1468                    if which_quantity == 'stage':
1469                        if elevations[0,k,j] <= 0:
1470                            plot(model_time[0:n[j]-1,k,j],
1471                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
1472                            axis(stage_axis)
1473                        else:
1474                            plot(model_time[0:n[j]-1,k,j],
1475                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
1476                            #axis(depth_axis)                 
1477                        units = 'm'
1478                    if which_quantity == 'momentum':
1479                        plot(model_time[0:n[j]-1,k,j],
1480                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
1481                        axis(mom_axis)
1482                        units = 'm^2 / sec'
1483                    if which_quantity == 'xmomentum':
1484                        plot(model_time[0:n[j]-1,k,j],
1485                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
1486                        axis(xmom_axis)
1487                        units = 'm^2 / sec'
1488                    if which_quantity == 'ymomentum':
1489                        plot(model_time[0:n[j]-1,k,j],
1490                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
1491                        axis(ymom_axis)
1492                        units = 'm^2 / sec'
1493                    if which_quantity == 'speed':
1494                        plot(model_time[0:n[j]-1,k,j],
1495                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
1496                        axis(vel_axis)
1497                        units = 'm / sec'
1498                    if which_quantity == 'bearing':
1499                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
1500                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
1501                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
1502                        units = 'degrees from North'
1503                        #ax = axis([time_min, time_max, 0.0, 360.0])
1504                        legend(('Bearing','West','East'))
1505
1506                    if time_unit is 'mins': xlabel('time (mins)')
1507                    if time_unit is 'hours': xlabel('time (hours)')
1508                    #if which_quantity == 'stage' \
1509                    #   and elevations[0:n[j]-1,k,j] > 0:
1510                    #    ylabel('%s (%s)' %('depth', units))
1511                    #else:
1512                    #    ylabel('%s (%s)' %(which_quantity, units))
1513                        #ylabel('%s (%s)' %('wave height', units))
1514                    ylabel('%s (%s)' %(which_quantity, units))
1515                    if len(label_id) > 1: legend((leg_label),loc='upper right')
1516
1517                    #gaugeloc1 = gaugeloc.replace(' ','')
1518                    #gaugeloc2 = gaugeloc1.replace('_','')
1519                    gaugeloc2 = str(locations[k]).replace(' ','')
1520                    graphname = '%sgauge%s_%s' %(file_loc[j],
1521                                                 gaugeloc2,
1522                                                 which_quantity)
1523
1524                    if report == True and len(label_id) > 1:
1525                        figdir = getcwd()+sep+'report_figures'+sep
1526                        if access(figdir,F_OK) == 0 :
1527                            mkdir (figdir)
1528                        latex_file_loc = figdir.replace(sep,altsep) 
1529                        # storing files in production directory   
1530                        graphname_latex = '%sgauge%s%s' \
1531                                          % (latex_file_loc, gaugeloc2,
1532                                             which_quantity)
1533                        # giving location in latex output file
1534                        graphname_report_input = '%sgauge%s%s' % \
1535                                                 ('..' + altsep + 
1536                                                      'report_figures' + altsep,
1537                                                  gaugeloc2, which_quantity)
1538                        graphname_report.append(graphname_report_input)
1539                       
1540                        # save figures in production directory for report
1541                        savefig(graphname_latex)
1542
1543                    if report == True:
1544                        figdir = getcwd() + sep + 'report_figures' + sep
1545                        if access(figdir,F_OK) == 0:
1546                            mkdir(figdir)
1547                        latex_file_loc = figdir.replace(sep,altsep)   
1548
1549                        if len(label_id) == 1: 
1550                            # storing files in production directory 
1551                            graphname_latex = '%sgauge%s%s%s' % \
1552                                              (latex_file_loc, gaugeloc2,
1553                                               which_quantity, label_id2)
1554                            # giving location in latex output file
1555                            graphname_report = '%sgauge%s%s%s' % \
1556                                               ('..' + altsep +
1557                                                    'report_figures' + altsep,
1558                                                gaugeloc2, which_quantity,
1559                                                label_id2)
1560                            s = '\includegraphics' \
1561                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1562                                (graphname_report, '.png')
1563                            fid.write(s)
1564                            if where1 % 2 == 0:
1565                                s = '\\\\ \n'
1566                                where1 = 0
1567                            else:
1568                                s = '& \n'
1569                            fid.write(s)
1570                            savefig(graphname_latex)
1571                   
1572                    if title_on == True:
1573                        title('%s scenario: %s at %s gauge' % \
1574                              (label_id, which_quantity, gaugeloc2))
1575                        #title('Gauge %s (MOST elevation %.2f, ' \
1576                        #      'ANUGA elevation %.2f)' % \
1577                        #      (gaugeloc2, elevations[10,k,0],
1578                        #       elevations[10,k,1]))
1579
1580                    savefig(graphname) # save figures with sww file
1581
1582                if report == True and len(label_id) == 1:
1583                    for i in range(nn-1):
1584                        if nn > 2:
1585                            if plot_quantity[i] == 'stage' \
1586                               and elevations[0,k,j] > 0:
1587                                word_quantity += 'depth' + ', '
1588                            else:
1589                                word_quantity += plot_quantity[i] + ', '
1590                        else:
1591                            if plot_quantity[i] == 'stage' \
1592                               and elevations[0,k,j] > 0:
1593                                word_quantity += 'depth' + ', '
1594                            else:
1595                                word_quantity += plot_quantity[i]
1596                       
1597                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
1598                        word_quantity += ' and ' + 'depth'
1599                    else:
1600                        word_quantity += ' and ' + plot_quantity[nn-1]
1601                    caption = 'Time series for %s at %s location ' \
1602                              '(elevation %.2fm)' % \
1603                              (word_quantity, locations[k], elev[k])
1604                    if elev[k] == 0.0:
1605                        caption = 'Time series for %s at %s location ' \
1606                                  '(elevation %.2fm)' % \
1607                                  (word_quantity, locations[k],
1608                                   elevations[0,k,j])
1609                        east = gauges[0]
1610                        north = gauges[1]
1611                        elev_output.append([locations[k], east, north,
1612                                            elevations[0,k,j]])
1613                    label = '%sgauge%s' % (label_id2, gaugeloc2)
1614                    s = '\end{tabular} \n' \
1615                        '\\caption{%s} \n' \
1616                        '\label{fig:%s} \n' \
1617                        '\end{figure} \n \n' % (caption, label)
1618                    fid.write(s)
1619                    cc += 1
1620                    if cc % 6 == 0: fid.write('\\clearpage \n')
1621                    savefig(graphname_latex)               
1622                   
1623            if report == True and len(label_id) > 1:
1624                for i in range(nn-1):
1625                    if nn > 2:
1626                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1627                            word_quantity += 'depth' + ','
1628                        else:
1629                            word_quantity += plot_quantity[i] + ', '
1630                    else:
1631                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1632                            word_quantity += 'depth'
1633                        else:
1634                            word_quantity += plot_quantity[i]
1635                    where1 = 0
1636                    count1 += 1
1637                    index = j*len(plot_quantity)
1638                    for which_quantity in plot_quantity:
1639                        where1 += 1
1640                        s = '\includegraphics' \
1641                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1642                            (graphname_report[index], '.png')
1643                        index += 1
1644                        fid.write(s)
1645                        if where1 % 2 == 0:
1646                            s = '\\\\ \n'
1647                            where1 = 0
1648                        else:
1649                            s = '& \n'
1650                        fid.write(s)
1651                word_quantity += ' and ' + plot_quantity[nn-1]           
1652                label = 'gauge%s' %(gaugeloc2) 
1653                caption = 'Time series for %s at %s location ' \
1654                          '(elevation %.2fm)' % \
1655                          (word_quantity, locations[k], elev[k])
1656                if elev[k] == 0.0:
1657                        caption = 'Time series for %s at %s location ' \
1658                                  '(elevation %.2fm)' % \
1659                                  (word_quantity, locations[k],
1660                                   elevations[0,k,j])
1661                        thisgauge = gauges[k]
1662                        east = thisgauge[0]
1663                        north = thisgauge[1]
1664                        elev_output.append([locations[k], east, north,
1665                                            elevations[0,k,j]])
1666                       
1667                s = '\end{tabular} \n' \
1668                    '\\caption{%s} \n' \
1669                    '\label{fig:%s} \n' \
1670                    '\end{figure} \n \n' % (caption, label)
1671                fid.write(s)
1672                if float((k+1)/div - pp) == 0.:
1673                    fid.write('\\clearpage \n')
1674                    pp += 1
1675                #### finished generating figures ###
1676
1677            close('all')
1678       
1679    return texfile2, elev_output
1680
1681
1682# FIXME (DSG): Add unit test, make general, not just 2 files,
1683# but any number of files.
1684##
1685# @brief ??
1686# @param dir_name ??
1687# @param filename1 ??
1688# @param filename2 ??
1689# @return ??
1690# @note TEMP
1691def copy_code_files(dir_name, filename1, filename2):
1692    """Temporary Interface to new location"""
1693
1694    from anuga.shallow_water.data_manager import \
1695                    copy_code_files as dm_copy_code_files
1696    print 'copy_code_files has moved from util.py.  ',
1697    print 'Please use "from anuga.shallow_water.data_manager \
1698                                        import copy_code_files"'
1699   
1700    return dm_copy_code_files(dir_name, filename1, filename2)
1701
1702
1703##
1704# @brief Create a nested sub-directory path.
1705# @param root_directory The base diretory path.
1706# @param directories An iterable of sub-directory names.
1707# @return The final joined directory path.
1708# @note If each sub-directory doesn't exist, it will be created.
1709def add_directories(root_directory, directories):
1710    """
1711    Add the first sub-directory in 'directories' to root_directory.
1712    Then add the second sub-directory to the accumulating path and so on.
1713
1714    Return the path of the final directory.
1715
1716    This is handy for specifying and creating a directory where data will go.
1717    """
1718    dir = root_directory
1719    for new_dir in directories:
1720        dir = os.path.join(dir, new_dir)
1721        if not access(dir,F_OK):
1722            mkdir(dir)
1723    return dir
1724
1725
1726##
1727# @brief
1728# @param filename
1729# @param separator_value
1730# @return
1731# @note TEMP
1732def get_data_from_file(filename, separator_value=','):
1733    """Temporary Interface to new location"""
1734    from anuga.shallow_water.data_manager import \
1735                        get_data_from_file as dm_get_data_from_file
1736    print 'get_data_from_file has moved from util.py'
1737    print 'Please use "from anuga.shallow_water.data_manager \
1738                                     import get_data_from_file"'
1739   
1740    return dm_get_data_from_file(filename,separator_value = ',')
1741
1742
1743##
1744# @brief
1745# @param verbose
1746# @param kwargs
1747# @return
1748# @note TEMP
1749def store_parameters(verbose=False,**kwargs):
1750    """Temporary Interface to new location"""
1751   
1752    from anuga.shallow_water.data_manager \
1753                    import store_parameters as dm_store_parameters
1754    print 'store_parameters has moved from util.py.'
1755    print 'Please use "from anuga.shallow_water.data_manager \
1756                                     import store_parameters"'
1757   
1758    return dm_store_parameters(verbose=False,**kwargs)
1759
1760
1761##
1762# @brief Remove vertices that are not associated with any triangle.
1763# @param verts An iterable (or array) of points.
1764# @param triangles An iterable of 3 element tuples.
1765# @param number_of_full_nodes ??
1766# @return (verts, triangles) where 'verts' has been updated.
1767def remove_lone_verts(verts, triangles, number_of_full_nodes=None):
1768    """Removes vertices that are not associated with any triangles.
1769
1770    verts is a list/array of points.
1771    triangles is a list of 3 element tuples.  Each tuple represents a triangle.
1772    number_of_full_nodes relate to parallelism when a mesh has an
1773        extra layer of ghost points.
1774    """
1775
1776    verts = ensure_numeric(verts)
1777    triangles = ensure_numeric(triangles)
1778   
1779    N = len(verts)
1780   
1781    # initialise the array to easily find the index of the first loner
1782    # ie, if N=3 -> [6,5,4]
1783    loners=num.arange(2*N, N, -1)
1784    for t in triangles:
1785        for vert in t:
1786            loners[vert]= vert # all non-loners will have loners[i]=i
1787
1788    lone_start = 2*N - max(loners) # The index of the first loner
1789
1790    if lone_start-1 == N:
1791        # no loners
1792        pass
1793    elif min(loners[lone_start:N]) > N:
1794        # All the loners are at the end of the vert array
1795        verts = verts[0:lone_start]
1796    else:
1797        # change the loners list so it can be used to modify triangles
1798        # Remove the loners from verts
1799        # Could've used X=compress(less(loners,N),loners)
1800        # verts=num.take(verts,X)  to Remove the loners from verts
1801        # but I think it would use more memory
1802        new_i = lone_start      # point at first loner - 'shuffle down' target
1803        for i in range(lone_start, N):
1804            if loners[i] >= N:  # [i] is a loner, leave alone
1805                pass
1806            else:               # a non-loner, move down
1807                loners[i] = new_i
1808                verts[new_i] = verts[i]
1809                new_i += 1
1810        verts = verts[0:new_i]
1811
1812        # Modify the triangles
1813        #print "loners", loners
1814        #print "triangles before", triangles
1815        triangles = num.choose(triangles,loners)
1816        #print "triangles after", triangles
1817    return verts, triangles
1818
1819
1820##
1821# @brief Compute centroid values from vertex values
1822# @param x Values at vertices of triangular mesh.
1823# @param triangles Nx3 integer array pointing to vertex information.
1824# @return [N] array of centroid values.
1825def get_centroid_values(x, triangles):
1826    """Compute centroid values from vertex values
1827   
1828    x: Values at vertices of triangular mesh
1829    triangles: Nx3 integer array pointing to vertex information
1830    for each of the N triangels. Elements of triangles are
1831    indices into x
1832    """
1833       
1834    xc = num.zeros(triangles.shape[0], num.float) # Space for centroid info
1835   
1836    for k in range(triangles.shape[0]):
1837        # Indices of vertices
1838        i0 = triangles[k][0]
1839        i1 = triangles[k][1]
1840        i2 = triangles[k][2]       
1841       
1842        xc[k] = (x[i0] + x[i1] + x[i2])/3
1843
1844    return xc
1845
1846
1847# @note TEMP
1848def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]},
1849                                output_dir='',
1850                                base_name='',
1851                                plot_numbers=['3-5'],
1852                                quantities=['speed','stage','momentum'],
1853                                assess_all_csv_files=True,
1854                                extra_plot_name='test'):
1855
1856    msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs '
1857    msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import ' \
1858           'csv2timeseries_graphs"'
1859    raise Exception, msg
1860
1861    return csv2timeseries_graphs(directories_dic,
1862                                 output_dir,
1863                                 base_name,
1864                                 plot_numbers,
1865                                 quantities,
1866                                 extra_plot_name,
1867                                 assess_all_csv_files)
1868
1869
1870##
1871# @brief Plot time series from CSV files.
1872# @param directories_dic
1873# @param output_dir
1874# @param base_name
1875# @param plot_numbers
1876# @param quantities
1877# @param extra_plot_name
1878# @param assess_all_csv_files
1879# @param create_latex
1880# @param verbose
1881# @note Assumes that 'elevation' is in the CSV file(s).
1882def csv2timeseries_graphs(directories_dic={},
1883                          output_dir='',
1884                          base_name=None,
1885                          plot_numbers='',
1886                          quantities=['stage'],
1887                          extra_plot_name='',
1888                          assess_all_csv_files=True,
1889                          create_latex=False,
1890                          verbose=False):
1891                               
1892    """
1893    Read in csv files that have the right header information and
1894    plot time series such as Stage, Speed, etc. Will also plot several
1895    time series on one plot. Filenames must follow this convention,
1896    <base_name><plot_number>.csv eg gauge_timeseries3.csv
1897   
1898    NOTE: relies that 'elevation' is in the csv file!
1899
1900    Each file represents a location and within each file there are
1901    time, quantity columns.
1902   
1903    For example:   
1904    if "directories_dic" defines 4 directories and in each directories
1905    there is a csv files corresponding to the right "plot_numbers",
1906    this will create a plot with 4 lines one for each directory AND
1907    one plot for each "quantities".  ??? FIXME: unclear.
1908   
1909    Usage:
1910        csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1911                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1912                            output_dir='fixed_wave'+sep,
1913                            base_name='gauge_timeseries_',
1914                            plot_numbers='',
1915                            quantities=['stage','speed'],
1916                            extra_plot_name='',
1917                            assess_all_csv_files=True,                           
1918                            create_latex=False,
1919                            verbose=True)
1920            this will create one plot for stage with both 'slide' and
1921            'fixed_wave' lines on it for stage and speed for each csv
1922            file with 'gauge_timeseries_' as the prefix. The graghs
1923            will be in the output directory 'fixed_wave' and the graph
1924            axis will be determined by assessing all the
1925   
1926    ANOTHER EXAMPLE
1927        new_csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1928                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1929                            output_dir='fixed_wave'+sep,
1930                            base_name='gauge_timeseries_',
1931                            plot_numbers=['1-3'],
1932                            quantities=['stage','speed'],
1933                            extra_plot_name='',
1934                            assess_all_csv_files=False,                           
1935                            create_latex=False,
1936                            verbose=True)
1937        This will plot csv files called gauge_timeseries_1.csv and
1938        gauge_timeseries3.csv from both 'slide' and 'fixed_wave' directories
1939        to 'fixed_wave'. There will be 4 plots created two speed and two stage
1940        one for each csv file. There will be two lines on each of these plots.
1941        And the axis will have been determined from only these files, had
1942        assess_all_csv_files = True all csv file with 'gauges_timeseries_' prefix
1943        would of been assessed.
1944   
1945    ANOTHER EXAMPLE   
1946         csv2timeseries_graphs({'J:'+sep+'anuga_validation'+sep:['new',20,-.1],
1947                                   'J:'+sep+'conical_island'+sep:['test',0,0]},
1948                                   output_dir='',
1949                                   plot_numbers=['1','3'],
1950                                   quantities=['stage','depth','bearing'],
1951                                   base_name='gauge_b',
1952                                   assess_all_csv_files=True,
1953                                  verbose=True)   
1954       
1955            This will produce one plot for each quantity (therefore 3) in the
1956            current directory, each plot will have 2 lines on them. The first
1957            plot named 'new' will have the time offseted by 20secs and the stage
1958            height adjusted by -0.1m
1959       
1960    Inputs:
1961        directories_dic: dictionary of directory with values (plot
1962                         legend name for directory), (start time of
1963                         the time series) and the (value to add to
1964                         stage if needed). For example
1965                         {dir1:['Anuga_ons',5000, 0],
1966                          dir2:['b_emoth',5000,1.5],
1967                          dir3:['b_ons',5000,1.5]}
1968                         Having multiple directories defined will plot them on
1969                         one plot, therefore there will be 3 lines on each of
1970                         these plot. If you only want one line per plot call
1971                         csv2timeseries_graph separately for each directory,
1972                         eg only have one directory in the 'directories_dic' in
1973                         each call.
1974                         
1975        output_dir: directory for the plot outputs. Only important to define when
1976                    you have more than one directory in your directories_dic, if
1977                    you have not defined it and you have multiple directories in
1978                    'directories_dic' there will be plots in each directory,
1979                    however only one directory will contain the complete
1980                    plot/graphs.
1981       
1982        base_name: Is used a couple of times.
1983                   1) to find the csv files to be plotted if there is no
1984                      'plot_numbers' then csv files with 'base_name' are plotted
1985                   2) in the title of the plots, the length of base_name is
1986                      removed from the front of the filename to be used in the
1987                      title.
1988                   This could be changed if needed.
1989                   Note is ignored if assess_all_csv_files=True
1990       
1991        plot_numbers: a String list of numbers to plot. For example
1992                      [0-4,10,15-17] will read and attempt to plot
1993                      the follow 0,1,2,3,4,10,15,16,17
1994                      NOTE: if no plot numbers this will create one plot per
1995                            quantity, per gauge
1996
1997        quantities: Will get available quantities from the header in the csv
1998                    file.  Quantities must be one of these.
1999                    NOTE: ALL QUANTITY NAMES MUST BE lower case!
2000                   
2001        extra_plot_name: A string that is appended to the end of the
2002                         output filename.
2003                   
2004        assess_all_csv_files: if true it will read ALL csv file with
2005                             "base_name", regardless of 'plot_numbers'
2006                              and determine a uniform set of axes for
2007                              Stage, Speed and Momentum. IF FALSE it
2008                              will only read the csv file within the
2009                             'plot_numbers'
2010                             
2011        create_latex: NOT IMPLEMENTED YET!! sorry Jane....
2012       
2013    OUTPUTS: saves the plots to
2014              <output_dir><base_name><plot_number><extra_plot_name>.png
2015    """
2016
2017    try: 
2018        import pylab
2019    except ImportError:
2020        msg='csv2timeseries_graphs needs pylab to be installed correctly'
2021        raise msg
2022            #ANUGA don't need pylab to work so the system doesn't
2023            #rely on pylab being installed
2024        return
2025
2026    from os import sep
2027    from anuga.shallow_water.data_manager import \
2028                               get_all_files_with_extension, csv2dict
2029
2030    seconds_in_hour = 3600
2031    seconds_in_minutes = 60
2032   
2033    quantities_label={}
2034#    quantities_label['time'] = 'time (hours)'
2035    quantities_label['time'] = 'time (minutes)'
2036    quantities_label['stage'] = 'wave height (m)'
2037    quantities_label['speed'] = 'speed (m/s)'
2038    quantities_label['momentum'] = 'momentum (m^2/sec)'
2039    quantities_label['depth'] = 'water depth (m)'
2040    quantities_label['xmomentum'] = 'momentum (m^2/sec)'
2041    quantities_label['ymomentum'] = 'momentum (m^2/sec)'
2042    quantities_label['bearing'] = 'degrees (o)'
2043    quantities_label['elevation'] = 'elevation (m)'
2044   
2045    if extra_plot_name != '':
2046        extra_plot_name = '_' + extra_plot_name
2047
2048    new_plot_numbers=[]
2049    #change plot_numbers to list, eg ['0-4','10']
2050    #to ['0','1','2','3','4','10']
2051    for i, num_string in enumerate(plot_numbers):
2052        if '-' in num_string: 
2053            start = int(num_string[:num_string.rfind('-')])
2054            end = int(num_string[num_string.rfind('-') + 1:]) + 1
2055            for x in range(start, end):
2056                new_plot_numbers.append(str(x))
2057        else:
2058            new_plot_numbers.append(num_string)
2059
2060    #finds all the files that fit the specs provided and return a list of them
2061    #so to help find a uniform max and min for the plots...
2062    list_filenames=[]
2063    all_csv_filenames=[]
2064    if verbose: print 'Determining files to access for axes ranges.'
2065   
2066    for i,directory in enumerate(directories_dic.keys()):
2067        all_csv_filenames.append(get_all_files_with_extension(directory,
2068                                                              base_name, '.csv'))
2069
2070        filenames=[]
2071        if plot_numbers == '': 
2072            list_filenames.append(get_all_files_with_extension(directory,
2073                                                               base_name,'.csv'))
2074        else:
2075            for number in new_plot_numbers:
2076                filenames.append(base_name + number)
2077            list_filenames.append(filenames)
2078
2079    #use all the files to get the values for the plot axis
2080    max_start_time= -1000.
2081    min_start_time = 100000 
2082   
2083    if verbose: print 'Determining uniform axes' 
2084
2085    #this entire loop is to determine the min and max range for the
2086    #axes of the plots
2087
2088#    quantities.insert(0,'elevation')
2089    quantities.insert(0,'time')
2090
2091    directory_quantity_value={}
2092#    quantity_value={}
2093    min_quantity_value={}
2094    max_quantity_value={}
2095
2096    for i, directory in enumerate(directories_dic.keys()):
2097        filename_quantity_value = {}
2098        if assess_all_csv_files == False:
2099            which_csv_to_assess = list_filenames[i]
2100        else:
2101            #gets list of filenames for directory "i"
2102            which_csv_to_assess = all_csv_filenames[i]
2103       
2104        for j, filename in enumerate(which_csv_to_assess):
2105            quantity_value = {}
2106
2107            dir_filename = join(directory,filename)
2108            attribute_dic, title_index_dic = csv2dict(dir_filename + '.csv')
2109            directory_start_time = directories_dic[directory][1]
2110            directory_add_tide = directories_dic[directory][2]
2111
2112            if verbose: print 'reading: %s.csv' %dir_filename
2113
2114            #add time to get values
2115            for k, quantity in enumerate(quantities):
2116                quantity_value[quantity] = [float(x) for
2117                                                x in attribute_dic[quantity]]
2118
2119                #add tide to stage if provided
2120                if quantity == 'stage':
2121                     quantity_value[quantity] = num.array(quantity_value[quantity], num.float) \
2122                                                          + directory_add_tide
2123
2124                #condition to find max and mins for all the plots
2125                # populate the list with something when i=0 and j=0 and
2126                # then compare to the other values to determine abs max and min
2127                if i==0 and j==0:
2128                    min_quantity_value[quantity], \
2129                        max_quantity_value[quantity] = \
2130                            get_min_max_values(quantity_value[quantity])
2131
2132                    if quantity != 'time':
2133                        min_quantity_value[quantity] = \
2134                            min_quantity_value[quantity] *1.1
2135                        max_quantity_value[quantity] = \
2136                            max_quantity_value[quantity] *1.1
2137                else:
2138                    min, max = get_min_max_values(quantity_value[quantity])
2139               
2140                    # min and max are multipled by "1+increase_axis" to get axes
2141                    # that are slighty bigger than the max and mins
2142                    # so the plots look good.
2143
2144                    increase_axis = (max-min)*0.05
2145                    if min <= min_quantity_value[quantity]:
2146                        if quantity == 'time': 
2147                            min_quantity_value[quantity] = min
2148                        else:
2149                            if round(min,2) == 0.00:
2150                                min_quantity_value[quantity] = -increase_axis
2151#                                min_quantity_value[quantity] = -2.
2152                                #min_quantity_value[quantity] = \
2153                                #    -max_quantity_value[quantity]*increase_axis
2154                            else:
2155#                                min_quantity_value[quantity] = \
2156#                                    min*(1+increase_axis)
2157                                min_quantity_value[quantity]=min-increase_axis
2158                   
2159                    if max > max_quantity_value[quantity]: 
2160                        if quantity == 'time': 
2161                            max_quantity_value[quantity] = max
2162                        else:
2163                            max_quantity_value[quantity] = max + increase_axis
2164#                            max_quantity_value[quantity]=max*(1+increase_axis)
2165
2166            #set the time... ???
2167            if min_start_time > directory_start_time: 
2168                min_start_time = directory_start_time
2169            if max_start_time < directory_start_time: 
2170                max_start_time = directory_start_time
2171           
2172            filename_quantity_value[filename]=quantity_value
2173           
2174        directory_quantity_value[directory]=filename_quantity_value
2175   
2176    #final step to unifrom axis for the graphs
2177    quantities_axis={}
2178   
2179    for i, quantity in enumerate(quantities):
2180        quantities_axis[quantity] = (float(min_start_time) \
2181                                         / float(seconds_in_minutes),
2182                                     (float(max_quantity_value['time']) \
2183                                          + float(max_start_time)) \
2184                                              / float(seconds_in_minutes),
2185                                     min_quantity_value[quantity],
2186                                     max_quantity_value[quantity])
2187
2188        if verbose and (quantity != 'time' and quantity != 'elevation'): 
2189            print 'axis for quantity %s are x:(%s to %s)%s and y:(%s to %s)%s' \
2190                  % (quantity, 
2191                     quantities_axis[quantity][0],
2192                     quantities_axis[quantity][1],
2193                     quantities_label['time'],
2194                     quantities_axis[quantity][2],
2195                     quantities_axis[quantity][3],
2196                     quantities_label[quantity])
2197
2198    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
2199
2200    if verbose: print 'Now start to plot \n'
2201   
2202    i_max = len(directories_dic.keys())
2203    legend_list_dic = {}
2204    legend_list = []
2205    for i, directory in enumerate(directories_dic.keys()):
2206        if verbose: print 'Plotting in %s %s' % (directory, new_plot_numbers)
2207
2208        # FIXME THIS SORT IS VERY IMPORTANT
2209        # Without it the assigned plot numbers may not work correctly
2210        # there must be a better way
2211        list_filenames[i].sort()
2212        for j, filename in enumerate(list_filenames[i]):
2213            if verbose: print 'Starting %s' % filename 
2214
2215            directory_name = directories_dic[directory][0]
2216            directory_start_time = directories_dic[directory][1]
2217            directory_add_tide = directories_dic[directory][2]
2218           
2219            # create an if about the start time and tide height if don't exist
2220            attribute_dic, title_index_dic = csv2dict(directory + sep
2221                                                      + filename + '.csv')
2222            #get data from dict in to list
2223            #do maths to list by changing to array
2224            t = (num.array(directory_quantity_value[directory][filename]['time'])
2225                     + directory_start_time) / seconds_in_minutes
2226
2227            #finds the maximum elevation, used only as a test
2228            # and as info in the graphs
2229            max_ele=-100000
2230            min_ele=100000
2231            elevation = [float(x) for x in attribute_dic["elevation"]]
2232           
2233            min_ele, max_ele = get_min_max_values(elevation)
2234           
2235            if min_ele != max_ele:
2236                print "Note! Elevation changes in %s" %dir_filename
2237
2238            # creates a dictionary with keys that is the filename and attributes
2239            # are a list of lists containing 'directory_name' and 'elevation'.
2240            # This is used to make the contents for the legends in the graphs,
2241            # this is the name of the model and the elevation.  All in this
2242            # great one liner from DG. If the key 'filename' doesn't exist it
2243            # creates the entry if the entry exist it appends to the key.
2244
2245            legend_list_dic.setdefault(filename,[]) \
2246                .append([directory_name, round(max_ele, 3)])
2247
2248            # creates a LIST for the legend on the last iteration of the
2249            # directories which is when "legend_list_dic" has been fully
2250            # populated. Creates a list of strings which is used in the legend
2251            # only runs on the last iteration for all the gauges(csv) files
2252            # empties the list before creating it
2253
2254            if i == i_max - 1:
2255                legend_list = []
2256   
2257                for name_and_elevation in legend_list_dic[filename]:
2258                    legend_list.append('%s (elevation = %sm)'\
2259                                       % (name_and_elevation[0],
2260                                          name_and_elevation[1]))
2261           
2262            #skip time and elevation so it is not plotted!
2263            for k, quantity in enumerate(quantities):
2264                if quantity != 'time' and quantity != 'elevation':
2265                    pylab.figure(int(k*100+j))
2266                    pylab.ylabel(quantities_label[quantity])
2267                    pylab.plot(t,
2268                               directory_quantity_value[directory]\
2269                                                       [filename][quantity],
2270                               c = cstr[i], linewidth=1)
2271                    pylab.xlabel(quantities_label['time'])
2272                    pylab.axis(quantities_axis[quantity])
2273                    pylab.legend(legend_list,loc='upper right')
2274                   
2275                    pylab.title('%s at %s gauge'
2276                                % (quantity, filename[len(base_name):]))
2277
2278                    if output_dir == '':
2279                        figname = '%s%s%s_%s%s.png' \
2280                                  % (directory, sep, filename, quantity,
2281                                     extra_plot_name)
2282                    else:
2283                        figname = '%s%s%s_%s%s.png' \
2284                                  % (output_dir, sep, filename, quantity,
2285                                     extra_plot_name)
2286
2287                    if verbose: print 'saving figure here %s' %figname
2288
2289                    pylab.savefig(figname)
2290           
2291    if verbose: print 'Closing all plots'
2292
2293    pylab.close('all')
2294    del pylab
2295
2296    if verbose: print 'Finished closing plots'
2297
2298##
2299# @brief Return min and max of an iterable.
2300# @param list The iterable to return min & max of.
2301# @return (min, max) of 'list'.
2302def get_min_max_values(list=None):
2303    """
2304    Returns the min and max of the list it was provided.
2305    """
2306
2307    if list == None: print 'List must be provided'
2308       
2309    return min(list), max(list)
2310
2311
2312##
2313# @brief Get runup around a point in a CSV file.
2314# @param gauge_filename gauge file name.
2315# @param sww_filename SWW file name.
2316# @param runup_filename Name of file to report into.
2317# @param size ??
2318# @param verbose ??
2319def get_runup_data_for_locations_from_file(gauge_filename,
2320                                           sww_filename,
2321                                           runup_filename,
2322                                           size=10,
2323                                           verbose=False):
2324    """this will read a csv file with the header x,y. Then look in a square
2325    'size'x2 around this position for the 'max_inundaiton_height' in the
2326    'sww_filename' and report the findings in the 'runup_filename'.
2327   
2328    WARNING: NO TESTS!
2329    """
2330
2331    from anuga.shallow_water.data_manager import get_all_directories_with_name,\
2332                                                 get_maximum_inundation_data,\
2333                                                 csv2dict
2334                                                 
2335    file = open(runup_filename, "w")
2336    file.write("easting,northing,runup \n ")
2337    file.close()
2338   
2339    #read gauge csv file to dictionary
2340    attribute_dic, title_index_dic = csv2dict(gauge_filename)
2341    northing = [float(x) for x in attribute_dic["y"]]
2342    easting = [float(x) for x in attribute_dic["x"]]
2343
2344    print 'Reading %s' %sww_filename
2345
2346    runup_locations=[]
2347    for i, x in enumerate(northing):
2348        poly = [[int(easting[i]+size),int(northing[i]+size)],
2349                [int(easting[i]+size),int(northing[i]-size)],
2350                [int(easting[i]-size),int(northing[i]-size)],
2351                [int(easting[i]-size),int(northing[i]+size)]]
2352       
2353        run_up, x_y = get_maximum_inundation_data(filename=sww_filename,
2354                                                  polygon=poly,
2355                                                  verbose=False) 
2356
2357        #if no runup will return 0 instead of NONE
2358        if run_up==None: run_up=0
2359        if x_y==None: x_y=[0,0]
2360       
2361        if verbose:
2362            print 'maximum inundation runup near %s is %s meters' %(x_y,run_up)
2363       
2364        #writes to file
2365        file = open(runup_filename, "a")
2366        temp = '%s,%s,%s \n' % (x_y[0], x_y[1], run_up)
2367        file.write(temp)
2368        file.close()
2369
2370
2371##
2372# @brief ??
2373# @param sww_file ??
2374# @param gauge_file ??
2375# @param out_name ??
2376# @param quantities ??
2377# @param verbose ??
2378# @param use_cache ??
2379def sww2csv_gauges(sww_file,
2380                   gauge_file,
2381                   out_name='gauge_',
2382                   quantities=['stage', 'depth', 'elevation',
2383                               'xmomentum', 'ymomentum'],
2384                   verbose=False,
2385                   use_cache = True):
2386    """
2387   
2388    Inputs:
2389       
2390        NOTE: if using csv2timeseries_graphs after creating csv file,
2391        it is essential to export quantities 'depth' and 'elevation'.
2392        'depth' is good to analyse gauges on land and elevation is used
2393        automatically by csv2timeseries_graphs in the legend.
2394       
2395        sww_file: path to any sww file
2396       
2397        gauge_file: Assumes that it follows this format
2398            name, easting, northing, elevation
2399            point1, 100.3, 50.2, 10.0
2400            point2, 10.3, 70.3, 78.0
2401       
2402        NOTE: order of column can change but names eg 'easting', 'elevation'
2403        must be the same! ALL lowercaps!
2404
2405        out_name: prefix for output file name (default is 'gauge_')
2406       
2407    Outputs:
2408        one file for each gauge/point location in the points file. They
2409        will be named with this format in the same directory as the 'sww_file'
2410            <out_name><name>.csv
2411        eg gauge_point1.csv if <out_name> not supplied
2412           myfile_2_point1.csv if <out_name> ='myfile_2_'
2413           
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#    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
2440    #print "points",points
2441
2442    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
2443    assert type(out_name) == type(''), 'Output filename prefix must be a string'
2444   
2445    try:
2446        point_reader = reader(file(gauge_file))
2447    except Exception, e:
2448        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e)
2449        raise msg
2450
2451    if verbose: print '\n Gauges obtained from: %s \n' %gauge_file
2452   
2453    point_reader = reader(file(gauge_file))
2454    points = []
2455    point_name = []
2456   
2457    #read point info from file
2458    for i,row in enumerate(point_reader):
2459        #read header and determine the column numbers to read correcty.
2460        if i==0:
2461            for j,value in enumerate(row):
2462                if value.strip()=='easting':easting=j
2463                if value.strip()=='northing':northing=j
2464                if value.strip()=='name':name=j
2465                if value.strip()=='elevation':elevation=j
2466        else:
2467            points.append([float(row[easting]),float(row[northing])])
2468            point_name.append(row[name])
2469       
2470    #convert to array for file_function
2471    points_array = num.array(points,num.float)
2472       
2473    points_array = ensure_absolute(points_array)
2474
2475    dir_name, base = os.path.split(sww_file)   
2476
2477    #need to get current directory so when path and file
2478    #are "joined" below the directory is correct
2479    if dir_name == '':
2480        dir_name =getcwd()
2481       
2482    if access(sww_file,R_OK):
2483        if verbose: print 'File %s exists' %(sww_file)
2484    else:
2485        msg = 'File "%s" could not be opened: no read permission' % sww_file
2486        raise msg
2487
2488    sww_files = get_all_swwfiles(look_in_dir=dir_name,
2489                                 base_name=base,
2490                                 verbose=verbose)
2491   
2492    #to make all the quantities lower case for file_function
2493    quantities = [quantity.lower() for quantity in quantities]
2494
2495    # what is quantities are needed from sww file to calculate output quantities
2496    # also
2497
2498    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
2499   
2500    for sww_file in sww_files:
2501        sww_file = join(dir_name, sww_file+'.sww')
2502        callable_sww = file_function(sww_file,
2503                                     quantities=core_quantities,
2504                                     interpolation_points=points_array,
2505                                     verbose=verbose,
2506                                     use_cache=use_cache)
2507    gauge_file = out_name
2508
2509    heading = [quantity for quantity in quantities]
2510    heading.insert(0,'time')
2511
2512    #create a list of csv writers for all the points and write header
2513    points_writer = []
2514    for i,point in enumerate(points):
2515        points_writer.append(writer(file(dir_name + sep + gauge_file
2516                                         + point_name[i] + '.csv', "wb")))
2517        points_writer[i].writerow(heading)
2518   
2519    if verbose: print 'Writing csv files'
2520
2521    for time in callable_sww.get_time():
2522        for point_i, point in enumerate(points_array):
2523            #add domain starttime to relative time.
2524            points_list = [time + callable_sww.starttime]
2525            point_quantities = callable_sww(time,point_i)
2526           
2527            for quantity in quantities:
2528                if quantity == NAN:
2529                    print 'quantity does not exist in' % callable_sww.get_name
2530                else:
2531                    if quantity == 'stage':
2532                        points_list.append(point_quantities[0])
2533                       
2534                    if quantity == 'elevation':
2535                        points_list.append(point_quantities[1])
2536                       
2537                    if quantity == 'xmomentum':
2538                        points_list.append(point_quantities[2])
2539                       
2540                    if quantity == 'ymomentum':
2541                        points_list.append(point_quantities[3])
2542                       
2543                    if quantity == 'depth':
2544                        points_list.append(point_quantities[0] 
2545                                           - point_quantities[1])
2546
2547                    if quantity == 'momentum':
2548                        momentum = sqrt(point_quantities[2]**2 
2549                                        + point_quantities[3]**2)
2550                        points_list.append(momentum)
2551                       
2552                    if quantity == 'speed':
2553                        #if depth is less than 0.001 then speed = 0.0
2554                        if point_quantities[0] - point_quantities[1] < 0.001:
2555                            vel = 0.0
2556                        else:
2557                            if point_quantities[2] < 1.0e6:
2558                                momentum = sqrt(point_quantities[2]**2
2559                                                + point_quantities[3]**2)
2560    #                            vel = momentum/depth             
2561                                vel = momentum / (point_quantities[0] 
2562                                                  - point_quantities[1])
2563    #                            vel = momentum/(depth + 1.e-6/depth)
2564                            else:
2565                                momentum = 0
2566                                vel = 0
2567                           
2568                        points_list.append(vel)
2569                       
2570                    if quantity == 'bearing':
2571                        points_list.append(calc_bearing(point_quantities[2],
2572                                                        point_quantities[3]))
2573
2574            points_writer[point_i].writerow(points_list)
2575       
2576
2577##
2578# @brief Get a wave height at a certain depth given wave height at another depth.
2579# @param d1 The first depth.
2580# @param d2 The second depth.
2581# @param h1 Wave ampitude at d1
2582# @param verbose True if this function is to be verbose.
2583# @return The wave height at d2.
2584def greens_law(d1, d2, h1, verbose=False):
2585    """Green's Law
2586
2587    Green's Law allows an approximation of wave amplitude at
2588    a given depth based on the fourh root of the ratio of two depths
2589    and the amplitude at another given depth.
2590
2591    Note, wave amplitude is equal to stage.
2592   
2593    Inputs:
2594
2595    d1, d2 - the two depths
2596    h1     - the wave amplitude at d1
2597    h2     - the derived amplitude at d2
2598
2599    h2 = h1 (d1/d2)^(1/4), where d2 cannot equal 0.
2600   
2601    """
2602
2603    d1 = ensure_numeric(d1)
2604    d2 = ensure_numeric(d2)
2605    h1 = ensure_numeric(h1)
2606
2607    if d1 <= 0.0:
2608        msg = 'the first depth, d1 (%f), must be strictly positive' % (d1)
2609        raise Exception(msg)
2610
2611    if d2 <= 0.0:
2612        msg = 'the second depth, d2 (%f), must be strictly positive' % (d2)
2613        raise Exception(msg)
2614   
2615    if h1 <= 0.0:
2616        msg = 'the wave amplitude, h1 (%f), must be strictly positive' % (h1)
2617        raise Exception(msg)
2618   
2619    h2 = h1*(d1/d2)**0.25
2620
2621    assert h2 > 0
2622   
2623    return h2
2624       
2625
2626##
2627# @brief Get the square-root of a value.
2628# @param s The value to get the square-root of.
2629# @return The square-root of 's'.
2630def square_root(s):
2631    return sqrt(s)
2632
2633
Note: See TracBrowser for help on using the repository browser.