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

Last change on this file since 5016 was 5016, checked in by ole, 16 years ago

Added IP tracking to test_all and updated create_distribution

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