source: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py @ 7780

Last change on this file since 7780 was 7780, checked in by hudson, 13 years ago

Almost all failing tests fixed.
File size: 46.5 KB
Line 
1"""Gauge functions
2   
3   High-level functions for converting gauge and sww files into timeseries plots.
4
5
6   Copyright 2010
7   Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou, James Hudson
8   Geoscience Australia
9"""
10
11import numpy as num
12
13from anuga.geospatial_data.geospatial_data import ensure_absolute
14from util import check_list, calc_bearing
15from file_function import file_function
16
17import os
18
19from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep, getcwd
20from os.path import exists, split, join
21import anuga.utilities.log as log
22
23from math import sqrt
24
25def _quantities2csv(quantities, point_quantities, centroids, point_i):
26    points_list = []
27   
28    for quantity in quantities:
29        #core quantities that are exported from the interpolator     
30        if quantity == 'stage':
31            points_list.append(point_quantities[0])
32           
33        if quantity == 'elevation':
34            points_list.append(point_quantities[1])
35           
36        if quantity == 'xmomentum':
37            points_list.append(point_quantities[2])
38           
39        if quantity == 'ymomentum':
40            points_list.append(point_quantities[3])
41
42        #derived quantities that are calculated from the core ones
43        if quantity == 'depth':
44            points_list.append(point_quantities[0] 
45                               - point_quantities[1])
46
47        if quantity == 'momentum':
48            momentum = sqrt(point_quantities[2]**2 
49                            + point_quantities[3]**2)
50            points_list.append(momentum)
51           
52        if quantity == 'speed':
53            #if depth is less than 0.001 then speed = 0.0
54            if point_quantities[0] - point_quantities[1] < 0.001:
55                vel = 0.0
56            else:
57                if point_quantities[2] < 1.0e6:
58                    momentum = sqrt(point_quantities[2]**2
59                                    + point_quantities[3]**2)
60                    vel = momentum / (point_quantities[0] 
61                                      - point_quantities[1])
62                else:
63                    momentum = 0
64                    vel = 0
65               
66            points_list.append(vel)
67           
68        if quantity == 'bearing':
69            points_list.append(calc_bearing(point_quantities[2],
70                                            point_quantities[3]))
71        if quantity == 'xcentroid':
72            points_list.append(centroids[point_i][0])
73
74        if quantity == 'ycentroid':
75            points_list.append(centroids[point_i][1])
76
77    return points_list
78   
79   
80##
81# @brief Take gauge readings, given a gauge file and a sww mesh
82#
83#        Use this function to take a timeseries sample, given a list of gauge points
84# @param  sww_file sww file to use as input
85# @param gauge_file gauge file as input, containing list of gauge points to sample
86# @param out_name output file prefix
87# @param quantities which quantities in the sww file we want to export
88# @param verbose show extra logging information for debug purposes
89# @param use_cache cache requests if possible, for speed
90# @param output_centroids Set to true to output the values at the centroid of the mesh triangle
91def sww2csv_gauges(sww_file,
92                   gauge_file,
93                   out_name='gauge_',
94                   quantities=['stage', 'depth', 'elevation',
95                               'xmomentum', 'ymomentum'],
96                   verbose=False,
97                   use_cache=True,
98                   output_centroids=False):
99    """
100   
101    Inputs:
102        NOTE: if using csv2timeseries_graphs after creating csv file,
103        it is essential to export quantities 'depth' and 'elevation'.
104        'depth' is good to analyse gauges on land and elevation is used
105        automatically by csv2timeseries_graphs in the legend.
106       
107        sww_file: path to any sww file
108       
109        gauge_file: Assumes that it follows this format
110            name, easting, northing, elevation
111            point1, 100.3, 50.2, 10.0
112            point2, 10.3, 70.3, 78.0
113       
114        NOTE: order of column can change but names eg 'easting', 'elevation'
115        must be the same! ALL lowercaps!
116
117        out_name: prefix for output file name (default is 'gauge_')
118       
119    Outputs:
120        one file for each gauge/point location in the points file. They
121        will be named with this format in the same directory as the 'sww_file'
122            <out_name><name>.csv
123        eg gauge_point1.csv if <out_name> not supplied
124           myfile_2_point1.csv if <out_name> ='myfile_2_'
125           
126        They will all have a header
127   
128    Usage: sww2csv_gauges(sww_file='test1.sww',
129                          quantities = ['stage', 'elevation','depth','bearing'],
130                          gauge_file='gauge.txt')   
131   
132    Interpolate the quantities at a given set of locations, given
133    an sww file.
134    The results are written to a csv file.
135
136    In the future let points be a points file.
137    And the user choose the quantities.
138
139    This is currently quite specific.
140    If it needs to be more general, change things.
141
142    This is really returning speed, not velocity.
143    """
144   
145    from csv import reader,writer
146    from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
147    import string
148    from anuga.utilities.file_utils import get_all_swwfiles
149    from anuga.abstract_2d_finite_volumes.util import file_function   
150
151    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
152    assert type(out_name) == type(''), 'Output filename prefix must be a string'
153   
154    try:
155        point_reader = reader(file(gauge_file))
156    except Exception, e:
157        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e)
158        raise msg
159
160    if verbose: log.critical('Gauges obtained from: %s' % gauge_file)
161   
162    point_reader = reader(file(gauge_file))
163    points = []
164    point_name = []
165   
166    # read point info from file
167    for i,row in enumerate(point_reader):
168        # read header and determine the column numbers to read correctly.
169        if i==0:
170            for j,value in enumerate(row):
171                if value.strip()=='easting':easting=j
172                if value.strip()=='northing':northing=j
173                if value.strip()=='name':name=j
174                if value.strip()=='elevation':elevation=j
175        else:
176            points.append([float(row[easting]),float(row[northing])])
177            point_name.append(row[name])
178       
179    #convert to array for file_function
180    points_array = num.array(points,num.float)
181       
182    points_array = ensure_absolute(points_array)
183
184    dir_name, base = os.path.split(sww_file)   
185
186    #need to get current directory so when path and file
187    #are "joined" below the directory is correct
188    if dir_name == '':
189        dir_name =getcwd()
190       
191    if access(sww_file,R_OK):
192        if verbose: log.critical('File %s exists' % sww_file)
193    else:
194        msg = 'File "%s" could not be opened: no read permission' % sww_file
195        raise msg
196
197    sww_files = get_all_swwfiles(look_in_dir=dir_name,
198                                 base_name=base,
199                                 verbose=verbose)
200
201    # fudge to get SWW files in 'correct' order, oldest on the left
202    sww_files.sort()
203
204    if verbose:
205        log.critical('sww files=%s' % sww_files)
206   
207    #to make all the quantities lower case for file_function
208    quantities = [quantity.lower() for quantity in quantities]
209
210    # what is quantities are needed from sww file to calculate output quantities
211    # also
212
213    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
214
215    gauge_file = out_name
216
217    heading = [quantity for quantity in quantities]
218    heading.insert(0,'time')
219    heading.insert(1,'hours')
220   
221    if verbose: log.critical('Writing csv files')
222
223    quake_offset_time = None
224
225    for sww_file in sww_files:
226        sww_file = join(dir_name, sww_file+'.sww')
227        callable_sww = file_function(sww_file,
228                                     quantities=core_quantities,
229                                     interpolation_points=points_array,
230                                     verbose=verbose,
231                                     use_cache=use_cache,
232                                     output_centroids = output_centroids)
233
234        if quake_offset_time is None:
235            quake_offset_time = callable_sww.starttime
236
237    for point_i, point in enumerate(points_array):
238        is_opened = False   
239        for time in callable_sww.get_time():
240            #add domain starttime to relative time.
241            quake_time = time + quake_offset_time
242            point_quantities = callable_sww(time, point_i) # __call__ is overridden
243
244            if point_quantities[0] != NAN:
245                if is_opened == False:
246                    points_writer = writer(file(dir_name + sep + gauge_file
247                                                        + point_name[point_i] + '.csv', "wb"))
248                    points_writer.writerow(heading)
249                    is_opened = True
250                points_list = [quake_time, quake_time/3600.] +  _quantities2csv(quantities, point_quantities, callable_sww.centroids, point_i)
251                points_writer.writerow(points_list)
252            else:
253                if verbose:
254                    msg = 'gauge' + point_name[point_i] + 'falls off the mesh in file ' + sww_file + '.'
255                    log.warning(msg)
256##
257# @brief Read a .sww file and plot the time series.
258# @param swwfiles Dictionary of .sww files.
259# @param gauge_filename Name of gauge file.
260# @param production_dirs ??
261# @param report If True, write figures to report directory.
262# @param reportname Name of generated report (if any).
263# @param plot_quantity List containing quantities to plot.
264# @param generate_fig If True, generate figures as well as CSV files.
265# @param surface If True, then generate solution surface with 3d plot.
266# @param time_min Beginning of user defined time range for plotting purposes.
267# @param time_max End of user defined time range for plotting purposes.
268# @param time_thinning ??
269# @param time_unit ??
270# @param title_on If True, export standard graphics with title.
271# @param use_cache If True, use caching.
272# @param verbose If True, this function is verbose.
273def sww2timeseries(swwfiles,
274                   gauge_filename,
275                   production_dirs,
276                   report=None,
277                   reportname=None,
278                   plot_quantity=None,
279                   generate_fig=False,
280                   surface=None,
281                   time_min=None,
282                   time_max=None,
283                   time_thinning=1,                   
284                   time_unit=None,
285                   title_on=None,
286                   use_cache=False,
287                   verbose=False,
288                   output_centroids=False):
289    """ Read sww file and plot the time series for the
290    prescribed quantities at defined gauge locations and
291    prescribed time range.
292
293    Input variables:
294
295    swwfiles        - dictionary of sww files with label_ids (used in
296                      generating latex output. It will be part of
297                      the directory name of file_loc (typically the timestamp).
298                      Helps to differentiate latex files for different
299                      simulations for a particular scenario. 
300                    - assume that all conserved quantities have been stored
301                    - assume each sww file has been simulated with same timestep
302   
303    gauge_filename  - name of file containing gauge data
304                        - easting, northing, name , elevation?
305                    - OR (this is not yet done)
306                        - structure which can be converted to a numeric array,
307                          such as a geospatial data object
308                     
309    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
310                                                'boundaries': 'urs boundary'}
311                      this will use the second part as the label and the
312                      first part as the ?
313                      #FIXME: Is it a list or a dictionary
314                      # This is probably obsolete by now
315                     
316    report          - if True, then write figures to report_figures directory in
317                      relevant production directory
318                    - if False, figures are already stored with sww file
319                    - default to False
320
321    reportname      - name for report if wishing to generate report
322   
323    plot_quantity   - list containing quantities to plot, they must
324                      be the name of an existing quantity or one of
325                      the following possibilities
326                    - possibilities:
327                        - stage; 'stage'
328                        - depth; 'depth'
329                        - speed; calculated as absolute momentum
330                         (pointwise) divided by depth; 'speed'
331                        - bearing; calculated as the angle of the momentum
332                          vector (xmomentum, ymomentum) from the North; 'bearing'
333                        - absolute momentum; calculated as
334                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
335                        - x momentum; 'xmomentum'
336                        - y momentum; 'ymomentum'
337                    - default will be ['stage', 'speed', 'bearing']
338
339    generate_fig     - if True, generate figures as well as csv file
340                     - if False, csv files created only
341                     
342    surface          - if True, then generate solution surface with 3d plot
343                       and save to current working directory
344                     - default = False
345   
346    time_min         - beginning of user defined time range for plotting purposes
347                        - default will be first available time found in swwfile
348                       
349    time_max         - end of user defined time range for plotting purposes
350                        - default will be last available time found in swwfile
351                       
352    title_on        - if True, export standard graphics with title
353                    - if False, export standard graphics without title
354
355
356    Output:
357   
358    - time series data stored in .csv for later use if required.
359      Name = gauges_timeseries followed by gauge name
360    - latex file will be generated in same directory as where script is
361      run (usually production scenario directory.
362      Name = latexoutputlabel_id.tex
363
364    Other important information:
365   
366    It is assumed that the used has stored all the conserved quantities
367    and elevation during the scenario run, i.e.
368    ['stage', 'elevation', 'xmomentum', 'ymomentum']
369    If this has not occurred then sww2timeseries will not work.
370
371
372    Usage example
373    texname = sww2timeseries({project.boundary_name + '.sww': ''},
374                             project.polygons_dir + sep + 'boundary_extent.csv',
375                             project.anuga_dir,
376                             report = False,
377                             plot_quantity = ['stage', 'speed', 'bearing'],
378                             time_min = None,
379                             time_max = None,
380                             title_on = True,   
381                             verbose = True)
382   
383    """
384
385    msg = 'NOTE: A new function is available to create csv files from sww '
386    msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
387    msg += ' PLUS another new function to create graphs from csv files called '
388    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
389    log.critical(msg)
390   
391    k = _sww2timeseries(swwfiles,
392                        gauge_filename,
393                        production_dirs,
394                        report,
395                        reportname,
396                        plot_quantity,
397                        generate_fig,
398                        surface,
399                        time_min,
400                        time_max,
401                        time_thinning,                       
402                        time_unit,
403                        title_on,
404                        use_cache,
405                        verbose,
406                        output_centroids = output_centroids)
407    return k
408
409
410##
411# @brief Read a .sww file and plot the time series.
412# @param swwfiles Dictionary of .sww files.
413# @param gauge_filename Name of gauge file.
414# @param production_dirs ??
415# @param report If True, write figures to report directory.
416# @param reportname Name of generated report (if any).
417# @param plot_quantity List containing quantities to plot.
418# @param generate_fig If True, generate figures as well as CSV files.
419# @param surface If True, then generate solution surface with 3d plot.
420# @param time_min Beginning of user defined time range for plotting purposes.
421# @param time_max End of user defined time range for plotting purposes.
422# @param time_thinning ??
423# @param time_unit ??
424# @param title_on If True, export standard graphics with title.
425# @param use_cache If True, use caching.
426# @param verbose If True, this function is verbose.
427def _sww2timeseries(swwfiles,
428                    gauge_filename,
429                    production_dirs,
430                    report = None,
431                    reportname = None,
432                    plot_quantity = None,
433                    generate_fig = False,
434                    surface = None,
435                    time_min = None,
436                    time_max = None,
437                    time_thinning = 1,                   
438                    time_unit = None,
439                    title_on = None,
440                    use_cache = False,
441                    verbose = False,
442                    output_centroids = False):   
443       
444    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
445    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
446   
447    try:
448        fid = open(gauge_filename)
449    except Exception, e:
450        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
451        raise msg
452
453    if report is None:
454        report = False
455       
456    if plot_quantity is None:
457        plot_quantity = ['depth', 'speed']
458    else:
459        assert type(plot_quantity) == list, 'plot_quantity must be a list'
460        check_list(plot_quantity)
461
462    if surface is None:
463        surface = False
464
465    if time_unit is None:
466        time_unit = 'hours'
467   
468    if title_on is None:
469        title_on = True
470   
471    if verbose: log.critical('Gauges obtained from: %s' % gauge_filename)
472
473    gauges, locations, elev = gauge_get_from_file(gauge_filename)
474
475    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
476
477    file_loc = []
478    f_list = []
479    label_id = []
480    leg_label = []
481    themaxT = 0.0
482    theminT = 0.0
483
484    for swwfile in swwfiles.keys():
485        try:
486            fid = open(swwfile)
487        except Exception, e:
488            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
489            raise msg
490
491        if verbose:
492            log.critical('swwfile = %s' % swwfile)
493
494        # Extract parent dir name and use as label
495        path, _ = os.path.split(swwfile)
496        _, label = os.path.split(path)       
497       
498        leg_label.append(label)
499
500        f = file_function(swwfile,
501                          quantities = sww_quantity,
502                          interpolation_points = gauges,
503                          time_thinning = time_thinning,
504                          verbose = verbose,
505                          use_cache = use_cache,
506                          output_centroids = output_centroids)
507
508        # determine which gauges are contained in sww file
509        count = 0
510        gauge_index = []
511        for k, g in enumerate(gauges):
512            if f(0.0, point_id = k)[2] > 1.0e6:
513                count += 1
514                if count == 1: log.critical('Gauges not contained here:')
515                log.critical(locations[k])
516            else:
517                gauge_index.append(k)
518
519        if len(gauge_index) > 0:
520            log.critical('Gauges contained here:')
521        else:
522            log.critical('No gauges contained here.')
523        for i in range(len(gauge_index)):
524             log.critical(locations[gauge_index[i]])
525             
526        index = swwfile.rfind(sep)
527        file_loc.append(swwfile[:index+1])
528        label_id.append(swwfiles[swwfile])
529       
530        f_list.append(f)
531        maxT = max(f.get_time())
532        minT = min(f.get_time())
533        if maxT > themaxT: themaxT = maxT
534        if minT > theminT: theminT = minT
535
536    if time_min is None:
537        time_min = theminT # min(T)
538    else:
539        if time_min < theminT: # min(T):
540            msg = 'Minimum time entered not correct - please try again'
541            raise Exception, msg
542
543    if time_max is None:
544        time_max = themaxT # max(T)
545    else:
546        if time_max > themaxT: # max(T):
547            msg = 'Maximum time entered not correct - please try again'
548            raise Exception, msg
549
550    if verbose and len(gauge_index) > 0:
551         log.critical('Inputs OK - going to generate figures')
552
553    if len(gauge_index) <> 0:
554        texfile, elev_output = \
555            _generate_figures(plot_quantity, file_loc, report, reportname,
556                             surface, leg_label, f_list, gauges, locations,
557                             elev, gauge_index, production_dirs, time_min,
558                             time_max, time_unit, title_on, label_id,
559                             generate_fig, verbose)
560    else:
561        texfile = ''
562        elev_output = []
563
564    return texfile, elev_output
565
566
567##
568# @brief Read gauge info from a file.
569# @param filename The name of the file to read.
570# @return A (gauges, gaugelocation, elev) tuple.
571def gauge_get_from_file(filename):
572    """ Read in gauge information from file
573    """
574
575    from os import sep, getcwd, access, F_OK, mkdir
576
577    # Get data from the gauge file
578    fid = open(filename)
579    lines = fid.readlines()
580    fid.close()
581   
582    gauges = []
583    gaugelocation = []
584    elev = []
585
586    # Check header information   
587    line1 = lines[0]
588    line11 = line1.split(',')
589
590    if isinstance(line11[0], str) is True:
591        # We have found text in the first line
592        east_index = None
593        north_index = None
594        name_index = None
595        elev_index = None
596
597        for i in range(len(line11)):
598            if line11[i].strip().lower() == 'easting':   east_index = i
599            if line11[i].strip().lower() == 'northing':  north_index = i
600            if line11[i].strip().lower() == 'name':      name_index = i
601            if line11[i].strip().lower() == 'elevation': elev_index = i
602
603        if east_index < len(line11) and north_index < len(line11):
604            pass
605        else:
606            msg = 'WARNING: %s does not contain correct header information' \
607                  % filename
608            msg += 'The header must be: easting, northing, name, elevation'
609            raise Exception, msg
610
611        if elev_index is None: 
612            raise Exception
613   
614        if name_index is None: 
615            raise Exception
616
617        lines = lines[1:] # Remove header from data
618    else:
619        # No header, assume that this is a simple easting, northing file
620
621        msg = 'There was no header in file %s and the number of columns is %d' \
622              % (filename, len(line11))
623        msg += '- was assuming two columns corresponding to Easting and Northing'
624        assert len(line11) == 2, msg
625
626        east_index = 0
627        north_index = 1
628
629        N = len(lines)
630        elev = [-9999]*N
631        gaugelocation = range(N)
632       
633    # Read in gauge data
634    for line in lines:
635        fields = line.split(',')
636
637        gauges.append([float(fields[east_index]), float(fields[north_index])])
638
639        if len(fields) > 2:
640            elev.append(float(fields[elev_index]))
641            loc = fields[name_index]
642            gaugelocation.append(loc.strip('\n'))
643
644    return gauges, gaugelocation, elev
645
646   
647##
648# @brief Generate figures from quantities and gauges for each sww file.
649#        This is a legacy function, used only by sww2timeseries
650# @param plot_quantity  ??
651# @param file_loc ??
652# @param report ??
653# @param reportname ??
654# @param surface ??
655# @param leg_label ??
656# @param f_list ??
657# @param gauges ??
658# @param locations ??
659# @param elev ??
660# @param gauge_index ??
661# @param production_dirs ??
662# @param time_min ??
663# @param time_max ??
664# @param time_unit ??
665# @param title_on ??
666# @param label_id ??
667# @param generate_fig ??
668# @param verbose??
669# @return (texfile2, elev_output)
670def _generate_figures(plot_quantity, file_loc, report, reportname, surface,
671                     leg_label, f_list, gauges, locations, elev, gauge_index,
672                     production_dirs, time_min, time_max, time_unit,
673                     title_on, label_id, generate_fig, verbose):
674    """ Generate figures based on required quantities and gauges for
675    each sww file
676    """
677    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
678
679    if generate_fig is True:
680        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
681             xlabel, ylabel, title, close, subplot
682   
683        if surface is True:
684            import pylab as p1
685            import mpl3d.mplot3d as p3
686       
687    if report == True:   
688        texdir = getcwd()+sep+'report'+sep
689        if access(texdir,F_OK) == 0:
690            mkdir (texdir)
691        if len(label_id) == 1:
692            label_id1 = label_id[0].replace(sep,'')
693            label_id2 = label_id1.replace('_','')
694            texfile = texdir + reportname + '%s' % label_id2
695            texfile2 = reportname + '%s' % label_id2
696            texfilename = texfile + '.tex'
697            fid = open(texfilename, 'w')
698
699            if verbose: log.critical('Latex output printed to %s' % texfilename)
700        else:
701            texfile = texdir+reportname
702            texfile2 = reportname
703            texfilename = texfile + '.tex' 
704            fid = open(texfilename, 'w')
705
706            if verbose: log.critical('Latex output printed to %s' % texfilename)
707    else:
708        texfile = ''
709        texfile2 = ''
710
711    p = len(f_list)
712    n = []
713    n0 = 0
714    for i in range(len(f_list)):
715        n.append(len(f_list[i].get_time()))
716        if n[i] > n0: n0 = n[i] 
717    n0 = int(n0)
718    m = len(locations)
719    model_time = num.zeros((n0, m, p), num.float) 
720    stages = num.zeros((n0, m, p), num.float)
721    elevations = num.zeros((n0, m, p), num.float) 
722    momenta = num.zeros((n0, m, p), num.float)
723    xmom = num.zeros((n0, m, p), num.float)
724    ymom = num.zeros((n0, m, p), num.float)
725    speed = num.zeros((n0, m, p), num.float)
726    bearings = num.zeros((n0, m, p), num.float)
727    due_east = 90.0*num.ones((n0, 1), num.float)
728    due_west = 270.0*num.ones((n0, 1), num.float)
729    depths = num.zeros((n0, m, p), num.float)
730    eastings = num.zeros((n0, m, p), num.float)
731    min_stages = []
732    max_stages = []
733    min_momentums = []   
734    max_momentums = []
735    max_xmomentums = []
736    max_ymomentums = []
737    min_xmomentums = []
738    min_ymomentums = []
739    max_speeds = []
740    min_speeds = []   
741    max_depths = []
742    model_time_plot3d = num.zeros((n0, m), num.float)
743    stages_plot3d = num.zeros((n0, m), num.float)
744    eastings_plot3d = num.zeros((n0, m),num.float)
745    if time_unit is 'mins': scale = 60.0
746    if time_unit is 'hours': scale = 3600.0
747
748    ##### loop over each swwfile #####
749    for j, f in enumerate(f_list):
750        if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
751
752        starttime = f.starttime
753        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
754        fid_compare = open(comparefile, 'w')
755        file0 = file_loc[j] + 'gauges_t0.csv'
756        fid_0 = open(file0, 'w')
757
758        ##### loop over each gauge #####
759        for k in gauge_index:
760            if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
761
762            g = gauges[k]
763            min_stage = 10
764            max_stage = 0
765            max_momentum = max_xmomentum = max_ymomentum = 0
766            min_momentum = min_xmomentum = min_ymomentum = 100
767            max_speed = 0
768            min_speed = 0           
769            max_depth = 0           
770            gaugeloc = str(locations[k])
771            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
772                       + gaugeloc + '.csv'
773            if j == 0:
774                fid_out = open(thisfile, 'w')
775                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
776                fid_out.write(s)           
777
778            #### generate quantities #######
779            for i, t in enumerate(f.get_time()):
780                if time_min <= t <= time_max:
781                    w = f(t, point_id = k)[0]
782                    z = f(t, point_id = k)[1]
783                    uh = f(t, point_id = k)[2]
784                    vh = f(t, point_id = k)[3]
785                    depth = w-z     
786                    m = sqrt(uh*uh + vh*vh)
787                    if depth < 0.001:
788                        vel = 0.0
789                    else:
790                        vel = m / (depth + 1.e-6/depth) 
791                    bearing = calc_bearing(uh, vh)                   
792                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
793                    stages[i,k,j] = w
794                    elevations[i,k,j] = z
795                    xmom[i,k,j] = uh
796                    ymom[i,k,j] = vh
797                    momenta[i,k,j] = m
798                    speed[i,k,j] = vel
799                    bearings[i,k,j] = bearing
800                    depths[i,k,j] = depth
801                    thisgauge = gauges[k]
802                    eastings[i,k,j] = thisgauge[0]
803                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
804                            % (t, w, m, vel, z, uh, vh, bearing)
805                    fid_out.write(s)
806                    if t == 0:
807                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
808                        fid_0.write(s)
809                    if t/60.0 <= 13920: tindex = i
810                    if w > max_stage: max_stage = w
811                    if w < min_stage: min_stage = w
812                    if m > max_momentum: max_momentum = m
813                    if m < min_momentum: min_momentum = m                   
814                    if uh > max_xmomentum: max_xmomentum = uh
815                    if vh > max_ymomentum: max_ymomentum = vh
816                    if uh < min_xmomentum: min_xmomentum = uh
817                    if vh < min_ymomentum: min_ymomentum = vh
818                    if vel > max_speed: max_speed = vel
819                    if vel < min_speed: min_speed = vel                   
820                    if z > 0 and depth > max_depth: max_depth = depth
821                   
822                   
823            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
824                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
825            fid_compare.write(s)
826            max_stages.append(max_stage)
827            min_stages.append(min_stage)
828            max_momentums.append(max_momentum)
829            max_xmomentums.append(max_xmomentum)
830            max_ymomentums.append(max_ymomentum)
831            min_xmomentums.append(min_xmomentum)
832            min_ymomentums.append(min_ymomentum)
833            min_momentums.append(min_momentum)           
834            max_depths.append(max_depth)
835            max_speeds.append(max_speed)
836            min_speeds.append(min_speed)           
837            #### finished generating quantities for each swwfile #####
838       
839        model_time_plot3d[:,:] = model_time[:,:,j]
840        stages_plot3d[:,:] = stages[:,:,j]
841        eastings_plot3d[:,] = eastings[:,:,j]
842           
843        if surface is True:
844            log.critical('Printing surface figure')
845            for i in range(2):
846                fig = p1.figure(10)
847                ax = p3.Axes3D(fig)
848                if len(gauges) > 80:
849                    ax.plot_surface(model_time[:,:,j],
850                                    eastings[:,:,j],
851                                    stages[:,:,j])
852                else:
853                    ax.plot3D(num.ravel(eastings[:,:,j]),
854                              num.ravel(model_time[:,:,j]),
855                              num.ravel(stages[:,:,j]))
856                ax.set_xlabel('time')
857                ax.set_ylabel('x')
858                ax.set_zlabel('stage')
859                fig.add_axes(ax)
860                p1.show()
861                surfacefig = 'solution_surface%s' % leg_label[j]
862                p1.savefig(surfacefig)
863                p1.close()
864           
865    #### finished generating quantities for all swwfiles #####
866
867    # x profile for given time
868    if surface is True:
869        figure(11)
870        plot(eastings[tindex,:,j], stages[tindex,:,j])
871        xlabel('x')
872        ylabel('stage')
873        profilefig = 'solution_xprofile' 
874        savefig('profilefig')
875
876    elev_output = []
877    if generate_fig is True:
878        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
879                           max(max_depths)*1.1])
880        stage_axis = axis([starttime/scale, time_max/scale,
881                           min(min_stages), max(max_stages)*1.1])
882        vel_axis = axis([starttime/scale, time_max/scale,
883                         min(min_speeds), max(max_speeds)*1.1])
884        mom_axis = axis([starttime/scale, time_max/scale,
885                         min(min_momentums), max(max_momentums)*1.1])
886        xmom_axis = axis([starttime/scale, time_max/scale,
887                          min(min_xmomentums), max(max_xmomentums)*1.1])
888        ymom_axis = axis([starttime/scale, time_max/scale,
889                          min(min_ymomentums), max(max_ymomentums)*1.1])
890        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
891        nn = len(plot_quantity)
892        no_cols = 2
893       
894        if len(label_id) > 1: graphname_report = []
895        pp = 1
896        div = 11.
897        cc = 0
898        for k in gauge_index:
899            g = gauges[k]
900            count1 = 0
901            if report == True and len(label_id) > 1:
902                s = '\\begin{figure}[ht] \n' \
903                    '\\centering \n' \
904                    '\\begin{tabular}{cc} \n'
905                fid.write(s)
906            if len(label_id) > 1: graphname_report = []
907
908            #### generate figures for each gauge ####
909            for j, f in enumerate(f_list):
910                ion()
911                hold(True)
912                count = 0
913                where1 = 0
914                where2 = 0
915                word_quantity = ''
916                if report == True and len(label_id) == 1:
917                    s = '\\begin{figure}[hbt] \n' \
918                        '\\centering \n' \
919                        '\\begin{tabular}{cc} \n'
920                    fid.write(s)
921                   
922                for which_quantity in plot_quantity:
923                    count += 1
924                    where1 += 1
925                    figure(count, frameon = False)
926                    if which_quantity == 'depth':
927                        plot(model_time[0:n[j]-1,k,j],
928                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
929                        units = 'm'
930                        axis(depth_axis)
931                    if which_quantity == 'stage':
932                        if elevations[0,k,j] <= 0:
933                            plot(model_time[0:n[j]-1,k,j],
934                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
935                            axis(stage_axis)
936                        else:
937                            plot(model_time[0:n[j]-1,k,j],
938                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
939                            #axis(depth_axis)                 
940                        units = 'm'
941                    if which_quantity == 'momentum':
942                        plot(model_time[0:n[j]-1,k,j],
943                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
944                        axis(mom_axis)
945                        units = 'm^2 / sec'
946                    if which_quantity == 'xmomentum':
947                        plot(model_time[0:n[j]-1,k,j],
948                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
949                        axis(xmom_axis)
950                        units = 'm^2 / sec'
951                    if which_quantity == 'ymomentum':
952                        plot(model_time[0:n[j]-1,k,j],
953                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
954                        axis(ymom_axis)
955                        units = 'm^2 / sec'
956                    if which_quantity == 'speed':
957                        plot(model_time[0:n[j]-1,k,j],
958                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
959                        axis(vel_axis)
960                        units = 'm / sec'
961                    if which_quantity == 'bearing':
962                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
963                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
964                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
965                        units = 'degrees from North'
966                        #ax = axis([time_min, time_max, 0.0, 360.0])
967                        legend(('Bearing','West','East'))
968
969                    if time_unit is 'mins': xlabel('time (mins)')
970                    if time_unit is 'hours': xlabel('time (hours)')
971                    #if which_quantity == 'stage' \
972                    #   and elevations[0:n[j]-1,k,j] > 0:
973                    #    ylabel('%s (%s)' %('depth', units))
974                    #else:
975                    #    ylabel('%s (%s)' %(which_quantity, units))
976                        #ylabel('%s (%s)' %('wave height', units))
977                    ylabel('%s (%s)' %(which_quantity, units))
978                    if len(label_id) > 1: legend((leg_label),loc='upper right')
979
980                    #gaugeloc1 = gaugeloc.replace(' ','')
981                    #gaugeloc2 = gaugeloc1.replace('_','')
982                    gaugeloc2 = str(locations[k]).replace(' ','')
983                    graphname = '%sgauge%s_%s' %(file_loc[j],
984                                                 gaugeloc2,
985                                                 which_quantity)
986
987                    if report == True and len(label_id) > 1:
988                        figdir = getcwd()+sep+'report_figures'+sep
989                        if access(figdir,F_OK) == 0 :
990                            mkdir (figdir)
991                        latex_file_loc = figdir.replace(sep,altsep) 
992                        # storing files in production directory   
993                        graphname_latex = '%sgauge%s%s' \
994                                          % (latex_file_loc, gaugeloc2,
995                                             which_quantity)
996                        # giving location in latex output file
997                        graphname_report_input = '%sgauge%s%s' % \
998                                                 ('..' + altsep + 
999                                                      'report_figures' + altsep,
1000                                                  gaugeloc2, which_quantity)
1001                        graphname_report.append(graphname_report_input)
1002                       
1003                        # save figures in production directory for report
1004                        savefig(graphname_latex)
1005
1006                    if report == True:
1007                        figdir = getcwd() + sep + 'report_figures' + sep
1008                        if access(figdir,F_OK) == 0:
1009                            mkdir(figdir)
1010                        latex_file_loc = figdir.replace(sep,altsep)   
1011
1012                        if len(label_id) == 1: 
1013                            # storing files in production directory 
1014                            graphname_latex = '%sgauge%s%s%s' % \
1015                                              (latex_file_loc, gaugeloc2,
1016                                               which_quantity, label_id2)
1017                            # giving location in latex output file
1018                            graphname_report = '%sgauge%s%s%s' % \
1019                                               ('..' + altsep +
1020                                                    'report_figures' + altsep,
1021                                                gaugeloc2, which_quantity,
1022                                                label_id2)
1023                            s = '\includegraphics' \
1024                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1025                                (graphname_report, '.png')
1026                            fid.write(s)
1027                            if where1 % 2 == 0:
1028                                s = '\\\\ \n'
1029                                where1 = 0
1030                            else:
1031                                s = '& \n'
1032                            fid.write(s)
1033                            savefig(graphname_latex)
1034                   
1035                    if title_on == True:
1036                        title('%s scenario: %s at %s gauge' % \
1037                              (label_id, which_quantity, gaugeloc2))
1038                        #title('Gauge %s (MOST elevation %.2f, ' \
1039                        #      'ANUGA elevation %.2f)' % \
1040                        #      (gaugeloc2, elevations[10,k,0],
1041                        #       elevations[10,k,1]))
1042
1043                    savefig(graphname) # save figures with sww file
1044
1045                if report == True and len(label_id) == 1:
1046                    for i in range(nn-1):
1047                        if nn > 2:
1048                            if plot_quantity[i] == 'stage' \
1049                               and elevations[0,k,j] > 0:
1050                                word_quantity += 'depth' + ', '
1051                            else:
1052                                word_quantity += plot_quantity[i] + ', '
1053                        else:
1054                            if plot_quantity[i] == 'stage' \
1055                               and elevations[0,k,j] > 0:
1056                                word_quantity += 'depth' + ', '
1057                            else:
1058                                word_quantity += plot_quantity[i]
1059                       
1060                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
1061                        word_quantity += ' and ' + 'depth'
1062                    else:
1063                        word_quantity += ' and ' + plot_quantity[nn-1]
1064                    caption = 'Time series for %s at %s location ' \
1065                              '(elevation %.2fm)' % \
1066                              (word_quantity, locations[k], elev[k])
1067                    if elev[k] == 0.0:
1068                        caption = 'Time series for %s at %s location ' \
1069                                  '(elevation %.2fm)' % \
1070                                  (word_quantity, locations[k],
1071                                   elevations[0,k,j])
1072                        east = gauges[0]
1073                        north = gauges[1]
1074                        elev_output.append([locations[k], east, north,
1075                                            elevations[0,k,j]])
1076                    label = '%sgauge%s' % (label_id2, gaugeloc2)
1077                    s = '\end{tabular} \n' \
1078                        '\\caption{%s} \n' \
1079                        '\label{fig:%s} \n' \
1080                        '\end{figure} \n \n' % (caption, label)
1081                    fid.write(s)
1082                    cc += 1
1083                    if cc % 6 == 0: fid.write('\\clearpage \n')
1084                    savefig(graphname_latex)               
1085                   
1086            if report == True and len(label_id) > 1:
1087                for i in range(nn-1):
1088                    if nn > 2:
1089                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1090                            word_quantity += 'depth' + ','
1091                        else:
1092                            word_quantity += plot_quantity[i] + ', '
1093                    else:
1094                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1095                            word_quantity += 'depth'
1096                        else:
1097                            word_quantity += plot_quantity[i]
1098                    where1 = 0
1099                    count1 += 1
1100                    index = j*len(plot_quantity)
1101                    for which_quantity in plot_quantity:
1102                        where1 += 1
1103                        s = '\includegraphics' \
1104                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1105                            (graphname_report[index], '.png')
1106                        index += 1
1107                        fid.write(s)
1108                        if where1 % 2 == 0:
1109                            s = '\\\\ \n'
1110                            where1 = 0
1111                        else:
1112                            s = '& \n'
1113                        fid.write(s)
1114                word_quantity += ' and ' + plot_quantity[nn-1]           
1115                label = 'gauge%s' %(gaugeloc2) 
1116                caption = 'Time series for %s at %s location ' \
1117                          '(elevation %.2fm)' % \
1118                          (word_quantity, locations[k], elev[k])
1119                if elev[k] == 0.0:
1120                        caption = 'Time series for %s at %s location ' \
1121                                  '(elevation %.2fm)' % \
1122                                  (word_quantity, locations[k],
1123                                   elevations[0,k,j])
1124                        thisgauge = gauges[k]
1125                        east = thisgauge[0]
1126                        north = thisgauge[1]
1127                        elev_output.append([locations[k], east, north,
1128                                            elevations[0,k,j]])
1129                       
1130                s = '\end{tabular} \n' \
1131                    '\\caption{%s} \n' \
1132                    '\label{fig:%s} \n' \
1133                    '\end{figure} \n \n' % (caption, label)
1134                fid.write(s)
1135                if float((k+1)/div - pp) == 0.:
1136                    fid.write('\\clearpage \n')
1137                    pp += 1
1138                #### finished generating figures ###
1139
1140            close('all')
1141       
1142    return texfile2, elev_output
Note: See TracBrowser for help on using the repository browser.