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

Last change on this file since 7819 was 7819, checked in by hudson, 12 years ago

Semi-working code for SWW merge util.

File size: 4.3 KB
Line 
1import numpy as num
2from anuga.utilities.numerical_tools import ensure_numeric
3
4from Scientific.IO.NetCDF import NetCDFFile
5from anuga.config import netcdf_mode_r, netcdf_mode_w, \
6                            netcdf_mode_a, netcdf_float
7from anuga.file.sww import SWW_file, Write_sww
8
9
10def sww_merge(swwfiles, output, verbose = False):
11    """
12        Merge a list of sww files into a single file.
13    """
14   
15    static_quantities = ['elevation']
16    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
17   
18    first_file = True
19    tri_offset = 0
20    for filename in swwfiles:
21        if verbose:
22            print 'Reading file ', filename, ':'   
23   
24        fid = NetCDFFile(filename, netcdf_mode_r)
25        tris = fid.variables['volumes'][:]       
26         
27        if first_file:
28            times = fid.variables['time'][:]
29            x = []
30            y = []
31            out_tris = list(tris) 
32            out_s_quantities = {}
33            out_d_quantities = {}
34           
35            for quantity in static_quantities:
36                out_s_quantities[quantity] = []
37
38            for quantity in dynamic_quantities:
39                out_d_quantities[quantity] = [ [] for _ in range(len(times))]
40                 
41            description = 'merged:' + getattr(fid, 'description')         
42            first_file = False
43        else:
44            for tri in tris:
45                verts = [vertex+tri_offset for vertex in tri]
46                out_tris.append(verts)
47
48        num_pts = fid.dimensions['number_of_points']
49        tri_offset += num_pts
50       
51        if verbose:
52            print '  new triangle index offset is ', tri_offset
53           
54        x.extend(list(fid.variables['x'][:]))
55        y.extend(list(fid.variables['y'][:]))
56       
57        for quantity in static_quantities:
58            out_s_quantities[quantity].extend(fid.variables[quantity][:])
59           
60        for quantity in dynamic_quantities:
61            time_chunks = fid.variables[quantity][:]
62            for i, time_chunk in enumerate(time_chunks):
63                out_d_quantities[quantity][i].extend(time_chunk)           
64       
65    points = [[xx, yy] for xx, yy in zip(x, y)]
66    fid.close()
67   
68    # NetCDF file definition
69    fido = NetCDFFile(output, netcdf_mode_w)
70    sww = Write_sww(static_quantities, dynamic_quantities)
71    sww.store_header(fido, times,
72                             len(out_tris),
73                             len(points),
74                             description=description,
75                             sww_precision=netcdf_float)
76
77
78    sww.store_triangulation(fido, points, out_tris)
79       
80    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
81
82    for q in dynamic_quantities:
83        q_values = ensure_numeric(out_d_quantities[quantity])
84        x = q_values.astype(netcdf_float)
85        fido.variables[q] = x
86
87        # This updates the _range values
88        q_range = fido.variables[q + Write_sww.RANGE][:]
89        q_values_min = num.min(q_values)
90        if q_values_min < q_range[0]:
91            fido.variables[q + Write_sww.RANGE][0] = q_values_min
92        q_values_max = num.max(q_values)
93        if q_values_max > q_range[1]:
94            fido.variables[q + Write_sww.RANGE][1] = q_values_max       
95
96                                       
97    fido.close()
98   
99
100from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular, \
101                                                            rectangular_cross
102from anuga.shallow_water.shallow_water_domain import Domain
103from anuga.file.sww import SWW_file
104from anuga.abstract_2d_finite_volumes.generic_boundary_conditions import \
105    Dirichlet_boundary
106
107Bd = Dirichlet_boundary([0.5, 0., 0.])
108
109# Create shallow water domain
110domain = Domain(*rectangular_cross(2, 2))
111domain.set_name('test1')
112domain.set_quantity('elevation', 2)
113domain.set_quantity('stage', 5)
114domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
115for t in domain.evolve(yieldstep=0.5, finaltime=1):
116    pass
117   
118domain = Domain(*rectangular(3, 3))
119domain.set_name('test2')
120domain.set_quantity('elevation', 3)
121domain.set_quantity('stage', 50)
122domain.set_boundary({'left': Bd, 'right': Bd, 'top': Bd, 'bottom': Bd})
123for t in domain.evolve(yieldstep=0.5, finaltime=1):
124    pass
125       
126sww_merge(['test1.sww', 'test2.sww'], 'test_out.sww', verbose = True)
Note: See TracBrowser for help on using the repository browser.