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

Last change on this file since 4634 was 4634, checked in by nick, 17 years ago

New "make_plots_from_csv_files" has been added, read doc string for info on usage. Plus added "get_min_max_values" from a list and unit tests for both.

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