source: anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py @ 7673

Last change on this file since 7673 was 7673, checked in by hudson, 14 years ago

Ticket 113 - added an output_centroids parameter to allow centroid data to be written to gauges, rather than the data at the sampled point.

Added 2 new unit tests for this functionality.

File size: 24.2 KB
Line 
1"""Gauge functions
2   
3   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
14
15import os
16
17from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep, getcwd
18from os.path import exists, split, join
19import anuga.utilities.log as log
20
21from math import sqrt
22
23##
24# @brief Take gauge readings, given a gauge file and a sww mesh
25#
26#        Use this function to take a timeseries sample, given a list of gauge points
27# @param  sww_file sww file to use as input
28# @param gauge_file gauge file as input, containing list of gauge points to sample
29# @param out_name output file prefix
30# @param quantities which quantities in the sww file we want to export
31# @param verbose show extra logging information for debug purposes
32# @param use_cache cache requests if possible, for speed
33# @param output_centroids Set to true to output the values at the centroid of the mesh triangle
34def gauge_sww2csv(sww_file,
35                   gauge_file,
36                   out_name='gauge_',
37                   quantities=['stage', 'depth', 'elevation',
38                               'xmomentum', 'ymomentum'],
39                   verbose=False,
40                   use_cache=True,
41                   output_centroids=False):
42    """
43   
44    Inputs:
45        NOTE: if using csv2timeseries_graphs after creating csv file,
46        it is essential to export quantities 'depth' and 'elevation'.
47        'depth' is good to analyse gauges on land and elevation is used
48        automatically by csv2timeseries_graphs in the legend.
49       
50        sww_file: path to any sww file
51       
52        gauge_file: Assumes that it follows this format
53            name, easting, northing, elevation
54            point1, 100.3, 50.2, 10.0
55            point2, 10.3, 70.3, 78.0
56       
57        NOTE: order of column can change but names eg 'easting', 'elevation'
58        must be the same! ALL lowercaps!
59
60        out_name: prefix for output file name (default is 'gauge_')
61       
62    Outputs:
63        one file for each gauge/point location in the points file. They
64        will be named with this format in the same directory as the 'sww_file'
65            <out_name><name>.csv
66        eg gauge_point1.csv if <out_name> not supplied
67           myfile_2_point1.csv if <out_name> ='myfile_2_'
68           
69        They will all have a header
70   
71    Usage: sww2csv_gauges(sww_file='test1.sww',
72                          quantities = ['stage', 'elevation','depth','bearing'],
73                          gauge_file='gauge.txt')   
74   
75    Interpolate the quantities at a given set of locations, given
76    an sww file.
77    The results are written to a csv file.
78
79    In the future let points be a points file.
80    And the user choose the quantities.
81
82    This is currently quite specific.
83    If it needs to be more general, change things.
84
85    This is really returning speed, not velocity.
86    """
87   
88    from csv import reader,writer
89    from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN
90    import string
91    from anuga.shallow_water.data_manager import get_all_swwfiles
92    from anuga.abstract_2d_finite_volumes.util import file_function     
93
94    assert type(gauge_file) == type(''), 'Gauge filename must be a string'
95    assert type(out_name) == type(''), 'Output filename prefix must be a string'
96   
97    try:
98        point_reader = reader(file(gauge_file))
99    except Exception, e:
100        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_file, e)
101        raise msg
102
103    if verbose: log.critical('Gauges obtained from: %s' % gauge_file)
104   
105    point_reader = reader(file(gauge_file))
106    points = []
107    point_name = []
108   
109    # read point info from file
110    for i,row in enumerate(point_reader):
111        # read header and determine the column numbers to read correctly.
112        if i==0:
113            for j,value in enumerate(row):
114                if value.strip()=='easting':easting=j
115                if value.strip()=='northing':northing=j
116                if value.strip()=='name':name=j
117                if value.strip()=='elevation':elevation=j
118        else:
119            points.append([float(row[easting]),float(row[northing])])
120            point_name.append(row[name])
121       
122    #convert to array for file_function
123    points_array = num.array(points,num.float)
124       
125    points_array = ensure_absolute(points_array)
126
127    dir_name, base = os.path.split(sww_file)   
128
129    #need to get current directory so when path and file
130    #are "joined" below the directory is correct
131    if dir_name == '':
132        dir_name =getcwd()
133       
134    if access(sww_file,R_OK):
135        if verbose: log.critical('File %s exists' % sww_file)
136    else:
137        msg = 'File "%s" could not be opened: no read permission' % sww_file
138        raise msg
139
140    sww_files = get_all_swwfiles(look_in_dir=dir_name,
141                                 base_name=base,
142                                 verbose=verbose)
143
144    # fudge to get SWW files in 'correct' order, oldest on the left
145    sww_files.sort()
146
147    if verbose:
148        log.critical('sww files=%s' % sww_files)
149   
150    #to make all the quantities lower case for file_function
151    quantities = [quantity.lower() for quantity in quantities]
152
153    # what is quantities are needed from sww file to calculate output quantities
154    # also
155
156    core_quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
157
158    gauge_file = out_name
159
160    heading = [quantity for quantity in quantities]
161    heading.insert(0,'time')
162    heading.insert(1,'hours')
163
164    #create a list of csv writers for all the points and write header
165    points_writer = []
166    for point_i,point in enumerate(points):
167        points_writer.append(writer(file(dir_name + sep + gauge_file
168                                         + point_name[point_i] + '.csv', "wb")))
169        points_writer[point_i].writerow(heading)
170   
171    if verbose: log.critical('Writing csv files')
172
173    quake_offset_time = None
174
175    for sww_file in sww_files:
176        sww_file = join(dir_name, sww_file+'.sww')
177        callable_sww = file_function(sww_file,
178                                     quantities=core_quantities,
179                                     interpolation_points=points_array,
180                                     verbose=verbose,
181                                     use_cache=use_cache,
182                                     output_centroids = output_centroids)
183
184        if quake_offset_time is None:
185            quake_offset_time = callable_sww.starttime
186
187        for time in callable_sww.get_time():
188            for point_i, point in enumerate(points_array):
189                #add domain starttime to relative time.
190                quake_time = time + quake_offset_time
191                points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
192                point_quantities = callable_sww(time,point_i) # __call__ is overridden
193                               
194                for quantity in quantities:
195                    if quantity == NAN:
196                        log.critical('quantity does not exist in %s'
197                                     % callable_sww.get_name)
198                    else:
199                        if quantity == 'stage':
200                            points_list.append(point_quantities[0])
201                           
202                        if quantity == 'elevation':
203                            points_list.append(point_quantities[1])
204                           
205                        if quantity == 'xmomentum':
206                            points_list.append(point_quantities[2])
207                           
208                        if quantity == 'ymomentum':
209                            points_list.append(point_quantities[3])
210                           
211                        if quantity == 'depth':
212                            points_list.append(point_quantities[0] 
213                                               - point_quantities[1])
214
215                        if quantity == 'momentum':
216                            momentum = sqrt(point_quantities[2]**2 
217                                            + point_quantities[3]**2)
218                            points_list.append(momentum)
219                           
220                        if quantity == 'speed':
221                            #if depth is less than 0.001 then speed = 0.0
222                            if point_quantities[0] - point_quantities[1] < 0.001:
223                                vel = 0.0
224                            else:
225                                if point_quantities[2] < 1.0e6:
226                                    momentum = sqrt(point_quantities[2]**2
227                                                    + point_quantities[3]**2)
228                                    vel = momentum / (point_quantities[0] 
229                                                      - point_quantities[1])
230                                else:
231                                    momentum = 0
232                                    vel = 0
233                               
234                            points_list.append(vel)
235                           
236                        if quantity == 'bearing':
237                            points_list.append(calc_bearing(point_quantities[2],
238                                                            point_quantities[3]))
239
240                points_writer[point_i].writerow(points_list)
241
242##
243# @brief Read a .sww file and plot the time series.
244# @param swwfiles Dictionary of .sww files.
245# @param gauge_filename Name of gauge file.
246# @param production_dirs ??
247# @param report If True, write figures to report directory.
248# @param reportname Name of generated report (if any).
249# @param plot_quantity List containing quantities to plot.
250# @param generate_fig If True, generate figures as well as CSV files.
251# @param surface If True, then generate solution surface with 3d plot.
252# @param time_min Beginning of user defined time range for plotting purposes.
253# @param time_max End of user defined time range for plotting purposes.
254# @param time_thinning ??
255# @param time_unit ??
256# @param title_on If True, export standard graphics with title.
257# @param use_cache If True, use caching.
258# @param verbose If True, this function is verbose.
259def gauge_sww2timeseries(swwfiles,
260                   gauge_filename,
261                   production_dirs,
262                   report=None,
263                   reportname=None,
264                   plot_quantity=None,
265                   generate_fig=False,
266                   surface=None,
267                   time_min=None,
268                   time_max=None,
269                   time_thinning=1,                   
270                   time_unit=None,
271                   title_on=None,
272                   use_cache=False,
273                   verbose=False,
274                                   output_centroids=False):
275    """ Read sww file and plot the time series for the
276    prescribed quantities at defined gauge locations and
277    prescribed time range.
278
279    Input variables:
280
281    swwfiles        - dictionary of sww files with label_ids (used in
282                      generating latex output. It will be part of
283                      the directory name of file_loc (typically the timestamp).
284                      Helps to differentiate latex files for different
285                      simulations for a particular scenario. 
286                    - assume that all conserved quantities have been stored
287                    - assume each sww file has been simulated with same timestep
288   
289    gauge_filename  - name of file containing gauge data
290                        - easting, northing, name , elevation?
291                    - OR (this is not yet done)
292                        - structure which can be converted to a numeric array,
293                          such as a geospatial data object
294                     
295    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
296                                                'boundaries': 'urs boundary'}
297                      this will use the second part as the label and the
298                      first part as the ?
299                      #FIXME: Is it a list or a dictionary
300                      # This is probably obsolete by now
301                     
302    report          - if True, then write figures to report_figures directory in
303                      relevant production directory
304                    - if False, figures are already stored with sww file
305                    - default to False
306
307    reportname      - name for report if wishing to generate report
308   
309    plot_quantity   - list containing quantities to plot, they must
310                      be the name of an existing quantity or one of
311                      the following possibilities
312                    - possibilities:
313                        - stage; 'stage'
314                        - depth; 'depth'
315                        - speed; calculated as absolute momentum
316                         (pointwise) divided by depth; 'speed'
317                        - bearing; calculated as the angle of the momentum
318                          vector (xmomentum, ymomentum) from the North; 'bearing'
319                        - absolute momentum; calculated as
320                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
321                        - x momentum; 'xmomentum'
322                        - y momentum; 'ymomentum'
323                    - default will be ['stage', 'speed', 'bearing']
324
325    generate_fig     - if True, generate figures as well as csv file
326                     - if False, csv files created only
327                     
328    surface          - if True, then generate solution surface with 3d plot
329                       and save to current working directory
330                     - default = False
331   
332    time_min         - beginning of user defined time range for plotting purposes
333                        - default will be first available time found in swwfile
334                       
335    time_max         - end of user defined time range for plotting purposes
336                        - default will be last available time found in swwfile
337                       
338    title_on        - if True, export standard graphics with title
339                    - if False, export standard graphics without title
340
341
342    Output:
343   
344    - time series data stored in .csv for later use if required.
345      Name = gauges_timeseries followed by gauge name
346    - latex file will be generated in same directory as where script is
347      run (usually production scenario directory.
348      Name = latexoutputlabel_id.tex
349
350    Other important information:
351   
352    It is assumed that the used has stored all the conserved quantities
353    and elevation during the scenario run, i.e.
354    ['stage', 'elevation', 'xmomentum', 'ymomentum']
355    If this has not occurred then sww2timeseries will not work.
356
357
358    Usage example
359    texname = sww2timeseries({project.boundary_name + '.sww': ''},
360                             project.polygons_dir + sep + 'boundary_extent.csv',
361                             project.anuga_dir,
362                             report = False,
363                             plot_quantity = ['stage', 'speed', 'bearing'],
364                             time_min = None,
365                             time_max = None,
366                             title_on = True,   
367                             verbose = True)
368   
369    """
370
371    msg = 'NOTE: A new function is available to create csv files from sww '
372    msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
373    msg += ' PLUS another new function to create graphs from csv files called '
374    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
375    log.critical(msg)
376   
377    k = _sww2timeseries(swwfiles,
378                        gauge_filename,
379                        production_dirs,
380                        report,
381                        reportname,
382                        plot_quantity,
383                        generate_fig,
384                        surface,
385                        time_min,
386                        time_max,
387                        time_thinning,                       
388                        time_unit,
389                        title_on,
390                        use_cache,
391                        verbose,
392                                                output_centroids = output_centroids)
393    return k
394
395
396##
397# @brief Read a .sww file and plot the time series.
398# @param swwfiles Dictionary of .sww files.
399# @param gauge_filename Name of gauge file.
400# @param production_dirs ??
401# @param report If True, write figures to report directory.
402# @param reportname Name of generated report (if any).
403# @param plot_quantity List containing quantities to plot.
404# @param generate_fig If True, generate figures as well as CSV files.
405# @param surface If True, then generate solution surface with 3d plot.
406# @param time_min Beginning of user defined time range for plotting purposes.
407# @param time_max End of user defined time range for plotting purposes.
408# @param time_thinning ??
409# @param time_unit ??
410# @param title_on If True, export standard graphics with title.
411# @param use_cache If True, use caching.
412# @param verbose If True, this function is verbose.
413def _sww2timeseries(swwfiles,
414                    gauge_filename,
415                    production_dirs,
416                    report = None,
417                    reportname = None,
418                    plot_quantity = None,
419                    generate_fig = False,
420                    surface = None,
421                    time_min = None,
422                    time_max = None,
423                    time_thinning = 1,                   
424                    time_unit = None,
425                    title_on = None,
426                    use_cache = False,
427                    verbose = False,
428                                        output_centroids = False):   
429       
430    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
431    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
432   
433    try:
434        fid = open(gauge_filename)
435    except Exception, e:
436        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
437        raise msg
438
439    if report is None:
440        report = False
441       
442    if plot_quantity is None:
443        plot_quantity = ['depth', 'speed']
444    else:
445        assert type(plot_quantity) == list, 'plot_quantity must be a list'
446        check_list(plot_quantity)
447
448    if surface is None:
449        surface = False
450
451    if time_unit is None:
452        time_unit = 'hours'
453   
454    if title_on is None:
455        title_on = True
456   
457    if verbose: log.critical('Gauges obtained from: %s' % gauge_filename)
458
459    gauges, locations, elev = get_gauges_from_file(gauge_filename)
460
461    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
462
463    file_loc = []
464    f_list = []
465    label_id = []
466    leg_label = []
467    themaxT = 0.0
468    theminT = 0.0
469
470    for swwfile in swwfiles.keys():
471        try:
472            fid = open(swwfile)
473        except Exception, e:
474            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
475            raise msg
476
477        if verbose:
478            log.critical('swwfile = %s' % swwfile)
479
480        # Extract parent dir name and use as label
481        path, _ = os.path.split(swwfile)
482        _, label = os.path.split(path)       
483       
484        leg_label.append(label)
485
486        f = file_function(swwfile,
487                          quantities = sww_quantity,
488                          interpolation_points = gauges,
489                          time_thinning = time_thinning,
490                          verbose = verbose,
491                          use_cache = use_cache,
492                                                  output_centroids = output_centroids)
493
494        # determine which gauges are contained in sww file
495        count = 0
496        gauge_index = []
497        for k, g in enumerate(gauges):
498            if f(0.0, point_id = k)[2] > 1.0e6:
499                count += 1
500                if count == 1: log.critical('Gauges not contained here:')
501                log.critical(locations[k])
502            else:
503                gauge_index.append(k)
504
505        if len(gauge_index) > 0:
506            log.critical('Gauges contained here:')
507        else:
508            log.critical('No gauges contained here.')
509        for i in range(len(gauge_index)):
510             log.critical(locations[gauge_index[i]])
511             
512        index = swwfile.rfind(sep)
513        file_loc.append(swwfile[:index+1])
514        label_id.append(swwfiles[swwfile])
515       
516        f_list.append(f)
517        maxT = max(f.get_time())
518        minT = min(f.get_time())
519        if maxT > themaxT: themaxT = maxT
520        if minT > theminT: theminT = minT
521
522    if time_min is None:
523        time_min = theminT # min(T)
524    else:
525        if time_min < theminT: # min(T):
526            msg = 'Minimum time entered not correct - please try again'
527            raise Exception, msg
528
529    if time_max is None:
530        time_max = themaxT # max(T)
531    else:
532        if time_max > themaxT: # max(T):
533            msg = 'Maximum time entered not correct - please try again'
534            raise Exception, msg
535
536    if verbose and len(gauge_index) > 0:
537         log.critical('Inputs OK - going to generate figures')
538
539    if len(gauge_index) <> 0:
540        texfile, elev_output = \
541            generate_figures(plot_quantity, file_loc, report, reportname,
542                             surface, leg_label, f_list, gauges, locations,
543                             elev, gauge_index, production_dirs, time_min,
544                             time_max, time_unit, title_on, label_id,
545                             generate_fig, verbose)
546    else:
547        texfile = ''
548        elev_output = []
549
550    return texfile, elev_output
551
552
553##
554# @brief Read gauge info from a file.
555# @param filename The name of the file to read.
556# @return A (gauges, gaugelocation, elev) tuple.
557def get_from_file(filename):
558    """ Read in gauge information from file
559    """
560
561    from os import sep, getcwd, access, F_OK, mkdir
562
563    # Get data from the gauge file
564    fid = open(filename)
565    lines = fid.readlines()
566    fid.close()
567   
568    gauges = []
569    gaugelocation = []
570    elev = []
571
572    # Check header information   
573    line1 = lines[0]
574    line11 = line1.split(',')
575
576    if isinstance(line11[0], str) is True:
577        # We have found text in the first line
578        east_index = None
579        north_index = None
580        name_index = None
581        elev_index = None
582
583        for i in range(len(line11)):
584            if line11[i].strip().lower() == 'easting':   east_index = i
585            if line11[i].strip().lower() == 'northing':  north_index = i
586            if line11[i].strip().lower() == 'name':      name_index = i
587            if line11[i].strip().lower() == 'elevation': elev_index = i
588
589        if east_index < len(line11) and north_index < len(line11):
590            pass
591        else:
592            msg = 'WARNING: %s does not contain correct header information' \
593                  % filename
594            msg += 'The header must be: easting, northing, name, elevation'
595            raise Exception, msg
596
597        if elev_index is None: 
598            raise Exception
599   
600        if name_index is None: 
601            raise Exception
602
603        lines = lines[1:] # Remove header from data
604    else:
605        # No header, assume that this is a simple easting, northing file
606
607        msg = 'There was no header in file %s and the number of columns is %d' \
608              % (filename, len(line11))
609        msg += '- was assuming two columns corresponding to Easting and Northing'
610        assert len(line11) == 2, msg
611
612        east_index = 0
613        north_index = 1
614
615        N = len(lines)
616        elev = [-9999]*N
617        gaugelocation = range(N)
618       
619    # Read in gauge data
620    for line in lines:
621        fields = line.split(',')
622
623        gauges.append([float(fields[east_index]), float(fields[north_index])])
624
625        if len(fields) > 2:
626            elev.append(float(fields[elev_index]))
627            loc = fields[name_index]
628            gaugelocation.append(loc.strip('\n'))
629
630    return gauges, gaugelocation, elev
Note: See TracBrowser for help on using the repository browser.