source: trunk/anuga_core/source/anuga/utilities/sww_merge.py @ 7820

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

Added unit test for sww_merge.

File size: 3.5 KB
Line 
1"""
2    Merge a list of .sww files together into a single file.
3"""
4
5import numpy as num
6from anuga.utilities.numerical_tools import ensure_numeric
7
8from Scientific.IO.NetCDF import NetCDFFile
9from anuga.config import netcdf_mode_r, netcdf_mode_w, \
10                            netcdf_mode_a, netcdf_float
11from anuga.file.sww import SWW_file, Write_sww
12
13
14def sww_merge(swwfiles, output, verbose = False):
15    """
16        Merge a list of sww files into a single file.
17       
18        May be useful for parallel runs. Note that some advanced information
19        and custom quantities may not be exported.
20    """
21   
22    static_quantities = ['elevation']
23    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
24   
25    first_file = True
26    tri_offset = 0
27    for filename in swwfiles:
28        if verbose:
29            print 'Reading file ', filename, ':'   
30   
31        fid = NetCDFFile(filename, netcdf_mode_r)
32        tris = fid.variables['volumes'][:]       
33         
34        if first_file:
35            times = fid.variables['time'][:]
36            x = []
37            y = []
38            out_tris = list(tris) 
39            out_s_quantities = {}
40            out_d_quantities = {}
41           
42            for quantity in static_quantities:
43                out_s_quantities[quantity] = []
44
45            for quantity in dynamic_quantities:
46                out_d_quantities[quantity] = [ [] for _ in range(len(times))]
47                 
48            description = 'merged:' + getattr(fid, 'description')         
49            first_file = False
50        else:
51            for tri in tris:
52                verts = [vertex+tri_offset for vertex in tri]
53                out_tris.append(verts)
54
55        num_pts = fid.dimensions['number_of_points']
56        tri_offset += num_pts
57       
58        if verbose:
59            print '  new triangle index offset is ', tri_offset
60           
61        x.extend(list(fid.variables['x'][:]))
62        y.extend(list(fid.variables['y'][:]))
63       
64        for quantity in static_quantities:
65            out_s_quantities[quantity].extend(fid.variables[quantity][:])
66           
67        for quantity in dynamic_quantities:
68            time_chunks = fid.variables[quantity][:]
69            for i, time_chunk in enumerate(time_chunks):
70                out_d_quantities[quantity][i].extend(time_chunk)           
71       
72    points = [[xx, yy] for xx, yy in zip(x, y)]
73    fid.close()
74   
75    # NetCDF file definition
76    fido = NetCDFFile(output, netcdf_mode_w)
77    sww = Write_sww(static_quantities, dynamic_quantities)
78    sww.store_header(fido, times,
79                             len(out_tris),
80                             len(points),
81                             description=description,
82                             sww_precision=netcdf_float)
83
84
85    sww.store_triangulation(fido, points, out_tris)
86       
87    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
88
89    for q in dynamic_quantities:
90        q_values = out_d_quantities[q]
91        for i, time_slice in enumerate(q_values):
92            fido.variables[q][i] = num.array(time_slice, netcdf_float)
93       
94        # This updates the _range values
95        q_range = fido.variables[q + Write_sww.RANGE][:]
96        q_values_min = num.min(q_values)
97        if q_values_min < q_range[0]:
98            fido.variables[q + Write_sww.RANGE][0] = q_values_min
99        q_values_max = num.max(q_values)
100        if q_values_max > q_range[1]:
101            fido.variables[q + Write_sww.RANGE][1] = q_values_max       
102
103                                       
104    fido.close()
105   
106
Note: See TracBrowser for help on using the repository browser.