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

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

Ticket 344 completed and unit tests pass.
Stripped out a heap of tabs and replaced them with 4 spaces.
Removed some redundant imports and deprecated functions.

File size: 64.6 KB
Line 
1"""This module contains various auxiliary function used by pyvolution.
2
3It is also a clearing house for functions that may later earn a module
4of their own.
5"""
6
7import anuga.utilities.polygon
8import sys
9import os
10
11from os import remove, mkdir, access, F_OK, R_OK, W_OK, sep,getcwd
12from os.path import exists, basename, split,join
13from warnings import warn
14from shutil import copy
15
16from anuga.utilities.numerical_tools import ensure_numeric
17
18from math import sqrt, atan, degrees
19
20# FIXME (Ole): Temporary short cuts -
21# FIXME (Ole): remove and update scripts where they are used
22from anuga.utilities.system_tools import get_revision_number
23from anuga.utilities.system_tools import store_version_info
24
25import anuga.utilities.log as log
26
27import numpy as num
28
29
30def file_function(filename,
31                  domain=None,
32                  quantities=None,
33                  interpolation_points=None,
34                  time_thinning=1,
35                  time_limit=None,
36                  verbose=False,
37                  use_cache=False,
38                  boundary_polygon=None,
39                  output_centroids=False):
40    from file_function import file_function as file_function_new
41    return file_function_new(filename, domain, quantities, interpolation_points,
42                      time_thinning, time_limit, verbose, use_cache,
43                      boundary_polygon, output_centroids)
44
45##
46# @brief Replace multiple substrings in a string.
47# @param text The string to operate on.
48# @param dictionary A dict containing replacements, key->value.
49# @return The new string.
50def multiple_replace(text, dictionary):
51    """Multiple replace of words in text
52
53    text:       String to be modified
54    dictionary: Mapping of words that are to be substituted
55
56    Python Cookbook 3.14 page 88 and page 90
57    http://code.activestate.com/recipes/81330/
58    """
59
60    import re
61   
62    #Create a regular expression from all of the dictionary keys
63    #matching only entire words
64    regex = re.compile(r'\b'+ \
65                       r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \
66                       r'\b' )
67
68    #For each match, lookup the corresponding value in the dictionary
69    return regex.sub(lambda match: dictionary[match.group(0)], text)
70
71
72##
73# @brief Apply arbitrary expressions to the values of a dict.
74# @param expression A string expression to apply.
75# @param dictionary The dictionary to apply the expression to.
76def apply_expression_to_dictionary(expression, dictionary):
77    """Apply arbitrary expression to values of dictionary
78
79    Given an expression in terms of the keys, replace key by the
80    corresponding values and evaluate.
81
82    expression: Arbitrary, e.g. arithmetric, expression relating keys
83                from dictionary.
84
85    dictionary: Mapping between symbols in expression and objects that
86                will be evaluated by expression.
87                Values in dictionary must support operators given in
88                expression e.g. by overloading
89
90    Due to a limitation with numeric, this can not evaluate 0/0
91    In general, the user can fix by adding 1e-30 to the numerator.
92    SciPy core can handle this situation.
93    """
94
95    import types
96    import re
97
98    assert isinstance(expression, basestring)
99    assert type(dictionary) == types.DictType
100
101    #Convert dictionary values to textual representations suitable for eval
102    D = {}
103    for key in dictionary:
104        D[key] = 'dictionary["%s"]' % key
105
106    #Perform substitution of variables   
107    expression = multiple_replace(expression, D)
108
109    #Evaluate and return
110    try:
111        return eval(expression)
112    except NameError, e:
113        msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)
114        raise NameError, msg
115    except ValueError, e:
116        msg = 'Expression "%s" could not be evaluated: %s' % (expression, e)
117        raise ValueError, msg
118
119
120##
121# @brief Format a float into a string.
122# @param value Float value to format.
123# @param format The format to use (%.2f is default).
124# @return The formatted float as a string.
125def get_textual_float(value, format = '%.2f'):
126    """Get textual representation of floating point numbers
127    and accept None as valid entry
128
129    format is a string - default = '%.2f'
130    """
131
132    if value is None:
133        return 'None'
134    else:
135        try:
136            float(value)
137        except:
138            # May this is a vector
139            if len(value) > 1:
140                s = '('
141                for v in value:
142                    s += get_textual_float(v, format) + ', '
143                   
144                s = s[:-2] + ')' # Strip trailing comma and close
145                return s
146            else:
147                raise 'Illegal input to get_textual_float:', value
148        else:
149            return format % float(value)
150
151
152#################################################################################
153# OBSOLETE STUFF
154#################################################################################
155
156# @note TEMP
157def angle(v1, v2):
158    """Temporary Interface to new location"""
159
160    import anuga.utilities.numerical_tools as NT   
161   
162    msg = 'angle has moved from util.py.  '
163    msg += 'Please use "from anuga.utilities.numerical_tools import angle"'
164    warn(msg, DeprecationWarning) 
165
166    return NT.angle(v1, v2)
167
168   
169# @note TEMP
170def anglediff(v0, v1):
171    """Temporary Interface to new location"""
172
173    import anuga.utilities.numerical_tools as NT
174   
175    msg = 'anglediff has moved from util.py.  '
176    msg += 'Please use "from anuga.utilities.numerical_tools import anglediff"'
177    warn(msg, DeprecationWarning) 
178
179    return NT.anglediff(v0, v1)   
180
181# @note TEMP
182def point_on_line(*args, **kwargs):
183    """Temporary Interface to new location"""
184
185    msg = 'point_on_line has moved from util.py.  '
186    msg += 'Please use "from anuga.utilities.polygon import point_on_line"'
187    warn(msg, DeprecationWarning) 
188
189    return utilities.polygon.point_on_line(*args, **kwargs)   
190
191   
192# @note TEMP
193def inside_polygon(*args, **kwargs):
194    """Temporary Interface to new location"""
195
196    log.critical('inside_polygon has moved from util.py.')
197    log.critical('Please use '
198                 '"from anuga.utilities.polygon import inside_polygon"')
199
200    return utilities.polygon.inside_polygon(*args, **kwargs)   
201
202   
203# @note TEMP
204def outside_polygon(*args, **kwargs):
205    """Temporary Interface to new location"""
206
207    log.critical('outside_polygon has moved from util.py.')
208    log.critical('Please use '
209                 '"from anuga.utilities.polygon import outside_polygon"')
210
211    return utilities.polygon.outside_polygon(*args, **kwargs)   
212
213   
214# @note TEMP
215def read_polygon(*args, **kwargs):
216    """Temporary Interface to new location"""
217
218    log.critical('read_polygon has moved from util.py.')
219    log.critical('Please use '
220                 '"from anuga.utilities.polygon import read_polygon"')
221
222    return utilities.polygon.read_polygon(*args, **kwargs)   
223
224
225# @note TEMP
226def populate_polygon(*args, **kwargs):
227    """Temporary Interface to new location"""
228
229    log.critical('populate_polygon has moved from util.py.')
230    log.critical('Please use '
231                 '"from anuga.utilities.polygon import populate_polygon"')
232
233    return utilities.polygon.populate_polygon(*args, **kwargs)   
234
235
236#################################################################################
237# End of obsolete stuff ?
238#################################################################################
239
240# @note TEMP
241def start_screen_catcher(dir_name, myid='', numprocs='', extra_info='',
242                         verbose=False):
243    """Temporary Interface to new location"""
244    from anuga.shallow_water.data_manager import start_screen_catcher \
245         as dm_start_screen_catcher
246
247    log.critical('start_screen_catcher has moved from util.py.')
248    log.critical('Please use "from anuga.shallow_water.data_manager import '
249                 'start_screen_catcher"')
250   
251    return dm_start_screen_catcher(dir_name, myid='', numprocs='',
252                                   extra_info='', verbose=False)
253
254
255
256##
257# @brief Read gauge info from a file.
258# @param filename The name of the file to read.
259# @return A (gauges, gaugelocation, elev) tuple.
260def get_gauges_from_file(filename):
261    return gauge_get_from_file(filename)
262
263
264##
265# @brief Check that input quantities in quantity list are legal.
266# @param quantity Quantity list to check.
267# @note Raises an exception of list is not legal.
268def check_list(quantity):
269    """ Check that input quantities in quantity list are possible
270    """
271    import sys
272
273    all_quantity = ['stage', 'depth', 'momentum', 'xmomentum',
274                    'ymomentum', 'speed', 'bearing', 'elevation']
275
276    # convert all quanitiy names to lowercase
277    for i,j in enumerate(quantity):
278        quantity[i] = quantity[i].lower()
279
280    # check that all names in 'quantity' appear in 'all_quantity'
281    p = list(set(quantity).difference(set(all_quantity)))
282    if len(p) != 0:
283        msg = 'Quantities %s do not exist - please try again' %p
284        raise Exception, msg
285
286
287##
288# @brief Calculate velocity bearing from North.
289# @param uh ??
290# @param vh ??
291# @return The calculated bearing.
292def calc_bearing(uh, vh):
293    """ Calculate velocity bearing from North
294    """
295    #FIXME (Ole): I reckon we should refactor this one to use
296    #             the function angle() in utilities/numerical_tools
297    #
298    #             It will be a simple matter of
299    # * converting from radians to degrees
300    # * moving the reference direction from [1,0] to North
301    # * changing from counter clockwise to clocwise.
302       
303    angle = degrees(atan(vh/(uh+1.e-15)))
304
305    if (0 < angle < 90.0):
306        if vh > 0:
307            bearing = 90.0 - abs(angle)
308        if vh < 0:
309            bearing = 270.0 - abs(angle)
310   
311    if (-90 < angle < 0):
312        if vh < 0:
313            bearing = 90.0 - (angle)
314        if vh > 0:
315            bearing = 270.0 - (angle)
316    if angle == 0: bearing = 0.0
317
318    return bearing
319
320
321##
322# @brief Generate figures from quantities and gauges for each sww file.
323# @param plot_quantity  ??
324# @param file_loc ??
325# @param report ??
326# @param reportname ??
327# @param surface ??
328# @param leg_label ??
329# @param f_list ??
330# @param gauges ??
331# @param locations ??
332# @param elev ??
333# @param gauge_index ??
334# @param production_dirs ??
335# @param time_min ??
336# @param time_max ??
337# @param time_unit ??
338# @param title_on ??
339# @param label_id ??
340# @param generate_fig ??
341# @param verbose??
342# @return (texfile2, elev_output)
343def generate_figures(plot_quantity, file_loc, report, reportname, surface,
344                     leg_label, f_list, gauges, locations, elev, gauge_index,
345                     production_dirs, time_min, time_max, time_unit,
346                     title_on, label_id, generate_fig, verbose):
347    """ Generate figures based on required quantities and gauges for
348    each sww file
349    """
350    from os import sep, altsep, getcwd, mkdir, access, F_OK, environ
351
352    if generate_fig is True:
353        from pylab import ion, hold, plot, axis, figure, legend, savefig, \
354             xlabel, ylabel, title, close, subplot
355   
356        if surface is True:
357            import pylab as p1
358            import mpl3d.mplot3d as p3
359       
360    if report == True:   
361        texdir = getcwd()+sep+'report'+sep
362        if access(texdir,F_OK) == 0:
363            mkdir (texdir)
364        if len(label_id) == 1:
365            label_id1 = label_id[0].replace(sep,'')
366            label_id2 = label_id1.replace('_','')
367            texfile = texdir + reportname + '%s' % label_id2
368            texfile2 = reportname + '%s' % label_id2
369            texfilename = texfile + '.tex'
370            fid = open(texfilename, 'w')
371
372            if verbose: log.critical('Latex output printed to %s' % texfilename)
373        else:
374            texfile = texdir+reportname
375            texfile2 = reportname
376            texfilename = texfile + '.tex' 
377            fid = open(texfilename, 'w')
378
379            if verbose: log.critical('Latex output printed to %s' % texfilename)
380    else:
381        texfile = ''
382        texfile2 = ''
383
384    p = len(f_list)
385    n = []
386    n0 = 0
387    for i in range(len(f_list)):
388        n.append(len(f_list[i].get_time()))
389        if n[i] > n0: n0 = n[i] 
390    n0 = int(n0)
391    m = len(locations)
392    model_time = num.zeros((n0, m, p), num.float) 
393    stages = num.zeros((n0, m, p), num.float)
394    elevations = num.zeros((n0, m, p), num.float) 
395    momenta = num.zeros((n0, m, p), num.float)
396    xmom = num.zeros((n0, m, p), num.float)
397    ymom = num.zeros((n0, m, p), num.float)
398    speed = num.zeros((n0, m, p), num.float)
399    bearings = num.zeros((n0, m, p), num.float)
400    due_east = 90.0*num.ones((n0, 1), num.float)
401    due_west = 270.0*num.ones((n0, 1), num.float)
402    depths = num.zeros((n0, m, p), num.float)
403    eastings = num.zeros((n0, m, p), num.float)
404    min_stages = []
405    max_stages = []
406    min_momentums = []   
407    max_momentums = []
408    max_xmomentums = []
409    max_ymomentums = []
410    min_xmomentums = []
411    min_ymomentums = []
412    max_speeds = []
413    min_speeds = []   
414    max_depths = []
415    model_time_plot3d = num.zeros((n0, m), num.float)
416    stages_plot3d = num.zeros((n0, m), num.float)
417    eastings_plot3d = num.zeros((n0, m),num.float)
418    if time_unit is 'mins': scale = 60.0
419    if time_unit is 'hours': scale = 3600.0
420
421    ##### loop over each swwfile #####
422    for j, f in enumerate(f_list):
423        if verbose: log.critical('swwfile %d of %d' % (j, len(f_list)))
424
425        starttime = f.starttime
426        comparefile = file_loc[j] + sep + 'gauges_maxmins' + '.csv'
427        fid_compare = open(comparefile, 'w')
428        file0 = file_loc[j] + 'gauges_t0.csv'
429        fid_0 = open(file0, 'w')
430
431        ##### loop over each gauge #####
432        for k in gauge_index:
433            if verbose: log.critical('Gauge %d of %d' % (k, len(gauges)))
434
435            g = gauges[k]
436            min_stage = 10
437            max_stage = 0
438            max_momentum = max_xmomentum = max_ymomentum = 0
439            min_momentum = min_xmomentum = min_ymomentum = 100
440            max_speed = 0
441            min_speed = 0           
442            max_depth = 0           
443            gaugeloc = str(locations[k])
444            thisfile = file_loc[j] + sep + 'gauges_time_series' + '_' \
445                       + gaugeloc + '.csv'
446            if j == 0:
447                fid_out = open(thisfile, 'w')
448                s = 'Time, Stage, Momentum, Speed, Elevation, xmom, ymom, Bearing \n'
449                fid_out.write(s)           
450
451            #### generate quantities #######
452            for i, t in enumerate(f.get_time()):
453                if time_min <= t <= time_max:
454                    w = f(t, point_id = k)[0]
455                    z = f(t, point_id = k)[1]
456                    uh = f(t, point_id = k)[2]
457                    vh = f(t, point_id = k)[3]
458                    depth = w-z     
459                    m = sqrt(uh*uh + vh*vh)
460                    if depth < 0.001:
461                        vel = 0.0
462                    else:
463                        vel = m / (depth + 1.e-6/depth) 
464                    bearing = calc_bearing(uh, vh)                   
465                    model_time[i,k,j] = (t + starttime)/scale #t/60.0
466                    stages[i,k,j] = w
467                    elevations[i,k,j] = z
468                    xmom[i,k,j] = uh
469                    ymom[i,k,j] = vh
470                    momenta[i,k,j] = m
471                    speed[i,k,j] = vel
472                    bearings[i,k,j] = bearing
473                    depths[i,k,j] = depth
474                    thisgauge = gauges[k]
475                    eastings[i,k,j] = thisgauge[0]
476                    s = '%.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f, %.2f,\n' \
477                            % (t, w, m, vel, z, uh, vh, bearing)
478                    fid_out.write(s)
479                    if t == 0:
480                        s = '%.2f, %.2f, %.2f\n' % (g[0], g[1], w)
481                        fid_0.write(s)
482                    if t/60.0 <= 13920: tindex = i
483                    if w > max_stage: max_stage = w
484                    if w < min_stage: min_stage = w
485                    if m > max_momentum: max_momentum = m
486                    if m < min_momentum: min_momentum = m                   
487                    if uh > max_xmomentum: max_xmomentum = uh
488                    if vh > max_ymomentum: max_ymomentum = vh
489                    if uh < min_xmomentum: min_xmomentum = uh
490                    if vh < min_ymomentum: min_ymomentum = vh
491                    if vel > max_speed: max_speed = vel
492                    if vel < min_speed: min_speed = vel                   
493                    if z > 0 and depth > max_depth: max_depth = depth
494                   
495                   
496            s = '%.2f, %.2f, %.2f, %.2f, %s\n' \
497                    % (max_stage, min_stage, z, thisgauge[0], leg_label[j])
498            fid_compare.write(s)
499            max_stages.append(max_stage)
500            min_stages.append(min_stage)
501            max_momentums.append(max_momentum)
502            max_xmomentums.append(max_xmomentum)
503            max_ymomentums.append(max_ymomentum)
504            min_xmomentums.append(min_xmomentum)
505            min_ymomentums.append(min_ymomentum)
506            min_momentums.append(min_momentum)           
507            max_depths.append(max_depth)
508            max_speeds.append(max_speed)
509            min_speeds.append(min_speed)           
510            #### finished generating quantities for each swwfile #####
511       
512        model_time_plot3d[:,:] = model_time[:,:,j]
513        stages_plot3d[:,:] = stages[:,:,j]
514        eastings_plot3d[:,] = eastings[:,:,j]
515           
516        if surface is True:
517            log.critical('Printing surface figure')
518            for i in range(2):
519                fig = p1.figure(10)
520                ax = p3.Axes3D(fig)
521                if len(gauges) > 80:
522                    ax.plot_surface(model_time[:,:,j],
523                                    eastings[:,:,j],
524                                    stages[:,:,j])
525                else:
526                    ax.plot3D(num.ravel(eastings[:,:,j]),
527                              num.ravel(model_time[:,:,j]),
528                              num.ravel(stages[:,:,j]))
529                ax.set_xlabel('time')
530                ax.set_ylabel('x')
531                ax.set_zlabel('stage')
532                fig.add_axes(ax)
533                p1.show()
534                surfacefig = 'solution_surface%s' % leg_label[j]
535                p1.savefig(surfacefig)
536                p1.close()
537           
538    #### finished generating quantities for all swwfiles #####
539
540    # x profile for given time
541    if surface is True:
542        figure(11)
543        plot(eastings[tindex,:,j], stages[tindex,:,j])
544        xlabel('x')
545        ylabel('stage')
546        profilefig = 'solution_xprofile' 
547        savefig('profilefig')
548
549    elev_output = []
550    if generate_fig is True:
551        depth_axis = axis([starttime/scale, time_max/scale, -0.1,
552                           max(max_depths)*1.1])
553        stage_axis = axis([starttime/scale, time_max/scale,
554                           min(min_stages), max(max_stages)*1.1])
555        vel_axis = axis([starttime/scale, time_max/scale,
556                         min(min_speeds), max(max_speeds)*1.1])
557        mom_axis = axis([starttime/scale, time_max/scale,
558                         min(min_momentums), max(max_momentums)*1.1])
559        xmom_axis = axis([starttime/scale, time_max/scale,
560                          min(min_xmomentums), max(max_xmomentums)*1.1])
561        ymom_axis = axis([starttime/scale, time_max/scale,
562                          min(min_ymomentums), max(max_ymomentums)*1.1])
563        cstr = ['g', 'r', 'b', 'c', 'm', 'y', 'k']
564        nn = len(plot_quantity)
565        no_cols = 2
566       
567        if len(label_id) > 1: graphname_report = []
568        pp = 1
569        div = 11.
570        cc = 0
571        for k in gauge_index:
572            g = gauges[k]
573            count1 = 0
574            if report == True and len(label_id) > 1:
575                s = '\\begin{figure}[ht] \n' \
576                    '\\centering \n' \
577                    '\\begin{tabular}{cc} \n'
578                fid.write(s)
579            if len(label_id) > 1: graphname_report = []
580
581            #### generate figures for each gauge ####
582            for j, f in enumerate(f_list):
583                ion()
584                hold(True)
585                count = 0
586                where1 = 0
587                where2 = 0
588                word_quantity = ''
589                if report == True and len(label_id) == 1:
590                    s = '\\begin{figure}[hbt] \n' \
591                        '\\centering \n' \
592                        '\\begin{tabular}{cc} \n'
593                    fid.write(s)
594                   
595                for which_quantity in plot_quantity:
596                    count += 1
597                    where1 += 1
598                    figure(count, frameon = False)
599                    if which_quantity == 'depth':
600                        plot(model_time[0:n[j]-1,k,j],
601                             depths[0:n[j]-1,k,j], '-', c = cstr[j])
602                        units = 'm'
603                        axis(depth_axis)
604                    if which_quantity == 'stage':
605                        if elevations[0,k,j] <= 0:
606                            plot(model_time[0:n[j]-1,k,j],
607                                 stages[0:n[j]-1,k,j], '-', c = cstr[j])
608                            axis(stage_axis)
609                        else:
610                            plot(model_time[0:n[j]-1,k,j],
611                                 depths[0:n[j]-1,k,j], '-', c = cstr[j])
612                            #axis(depth_axis)                 
613                        units = 'm'
614                    if which_quantity == 'momentum':
615                        plot(model_time[0:n[j]-1,k,j],
616                             momenta[0:n[j]-1,k,j], '-', c = cstr[j])
617                        axis(mom_axis)
618                        units = 'm^2 / sec'
619                    if which_quantity == 'xmomentum':
620                        plot(model_time[0:n[j]-1,k,j],
621                             xmom[0:n[j]-1,k,j], '-', c = cstr[j])
622                        axis(xmom_axis)
623                        units = 'm^2 / sec'
624                    if which_quantity == 'ymomentum':
625                        plot(model_time[0:n[j]-1,k,j],
626                             ymom[0:n[j]-1,k,j], '-', c = cstr[j])
627                        axis(ymom_axis)
628                        units = 'm^2 / sec'
629                    if which_quantity == 'speed':
630                        plot(model_time[0:n[j]-1,k,j],
631                             speed[0:n[j]-1,k,j], '-', c = cstr[j])
632                        axis(vel_axis)
633                        units = 'm / sec'
634                    if which_quantity == 'bearing':
635                        plot(model_time[0:n[j]-1,k,j],bearings[0:n[j]-1,k,j],'-',
636                             model_time[0:n[j]-1,k,j], due_west[0:n[j]-1], '-.', 
637                             model_time[0:n[j]-1,k,j], due_east[0:n[j]-1], '-.')
638                        units = 'degrees from North'
639                        #ax = axis([time_min, time_max, 0.0, 360.0])
640                        legend(('Bearing','West','East'))
641
642                    if time_unit is 'mins': xlabel('time (mins)')
643                    if time_unit is 'hours': xlabel('time (hours)')
644                    #if which_quantity == 'stage' \
645                    #   and elevations[0:n[j]-1,k,j] > 0:
646                    #    ylabel('%s (%s)' %('depth', units))
647                    #else:
648                    #    ylabel('%s (%s)' %(which_quantity, units))
649                        #ylabel('%s (%s)' %('wave height', units))
650                    ylabel('%s (%s)' %(which_quantity, units))
651                    if len(label_id) > 1: legend((leg_label),loc='upper right')
652
653                    #gaugeloc1 = gaugeloc.replace(' ','')
654                    #gaugeloc2 = gaugeloc1.replace('_','')
655                    gaugeloc2 = str(locations[k]).replace(' ','')
656                    graphname = '%sgauge%s_%s' %(file_loc[j],
657                                                 gaugeloc2,
658                                                 which_quantity)
659
660                    if report == True and len(label_id) > 1:
661                        figdir = getcwd()+sep+'report_figures'+sep
662                        if access(figdir,F_OK) == 0 :
663                            mkdir (figdir)
664                        latex_file_loc = figdir.replace(sep,altsep) 
665                        # storing files in production directory   
666                        graphname_latex = '%sgauge%s%s' \
667                                          % (latex_file_loc, gaugeloc2,
668                                             which_quantity)
669                        # giving location in latex output file
670                        graphname_report_input = '%sgauge%s%s' % \
671                                                 ('..' + altsep + 
672                                                      'report_figures' + altsep,
673                                                  gaugeloc2, which_quantity)
674                        graphname_report.append(graphname_report_input)
675                       
676                        # save figures in production directory for report
677                        savefig(graphname_latex)
678
679                    if report == True:
680                        figdir = getcwd() + sep + 'report_figures' + sep
681                        if access(figdir,F_OK) == 0:
682                            mkdir(figdir)
683                        latex_file_loc = figdir.replace(sep,altsep)   
684
685                        if len(label_id) == 1: 
686                            # storing files in production directory 
687                            graphname_latex = '%sgauge%s%s%s' % \
688                                              (latex_file_loc, gaugeloc2,
689                                               which_quantity, label_id2)
690                            # giving location in latex output file
691                            graphname_report = '%sgauge%s%s%s' % \
692                                               ('..' + altsep +
693                                                    'report_figures' + altsep,
694                                                gaugeloc2, which_quantity,
695                                                label_id2)
696                            s = '\includegraphics' \
697                                '[width=0.49\linewidth, height=50mm]{%s%s}' % \
698                                (graphname_report, '.png')
699                            fid.write(s)
700                            if where1 % 2 == 0:
701                                s = '\\\\ \n'
702                                where1 = 0
703                            else:
704                                s = '& \n'
705                            fid.write(s)
706                            savefig(graphname_latex)
707                   
708                    if title_on == True:
709                        title('%s scenario: %s at %s gauge' % \
710                              (label_id, which_quantity, gaugeloc2))
711                        #title('Gauge %s (MOST elevation %.2f, ' \
712                        #      'ANUGA elevation %.2f)' % \
713                        #      (gaugeloc2, elevations[10,k,0],
714                        #       elevations[10,k,1]))
715
716                    savefig(graphname) # save figures with sww file
717
718                if report == True and len(label_id) == 1:
719                    for i in range(nn-1):
720                        if nn > 2:
721                            if plot_quantity[i] == 'stage' \
722                               and elevations[0,k,j] > 0:
723                                word_quantity += 'depth' + ', '
724                            else:
725                                word_quantity += plot_quantity[i] + ', '
726                        else:
727                            if plot_quantity[i] == 'stage' \
728                               and elevations[0,k,j] > 0:
729                                word_quantity += 'depth' + ', '
730                            else:
731                                word_quantity += plot_quantity[i]
732                       
733                    if plot_quantity[nn-1] == 'stage' and elevations[0,k,j] > 0:
734                        word_quantity += ' and ' + 'depth'
735                    else:
736                        word_quantity += ' and ' + plot_quantity[nn-1]
737                    caption = 'Time series for %s at %s location ' \
738                              '(elevation %.2fm)' % \
739                              (word_quantity, locations[k], elev[k])
740                    if elev[k] == 0.0:
741                        caption = 'Time series for %s at %s location ' \
742                                  '(elevation %.2fm)' % \
743                                  (word_quantity, locations[k],
744                                   elevations[0,k,j])
745                        east = gauges[0]
746                        north = gauges[1]
747                        elev_output.append([locations[k], east, north,
748                                            elevations[0,k,j]])
749                    label = '%sgauge%s' % (label_id2, gaugeloc2)
750                    s = '\end{tabular} \n' \
751                        '\\caption{%s} \n' \
752                        '\label{fig:%s} \n' \
753                        '\end{figure} \n \n' % (caption, label)
754                    fid.write(s)
755                    cc += 1
756                    if cc % 6 == 0: fid.write('\\clearpage \n')
757                    savefig(graphname_latex)               
758                   
759            if report == True and len(label_id) > 1:
760                for i in range(nn-1):
761                    if nn > 2:
762                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
763                            word_quantity += 'depth' + ','
764                        else:
765                            word_quantity += plot_quantity[i] + ', '
766                    else:
767                        if plot_quantity[i] == 'stage' and elevations[0,k,j] > 0:
768                            word_quantity += 'depth'
769                        else:
770                            word_quantity += plot_quantity[i]
771                    where1 = 0
772                    count1 += 1
773                    index = j*len(plot_quantity)
774                    for which_quantity in plot_quantity:
775                        where1 += 1
776                        s = '\includegraphics' \
777                            '[width=0.49\linewidth, height=50mm]{%s%s}' % \
778                            (graphname_report[index], '.png')
779                        index += 1
780                        fid.write(s)
781                        if where1 % 2 == 0:
782                            s = '\\\\ \n'
783                            where1 = 0
784                        else:
785                            s = '& \n'
786                        fid.write(s)
787                word_quantity += ' and ' + plot_quantity[nn-1]           
788                label = 'gauge%s' %(gaugeloc2) 
789                caption = 'Time series for %s at %s location ' \
790                          '(elevation %.2fm)' % \
791                          (word_quantity, locations[k], elev[k])
792                if elev[k] == 0.0:
793                        caption = 'Time series for %s at %s location ' \
794                                  '(elevation %.2fm)' % \
795                                  (word_quantity, locations[k],
796                                   elevations[0,k,j])
797                        thisgauge = gauges[k]
798                        east = thisgauge[0]
799                        north = thisgauge[1]
800                        elev_output.append([locations[k], east, north,
801                                            elevations[0,k,j]])
802                       
803                s = '\end{tabular} \n' \
804                    '\\caption{%s} \n' \
805                    '\label{fig:%s} \n' \
806                    '\end{figure} \n \n' % (caption, label)
807                fid.write(s)
808                if float((k+1)/div - pp) == 0.:
809                    fid.write('\\clearpage \n')
810                    pp += 1
811                #### finished generating figures ###
812
813            close('all')
814       
815    return texfile2, elev_output
816
817
818# FIXME (DSG): Add unit test, make general, not just 2 files,
819# but any number of files.
820##
821# @brief ??
822# @param dir_name ??
823# @param filename1 ??
824# @param filename2 ??
825# @return ??
826# @note TEMP
827def copy_code_files(dir_name, filename1, filename2):
828    """Temporary Interface to new location"""
829
830    from anuga.shallow_water.data_manager import \
831                    copy_code_files as dm_copy_code_files
832    log.critical('copy_code_files has moved from util.py.')
833    log.critical('Please use "from anuga.shallow_water.data_manager import '
834                 'copy_code_files"')
835   
836    return dm_copy_code_files(dir_name, filename1, filename2)
837
838
839##
840# @brief Create a nested sub-directory path.
841# @param root_directory The base diretory path.
842# @param directories An iterable of sub-directory names.
843# @return The final joined directory path.
844# @note If each sub-directory doesn't exist, it will be created.
845def add_directories(root_directory, directories):
846    """
847    Add the first sub-directory in 'directories' to root_directory.
848    Then add the second sub-directory to the accumulating path and so on.
849
850    Return the path of the final directory.
851
852    This is handy for specifying and creating a directory where data will go.
853    """
854    dir = root_directory
855    for new_dir in directories:
856        dir = os.path.join(dir, new_dir)
857        if not access(dir,F_OK):
858            mkdir(dir)
859    return dir
860
861
862##
863# @brief
864# @param filename
865# @param separator_value
866# @return
867# @note TEMP
868def get_data_from_file(filename, separator_value=','):
869    """Temporary Interface to new location"""
870    from anuga.shallow_water.data_manager import \
871                        get_data_from_file as dm_get_data_from_file
872    log.critical('get_data_from_file has moved from util.py')
873    log.critical('Please use "from anuga.shallow_water.data_manager import '
874                 'get_data_from_file"')
875   
876    return dm_get_data_from_file(filename,separator_value = ',')
877
878
879##
880# @brief
881# @param verbose
882# @param kwargs
883# @return
884# @note TEMP
885def store_parameters(verbose=False,**kwargs):
886    """Temporary Interface to new location"""
887   
888    from anuga.shallow_water.data_manager \
889                    import store_parameters as dm_store_parameters
890    log.critical('store_parameters has moved from util.py.')
891    log.critical('Please use "from anuga.shallow_water.data_manager '
892                 'import store_parameters"')
893   
894    return dm_store_parameters(verbose=False,**kwargs)
895
896
897##
898# @brief Remove vertices that are not associated with any triangle.
899# @param verts An iterable (or array) of points.
900# @param triangles An iterable of 3 element tuples.
901# @param number_of_full_nodes ??
902# @return (verts, triangles) where 'verts' has been updated.
903def remove_lone_verts(verts, triangles, number_of_full_nodes=None):
904    """Removes vertices that are not associated with any triangles.
905
906    verts is a list/array of points.
907    triangles is a list of 3 element tuples.  Each tuple represents a triangle.
908    number_of_full_nodes relate to parallelism when a mesh has an
909        extra layer of ghost points.
910    """
911
912    verts = ensure_numeric(verts)
913    triangles = ensure_numeric(triangles)
914   
915    N = len(verts)
916   
917    # initialise the array to easily find the index of the first loner
918    # ie, if N=3 -> [6,5,4]
919    loners=num.arange(2*N, N, -1)
920    for t in triangles:
921        for vert in t:
922            loners[vert]= vert # all non-loners will have loners[i]=i
923
924    lone_start = 2*N - max(loners) # The index of the first loner
925
926    if lone_start-1 == N:
927        # no loners
928        pass
929    elif min(loners[lone_start:N]) > N:
930        # All the loners are at the end of the vert array
931        verts = verts[0:lone_start]
932    else:
933        # change the loners list so it can be used to modify triangles
934        # Remove the loners from verts
935        # Could've used X=compress(less(loners,N),loners)
936        # verts=num.take(verts,X,axis=0)  to Remove the loners from verts
937        # but I think it would use more memory
938        new_i = lone_start    # point at first loner - 'shuffle down' target
939        for i in range(lone_start, N):
940            if loners[i] >= N:    # [i] is a loner, leave alone
941                pass
942            else:        # a non-loner, move down
943                loners[i] = new_i
944                verts[new_i] = verts[i]
945                new_i += 1
946        verts = verts[0:new_i]
947
948        # Modify the triangles
949        triangles = num.choose(triangles,loners)
950    return verts, triangles
951
952
953##
954# @brief Compute centroid values from vertex values
955# @param x Values at vertices of triangular mesh.
956# @param triangles Nx3 integer array pointing to vertex information.
957# @return [N] array of centroid values.
958def get_centroid_values(x, triangles):
959    """Compute centroid values from vertex values
960   
961    x: Values at vertices of triangular mesh
962    triangles: Nx3 integer array pointing to vertex information
963    for each of the N triangels. Elements of triangles are
964    indices into x
965    """
966       
967    xc = num.zeros(triangles.shape[0], num.float) # Space for centroid info
968   
969    for k in range(triangles.shape[0]):
970        # Indices of vertices
971        i0 = triangles[k][0]
972        i1 = triangles[k][1]
973        i2 = triangles[k][2]       
974       
975        xc[k] = (x[i0] + x[i1] + x[i2])/3
976
977    return xc
978
979
980# @note TEMP
981def make_plots_from_csv_file(directories_dic={dir:['gauge', 0, 0]},
982                                output_dir='',
983                                base_name='',
984                                plot_numbers=['3-5'],
985                                quantities=['speed','stage','momentum'],
986                                assess_all_csv_files=True,
987                                extra_plot_name='test'):
988
989    msg = 'make_plots_from_csv_file has been replaced by csv2timeseries_graphs '
990    msg += 'Please use "from anuga.abstract_2d_finite_volumes.util import ' \
991           'csv2timeseries_graphs"'
992    raise Exception, msg
993
994    return csv2timeseries_graphs(directories_dic,
995                                 output_dir,
996                                 base_name,
997                                 plot_numbers,
998                                 quantities,
999                                 extra_plot_name,
1000                                 assess_all_csv_files)
1001
1002
1003##
1004# @brief Plot time series from CSV files.
1005# @param directories_dic
1006# @param output_dir
1007# @param base_name
1008# @param plot_numbers
1009# @param quantities
1010# @param extra_plot_name
1011# @param assess_all_csv_files
1012# @param create_latex
1013# @param verbose
1014# @note Assumes that 'elevation' is in the CSV file(s).
1015def csv2timeseries_graphs(directories_dic={},
1016                          output_dir='',
1017                          base_name=None,
1018                          plot_numbers='',
1019                          quantities=['stage'],
1020                          extra_plot_name='',
1021                          assess_all_csv_files=True,
1022                          create_latex=False,
1023                          verbose=False):
1024                               
1025    """
1026    Read in csv files that have the right header information and
1027    plot time series such as Stage, Speed, etc. Will also plot several
1028    time series on one plot. Filenames must follow this convention,
1029    <base_name><plot_number>.csv eg gauge_timeseries3.csv
1030   
1031    NOTE: relies that 'elevation' is in the csv file!
1032
1033    Each file represents a location and within each file there are
1034    time, quantity columns.
1035   
1036    For example:   
1037    if "directories_dic" defines 4 directories and in each directories
1038    there is a csv files corresponding to the right "plot_numbers",
1039    this will create a plot with 4 lines one for each directory AND
1040    one plot for each "quantities".  ??? FIXME: unclear.
1041   
1042    Usage:
1043        csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1044                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1045                            output_dir='fixed_wave'+sep,
1046                            base_name='gauge_timeseries_',
1047                            plot_numbers='',
1048                            quantities=['stage','speed'],
1049                            extra_plot_name='',
1050                            assess_all_csv_files=True,                           
1051                            create_latex=False,
1052                            verbose=True)
1053            this will create one plot for stage with both 'slide' and
1054            'fixed_wave' lines on it for stage and speed for each csv
1055            file with 'gauge_timeseries_' as the prefix. The graghs
1056            will be in the output directory 'fixed_wave' and the graph
1057            axis will be determined by assessing all the
1058   
1059    ANOTHER EXAMPLE
1060        new_csv2timeseries_graphs(directories_dic={'slide'+sep:['Slide',0, 0],
1061                                       'fixed_wave'+sep:['Fixed Wave',0,0]},
1062                            output_dir='fixed_wave'+sep,
1063                            base_name='gauge_timeseries_',
1064                            plot_numbers=['1-3'],
1065                            quantities=['stage','speed'],
1066                            extra_plot_name='',
1067                            assess_all_csv_files=False,                           
1068                            create_latex=False,
1069                            verbose=True)
1070        This will plot csv files called gauge_timeseries_1.csv and
1071        gauge_timeseries3.csv from both 'slide' and 'fixed_wave' directories
1072        to 'fixed_wave'. There will be 4 plots created two speed and two stage
1073        one for each csv file. There will be two lines on each of these plots.
1074        And the axis will have been determined from only these files, had
1075        assess_all_csv_files = True all csv file with 'gauges_timeseries_' prefix
1076        would of been assessed.
1077   
1078    ANOTHER EXAMPLE   
1079         csv2timeseries_graphs({'J:'+sep+'anuga_validation'+sep:['new',20,-.1],
1080                                   'J:'+sep+'conical_island'+sep:['test',0,0]},
1081                                   output_dir='',
1082                                   plot_numbers=['1','3'],
1083                                   quantities=['stage','depth','bearing'],
1084                                   base_name='gauge_b',
1085                                   assess_all_csv_files=True,
1086                                  verbose=True)   
1087       
1088            This will produce one plot for each quantity (therefore 3) in the
1089            current directory, each plot will have 2 lines on them. The first
1090            plot named 'new' will have the time offseted by 20secs and the stage
1091            height adjusted by -0.1m
1092       
1093    Inputs:
1094        directories_dic: dictionary of directory with values (plot
1095                         legend name for directory), (start time of
1096                         the time series) and the (value to add to
1097                         stage if needed). For example
1098                         {dir1:['Anuga_ons',5000, 0],
1099                          dir2:['b_emoth',5000,1.5],
1100                          dir3:['b_ons',5000,1.5]}
1101                         Having multiple directories defined will plot them on
1102                         one plot, therefore there will be 3 lines on each of
1103                         these plot. If you only want one line per plot call
1104                         csv2timeseries_graph separately for each directory,
1105                         eg only have one directory in the 'directories_dic' in
1106                         each call.
1107                         
1108        output_dir: directory for the plot outputs. Only important to define when
1109                    you have more than one directory in your directories_dic, if
1110                    you have not defined it and you have multiple directories in
1111                    'directories_dic' there will be plots in each directory,
1112                    however only one directory will contain the complete
1113                    plot/graphs.
1114       
1115        base_name: Is used a couple of times.
1116                   1) to find the csv files to be plotted if there is no
1117                      'plot_numbers' then csv files with 'base_name' are plotted
1118                   2) in the title of the plots, the length of base_name is
1119                      removed from the front of the filename to be used in the
1120                      title.
1121                   This could be changed if needed.
1122                   Note is ignored if assess_all_csv_files=True
1123       
1124        plot_numbers: a String list of numbers to plot. For example
1125                      [0-4,10,15-17] will read and attempt to plot
1126                      the follow 0,1,2,3,4,10,15,16,17
1127                      NOTE: if no plot numbers this will create one plot per
1128                            quantity, per gauge
1129
1130        quantities: Will get available quantities from the header in the csv
1131                    file.  Quantities must be one of these.
1132                    NOTE: ALL QUANTITY NAMES MUST BE lower case!
1133                   
1134        extra_plot_name: A string that is appended to the end of the
1135                         output filename.
1136                   
1137        assess_all_csv_files: if true it will read ALL csv file with
1138                             "base_name", regardless of 'plot_numbers'
1139                              and determine a uniform set of axes for
1140                              Stage, Speed and Momentum. IF FALSE it
1141                              will only read the csv file within the
1142                             'plot_numbers'
1143                             
1144        create_latex: NOT IMPLEMENTED YET!! sorry Jane....
1145       
1146    OUTPUTS: saves the plots to
1147              <output_dir><base_name><plot_number><extra_plot_name>.png
1148    """
1149
1150    try: 
1151        import pylab
1152    except ImportError:
1153        msg='csv2timeseries_graphs needs pylab to be installed correctly'
1154        raise msg
1155            #ANUGA don't need pylab to work so the system doesn't
1156            #rely on pylab being installed
1157        return
1158
1159    from os import sep
1160    from anuga.shallow_water.data_manager import \
1161                               get_all_files_with_extension, csv2dict
1162
1163    seconds_in_hour = 3600
1164    seconds_in_minutes = 60
1165   
1166    quantities_label={}
1167#    quantities_label['time'] = 'time (hours)'
1168    quantities_label['time'] = 'time (minutes)'
1169    quantities_label['stage'] = 'wave height (m)'
1170    quantities_label['speed'] = 'speed (m/s)'
1171    quantities_label['momentum'] = 'momentum (m^2/sec)'
1172    quantities_label['depth'] = 'water depth (m)'
1173    quantities_label['xmomentum'] = 'momentum (m^2/sec)'
1174    quantities_label['ymomentum'] = 'momentum (m^2/sec)'
1175    quantities_label['bearing'] = 'degrees (o)'
1176    quantities_label['elevation'] = 'elevation (m)'
1177   
1178    if extra_plot_name != '':
1179        extra_plot_name = '_' + extra_plot_name
1180
1181    new_plot_numbers=[]
1182    #change plot_numbers to list, eg ['0-4','10']
1183    #to ['0','1','2','3','4','10']
1184    for i, num_string in enumerate(plot_numbers):
1185        if '-' in num_string: 
1186            start = int(num_string[:num_string.rfind('-')])
1187            end = int(num_string[num_string.rfind('-') + 1:]) + 1
1188            for x in range(start, end):
1189                new_plot_numbers.append(str(x))
1190        else:
1191            new_plot_numbers.append(num_string)
1192
1193    #finds all the files that fit the specs provided and return a list of them
1194    #so to help find a uniform max and min for the plots...
1195    list_filenames=[]
1196    all_csv_filenames=[]
1197    if verbose: log.critical('Determining files to access for axes ranges.')
1198   
1199    for i,directory in enumerate(directories_dic.keys()):
1200        all_csv_filenames.append(get_all_files_with_extension(directory,
1201                                                              base_name, '.csv'))
1202
1203        filenames=[]
1204        if plot_numbers == '': 
1205            list_filenames.append(get_all_files_with_extension(directory,
1206                                                               base_name,'.csv'))
1207        else:
1208            for number in new_plot_numbers:
1209                filenames.append(base_name + number)
1210            list_filenames.append(filenames)
1211
1212    #use all the files to get the values for the plot axis
1213    max_start_time= -1000.
1214    min_start_time = 100000 
1215   
1216    if verbose: log.critical('Determining uniform axes')
1217
1218    #this entire loop is to determine the min and max range for the
1219    #axes of the plots
1220
1221#    quantities.insert(0,'elevation')
1222    quantities.insert(0,'time')
1223
1224    directory_quantity_value={}
1225#    quantity_value={}
1226    min_quantity_value={}
1227    max_quantity_value={}
1228
1229    for i, directory in enumerate(directories_dic.keys()):
1230        filename_quantity_value = {}
1231        if assess_all_csv_files == False:
1232            which_csv_to_assess = list_filenames[i]
1233        else:
1234            #gets list of filenames for directory "i"
1235            which_csv_to_assess = all_csv_filenames[i]
1236       
1237        for j, filename in enumerate(which_csv_to_assess):
1238            quantity_value = {}
1239
1240            dir_filename = join(directory,filename)
1241            attribute_dic, title_index_dic = csv2dict(dir_filename + '.csv')
1242            directory_start_time = directories_dic[directory][1]
1243            directory_add_tide = directories_dic[directory][2]
1244
1245            if verbose: log.critical('reading: %s.csv' % dir_filename)
1246
1247            #add time to get values
1248            for k, quantity in enumerate(quantities):
1249                quantity_value[quantity] = [float(x) for
1250                                                x in attribute_dic[quantity]]
1251
1252                #add tide to stage if provided
1253                if quantity == 'stage':
1254                     quantity_value[quantity] = num.array(quantity_value[quantity],
1255                                                          num.float) + directory_add_tide
1256
1257                #condition to find max and mins for all the plots
1258                # populate the list with something when i=0 and j=0 and
1259                # then compare to the other values to determine abs max and min
1260                if i==0 and j==0:
1261                    min_quantity_value[quantity], \
1262                        max_quantity_value[quantity] = \
1263                            get_min_max_values(quantity_value[quantity])
1264
1265                    if quantity != 'time':
1266                        min_quantity_value[quantity] = \
1267                            min_quantity_value[quantity] *1.1
1268                        max_quantity_value[quantity] = \
1269                            max_quantity_value[quantity] *1.1
1270                else:
1271                    min, max = get_min_max_values(quantity_value[quantity])
1272               
1273                    # min and max are multipled by "1+increase_axis" to get axes
1274                    # that are slighty bigger than the max and mins
1275                    # so the plots look good.
1276
1277                    increase_axis = (max-min)*0.05
1278                    if min <= min_quantity_value[quantity]:
1279                        if quantity == 'time': 
1280                            min_quantity_value[quantity] = min
1281                        else:
1282                            if round(min,2) == 0.00:
1283                                min_quantity_value[quantity] = -increase_axis
1284#                                min_quantity_value[quantity] = -2.
1285                                #min_quantity_value[quantity] = \
1286                                #    -max_quantity_value[quantity]*increase_axis
1287                            else:
1288#                                min_quantity_value[quantity] = \
1289#                                    min*(1+increase_axis)
1290                                min_quantity_value[quantity]=min-increase_axis
1291                   
1292                    if max > max_quantity_value[quantity]: 
1293                        if quantity == 'time': 
1294                            max_quantity_value[quantity] = max
1295                        else:
1296                            max_quantity_value[quantity] = max + increase_axis
1297#                            max_quantity_value[quantity]=max*(1+increase_axis)
1298
1299            #set the time... ???
1300            if min_start_time > directory_start_time: 
1301                min_start_time = directory_start_time
1302            if max_start_time < directory_start_time: 
1303                max_start_time = directory_start_time
1304           
1305            filename_quantity_value[filename]=quantity_value
1306           
1307        directory_quantity_value[directory]=filename_quantity_value
1308   
1309    #final step to unifrom axis for the graphs
1310    quantities_axis={}
1311   
1312    for i, quantity in enumerate(quantities):
1313        quantities_axis[quantity] = (float(min_start_time) \
1314                                         / float(seconds_in_minutes),
1315                                     (float(max_quantity_value['time']) \
1316                                          + float(max_start_time)) \
1317                                              / float(seconds_in_minutes),
1318                                     min_quantity_value[quantity],
1319                                     max_quantity_value[quantity])
1320
1321        if verbose and (quantity != 'time' and quantity != 'elevation'): 
1322            log.critical('axis for quantity %s are x:(%s to %s)%s '
1323                         'and y:(%s to %s)%s' 
1324                         % (quantity, quantities_axis[quantity][0],
1325                            quantities_axis[quantity][1],
1326                            quantities_label['time'],
1327                            quantities_axis[quantity][2],
1328                            quantities_axis[quantity][3],
1329                            quantities_label[quantity]))
1330
1331    cstr = ['b', 'r', 'g', 'c', 'm', 'y', 'k']
1332
1333    if verbose: log.critical('Now start to plot')
1334   
1335    i_max = len(directories_dic.keys())
1336    legend_list_dic = {}
1337    legend_list = []
1338    for i, directory in enumerate(directories_dic.keys()):
1339        if verbose: log.critical('Plotting in %s %s'
1340                                 % (directory, new_plot_numbers))
1341
1342        # FIXME THIS SORT IS VERY IMPORTANT
1343        # Without it the assigned plot numbers may not work correctly
1344        # there must be a better way
1345        list_filenames[i].sort()
1346        for j, filename in enumerate(list_filenames[i]):
1347            if verbose: log.critical('Starting %s' % filename)
1348
1349            directory_name = directories_dic[directory][0]
1350            directory_start_time = directories_dic[directory][1]
1351            directory_add_tide = directories_dic[directory][2]
1352           
1353            # create an if about the start time and tide height if don't exist
1354            attribute_dic, title_index_dic = csv2dict(directory + sep
1355                                                      + filename + '.csv')
1356            #get data from dict in to list
1357            #do maths to list by changing to array
1358            t = (num.array(directory_quantity_value[directory][filename]['time'])
1359                     + directory_start_time) / seconds_in_minutes
1360
1361            #finds the maximum elevation, used only as a test
1362            # and as info in the graphs
1363            max_ele=-100000
1364            min_ele=100000
1365            elevation = [float(x) for x in attribute_dic["elevation"]]
1366           
1367            min_ele, max_ele = get_min_max_values(elevation)
1368           
1369            if min_ele != max_ele:
1370                log.critical("Note! Elevation changes in %s" % dir_filename)
1371
1372            # creates a dictionary with keys that is the filename and attributes
1373            # are a list of lists containing 'directory_name' and 'elevation'.
1374            # This is used to make the contents for the legends in the graphs,
1375            # this is the name of the model and the elevation.  All in this
1376            # great one liner from DG. If the key 'filename' doesn't exist it
1377            # creates the entry if the entry exist it appends to the key.
1378
1379            legend_list_dic.setdefault(filename,[]) \
1380                .append([directory_name, round(max_ele, 3)])
1381
1382            # creates a LIST for the legend on the last iteration of the
1383            # directories which is when "legend_list_dic" has been fully
1384            # populated. Creates a list of strings which is used in the legend
1385            # only runs on the last iteration for all the gauges(csv) files
1386            # empties the list before creating it
1387
1388            if i == i_max - 1:
1389                legend_list = []
1390   
1391                for name_and_elevation in legend_list_dic[filename]:
1392                    legend_list.append('%s (elevation = %sm)'\
1393                                       % (name_and_elevation[0],
1394                                          name_and_elevation[1]))
1395           
1396            #skip time and elevation so it is not plotted!
1397            for k, quantity in enumerate(quantities):
1398                if quantity != 'time' and quantity != 'elevation':
1399                    pylab.figure(int(k*100+j))
1400                    pylab.ylabel(quantities_label[quantity])
1401                    pylab.plot(t,
1402                               directory_quantity_value[directory]\
1403                                                       [filename][quantity],
1404                               c = cstr[i], linewidth=1)
1405                    pylab.xlabel(quantities_label['time'])
1406                    pylab.axis(quantities_axis[quantity])
1407                    pylab.legend(legend_list,loc='upper right')
1408                   
1409                    pylab.title('%s at %s gauge'
1410                                % (quantity, filename[len(base_name):]))
1411
1412                    if output_dir == '':
1413                        figname = '%s%s%s_%s%s.png' \
1414                                  % (directory, sep, filename, quantity,
1415                                     extra_plot_name)
1416                    else:
1417                        figname = '%s%s%s_%s%s.png' \
1418                                  % (output_dir, sep, filename, quantity,
1419                                     extra_plot_name)
1420
1421                    if verbose: log.critical('saving figure here %s' % figname)
1422
1423                    pylab.savefig(figname)
1424           
1425    if verbose: log.critical('Closing all plots')
1426
1427    pylab.close('all')
1428    del pylab
1429
1430    if verbose: log.critical('Finished closing plots')
1431
1432##
1433# @brief Return min and max of an iterable.
1434# @param list The iterable to return min & max of.
1435# @return (min, max) of 'list'.
1436def get_min_max_values(list=None):
1437    """
1438    Returns the min and max of the list it was provided.
1439    """
1440
1441    if list == None: log.critical('List must be provided')
1442       
1443    return min(list), max(list)
1444
1445
1446##
1447# @brief Get runup around a point in a CSV file.
1448# @param gauge_filename gauge file name.
1449# @param sww_filename SWW file name.
1450# @param runup_filename Name of file to report into.
1451# @param size ??
1452# @param verbose ??
1453def get_runup_data_for_locations_from_file(gauge_filename,
1454                                           sww_filename,
1455                                           runup_filename,
1456                                           size=10,
1457                                           verbose=False):
1458    """this will read a csv file with the header x,y. Then look in a square
1459    'size'x2 around this position for the 'max_inundaiton_height' in the
1460    'sww_filename' and report the findings in the 'runup_filename'.
1461   
1462    WARNING: NO TESTS!
1463    """
1464
1465    from anuga.shallow_water.data_manager import get_all_directories_with_name,\
1466                                                 get_maximum_inundation_data,\
1467                                                 csv2dict
1468                                                 
1469    file = open(runup_filename, "w")
1470    file.write("easting,northing,runup \n ")
1471    file.close()
1472   
1473    #read gauge csv file to dictionary
1474    attribute_dic, title_index_dic = csv2dict(gauge_filename)
1475    northing = [float(x) for x in attribute_dic["y"]]
1476    easting = [float(x) for x in attribute_dic["x"]]
1477
1478    log.critical('Reading %s' % sww_filename)
1479
1480    runup_locations=[]
1481    for i, x in enumerate(northing):
1482        poly = [[int(easting[i]+size),int(northing[i]+size)],
1483                [int(easting[i]+size),int(northing[i]-size)],
1484                [int(easting[i]-size),int(northing[i]-size)],
1485                [int(easting[i]-size),int(northing[i]+size)]]
1486       
1487        run_up, x_y = get_maximum_inundation_data(filename=sww_filename,
1488                                                  polygon=poly,
1489                                                  verbose=False) 
1490
1491        #if no runup will return 0 instead of NONE
1492        if run_up==None: run_up=0
1493        if x_y==None: x_y=[0,0]
1494       
1495        if verbose:
1496            log.critical('maximum inundation runup near %s is %s meters'
1497                         % (x_y, run_up))
1498       
1499        #writes to file
1500        file = open(runup_filename, "a")
1501        temp = '%s,%s,%s \n' % (x_y[0], x_y[1], run_up)
1502        file.write(temp)
1503        file.close()
1504
1505##
1506# @brief ??
1507# @param  ??
1508# @param gauge_file ??
1509# @param out_name ??
1510# @param quantities ??
1511# @param verbose ??
1512# @param use_cache ??
1513def sww2csv_gauges(sww_file,
1514                   gauge_file,
1515                   out_name='gauge_',
1516                   quantities=['stage', 'depth', 'elevation',
1517                               'xmomentum', 'ymomentum'],
1518                   verbose=False,
1519                   use_cache=True):
1520    """
1521   
1522    Inputs:
1523        NOTE: if using csv2timeseries_graphs after creating csv file,
1524        it is essential to export quantities 'depth' and 'elevation'.
1525        'depth' is good to analyse gauges on land and elevation is used
1526        automatically by csv2timeseries_graphs in the legend.
1527       
1528        sww_file: path to any sww file
1529       
1530        gauge_file: Assumes that it follows this format
1531            name, easting, northing, elevation
1532            point1, 100.3, 50.2, 10.0
1533            point2, 10.3, 70.3, 78.0
1534       
1535        NOTE: order of column can change but names eg 'easting', 'elevation'
1536        must be the same! ALL lowercaps!
1537
1538        out_name: prefix for output file name (default is 'gauge_')
1539       
1540    Outputs:
1541        one file for each gauge/point location in the points file. They
1542        will be named with this format in the same directory as the 'sww_file'
1543            <out_name><name>.csv
1544        eg gauge_point1.csv if <out_name> not supplied
1545           myfile_2_point1.csv if <out_name> ='myfile_2_'
1546           
1547        They will all have a header
1548   
1549    Usage: sww2csv_gauges(sww_file='test1.sww',
1550                          quantities = ['stage', 'elevation','depth','bearing'],
1551                          gauge_file='gauge.txt')   
1552   
1553    Interpolate the quantities at a given set of locations, given
1554    an sww file.
1555    The results are written to a csv file.
1556
1557    In the future let points be a points file.
1558    And the user choose the quantities.
1559
1560    This is currently quite specific.
1561    If it needs to be more general, change things.
1562
1563    This is really returning speed, not velocity.
1564    """
1565    from gauge import sww2csv_gauges as sww2csv
1566   
1567    return sww2csv(sww_file, gauge_file, out_name, quantities, verbose, use_cache)
1568   
1569def sww2timeseries(swwfiles,
1570                   gauge_filename,
1571                   production_dirs,
1572                   report=None,
1573                   reportname=None,
1574                   plot_quantity=None,
1575                   generate_fig=False,
1576                   surface=None,
1577                   time_min=None,
1578                   time_max=None,
1579                   time_thinning=1,                   
1580                   time_unit=None,
1581                   title_on=None,
1582                   use_cache=False,
1583                   verbose=False,
1584                   output_centroids=False):
1585    from gauge import sww2timeseries as sww2timeseries_new
1586    return sww2timeseries_new(swwfiles,
1587                   gauge_filename,
1588                   production_dirs,
1589                   report,
1590                   reportname,
1591                   plot_quantity,
1592                   generate_fig,
1593                   surface,
1594                   time_min,
1595                   time_max,
1596                   time_thinning,                   
1597                   time_unit,
1598                   title_on,
1599                   use_cache,
1600                   verbose,
1601                   output_centroids)                   
1602   
1603##
1604# @brief Get a wave height at a certain depth given wave height at another depth.
1605# @param d1 The first depth.
1606# @param d2 The second depth.
1607# @param h1 Wave ampitude at d1
1608# @param verbose True if this function is to be verbose.
1609# @return The wave height at d2.
1610def greens_law(d1, d2, h1, verbose=False):
1611    """Green's Law
1612
1613    Green's Law allows an approximation of wave amplitude at
1614    a given depth based on the fourh root of the ratio of two depths
1615    and the amplitude at another given depth.
1616
1617    Note, wave amplitude is equal to stage.
1618   
1619    Inputs:
1620
1621    d1, d2 - the two depths
1622    h1     - the wave amplitude at d1
1623    h2     - the derived amplitude at d2
1624
1625    h2 = h1 (d1/d2)^(1/4), where d2 cannot equal 0.
1626   
1627    """
1628
1629    d1 = ensure_numeric(d1)
1630    d2 = ensure_numeric(d2)
1631    h1 = ensure_numeric(h1)
1632
1633    if d1 <= 0.0:
1634        msg = 'the first depth, d1 (%f), must be strictly positive' % (d1)
1635        raise Exception(msg)
1636
1637    if d2 <= 0.0:
1638        msg = 'the second depth, d2 (%f), must be strictly positive' % (d2)
1639        raise Exception(msg)
1640   
1641    if h1 <= 0.0:
1642        msg = 'the wave amplitude, h1 (%f), must be strictly positive' % (h1)
1643        raise Exception(msg)
1644   
1645    h2 = h1*(d1/d2)**0.25
1646
1647    assert h2 > 0
1648   
1649    return h2
1650       
1651
1652##
1653# @brief Get the square-root of a value.
1654# @param s The value to get the square-root of.
1655# @return The square-root of 's'.
1656def square_root(s):
1657    return sqrt(s)
1658
1659
Note: See TracBrowser for help on using the repository browser.