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

Last change on this file since 4983 was 4983, checked in by nick, 16 years ago

cleaned up old code

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