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

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

Changed .array(A) to .array(A, num.Int) where appropriate. Helps when going to numpy.

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