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

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

added a function to easily get the runup data near a several certain locations

File size: 73.2 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    if isinstance(line11[0],str) is True:
973        # We have found text in the first line
974        east_index = None
975        north_index = None
976        name_index = None
977        elev_index = None
978        for i in range(len(line11)):
979            if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'easting': east_index = i
980            if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'northing': north_index = i
981            if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'name': name_index = i
982            if line11[i].strip('\n').strip('\r').strip(' ').lower() == 'elevation': elev_index = i
983
984        if east_index < len(line11) and north_index < len(line11):
985            pass
986        else:
987            msg = 'WARNING: %s does not contain correct header information' %(filename)
988            msg += 'The header must be: easting, northing, name, elevation'
989            raise Exception, msg
990
991        if elev_index is None: 
992            raise Exception
993   
994        if name_index is None: 
995            raise Exception
996
997        lines = lines[1:] # Remove header from data
998    else:
999        # No header, assume that this is a simple easting, northing file
1000
1001        msg = 'There was no header in file %s and the number of columns is %d' %(filename, len(line11))
1002        msg += '- I was assuming two columns corresponding to Easting and Northing'
1003        assert len(line11) == 2, msg
1004
1005        east_index = 0
1006        north_index = 1
1007
1008        N = len(lines)
1009        elev = [-9999]*N
1010        gaugelocation = range(N)
1011       
1012    # Read in gauge data
1013    for line in lines:
1014        fields = line.split(',')
1015
1016        gauges.append([float(fields[east_index]), float(fields[north_index])])
1017
1018        if len(fields) > 2:
1019            elev.append(float(fields[elev_index]))
1020            loc = fields[name_index]
1021            gaugelocation.append(loc.strip('\n'))
1022
1023    return gauges, gaugelocation, elev
1024
1025
1026
1027def check_list(quantity):
1028    """ Check that input quantities in quantity list are possible
1029    """
1030    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
1031                    'ymomentum', 'speed', 'bearing', 'elevation']
1032
1033       
1034                   
1035    import sys
1036    #if not sys.version.startswith('2.4'):
1037        # Backwards compatibility
1038    #   from sets import Set as set
1039
1040    from sets import Set as set
1041           
1042    for i,j in enumerate(quantity):
1043        quantity[i] = quantity[i].lower()
1044    p = list(set(quantity).difference(set(all_quantity)))
1045    if len(p) <> 0:
1046        msg = 'Quantities %s do not exist - please try again' %p
1047        raise Exception, msg
1048       
1049    return 
1050
1051def calc_bearing(uh, vh):
1052    """ Calculate velocity bearing from North
1053    """
1054    from math import atan, degrees
1055   
1056    angle = degrees(atan(vh/(uh+1.e-15)))
1057    if (0 < angle < 90.0):
1058        if vh > 0:
1059            bearing = 90.0 - abs(angle)
1060        if vh < 0:
1061            bearing = 270.0 - abs(angle)
1062    if (-90 < angle < 0):
1063        if vh < 0:
1064            bearing = 90.0 - abs(angle)
1065        if vh > 0:
1066            bearing = 270.0 - abs(angle)
1067    if angle == 0: bearing = 0.0
1068           
1069    return bearing
1070
1071def generate_figures(plot_quantity, file_loc, report, reportname, surface,
1072                     leg_label, f_list, gauges, locations, elev, gauge_index,
1073                     production_dirs, time_min, time_max, time_unit,
1074                     title_on, label_id, generate_fig, verbose):
1075    """ Generate figures based on required quantities and gauges for each sww file
1076    """
1077    from math import sqrt, atan, degrees
1078    from Numeric import ones, allclose, zeros, Float, ravel
1079    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
1080
1081    if generate_fig is True:
1082        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
1083             xlabel, ylabel, title, close, subplot
1084   
1085        if surface is True:
1086            import pylab as p1
1087            import mpl3d.mplot3d as p3
1088       
1089    if report == True:   
1090        texdir = getcwd()+sep+'report'+sep
1091        if access(texdir,F_OK) == 0:
1092            mkdir (texdir)
1093        if len(label_id) == 1:
1094            label_id1 = label_id[0].replace(sep,'')
1095            label_id2 = label_id1.replace('_','')
1096            texfile = texdir+reportname+'%s' %(label_id2)
1097            texfile2 = reportname+'%s' %(label_id2)
1098            texfilename = texfile + '.tex'
1099            if verbose: print '\n Latex output printed to %s \n' %texfilename
1100            fid = open(texfilename, 'w')
1101        else:
1102            texfile = texdir+reportname
1103            texfile2 = reportname
1104            texfilename = texfile + '.tex' 
1105            if verbose: print '\n Latex output printed to %s \n' %texfilename
1106            fid = open(texfilename, 'w')
1107    else:
1108        texfile = ''
1109        texfile2 = ''
1110
1111    p = len(f_list)
1112    n = []
1113    n0 = 0
1114    for i in range(len(f_list)):
1115        n.append(len(f_list[i].get_time()))
1116        if n[i] > n0: n0 = n[i] 
1117    n0 = int(n0)
1118    m = len(locations)
1119    model_time = zeros((n0,m,p), Float) 
1120    stages = zeros((n0,m,p), Float)
1121    elevations = zeros((n0,m,p), Float) 
1122    momenta = zeros((n0,m,p), Float)
1123    xmom = zeros((n0,m,p), Float)
1124    ymom = zeros((n0,m,p), Float)
1125    speed = zeros((n0,m,p), Float)
1126    bearings = zeros((n0,m,p), Float)
1127    depths = zeros((n0,m,p), Float)
1128    eastings = zeros((n0,m,p), Float)
1129    min_stages = []
1130    max_stages = []
1131    min_momentums = []   
1132    max_momentums = []
1133    max_xmomentums = []
1134    max_ymomentums = []
1135    min_xmomentums = []
1136    min_ymomentums = []
1137    max_speeds = []
1138    min_speeds = []   
1139    max_depths = []
1140    model_time_plot3d = zeros((n0,m), Float)
1141    stages_plot3d = zeros((n0,m), Float)
1142    eastings_plot3d = zeros((n0,m),Float)
1143    if time_unit is 'mins': scale = 60.0
1144    if time_unit is 'hours': scale = 3600.0
1145    ##### loop over each swwfile #####
1146    for j, f in enumerate(f_list):
1147        starttime = f.starttime
1148        if verbose: print 'swwfile %d of %d' %(j, len(f_list))
1149        comparefile = file_loc[j]+sep+'gauges_maxmins'+'.csv'
1150        fid_compare = open(comparefile, 'w')
1151        file0 = file_loc[j]+'gauges_t0.csv'
1152        fid_0 = open(file0, 'w')
1153        ##### loop over each gauge #####
1154        for k in gauge_index:
1155            g = gauges[k]
1156            if verbose: print 'Gauge %d of %d' %(k, len(gauges))
1157            min_stage = 10
1158            max_stage = 0
1159            max_momentum = max_xmomentum = max_ymomentum = 0
1160            min_momentum = min_xmomentum = min_ymomentum = 100
1161            max_speed = 0
1162            min_speed = 0           
1163            max_depth = 0           
1164            gaugeloc = str(locations[k])
1165            thisfile = file_loc[j]+sep+'gauges_time_series'+'_'\
1166                       +gaugeloc+'.csv'
1167            fid_out = open(thisfile, 'w')
1168            s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom \n'
1169            fid_out.write(s)
1170            #### generate quantities #######
1171            for i, t in enumerate(f.get_time()):
1172                if time_min <= t <= time_max:
1173                    w = f(t, point_id = k)[0]
1174                    z = f(t, point_id = k)[1]
1175                    uh = f(t, point_id = k)[2]
1176                    vh = f(t, point_id = k)[3]
1177                    depth = w-z     
1178                    m = sqrt(uh*uh + vh*vh)
1179                    if depth < 0.001:
1180                        vel = 0.0
1181                    else:
1182                        vel = m / (depth + 1.e-6/depth) 
1183                    #bearing = calc_bearing(uh, vh)                   
1184                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
1185                    stages[i,k,j] = w
1186                    elevations[i,k,j] = z
1187                    xmom[i,k,j] = uh
1188                    ymom[i,k,j] = vh
1189                    momenta[i,k,j] = m
1190                    speed[i,k,j] = vel
1191                    #bearings[i,k,j] = bearing
1192                    depths[i,k,j] = depth
1193                    thisgauge = gauges[k]
1194                    eastings[i,k,j] = thisgauge[0]
1195                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f\n' %(t, w, m, vel, z, uh, vh)
1196                    fid_out.write(s)
1197                    if t == 0:
1198                        s = '%.2f, %.2f, %.2f\n' %(g[0], g[1], w)
1199                        fid_0.write(s)
1200                    if t/60.0 <= 13920: tindex = i
1201                    if w > max_stage: max_stage = w
1202                    if w < min_stage: min_stage = w
1203                    if m > max_momentum: max_momentum = m
1204                    if m < min_momentum: min_momentum = m                   
1205                    if uh > max_xmomentum: max_xmomentum = uh
1206                    if vh > max_ymomentum: max_ymomentum = vh
1207                    if uh < min_xmomentum: min_xmomentum = uh
1208                    if vh < min_ymomentum: min_ymomentum = vh
1209                    if vel > max_speed: max_speed = vel
1210                    if vel < min_speed: min_speed = vel                   
1211                    if z > 0 and depth > max_depth: max_depth = depth
1212                   
1213                   
1214            s = '%.2f, %.2f, %.2f, %.2f, %s\n' %(max_stage, min_stage, z, thisgauge[0], leg_label[j])
1215            fid_compare.write(s)
1216            max_stages.append(max_stage)
1217            min_stages.append(min_stage)
1218            max_momentums.append(max_momentum)
1219            max_xmomentums.append(max_xmomentum)
1220            max_ymomentums.append(max_ymomentum)
1221            min_xmomentums.append(min_xmomentum)
1222            min_ymomentums.append(min_ymomentum)
1223            min_momentums.append(min_momentum)           
1224            max_depths.append(max_depth)
1225            max_speeds.append(max_speed)
1226            min_speeds.append(min_speed)           
1227            #### finished generating quantities for each swwfile #####
1228       
1229        model_time_plot3d[:,:] = model_time[:,:,j]
1230        stages_plot3d[:,:] = stages[:,:,j]
1231        eastings_plot3d[:,] = eastings[:,:,j]
1232           
1233        if surface is  True:
1234            print 'Printing surface figure'
1235            for i in range(2):
1236                fig = p1.figure(10)
1237                ax = p3.Axes3D(fig)
1238                if len(gauges) > 80:
1239                    ax.plot_surface(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
1240                else:
1241                    #ax.plot_wireframe(model_time[:,:,j],eastings[:,:,j],stages[:,:,j])
1242                    ax.plot3D(ravel(eastings[:,:,j]),ravel(model_time[:,:,j]),ravel(stages[:,:,j]))
1243                ax.set_xlabel('time')
1244                ax.set_ylabel('x')
1245                ax.set_zlabel('stage')
1246                fig.add_axes(ax)
1247                p1.show()
1248                surfacefig = 'solution_surface%s' %leg_label[j]
1249                p1.savefig(surfacefig)
1250                p1.close()
1251           
1252    #### finished generating quantities for all swwfiles #####
1253
1254    # x profile for given time
1255    if surface is True:
1256        figure(11)
1257        plot(eastings[tindex,:,j],stages[tindex,:,j])
1258        xlabel('x')
1259        ylabel('stage')
1260        profilefig = 'solution_xprofile' 
1261        savefig('profilefig')
1262
1263    elev_output = []
1264    if generate_fig is True:
1265        depth_axis = axis([starttime/scale, time_max/scale, -0.1, max(max_depths)*1.1])
1266        stage_axis = axis([starttime/scale, time_max/scale, min(min_stages), max(max_stages)*1.1])
1267        vel_axis = axis([starttime/scale, time_max/scale, min(min_speeds), max(max_speeds)*1.1])
1268        mom_axis = axis([starttime/scale, time_max/scale, min(min_momentums), max(max_momentums)*1.1])
1269        xmom_axis = axis([starttime/scale, time_max/scale, min(min_xmomentums), max(max_xmomentums)*1.1])
1270        ymom_axis = axis([starttime/scale, time_max/scale, min(min_ymomentums), max(max_ymomentums)*1.1])
1271        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
1272        nn = len(plot_quantity)
1273        no_cols = 2
1274       
1275        if len(label_id) > 1: graphname_report = []
1276        pp = 1
1277        div = 11.
1278        cc = 0
1279        for k in gauge_index:
1280            g = gauges[k]
1281            count1 = 0
1282            if report == True and len(label_id) > 1:
1283                s = '\\begin{figure}[ht] \n \\centering \n \\begin{tabular}{cc} \n'
1284                fid.write(s)
1285            if len(label_id) > 1: graphname_report = []
1286            #### generate figures for each gauge ####
1287            for j, f in enumerate(f_list):
1288                ion()
1289                hold(True)
1290                count = 0
1291                where1 = 0
1292                where2 = 0
1293                word_quantity = ''
1294                if report == True and len(label_id) == 1:
1295                    s = '\\begin{figure}[hbt] \n \\centering \n \\begin{tabular}{cc} \n'
1296                    fid.write(s)
1297                   
1298                for which_quantity in plot_quantity:
1299                    count += 1
1300                    where1 += 1
1301                    figure(count, frameon = False)
1302                    if which_quantity == 'depth':
1303                        plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j])
1304                        units = 'm'
1305                        axis(depth_axis)
1306                    if which_quantity == 'stage':
1307                        if elevations[0,k,j] <= 0:
1308                            plot(model_time[0:n[j]-1,k,j], stages[0:n[j]-1,k,j], '-', c = cstr[j])
1309                            axis(stage_axis)
1310                        else:
1311                            plot(model_time[0:n[j]-1,k,j], depths[0:n[j]-1,k,j], '-', c = cstr[j])
1312                            #axis(depth_axis)                 
1313                        units = 'm'
1314                    if which_quantity == 'momentum':
1315                        plot(model_time[0:n[j]-1,k,j], momenta[0:n[j]-1,k,j], '-', c = cstr[j])
1316                        axis(mom_axis)
1317                        units = 'm^2 / sec'
1318                    if which_quantity == 'xmomentum':
1319                        plot(model_time[0:n[j]-1,k,j], xmom[0:n[j]-1,k,j], '-', c = cstr[j])
1320                        axis(xmom_axis)
1321                        units = 'm^2 / sec'
1322                    if which_quantity == 'ymomentum':
1323                        plot(model_time[0:n[j]-1,k,j], ymom[0:n[j]-1,k,j], '-', c = cstr[j])
1324                        axis(ymom_axis)
1325                        units = 'm^2 / sec'
1326                    if which_quantity == 'speed':
1327                        plot(model_time[0:n[j]-1,k,j], speed[0:n[j]-1,k,j], '-', c = cstr[j])
1328                        axis(vel_axis)
1329                        units = 'm / sec'
1330                    if which_quantity == 'bearing':
1331                        due_east = 90.0*ones(shape(model_time[0:n[j]-1,k,j],Float))
1332                        due_west = 270.0*ones(shape(model_time[0:n[j]-1,k,j],Float))
1333                        plot(model_time[0:n[j]-1,k,j], bearings, '-', 
1334                             model_time[0:n[j]-1,k,j], due_west, '-.', 
1335                             model_time[0:n[j]-1,k,j], due_east, '-.')
1336                        units = 'degrees from North'
1337                        ax = axis([time_min, time_max, 0.0, 360.0])
1338                        legend(('Bearing','West','East'))
1339
1340                    if time_unit is 'mins': xlabel('time (mins)')
1341                    if time_unit is 'hours': xlabel('time (hours)')
1342                    #if which_quantity == 'stage' and elevations[0:n[j]-1,k,j] > 0:
1343                    #    ylabel('%s (%s)' %('depth', units))
1344                    #else:
1345                    #    ylabel('%s (%s)' %(which_quantity, units))
1346                        #ylabel('%s (%s)' %('wave height', units))
1347                    ylabel('%s (%s)' %(which_quantity, units))
1348                    if len(label_id) > 1: legend((leg_label),loc='upper right')
1349
1350                    #gaugeloc1 = gaugeloc.replace(' ','')
1351                    #gaugeloc2 = gaugeloc1.replace('_','')
1352                    gaugeloc2 = str(locations[k]).replace(' ','')
1353                    graphname = '%sgauge%s_%s' %(file_loc[j],
1354                                                 gaugeloc2,
1355                                                 which_quantity)
1356
1357                    if report == True and len(label_id) > 1:
1358                        figdir = getcwd()+sep+'report_figures'+sep
1359                        if access(figdir,F_OK) == 0 :
1360                            mkdir (figdir)
1361                        latex_file_loc = figdir.replace(sep,altsep) 
1362                        graphname_latex = '%sgauge%s%s' %(latex_file_loc, gaugeloc2, which_quantity) # storing files in production directory   
1363                        graphname_report_input = '%sgauge%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity) # giving location in latex output file
1364                        graphname_report.append(graphname_report_input)
1365                       
1366                        savefig(graphname_latex) # save figures in production directory for report generation
1367
1368                    if report == True:
1369
1370                        figdir = getcwd()+sep+'report_figures'+sep
1371                        if access(figdir,F_OK) == 0 :
1372                            mkdir (figdir)
1373                        latex_file_loc = figdir.replace(sep,altsep)   
1374
1375                        if len(label_id) == 1: 
1376                            graphname_latex = '%sgauge%s%s%s' %(latex_file_loc, gaugeloc2, which_quantity, label_id2) # storing files in production directory 
1377                            graphname_report = '%sgauge%s%s%s' %('..'+altsep+'report_figures'+altsep, gaugeloc2, which_quantity, label_id2) # giving location in latex output file
1378                            s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report, '.png')
1379                            fid.write(s)
1380                            if where1 % 2 == 0:
1381                                s = '\\\\ \n'
1382                                where1 = 0
1383                            else:
1384                                s = '& \n'
1385                            fid.write(s)
1386                            savefig(graphname_latex)
1387                   
1388                    if title_on == True:
1389                        title('%s scenario: %s at %s gauge' %(label_id, which_quantity, gaugeloc2))
1390                        #title('Gauge %s (MOST elevation %.2f, ANUGA elevation %.2f)' %(gaugeloc2, elevations[10,k,0], elevations[10,k,1] ))
1391
1392                    savefig(graphname) # save figures with sww file
1393
1394                if report == True and len(label_id) == 1:
1395                    for i in range(nn-1):
1396                        if nn > 2:
1397                            if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1398                                word_quantity += 'depth' + ', '
1399                            else:
1400                                word_quantity += plot_quantity[i] + ', '
1401                        else:
1402                            if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1403                                word_quantity += 'depth' + ', '
1404                            else:
1405                                word_quantity += plot_quantity[i]
1406                       
1407                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
1408                        word_quantity += ' and ' + 'depth'
1409                    else:
1410                        word_quantity += ' and ' + plot_quantity[nn-1]
1411                    caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k]) #gaugeloc.replace('_',' '))
1412                    if elev[k] == 0.0:
1413                        caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j])
1414                        east = gauges[0]
1415                        north = gauges[1]
1416                        elev_output.append([locations[k],east,north,elevations[0,k,j]])
1417                    label = '%sgauge%s' %(label_id2, gaugeloc2)
1418                    s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
1419                    fid.write(s)
1420                    cc += 1
1421                    if cc % 6 == 0: fid.write('\\clearpage \n')
1422                    savefig(graphname_latex)               
1423                   
1424            if report == True and len(label_id) > 1:
1425                for i in range(nn-1):
1426                    if nn > 2:
1427                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1428                            word_quantity += 'depth' + ','
1429                        else:
1430                            word_quantity += plot_quantity[i] + ', '
1431                    else:
1432                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1433                            word_quantity += 'depth'
1434                        else:
1435                            word_quantity += plot_quantity[i]
1436                    where1 = 0
1437                    count1 += 1
1438                    index = j*len(plot_quantity)
1439                    for which_quantity in plot_quantity:
1440                        where1 += 1
1441                        s = '\includegraphics[width=0.49\linewidth, height=50mm]{%s%s}' %(graphname_report[index], '.png')
1442                        index += 1
1443                        fid.write(s)
1444                        if where1 % 2 == 0:
1445                            s = '\\\\ \n'
1446                            where1 = 0
1447                        else:
1448                            s = '& \n'
1449                        fid.write(s)
1450                word_quantity += ' and ' + plot_quantity[nn-1]           
1451                label = 'gauge%s' %(gaugeloc2) 
1452                caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elev[k])
1453                if elev[k] == 0.0:
1454                        caption = 'Time series for %s at %s location (elevation %.2fm)' %(word_quantity, locations[k], elevations[0,k,j])
1455                        thisgauge = gauges[k]
1456                        east = thisgauge[0]
1457                        north = thisgauge[1]
1458                        elev_output.append([locations[k],east,north,elevations[0,k,j]])
1459                       
1460                s = '\end{tabular} \n \\caption{%s} \n \label{fig:%s} \n \end{figure} \n \n' %(caption, label)
1461                fid.write(s)
1462                if float((k+1)/div - pp) == 0.:
1463                    fid.write('\\clearpage \n')
1464                    pp += 1
1465               
1466                #### finished generating figures ###
1467
1468            close('all')
1469       
1470    return texfile2, elev_output
1471
1472# FIXME (DSG): Add unit test, make general, not just 2 files,
1473# but any number of files.
1474def copy_code_files(dir_name, filename1, filename2):
1475    """Temporary Interface to new location"""
1476
1477    from anuga.shallow_water.data_manager import \
1478                    copy_code_files as dm_copy_code_files
1479    print 'copy_code_files has moved from util.py.  ',
1480    print 'Please use "from anuga.shallow_water.data_manager \
1481                                        import copy_code_files"'
1482   
1483    return dm_copy_code_files(dir_name, filename1, filename2)
1484
1485
1486def add_directories(root_directory, directories):
1487    """
1488    Add the first directory in directories to root_directory.
1489    Then add the second
1490    directory to the first directory and so on.
1491
1492    Return the path of the final directory.
1493
1494    This is handy for specifying and creating a directory
1495    where data will go.
1496    """
1497    dir = root_directory
1498    for new_dir in directories:
1499        dir = os.path.join(dir, new_dir)
1500        if not access(dir,F_OK):
1501            mkdir (dir)
1502    return dir
1503
1504def get_data_from_file(filename,separator_value = ','):
1505    """Temporary Interface to new location"""
1506    from anuga.shallow_water.data_manager import \
1507                        get_data_from_file as dm_get_data_from_file
1508    print 'get_data_from_file has moved from util.py'
1509    print 'Please use "from anuga.shallow_water.data_manager \
1510                                     import get_data_from_file"'
1511   
1512    return dm_get_data_from_file(filename,separator_value = ',')
1513
1514def store_parameters(verbose=False,**kwargs):
1515    """Temporary Interface to new location"""
1516   
1517    from anuga.shallow_water.data_manager \
1518                    import store_parameters as dm_store_parameters
1519    print 'store_parameters has moved from util.py.'
1520    print 'Please use "from anuga.shallow_water.data_manager \
1521                                     import store_parameters"'
1522   
1523    return dm_store_parameters(verbose=False,**kwargs)
1524
1525def remove_lone_verts(verts, triangles, number_of_full_nodes=None):
1526    """
1527    Removes vertices that are not associated with any triangles.
1528
1529    verts is a list/array of points
1530    triangles is a list of 3 element tuples. 
1531    Each tuple represents a triangle.
1532
1533    number_of_full_nodes relate to parallelism when a mesh has an
1534        extra layer of ghost points.
1535       
1536    """
1537    verts = ensure_numeric(verts)
1538    triangles = ensure_numeric(triangles)
1539   
1540    N = len(verts)
1541   
1542    # initialise the array to easily find the index of the first loner
1543    loners=arange(2*N, N, -1) # if N=3 [6,5,4]
1544    for t in triangles:
1545        for vert in t:
1546            loners[vert]= vert # all non-loners will have loners[i]=i
1547    #print loners
1548
1549    lone_start = 2*N - max(loners) # The index of the first loner
1550
1551    if lone_start-1 == N:
1552        # no loners
1553        pass
1554    elif min(loners[lone_start:N]) > N:
1555        # All the loners are at the end of the vert array
1556        verts = verts[0:lone_start]
1557    else:
1558        # change the loners list so it can be used to modify triangles
1559        # Remove the loners from verts
1560        # Could've used X=compress(less(loners,N),loners)
1561        # verts=take(verts,X)  to Remove the loners from verts
1562        # but I think it would use more memory
1563        new_i = 0
1564        for i in range(lone_start, N):
1565            if loners[i] >= N:
1566                # Loner!
1567                pass
1568            else:
1569                loners[i] = new_i
1570                verts[new_i] = verts[i]
1571                new_i += 1
1572        verts = verts[0:new_i]
1573
1574        # Modify the triangles
1575        #print "loners", loners
1576        #print "triangles before", triangles
1577        triangles = choose(triangles,loners)
1578        #print "triangles after", triangles
1579    return verts, triangles
1580
1581 
1582def get_centroid_values(x, triangles):
1583    """Compute centroid values from vertex values
1584   
1585    x: Values at vertices of triangular mesh
1586    triangles: Nx3 integer array pointing to vertex information
1587    for each of the N triangels. Elements of triangles are
1588    indices into x
1589    """
1590
1591       
1592    xc = zeros(triangles.shape[0], Float) # Space for centroid info
1593   
1594    for k in range(triangles.shape[0]):
1595        # Indices of vertices
1596        i0 = triangles[k][0]
1597        i1 = triangles[k][1]
1598        i2 = triangles[k][2]       
1599       
1600        xc[k] = (x[i0] + x[i1] + x[i2])/3
1601
1602
1603    return xc
1604
1605def make_plots_from_csv_file(directories_dic={},
1606                            output_dir='',
1607                            base_name=None,
1608                            plot_numbers='0',
1609                            quantities=['Stage'],
1610                            extra_plot_name='',
1611                            assess_all_csv_files=True,                           
1612                            create_latex=False,
1613                            verbose=False):
1614                               
1615    """Read in csv files that have the right header information and
1616    plot time series such as Stage, Speed, etc. Will also plot several
1617    time series on one plot. Filenames must follow this convention,
1618    <base_name><plot_number>.csv eg gauge_timeseries3.csv
1619   
1620    For example:   
1621    if "directories_dic" defines 4 directories and in each directories
1622    there is a csv files corresponding to the right "plot_numbers",
1623    this will create a plot with 4 lines one for each directory AND
1624    one plot for each "quantities". 
1625   
1626    Inputs:
1627        directories_dic: dictionary of directory with values (plot
1628                         legend name for directory), (start time of
1629                         the time series) and the (value to add to
1630                         stage if needed). For example
1631                        {dir1:['Anuga_ons',5000, 0],
1632                         dir2:['b_emoth',5000,1.5],
1633                         dir3:['b_ons',5000,1.5]}
1634                         
1635        output_dir: directory for the plot outputs
1636       
1637        base_name: common name to all the csv files to be read
1638       
1639        plot_numbers: a String list of numbers to plot. For example
1640                      [0-4,10,15-17] will read and attempt to plot
1641                       the follow 0,1,2,3,4,10,15,16,17
1642        quantities: Currently limited to "Stage", "Speed", and
1643                    "Momentum", should be changed to incorporate
1644                    any quantity read from CSV file....
1645                   
1646        extra_plot_name: A string that is appended to the end of the
1647                         output filename.
1648                   
1649        assess_all_csv_files: if true it will read ALL csv file with
1650                             "base_name", regardless of 'plot_numbers'
1651                              and determine a uniform set of axes for
1652                              Stage, Speed and Momentum. IF FALSE it
1653                              will only read the csv file within the
1654                             'plot_numbers'
1655                             
1656        create_latex: NOT IMPLEMENTED YET!! sorry Jane....
1657       
1658        OUTPUTS: None, it saves the plots to
1659              <output_dir><base_name><plot_number><extra_plot_name>.png
1660   
1661    """
1662    import pylab# import plot, xlabel, ylabel, savefig, \
1663                #ion, hold, axis, close, figure, legend
1664    from os import sep
1665    from anuga.shallow_water.data_manager import \
1666                               get_all_files_with_extension, csv2dict
1667    #find all the files that meet the specs
1668    #FIXME remove all references to word that such as Stage
1669    #(with a Capital letter) could use the header info from
1670    #the dictionary that is returned from the csv2dict this could
1671    #be very good and very flexable.... but this works for now!
1672
1673    #FIXME plot all available data from the csv file, currently will
1674    #only plot Stage,Speed and Momentum and only if the header is
1675    #named like wise
1676
1677#    quantities_dic={'time':'time (hour)'}
1678#    if 'time' in quantities:
1679#        quantities_dic['time']='time (hour)'
1680#    if 'Time' in quantities:
1681#        quantities_dic['Time']='time (hour)'
1682#    if 'stage' in quantities:
1683#        quantities_dic['stage']='wave height (m)'
1684#    if 'Stage' in quantities:
1685#        quantities_dic['Stage']='wave height (m)'
1686#    if 'speed' in quantities:
1687#        quantities_dic['speed']='speed (m/s)'
1688#    if 'Speed' in quantities:
1689#        quantities_dic['Speed']='speed (m/s)'
1690#    if 'stage' or 'Stage' in quantities:
1691#        quantities_dic['stage']='speed (m/s)'
1692#    if 'stage' or 'Stage' in quantities:
1693#        quantities_dic['stage']='speed (m/s)'
1694
1695    seconds_in_hour = 3600
1696    time_label = 'time (hour)'
1697    stage_label = 'wave height (m)'
1698    speed_label = 'speed (m/s)'
1699    momentum_label = 'momentum (m^2/sec)'
1700   
1701    if extra_plot_name != '':
1702        extra_plot_name='_'+extra_plot_name
1703
1704    #finds all the files that fit the specs and return a list of them
1705    #so to help find a uniform max and min for the plots...
1706    list_filenames=[]
1707    if verbose: print 'Determining files to access for axes ranges \n'
1708    for i,directory in enumerate(directories_dic.keys()):
1709        list_filenames.append(get_all_files_with_extension(directory,
1710                              base_name,'.csv'))
1711#    print 'list_filenames',list_filenames
1712
1713    #use all the files to get the values for the plot axis
1714    max_st=max_sp=max_mom=min_st=min_sp=min_mom=max_t=min_t=0.
1715    max_start_time= 0.
1716    min_start_time = 100000 
1717
1718   
1719    new_plot_numbers=[]
1720    #change plot_numbers to list, eg ['0-4','10']
1721    #to ['0','1','2','3','4','10']
1722    for i, num_string in enumerate(plot_numbers):
1723        if '-' in num_string: 
1724            start = int(num_string[:num_string.rfind('-')])
1725            end = int(num_string[num_string.rfind('-')+1:])+1
1726            for x in range(start, end):
1727                new_plot_numbers.append(str(x))
1728        else:
1729            new_plot_numbers.append(num_string)
1730#    print 'new_plot_numbers',new_plot_numbers
1731   
1732    if verbose: print 'Determining uniform axes \n' 
1733    #this entire loop is to determine the min and max range for the
1734    #axes of the plot
1735    for i, directory in enumerate(directories_dic.keys()):
1736       
1737        if assess_all_csv_files==False:
1738            which_csv_to_assess = new_plot_numbers
1739        else:
1740            which_csv_to_assess = list_filenames[i]
1741
1742        for j, filename in enumerate(which_csv_to_assess):
1743            if assess_all_csv_files==False:
1744                dir_filename=directory+base_name+sep+filename
1745            else:
1746                dir_filename=join(directory,filename)
1747               
1748            attribute_dic, title_index_dic = csv2dict(dir_filename+
1749                                                       '.csv')
1750
1751            directory_start_time = directories_dic[directory][1]
1752            directory_add_tide = directories_dic[directory][2]
1753
1754            time = [float(x) for x in attribute_dic["Time"]]
1755            min_t, max_t = get_min_max_values(time,min_t,max_t)
1756           
1757            stage = [float(x) for x in attribute_dic["Stage"]]
1758            stage =array(stage)+directory_add_tide
1759            min_st, max_st = get_min_max_values(stage,min_st,max_st)
1760           
1761            speed = [float(x) for x in attribute_dic["Speed"]]
1762            min_sp, max_sp = get_min_max_values(speed,min_sp,max_sp)
1763           
1764            momentum = [float(x) for x in attribute_dic["Momentum"]]
1765            min_mom, max_mom = get_min_max_values(momentum,
1766                                                  min_mom,
1767                                                  max_mom)
1768                                                                     
1769#            print 'min_sp, max_sp',min_sp, max_sp,
1770            # print directory_start_time
1771            if min_start_time > directory_start_time: 
1772                min_start_time = directory_start_time
1773            if max_start_time < directory_start_time: 
1774                max_start_time = directory_start_time
1775#            print 'start_time' , max_start_time, min_start_time
1776   
1777    stage_axis = (min_start_time/seconds_in_hour,
1778                 (max_t+max_start_time)/seconds_in_hour,
1779                  min_st, max_st)
1780    speed_axis = (min_start_time/seconds_in_hour,
1781                 (max_t+max_start_time)/seconds_in_hour,
1782                 min_sp, max_sp)
1783    momentum_axis = (min_start_time/seconds_in_hour,
1784                    (max_t+max_start_time)/seconds_in_hour,
1785                     min_mom, max_mom)
1786#    ion()
1787    pylab.hold()
1788   
1789   
1790    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
1791   
1792#    if verbose: print 'Now start to plot \n'
1793   
1794    i_max = len(directories_dic.keys())
1795    legend_list_dic =[]
1796    legend_list =[]
1797    for i, directory in enumerate(directories_dic.keys()): 
1798        if verbose: print'Plotting in %s', directory
1799        for j, number in enumerate(new_plot_numbers):
1800            if verbose: print'Starting %s',base_name+number 
1801            directory_name = directories_dic[directory][0]
1802            directory_start_time = directories_dic[directory][1]
1803            directory_add_tide = directories_dic[directory][2]
1804           
1805            #create an if about the start time and tide hieght
1806            #if they don't exist
1807            file=directory+sep+base_name+number
1808            #print 'i %s,j %s, number %s, file %s' %(i,j,number,file)
1809            attribute_dic, title_index_dic = csv2dict(file+'.csv')
1810            #get data from dict in to list
1811            t = [float(x) for x in attribute_dic["Time"]]
1812
1813            #do maths to list by changing to array
1814            t=(array(t)+directory_start_time)/seconds_in_hour
1815            stage = [float(x) for x in attribute_dic["Stage"]]
1816            speed = [float(x) for x in attribute_dic["Speed"]]
1817            momentum = [float(x) for x in attribute_dic["Momentum"]]
1818                       
1819            #Can add tide so plots are all on the same tide height
1820            stage =array(stage)+directory_add_tide
1821           
1822            #finds the maximum elevation
1823            max_ele=-100000
1824            min_ele=100000
1825            elevation = [float(x) for x in attribute_dic["Elevation"]]
1826           
1827            min_ele, max_ele = get_min_max_values(elevation,
1828                                                  min_ele,
1829                                                  max_ele)
1830            if min_ele != max_ele:
1831                print "Note! Elevation changes in %s" %dir_filename
1832#            print 'min_ele, max_ele',min_ele, max_ele
1833           
1834
1835            #populates the legend_list_dic with dir_name and the elevation
1836            if i==0:
1837                legend_list_dic.append({directory_name:max_ele})
1838            else:
1839                legend_list_dic[j][directory_name]=max_ele
1840           
1841            # creates a list for the legend after at "legend_dic" has been fully populated
1842            # only runs on the last iteration for all the gauges(csv) files
1843            # empties the list before creating it
1844            if i==i_max-1:
1845                legend_list=[]
1846                for k, l in legend_list_dic[j].iteritems():
1847                    legend_list.append('%s (elevation = %sm)'%(k,l))
1848                    #print k,l, legend_list_dic[j]
1849         
1850            if "Stage" in quantities:
1851                pylab.figure(100+j)
1852                pylab.plot(t, stage, c = cstr[i], linewidth=1)
1853                pylab.xlabel(time_label)
1854                pylab.ylabel(stage_label)
1855                pylab.axis(stage_axis)
1856                pylab.legend(legend_list,loc='upper right')
1857                figname = '%sstage_%s%s.png' %(output_dir+sep,
1858                                               base_name+number,
1859                                               extra_plot_name)
1860                pylab.savefig(figname)
1861            if "Speed" in quantities:
1862                pylab.figure(200+j)
1863                pylab.plot(t, speed, c = cstr[i], linewidth=1)
1864                pylab.xlabel(time_label)
1865                pylab.ylabel(speed_label)
1866                pylab.axis(speed_axis)
1867                pylab.legend(legend_list,loc='upper right')
1868                figname = '%sspeed_%s%s.png' %(output_dir+sep,
1869                                               base_name+number,
1870                                               extra_plot_name)
1871                pylab.savefig(figname)
1872            if "Momentum" in quantities:
1873                pylab.figure(300+j)
1874                pylab.plot(t, momentum, c = cstr[i], linewidth=1)
1875                pylab.xlabel(time_label)
1876                pylab.ylabel(momentum_label)
1877                pylab.axis(momentum_axis)
1878                pylab.legend(legend_list,loc='upper right')
1879                figname = '%smomentum_%s%s.png' %(output_dir+sep,
1880                                                  base_name+number,
1881                                                  extra_plot_name)
1882                pylab.savefig(figname)
1883    if verbose: print 'Closing all plots'
1884    pylab.close('all')
1885    del pylab
1886    if verbose: print 'Finished closing plots'
1887
1888def get_min_max_values(list=None,min1=100,max1=-100):
1889    """ Returns the min and max of the list it was provided.
1890    NOTE: default min and max may need to change depeending on
1891    your list
1892    """
1893    if list == None: print 'List must be provided'
1894#    min = max_list = 0
1895    if max(list) > max1: 
1896        max1 = max(list)
1897    if min(list) < min1: 
1898        min1 = min(list)
1899       
1900    return min1, max1
1901
1902
1903def get_runup_data_for_locations_from_file(gauge_filename,
1904                                           sww_filename,
1905                                           runup_filename,
1906                                           size=10,
1907                                           verbose=False):
1908
1909    '''this will read a csv file with the header x,y. Then look in a square 'size'x2
1910    around this position for the 'max_inundaiton_height' in the 'sww_filename' and
1911    report the findings in the 'runup_filename
1912   
1913    WARNING: NO TESTS!
1914    '''
1915
1916    from anuga.shallow_water.data_manager import get_all_directories_with_name,\
1917                                                 get_maximum_inundation_data,\
1918                                                 csv2dict
1919                                                 
1920    file = open(runup_filename,"w")
1921    file.write("easting,northing,runup \n ")
1922    file.close()
1923   
1924    #read gauge csv file to dictionary
1925    attribute_dic, title_index_dic = csv2dict(gauge_filename)
1926    northing = [float(x) for x in attribute_dic["y"]]
1927    easting = [float(x) for x in attribute_dic["x"]]
1928
1929    print 'Reading %s' %sww_filename
1930
1931    runup_locations=[]
1932    for i, x in enumerate(northing):
1933#        print 'easting,northing',i,easting[i],northing[i]
1934        poly = [[int(easting[i]+size),int(northing[i]+size)],
1935                [int(easting[i]+size),int(northing[i]-size)],
1936                [int(easting[i]-size),int(northing[i]-size)],
1937                [int(easting[i]-size),int(northing[i]+size)]]
1938       
1939        run_up, x_y = get_maximum_inundation_data(filename=sww_filename,
1940                                              polygon=poly,
1941                                              verbose=False) 
1942        #if no runup will return 0 instead of NONE
1943        if run_up==None: run_up=0
1944        if x_y==None: x_y=[0,0]
1945       
1946        if verbose:print 'maximum inundation runup near %s is %s meters' %(x_y,run_up)
1947       
1948        #writes to file
1949        file = open(runup_filename,"a")
1950        temp = '%s,%s,%s \n' %(x_y[0], x_y[1], run_up)
1951        file.write(temp)
1952        file.close()
Note: See TracBrowser for help on using the repository browser.