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

Last change on this file since 8068 was 8063, checked in by jakeman, 13 years ago

updated sww2csv_gauges so that it can now read multiple temporaly consecutive sww files. This was previously broken. Unit test also included in test_gauges.py called test_sww2csv_multiple_gauges
File size: 46.8 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    is_opened = [False]*len(points_array)
226    for sww_file in sww_files:
227        sww_file = join(dir_name, sww_file+'.sww')
228        callable_sww = file_function(sww_file,
229                                     quantities=core_quantities,
230                                     interpolation_points=points_array,
231                                     verbose=verbose,
232                                     use_cache=use_cache,
233                                     output_centroids = output_centroids)
234
235        if quake_offset_time is None:
236            quake_offset_time = callable_sww.starttime
237
238        for point_i, point in enumerate(points_array):
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[point_i] == 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[point_i] = True
250                    else:
251                        points_writer = writer(file(dir_name + sep + gauge_file
252                                                    + point_name[point_i] + '.csv', "ab"))
253
254                    points_list = [quake_time, quake_time/3600.] +  _quantities2csv(quantities, point_quantities, callable_sww.centroids, point_i)
255                    points_writer.writerow(points_list)
256                else:
257                    if verbose:
258                        msg = 'gauge' + point_name[point_i] + 'falls off the mesh in file ' + sww_file + '.'
259                        log.warning(msg)
260##
261# @brief Read a .sww file and plot the time series.
262# @param swwfiles Dictionary of .sww files.
263# @param gauge_filename Name of gauge file.
264# @param production_dirs ??
265# @param report If True, write figures to report directory.
266# @param reportname Name of generated report (if any).
267# @param plot_quantity List containing quantities to plot.
268# @param generate_fig If True, generate figures as well as CSV files.
269# @param surface If True, then generate solution surface with 3d plot.
270# @param time_min Beginning of user defined time range for plotting purposes.
271# @param time_max End of user defined time range for plotting purposes.
272# @param time_thinning ??
273# @param time_unit ??
274# @param title_on If True, export standard graphics with title.
275# @param use_cache If True, use caching.
276# @param verbose If True, this function is verbose.
277def sww2timeseries(swwfiles,
278                   gauge_filename,
279                   production_dirs,
280                   report=None,
281                   reportname=None,
282                   plot_quantity=None,
283                   generate_fig=False,
284                   surface=None,
285                   time_min=None,
286                   time_max=None,
287                   time_thinning=1,                   
288                   time_unit=None,
289                   title_on=None,
290                   use_cache=False,
291                   verbose=False,
292                   output_centroids=False):
293    """ Read sww file and plot the time series for the
294    prescribed quantities at defined gauge locations and
295    prescribed time range.
296
297    Input variables:
298
299    swwfiles        - dictionary of sww files with label_ids (used in
300                      generating latex output. It will be part of
301                      the directory name of file_loc (typically the timestamp).
302                      Helps to differentiate latex files for different
303                      simulations for a particular scenario. 
304                    - assume that all conserved quantities have been stored
305                    - assume each sww file has been simulated with same timestep
306   
307    gauge_filename  - name of file containing gauge data
308                        - easting, northing, name , elevation?
309                    - OR (this is not yet done)
310                        - structure which can be converted to a numeric array,
311                          such as a geospatial data object
312                     
313    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
314                                                'boundaries': 'urs boundary'}
315                      this will use the second part as the label and the
316                      first part as the ?
317                      #FIXME: Is it a list or a dictionary
318                      # This is probably obsolete by now
319                     
320    report          - if True, then write figures to report_figures directory in
321                      relevant production directory
322                    - if False, figures are already stored with sww file
323                    - default to False
324
325    reportname      - name for report if wishing to generate report
326   
327    plot_quantity   - list containing quantities to plot, they must
328                      be the name of an existing quantity or one of
329                      the following possibilities
330                    - possibilities:
331                        - stage; 'stage'
332                        - depth; 'depth'
333                        - speed; calculated as absolute momentum
334                         (pointwise) divided by depth; 'speed'
335                        - bearing; calculated as the angle of the momentum
336                          vector (xmomentum, ymomentum) from the North; 'bearing'
337                        - absolute momentum; calculated as
338                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
339                        - x momentum; 'xmomentum'
340                        - y momentum; 'ymomentum'
341                    - default will be ['stage', 'speed', 'bearing']
342
343    generate_fig     - if True, generate figures as well as csv file
344                     - if False, csv files created only
345                     
346    surface          - if True, then generate solution surface with 3d plot
347                       and save to current working directory
348                     - default = False
349   
350    time_min         - beginning of user defined time range for plotting purposes
351                        - default will be first available time found in swwfile
352                       
353    time_max         - end of user defined time range for plotting purposes
354                        - default will be last available time found in swwfile
355                       
356    title_on        - if True, export standard graphics with title
357                    - if False, export standard graphics without title
358
359
360    Output:
361   
362    - time series data stored in .csv for later use if required.
363      Name = gauges_timeseries followed by gauge name
364    - latex file will be generated in same directory as where script is
365      run (usually production scenario directory.
366      Name = latexoutputlabel_id.tex
367
368    Other important information:
369   
370    It is assumed that the used has stored all the conserved quantities
371    and elevation during the scenario run, i.e.
372    ['stage', 'elevation', 'xmomentum', 'ymomentum']
373    If this has not occurred then sww2timeseries will not work.
374
375
376    Usage example
377    texname = sww2timeseries({project.boundary_name + '.sww': ''},
378                             project.polygons_dir + sep + 'boundary_extent.csv',
379                             project.anuga_dir,
380                             report = False,
381                             plot_quantity = ['stage', 'speed', 'bearing'],
382                             time_min = None,
383                             time_max = None,
384                             title_on = True,   
385                             verbose = True)
386   
387    """
388
389    msg = 'NOTE: A new function is available to create csv files from sww '
390    msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
391    msg += ' PLUS another new function to create graphs from csv files called '
392    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
393    log.critical(msg)
394   
395    k = _sww2timeseries(swwfiles,
396                        gauge_filename,
397                        production_dirs,
398                        report,
399                        reportname,
400                        plot_quantity,
401                        generate_fig,
402                        surface,
403                        time_min,
404                        time_max,
405                        time_thinning,                       
406                        time_unit,
407                        title_on,
408                        use_cache,
409                        verbose,
410                        output_centroids = output_centroids)
411    return k
412
413
414##
415# @brief Read a .sww file and plot the time series.
416# @param swwfiles Dictionary of .sww files.
417# @param gauge_filename Name of gauge file.
418# @param production_dirs ??
419# @param report If True, write figures to report directory.
420# @param reportname Name of generated report (if any).
421# @param plot_quantity List containing quantities to plot.
422# @param generate_fig If True, generate figures as well as CSV files.
423# @param surface If True, then generate solution surface with 3d plot.
424# @param time_min Beginning of user defined time range for plotting purposes.
425# @param time_max End of user defined time range for plotting purposes.
426# @param time_thinning ??
427# @param time_unit ??
428# @param title_on If True, export standard graphics with title.
429# @param use_cache If True, use caching.
430# @param verbose If True, this function is verbose.
431def _sww2timeseries(swwfiles,
432                    gauge_filename,
433                    production_dirs,
434                    report = None,
435                    reportname = None,
436                    plot_quantity = None,
437                    generate_fig = False,
438                    surface = None,
439                    time_min = None,
440                    time_max = None,
441                    time_thinning = 1,                   
442                    time_unit = None,
443                    title_on = None,
444                    use_cache = False,
445                    verbose = False,
446                    output_centroids = False):   
447       
448    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
449    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
450   
451    try:
452        fid = open(gauge_filename)
453    except Exception, e:
454        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
455        raise msg
456
457    if report is None:
458        report = False
459       
460    if plot_quantity is None:
461        plot_quantity = ['depth', 'speed']
462    else:
463        assert type(plot_quantity) == list, 'plot_quantity must be a list'
464        check_list(plot_quantity)
465
466    if surface is None:
467        surface = False
468
469    if time_unit is None:
470        time_unit = 'hours'
471   
472    if title_on is None:
473        title_on = True
474   
475    if verbose: log.critical('Gauges obtained from: %s' % gauge_filename)
476
477    gauges, locations, elev = gauge_get_from_file(gauge_filename)
478
479    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
480
481    file_loc = []
482    f_list = []
483    label_id = []
484    leg_label = []
485    themaxT = 0.0
486    theminT = 0.0
487
488    for swwfile in swwfiles.keys():
489        try:
490            fid = open(swwfile)
491        except Exception, e:
492            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
493            raise msg
494
495        if verbose:
496            log.critical('swwfile = %s' % swwfile)
497
498        # Extract parent dir name and use as label
499        path, _ = os.path.split(swwfile)
500        _, label = os.path.split(path)       
501       
502        leg_label.append(label)
503
504        f = file_function(swwfile,
505                          quantities = sww_quantity,
506                          interpolation_points = gauges,
507                          time_thinning = time_thinning,
508                          verbose = verbose,
509                          use_cache = use_cache,
510                          output_centroids = output_centroids)
511
512        # determine which gauges are contained in sww file
513        count = 0
514        gauge_index = []
515        for k, g in enumerate(gauges):
516            if f(0.0, point_id = k)[2] > 1.0e6:
517                count += 1
518                if count == 1: log.critical('Gauges not contained here:')
519                log.critical(locations[k])
520            else:
521                gauge_index.append(k)
522
523        if len(gauge_index) > 0:
524            log.critical('Gauges contained here:')
525        else:
526            log.critical('No gauges contained here.')
527        for i in range(len(gauge_index)):
528             log.critical(locations[gauge_index[i]])
529             
530        index = swwfile.rfind(sep)
531        file_loc.append(swwfile[:index+1])
532        label_id.append(swwfiles[swwfile])
533       
534        f_list.append(f)
535        maxT = max(f.get_time())
536        minT = min(f.get_time())
537        if maxT > themaxT: themaxT = maxT
538        if minT > theminT: theminT = minT
539
540    if time_min is None:
541        time_min = theminT # min(T)
542    else:
543        if time_min < theminT: # min(T):
544            msg = 'Minimum time entered not correct - please try again'
545            raise Exception, msg
546
547    if time_max is None:
548        time_max = themaxT # max(T)
549    else:
550        if time_max > themaxT: # max(T):
551            msg = 'Maximum time entered not correct - please try again'
552            raise Exception, msg
553
554    if verbose and len(gauge_index) > 0:
555         log.critical('Inputs OK - going to generate figures')
556
557    if len(gauge_index) <> 0:
558        texfile, elev_output = \
559            _generate_figures(plot_quantity, file_loc, report, reportname,
560                             surface, leg_label, f_list, gauges, locations,
561                             elev, gauge_index, production_dirs, time_min,
562                             time_max, time_unit, title_on, label_id,
563                             generate_fig, verbose)
564    else:
565        texfile = ''
566        elev_output = []
567
568    return texfile, elev_output
569
570
571##
572# @brief Read gauge info from a file.
573# @param filename The name of the file to read.
574# @return A (gauges, gaugelocation, elev) tuple.
575def gauge_get_from_file(filename):
576    """ Read in gauge information from file
577    """
578
579    from os import sep, getcwd, access, F_OK, mkdir
580
581    # Get data from the gauge file
582    fid = open(filename)
583    lines = fid.readlines()
584    fid.close()
585   
586    gauges = []
587    gaugelocation = []
588    elev = []
589
590    # Check header information   
591    line1 = lines[0]
592    line11 = line1.split(',')
593
594    if isinstance(line11[0], str) is True:
595        # We have found text in the first line
596        east_index = None
597        north_index = None
598        name_index = None
599        elev_index = None
600
601        for i in range(len(line11)):
602            if line11[i].strip().lower() == 'easting':   east_index = i
603            if line11[i].strip().lower() == 'northing':  north_index = i
604            if line11[i].strip().lower() == 'name':      name_index = i
605            if line11[i].strip().lower() == 'elevation': elev_index = i
606
607        if east_index < len(line11) and north_index < len(line11):
608            pass
609        else:
610            msg = 'WARNING: %s does not contain correct header information' \
611                  % filename
612            msg += 'The header must be: easting, northing, name, elevation'
613            raise Exception, msg
614
615        if elev_index is None: 
616            raise Exception
617   
618        if name_index is None: 
619            raise Exception
620
621        lines = lines[1:] # Remove header from data
622    else:
623        # No header, assume that this is a simple easting, northing file
624
625        msg = 'There was no header in file %s and the number of columns is %d' \
626              % (filename, len(line11))
627        msg += '- was assuming two columns corresponding to Easting and Northing'
628        assert len(line11) == 2, msg
629
630        east_index = 0
631        north_index = 1
632
633        N = len(lines)
634        elev = [-9999]*N
635        gaugelocation = range(N)
636       
637    # Read in gauge data
638    for line in lines:
639        fields = line.split(',')
640
641        gauges.append([float(fields[east_index]), float(fields[north_index])])
642
643        if len(fields) > 2:
644            elev.append(float(fields[elev_index]))
645            loc = fields[name_index]
646            gaugelocation.append(loc.strip('\n'))
647
648    return gauges, gaugelocation, elev
649
650   
651##
652# @brief Generate figures from quantities and gauges for each sww file.
653#        This is a legacy function, used only by sww2timeseries
654# @param plot_quantity  ??
655# @param file_loc ??
656# @param report ??
657# @param reportname ??
658# @param surface ??
659# @param leg_label ??
660# @param f_list ??
661# @param gauges ??
662# @param locations ??
663# @param elev ??
664# @param gauge_index ??
665# @param production_dirs ??
666# @param time_min ??
667# @param time_max ??
668# @param time_unit ??
669# @param title_on ??
670# @param label_id ??
671# @param generate_fig ??
672# @param verbose??
673# @return (texfile2, elev_output)
674def _generate_figures(plot_quantity, file_loc, report, reportname, surface,
675                     leg_label, f_list, gauges, locations, elev, gauge_index,
676                     production_dirs, time_min, time_max, time_unit,
677                     title_on, label_id, generate_fig, verbose):
678    """ Generate figures based on required quantities and gauges for
679    each sww file
680    """
681    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
682
683    if generate_fig is True:
684        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
685             xlabel, ylabel, title, close, subplot
686   
687        if surface is True:
688            import pylab as p1
689            import mpl3d.mplot3d as p3
690       
691    if report == True:   
692        texdir = getcwd()+sep+'report'+sep
693        if access(texdir,F_OK) == 0:
694            mkdir (texdir)
695        if len(label_id) == 1:
696            label_id1 = label_id[0].replace(sep,'')
697            label_id2 = label_id1.replace('_','')
698            texfile = texdir + reportname + '%s' % label_id2
699            texfile2 = reportname + '%s' % label_id2
700            texfilename = texfile + '.tex'
701            fid = open(texfilename, 'w')
702
703            if verbose: log.critical('Latex output printed to %s' % texfilename)
704        else:
705            texfile = texdir+reportname
706            texfile2 = reportname
707            texfilename = texfile + '.tex' 
708            fid = open(texfilename, 'w')
709
710            if verbose: log.critical('Latex output printed to %s' % texfilename)
711    else:
712        texfile = ''
713        texfile2 = ''
714
715    p = len(f_list)
716    n = []
717    n0 = 0
718    for i in range(len(f_list)):
719        n.append(len(f_list[i].get_time()))
720        if n[i] > n0: n0 = n[i] 
721    n0 = int(n0)
722    m = len(locations)
723    model_time = num.zeros((n0, m, p), num.float) 
724    stages = num.zeros((n0, m, p), num.float)
725    elevations = num.zeros((n0, m, p), num.float) 
726    momenta = num.zeros((n0, m, p), num.float)
727    xmom = num.zeros((n0, m, p), num.float)
728    ymom = num.zeros((n0, m, p), num.float)
729    speed = num.zeros((n0, m, p), num.float)
730    bearings = num.zeros((n0, m, p), num.float)
731    due_east = 90.0*num.ones((n0, 1), num.float)
732    due_west = 270.0*num.ones((n0, 1), num.float)
733    depths = num.zeros((n0, m, p), num.float)
734    eastings = num.zeros((n0, m, p), num.float)
735    min_stages = []
736    max_stages = []
737    min_momentums = []   
738    max_momentums = []
739    max_xmomentums = []
740    max_ymomentums = []
741    min_xmomentums = []
742    min_ymomentums = []
743    max_speeds = []
744    min_speeds = []   
745    max_depths = []
746    model_time_plot3d = num.zeros((n0, m), num.float)
747    stages_plot3d = num.zeros((n0, m), num.float)
748    eastings_plot3d = num.zeros((n0, m),num.float)
749    if time_unit is 'mins': scale = 60.0
750    if time_unit is 'hours': scale = 3600.0
751
752    ##### loop over each swwfile #####
753    for j, f in enumerate(f_list):
754        if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
755
756        starttime = f.starttime
757        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
758        fid_compare = open(comparefile, 'w')
759        file0 = file_loc[j] + 'gauges_t0.csv'
760        fid_0 = open(file0, 'w')
761
762        ##### loop over each gauge #####
763        for k in gauge_index:
764            if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
765
766            g = gauges[k]
767            min_stage = 10
768            max_stage = 0
769            max_momentum = max_xmomentum = max_ymomentum = 0
770            min_momentum = min_xmomentum = min_ymomentum = 100
771            max_speed = 0
772            min_speed = 0           
773            max_depth = 0           
774            gaugeloc = str(locations[k])
775            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
776                       + gaugeloc + '.csv'
777            if j == 0:
778                fid_out = open(thisfile, 'w')
779                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
780                fid_out.write(s)           
781
782            #### generate quantities #######
783            for i, t in enumerate(f.get_time()):
784                if time_min <= t <= time_max:
785                    w = f(t, point_id = k)[0]
786                    z = f(t, point_id = k)[1]
787                    uh = f(t, point_id = k)[2]
788                    vh = f(t, point_id = k)[3]
789                    depth = w-z     
790                    m = sqrt(uh*uh + vh*vh)
791                    if depth < 0.001:
792                        vel = 0.0
793                    else:
794                        vel = m / (depth + 1.e-6/depth) 
795                    bearing = calc_bearing(uh, vh)                   
796                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
797                    stages[i,k,j] = w
798                    elevations[i,k,j] = z
799                    xmom[i,k,j] = uh
800                    ymom[i,k,j] = vh
801                    momenta[i,k,j] = m
802                    speed[i,k,j] = vel
803                    bearings[i,k,j] = bearing
804                    depths[i,k,j] = depth
805                    thisgauge = gauges[k]
806                    eastings[i,k,j] = thisgauge[0]
807                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
808                            % (t, w, m, vel, z, uh, vh, bearing)
809                    fid_out.write(s)
810                    if t == 0:
811                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
812                        fid_0.write(s)
813                    if t/60.0 <= 13920: tindex = i
814                    if w > max_stage: max_stage = w
815                    if w < min_stage: min_stage = w
816                    if m > max_momentum: max_momentum = m
817                    if m < min_momentum: min_momentum = m                   
818                    if uh > max_xmomentum: max_xmomentum = uh
819                    if vh > max_ymomentum: max_ymomentum = vh
820                    if uh < min_xmomentum: min_xmomentum = uh
821                    if vh < min_ymomentum: min_ymomentum = vh
822                    if vel > max_speed: max_speed = vel
823                    if vel < min_speed: min_speed = vel                   
824                    if z > 0 and depth > max_depth: max_depth = depth
825                   
826                   
827            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
828                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
829            fid_compare.write(s)
830            max_stages.append(max_stage)
831            min_stages.append(min_stage)
832            max_momentums.append(max_momentum)
833            max_xmomentums.append(max_xmomentum)
834            max_ymomentums.append(max_ymomentum)
835            min_xmomentums.append(min_xmomentum)
836            min_ymomentums.append(min_ymomentum)
837            min_momentums.append(min_momentum)           
838            max_depths.append(max_depth)
839            max_speeds.append(max_speed)
840            min_speeds.append(min_speed)           
841            #### finished generating quantities for each swwfile #####
842       
843        model_time_plot3d[:,:] = model_time[:,:,j]
844        stages_plot3d[:,:] = stages[:,:,j]
845        eastings_plot3d[:,] = eastings[:,:,j]
846           
847        if surface is True:
848            log.critical('Printing surface figure')
849            for i in range(2):
850                fig = p1.figure(10)
851                ax = p3.Axes3D(fig)
852                if len(gauges) > 80:
853                    ax.plot_surface(model_time[:,:,j],
854                                    eastings[:,:,j],
855                                    stages[:,:,j])
856                else:
857                    ax.plot3D(num.ravel(eastings[:,:,j]),
858                              num.ravel(model_time[:,:,j]),
859                              num.ravel(stages[:,:,j]))
860                ax.set_xlabel('time')
861                ax.set_ylabel('x')
862                ax.set_zlabel('stage')
863                fig.add_axes(ax)
864                p1.show()
865                surfacefig = 'solution_surface%s' % leg_label[j]
866                p1.savefig(surfacefig)
867                p1.close()
868           
869    #### finished generating quantities for all swwfiles #####
870
871    # x profile for given time
872    if surface is True:
873        figure(11)
874        plot(eastings[tindex,:,j], stages[tindex,:,j])
875        xlabel('x')
876        ylabel('stage')
877        profilefig = 'solution_xprofile' 
878        savefig('profilefig')
879
880    elev_output = []
881    if generate_fig is True:
882        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
883                           max(max_depths)*1.1])
884        stage_axis = axis([starttime/scale, time_max/scale,
885                           min(min_stages), max(max_stages)*1.1])
886        vel_axis = axis([starttime/scale, time_max/scale,
887                         min(min_speeds), max(max_speeds)*1.1])
888        mom_axis = axis([starttime/scale, time_max/scale,
889                         min(min_momentums), max(max_momentums)*1.1])
890        xmom_axis = axis([starttime/scale, time_max/scale,
891                          min(min_xmomentums), max(max_xmomentums)*1.1])
892        ymom_axis = axis([starttime/scale, time_max/scale,
893                          min(min_ymomentums), max(max_ymomentums)*1.1])
894        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
895        nn = len(plot_quantity)
896        no_cols = 2
897       
898        if len(label_id) > 1: graphname_report = []
899        pp = 1
900        div = 11.
901        cc = 0
902        for k in gauge_index:
903            g = gauges[k]
904            count1 = 0
905            if report == True and len(label_id) > 1:
906                s = '\\begin{figure}[ht] \n' \
907                    '\\centering \n' \
908                    '\\begin{tabular}{cc} \n'
909                fid.write(s)
910            if len(label_id) > 1: graphname_report = []
911
912            #### generate figures for each gauge ####
913            for j, f in enumerate(f_list):
914                ion()
915                hold(True)
916                count = 0
917                where1 = 0
918                where2 = 0
919                word_quantity = ''
920                if report == True and len(label_id) == 1:
921                    s = '\\begin{figure}[hbt] \n' \
922                        '\\centering \n' \
923                        '\\begin{tabular}{cc} \n'
924                    fid.write(s)
925                   
926                for which_quantity in plot_quantity:
927                    count += 1
928                    where1 += 1
929                    figure(count, frameon = False)
930                    if which_quantity == 'depth':
931                        plot(model_time[0:n[j]-1,k,j],
932                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
933                        units = 'm'
934                        axis(depth_axis)
935                    if which_quantity == 'stage':
936                        if elevations[0,k,j] <= 0:
937                            plot(model_time[0:n[j]-1,k,j],
938                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
939                            axis(stage_axis)
940                        else:
941                            plot(model_time[0:n[j]-1,k,j],
942                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
943                            #axis(depth_axis)                 
944                        units = 'm'
945                    if which_quantity == 'momentum':
946                        plot(model_time[0:n[j]-1,k,j],
947                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
948                        axis(mom_axis)
949                        units = 'm^2 / sec'
950                    if which_quantity == 'xmomentum':
951                        plot(model_time[0:n[j]-1,k,j],
952                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
953                        axis(xmom_axis)
954                        units = 'm^2 / sec'
955                    if which_quantity == 'ymomentum':
956                        plot(model_time[0:n[j]-1,k,j],
957                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
958                        axis(ymom_axis)
959                        units = 'm^2 / sec'
960                    if which_quantity == 'speed':
961                        plot(model_time[0:n[j]-1,k,j],
962                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
963                        axis(vel_axis)
964                        units = 'm / sec'
965                    if which_quantity == 'bearing':
966                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
967                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
968                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
969                        units = 'degrees from North'
970                        #ax = axis([time_min, time_max, 0.0, 360.0])
971                        legend(('Bearing','West','East'))
972
973                    if time_unit is 'mins': xlabel('time (mins)')
974                    if time_unit is 'hours': xlabel('time (hours)')
975                    #if which_quantity == 'stage' \
976                    #   and elevations[0:n[j]-1,k,j] > 0:
977                    #    ylabel('%s (%s)' %('depth', units))
978                    #else:
979                    #    ylabel('%s (%s)' %(which_quantity, units))
980                        #ylabel('%s (%s)' %('wave height', units))
981                    ylabel('%s (%s)' %(which_quantity, units))
982                    if len(label_id) > 1: legend((leg_label),loc='upper right')
983
984                    #gaugeloc1 = gaugeloc.replace(' ','')
985                    #gaugeloc2 = gaugeloc1.replace('_','')
986                    gaugeloc2 = str(locations[k]).replace(' ','')
987                    graphname = '%sgauge%s_%s' %(file_loc[j],
988                                                 gaugeloc2,
989                                                 which_quantity)
990
991                    if report == True and len(label_id) > 1:
992                        figdir = getcwd()+sep+'report_figures'+sep
993                        if access(figdir,F_OK) == 0 :
994                            mkdir (figdir)
995                        latex_file_loc = figdir.replace(sep,altsep) 
996                        # storing files in production directory   
997                        graphname_latex = '%sgauge%s%s' \
998                                          % (latex_file_loc, gaugeloc2,
999                                             which_quantity)
1000                        # giving location in latex output file
1001                        graphname_report_input = '%sgauge%s%s' % \
1002                                                 ('..' + altsep + 
1003                                                      'report_figures' + altsep,
1004                                                  gaugeloc2, which_quantity)
1005                        graphname_report.append(graphname_report_input)
1006                       
1007                        # save figures in production directory for report
1008                        savefig(graphname_latex)
1009
1010                    if report == True:
1011                        figdir = getcwd() + sep + 'report_figures' + sep
1012                        if access(figdir,F_OK) == 0:
1013                            mkdir(figdir)
1014                        latex_file_loc = figdir.replace(sep,altsep)   
1015
1016                        if len(label_id) == 1: 
1017                            # storing files in production directory 
1018                            graphname_latex = '%sgauge%s%s%s' % \
1019                                              (latex_file_loc, gaugeloc2,
1020                                               which_quantity, label_id2)
1021                            # giving location in latex output file
1022                            graphname_report = '%sgauge%s%s%s' % \
1023                                               ('..' + altsep +
1024                                                    'report_figures' + altsep,
1025                                                gaugeloc2, which_quantity,
1026                                                label_id2)
1027                            s = '\includegraphics' \
1028                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1029                                (graphname_report, '.png')
1030                            fid.write(s)
1031                            if where1 % 2 == 0:
1032                                s = '\\\\ \n'
1033                                where1 = 0
1034                            else:
1035                                s = '& \n'
1036                            fid.write(s)
1037                            savefig(graphname_latex)
1038                   
1039                    if title_on == True:
1040                        title('%s scenario: %s at %s gauge' % \
1041                              (label_id, which_quantity, gaugeloc2))
1042                        #title('Gauge %s (MOST elevation %.2f, ' \
1043                        #      'ANUGA elevation %.2f)' % \
1044                        #      (gaugeloc2, elevations[10,k,0],
1045                        #       elevations[10,k,1]))
1046
1047                    savefig(graphname) # save figures with sww file
1048
1049                if report == True and len(label_id) == 1:
1050                    for i in range(nn-1):
1051                        if nn > 2:
1052                            if plot_quantity[i] == 'stage' \
1053                               and elevations[0,k,j] > 0:
1054                                word_quantity += 'depth' + ', '
1055                            else:
1056                                word_quantity += plot_quantity[i] + ', '
1057                        else:
1058                            if plot_quantity[i] == 'stage' \
1059                               and elevations[0,k,j] > 0:
1060                                word_quantity += 'depth' + ', '
1061                            else:
1062                                word_quantity += plot_quantity[i]
1063                       
1064                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
1065                        word_quantity += ' and ' + 'depth'
1066                    else:
1067                        word_quantity += ' and ' + plot_quantity[nn-1]
1068                    caption = 'Time series for %s at %s location ' \
1069                              '(elevation %.2fm)' % \
1070                              (word_quantity, locations[k], elev[k])
1071                    if elev[k] == 0.0:
1072                        caption = 'Time series for %s at %s location ' \
1073                                  '(elevation %.2fm)' % \
1074                                  (word_quantity, locations[k],
1075                                   elevations[0,k,j])
1076                        east = gauges[0]
1077                        north = gauges[1]
1078                        elev_output.append([locations[k], east, north,
1079                                            elevations[0,k,j]])
1080                    label = '%sgauge%s' % (label_id2, gaugeloc2)
1081                    s = '\end{tabular} \n' \
1082                        '\\caption{%s} \n' \
1083                        '\label{fig:%s} \n' \
1084                        '\end{figure} \n \n' % (caption, label)
1085                    fid.write(s)
1086                    cc += 1
1087                    if cc % 6 == 0: fid.write('\\clearpage \n')
1088                    savefig(graphname_latex)               
1089                   
1090            if report == True and len(label_id) > 1:
1091                for i in range(nn-1):
1092                    if nn > 2:
1093                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1094                            word_quantity += 'depth' + ','
1095                        else:
1096                            word_quantity += plot_quantity[i] + ', '
1097                    else:
1098                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
1099                            word_quantity += 'depth'
1100                        else:
1101                            word_quantity += plot_quantity[i]
1102                    where1 = 0
1103                    count1 += 1
1104                    index = j*len(plot_quantity)
1105                    for which_quantity in plot_quantity:
1106                        where1 += 1
1107                        s = '\includegraphics' \
1108                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
1109                            (graphname_report[index], '.png')
1110                        index += 1
1111                        fid.write(s)
1112                        if where1 % 2 == 0:
1113                            s = '\\\\ \n'
1114                            where1 = 0
1115                        else:
1116                            s = '& \n'
1117                        fid.write(s)
1118                word_quantity += ' and ' + plot_quantity[nn-1]           
1119                label = 'gauge%s' %(gaugeloc2) 
1120                caption = 'Time series for %s at %s location ' \
1121                          '(elevation %.2fm)' % \
1122                          (word_quantity, locations[k], elev[k])
1123                if elev[k] == 0.0:
1124                        caption = 'Time series for %s at %s location ' \
1125                                  '(elevation %.2fm)' % \
1126                                  (word_quantity, locations[k],
1127                                   elevations[0,k,j])
1128                        thisgauge = gauges[k]
1129                        east = thisgauge[0]
1130                        north = thisgauge[1]
1131                        elev_output.append([locations[k], east, north,
1132                                            elevations[0,k,j]])
1133                       
1134                s = '\end{tabular} \n' \
1135                    '\\caption{%s} \n' \
1136                    '\label{fig:%s} \n' \
1137                    '\end{figure} \n \n' % (caption, label)
1138                fid.write(s)
1139                if float((k+1)/div - pp) == 0.:
1140                    fid.write('\\clearpage \n')
1141                    pp += 1
1142                #### finished generating figures ###
1143
1144            close('all')
1145       
1146    return texfile2, elev_output
Note: See TracBrowser for help on using the repository browser.