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

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

Ticket 113 is complete, and all tests pass.
A centroid list is built by Interpolation_function as it calculates the interpolation matrix, and this is passed out as extra quantities which are output into the gauge.csv file.

File size: 24.7 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                        #core quantities that are exported from the interpolator     
200                        if quantity == 'stage':
201                            points_list.append(point_quantities[0])
202                           
203                        if quantity == 'elevation':
204                            points_list.append(point_quantities[1])
205                           
206                        if quantity == 'xmomentum':
207                            points_list.append(point_quantities[2])
208                           
209                        if quantity == 'ymomentum':
210                            points_list.append(point_quantities[3])
211
212                        #derived quantities that are calculated from the core ones
213                        if quantity == 'depth':
214                            points_list.append(point_quantities[0] 
215                                               - point_quantities[1])
216
217                        if quantity == 'momentum':
218                            momentum = sqrt(point_quantities[2]**2 
219                                            + point_quantities[3]**2)
220                            points_list.append(momentum)
221                           
222                        if quantity == 'speed':
223                            #if depth is less than 0.001 then speed = 0.0
224                            if point_quantities[0] - point_quantities[1] < 0.001:
225                                vel = 0.0
226                            else:
227                                if point_quantities[2] < 1.0e6:
228                                    momentum = sqrt(point_quantities[2]**2
229                                                    + point_quantities[3]**2)
230                                    vel = momentum / (point_quantities[0] 
231                                                      - point_quantities[1])
232                                else:
233                                    momentum = 0
234                                    vel = 0
235                               
236                            points_list.append(vel)
237                           
238                        if quantity == 'bearing':
239                            points_list.append(calc_bearing(point_quantities[2],
240                                                            point_quantities[3]))
241                        if quantity == 'xcentroid':
242                            points_list.append(callable_sww.centroids[point_i][0])
243
244                        if quantity == 'ycentroid':
245                            points_list.append(callable_sww.centroids[point_i][1])
246                                                       
247                points_writer[point_i].writerow(points_list)
248
249##
250# @brief Read a .sww file and plot the time series.
251# @param swwfiles Dictionary of .sww files.
252# @param gauge_filename Name of gauge file.
253# @param production_dirs ??
254# @param report If True, write figures to report directory.
255# @param reportname Name of generated report (if any).
256# @param plot_quantity List containing quantities to plot.
257# @param generate_fig If True, generate figures as well as CSV files.
258# @param surface If True, then generate solution surface with 3d plot.
259# @param time_min Beginning of user defined time range for plotting purposes.
260# @param time_max End of user defined time range for plotting purposes.
261# @param time_thinning ??
262# @param time_unit ??
263# @param title_on If True, export standard graphics with title.
264# @param use_cache If True, use caching.
265# @param verbose If True, this function is verbose.
266def gauge_sww2timeseries(swwfiles,
267                   gauge_filename,
268                   production_dirs,
269                   report=None,
270                   reportname=None,
271                   plot_quantity=None,
272                   generate_fig=False,
273                   surface=None,
274                   time_min=None,
275                   time_max=None,
276                   time_thinning=1,                   
277                   time_unit=None,
278                   title_on=None,
279                   use_cache=False,
280                   verbose=False,
281                                   output_centroids=False):
282    """ Read sww file and plot the time series for the
283    prescribed quantities at defined gauge locations and
284    prescribed time range.
285
286    Input variables:
287
288    swwfiles        - dictionary of sww files with label_ids (used in
289                      generating latex output. It will be part of
290                      the directory name of file_loc (typically the timestamp).
291                      Helps to differentiate latex files for different
292                      simulations for a particular scenario. 
293                    - assume that all conserved quantities have been stored
294                    - assume each sww file has been simulated with same timestep
295   
296    gauge_filename  - name of file containing gauge data
297                        - easting, northing, name , elevation?
298                    - OR (this is not yet done)
299                        - structure which can be converted to a numeric array,
300                          such as a geospatial data object
301                     
302    production_dirs -  A list of list, example {20061101_121212: '1 in 10000',
303                                                'boundaries': 'urs boundary'}
304                      this will use the second part as the label and the
305                      first part as the ?
306                      #FIXME: Is it a list or a dictionary
307                      # This is probably obsolete by now
308                     
309    report          - if True, then write figures to report_figures directory in
310                      relevant production directory
311                    - if False, figures are already stored with sww file
312                    - default to False
313
314    reportname      - name for report if wishing to generate report
315   
316    plot_quantity   - list containing quantities to plot, they must
317                      be the name of an existing quantity or one of
318                      the following possibilities
319                    - possibilities:
320                        - stage; 'stage'
321                        - depth; 'depth'
322                        - speed; calculated as absolute momentum
323                         (pointwise) divided by depth; 'speed'
324                        - bearing; calculated as the angle of the momentum
325                          vector (xmomentum, ymomentum) from the North; 'bearing'
326                        - absolute momentum; calculated as
327                          sqrt(xmomentum^2 + ymomentum^2); 'momentum'
328                        - x momentum; 'xmomentum'
329                        - y momentum; 'ymomentum'
330                    - default will be ['stage', 'speed', 'bearing']
331
332    generate_fig     - if True, generate figures as well as csv file
333                     - if False, csv files created only
334                     
335    surface          - if True, then generate solution surface with 3d plot
336                       and save to current working directory
337                     - default = False
338   
339    time_min         - beginning of user defined time range for plotting purposes
340                        - default will be first available time found in swwfile
341                       
342    time_max         - end of user defined time range for plotting purposes
343                        - default will be last available time found in swwfile
344                       
345    title_on        - if True, export standard graphics with title
346                    - if False, export standard graphics without title
347
348
349    Output:
350   
351    - time series data stored in .csv for later use if required.
352      Name = gauges_timeseries followed by gauge name
353    - latex file will be generated in same directory as where script is
354      run (usually production scenario directory.
355      Name = latexoutputlabel_id.tex
356
357    Other important information:
358   
359    It is assumed that the used has stored all the conserved quantities
360    and elevation during the scenario run, i.e.
361    ['stage', 'elevation', 'xmomentum', 'ymomentum']
362    If this has not occurred then sww2timeseries will not work.
363
364
365    Usage example
366    texname = sww2timeseries({project.boundary_name + '.sww': ''},
367                             project.polygons_dir + sep + 'boundary_extent.csv',
368                             project.anuga_dir,
369                             report = False,
370                             plot_quantity = ['stage', 'speed', 'bearing'],
371                             time_min = None,
372                             time_max = None,
373                             title_on = True,   
374                             verbose = True)
375   
376    """
377
378    msg = 'NOTE: A new function is available to create csv files from sww '
379    msg += 'files called sww2csv_gauges in anuga.abstract_2d_finite_volumes.util'
380    msg += ' PLUS another new function to create graphs from csv files called '
381    msg += 'csv2timeseries_graphs in anuga.abstract_2d_finite_volumes.util'
382    log.critical(msg)
383   
384    k = _sww2timeseries(swwfiles,
385                        gauge_filename,
386                        production_dirs,
387                        report,
388                        reportname,
389                        plot_quantity,
390                        generate_fig,
391                        surface,
392                        time_min,
393                        time_max,
394                        time_thinning,                       
395                        time_unit,
396                        title_on,
397                        use_cache,
398                        verbose,
399                                                output_centroids = output_centroids)
400    return k
401
402
403##
404# @brief Read a .sww file and plot the time series.
405# @param swwfiles Dictionary of .sww files.
406# @param gauge_filename Name of gauge file.
407# @param production_dirs ??
408# @param report If True, write figures to report directory.
409# @param reportname Name of generated report (if any).
410# @param plot_quantity List containing quantities to plot.
411# @param generate_fig If True, generate figures as well as CSV files.
412# @param surface If True, then generate solution surface with 3d plot.
413# @param time_min Beginning of user defined time range for plotting purposes.
414# @param time_max End of user defined time range for plotting purposes.
415# @param time_thinning ??
416# @param time_unit ??
417# @param title_on If True, export standard graphics with title.
418# @param use_cache If True, use caching.
419# @param verbose If True, this function is verbose.
420def _sww2timeseries(swwfiles,
421                    gauge_filename,
422                    production_dirs,
423                    report = None,
424                    reportname = None,
425                    plot_quantity = None,
426                    generate_fig = False,
427                    surface = None,
428                    time_min = None,
429                    time_max = None,
430                    time_thinning = 1,                   
431                    time_unit = None,
432                    title_on = None,
433                    use_cache = False,
434                    verbose = False,
435                                        output_centroids = False):   
436       
437    # FIXME(Ole): Shouldn't print statements here be governed by verbose?
438    assert type(gauge_filename) == type(''), 'Gauge filename must be a string'
439   
440    try:
441        fid = open(gauge_filename)
442    except Exception, e:
443        msg = 'File "%s" could not be opened: Error="%s"' % (gauge_filename, e)
444        raise msg
445
446    if report is None:
447        report = False
448       
449    if plot_quantity is None:
450        plot_quantity = ['depth', 'speed']
451    else:
452        assert type(plot_quantity) == list, 'plot_quantity must be a list'
453        check_list(plot_quantity)
454
455    if surface is None:
456        surface = False
457
458    if time_unit is None:
459        time_unit = 'hours'
460   
461    if title_on is None:
462        title_on = True
463   
464    if verbose: log.critical('Gauges obtained from: %s' % gauge_filename)
465
466    gauges, locations, elev = get_gauges_from_file(gauge_filename)
467
468    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']
469
470    file_loc = []
471    f_list = []
472    label_id = []
473    leg_label = []
474    themaxT = 0.0
475    theminT = 0.0
476
477    for swwfile in swwfiles.keys():
478        try:
479            fid = open(swwfile)
480        except Exception, e:
481            msg = 'File "%s" could not be opened: Error="%s"' % (swwfile, e)
482            raise msg
483
484        if verbose:
485            log.critical('swwfile = %s' % swwfile)
486
487        # Extract parent dir name and use as label
488        path, _ = os.path.split(swwfile)
489        _, label = os.path.split(path)       
490       
491        leg_label.append(label)
492
493        f = file_function(swwfile,
494                          quantities = sww_quantity,
495                          interpolation_points = gauges,
496                          time_thinning = time_thinning,
497                          verbose = verbose,
498                          use_cache = use_cache,
499                                                  output_centroids = output_centroids)
500
501        # determine which gauges are contained in sww file
502        count = 0
503        gauge_index = []
504        for k, g in enumerate(gauges):
505            if f(0.0, point_id = k)[2] > 1.0e6:
506                count += 1
507                if count == 1: log.critical('Gauges not contained here:')
508                log.critical(locations[k])
509            else:
510                gauge_index.append(k)
511
512        if len(gauge_index) > 0:
513            log.critical('Gauges contained here:')
514        else:
515            log.critical('No gauges contained here.')
516        for i in range(len(gauge_index)):
517             log.critical(locations[gauge_index[i]])
518             
519        index = swwfile.rfind(sep)
520        file_loc.append(swwfile[:index+1])
521        label_id.append(swwfiles[swwfile])
522       
523        f_list.append(f)
524        maxT = max(f.get_time())
525        minT = min(f.get_time())
526        if maxT > themaxT: themaxT = maxT
527        if minT > theminT: theminT = minT
528
529    if time_min is None:
530        time_min = theminT # min(T)
531    else:
532        if time_min < theminT: # min(T):
533            msg = 'Minimum time entered not correct - please try again'
534            raise Exception, msg
535
536    if time_max is None:
537        time_max = themaxT # max(T)
538    else:
539        if time_max > themaxT: # max(T):
540            msg = 'Maximum time entered not correct - please try again'
541            raise Exception, msg
542
543    if verbose and len(gauge_index) > 0:
544         log.critical('Inputs OK - going to generate figures')
545
546    if len(gauge_index) <> 0:
547        texfile, elev_output = \
548            generate_figures(plot_quantity, file_loc, report, reportname,
549                             surface, leg_label, f_list, gauges, locations,
550                             elev, gauge_index, production_dirs, time_min,
551                             time_max, time_unit, title_on, label_id,
552                             generate_fig, verbose)
553    else:
554        texfile = ''
555        elev_output = []
556
557    return texfile, elev_output
558
559
560##
561# @brief Read gauge info from a file.
562# @param filename The name of the file to read.
563# @return A (gauges, gaugelocation, elev) tuple.
564def get_from_file(filename):
565    """ Read in gauge information from file
566    """
567
568    from os import sep, getcwd, access, F_OK, mkdir
569
570    # Get data from the gauge file
571    fid = open(filename)
572    lines = fid.readlines()
573    fid.close()
574   
575    gauges = []
576    gaugelocation = []
577    elev = []
578
579    # Check header information   
580    line1 = lines[0]
581    line11 = line1.split(',')
582
583    if isinstance(line11[0], str) is True:
584        # We have found text in the first line
585        east_index = None
586        north_index = None
587        name_index = None
588        elev_index = None
589
590        for i in range(len(line11)):
591            if line11[i].strip().lower() == 'easting':   east_index = i
592            if line11[i].strip().lower() == 'northing':  north_index = i
593            if line11[i].strip().lower() == 'name':      name_index = i
594            if line11[i].strip().lower() == 'elevation': elev_index = i
595
596        if east_index < len(line11) and north_index < len(line11):
597            pass
598        else:
599            msg = 'WARNING: %s does not contain correct header information' \
600                  % filename
601            msg += 'The header must be: easting, northing, name, elevation'
602            raise Exception, msg
603
604        if elev_index is None: 
605            raise Exception
606   
607        if name_index is None: 
608            raise Exception
609
610        lines = lines[1:] # Remove header from data
611    else:
612        # No header, assume that this is a simple easting, northing file
613
614        msg = 'There was no header in file %s and the number of columns is %d' \
615              % (filename, len(line11))
616        msg += '- was assuming two columns corresponding to Easting and Northing'
617        assert len(line11) == 2, msg
618
619        east_index = 0
620        north_index = 1
621
622        N = len(lines)
623        elev = [-9999]*N
624        gaugelocation = range(N)
625       
626    # Read in gauge data
627    for line in lines:
628        fields = line.split(',')
629
630        gauges.append([float(fields[east_index]), float(fields[north_index])])
631
632        if len(fields) > 2:
633            elev.append(float(fields[elev_index]))
634            loc = fields[name_index]
635            gaugelocation.append(loc.strip('\n'))
636
637    return gauges, gaugelocation, elev
Note: See TracBrowser for help on using the repository browser.