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

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

Added sww_merge utility module for merging sww files into one.

File size: 2.8 KB
Line 
1import numpy as num
2
3from Scientific.IO.NetCDF import NetCDFFile
4from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
5from anuga.file.sww import SWW_file, Write_sww
6
7def sww_merge(swwfiles, output, verbose = False):
8    """
9        Merge a list of sww files into a single file.
10    """
11   
12    first_file = True
13    tri_offset = 0
14    for filename in swwfiles:
15        if verbose:
16            print 'Reading file ', filename, ':'   
17   
18        fid = NetCDFFile(filename, netcdf_mode_r)
19        tris = fid.variables['volumes'][:]       
20         
21        if first_file:
22            times = fid.variables['time'][:]
23            x = []
24            y = []
25            out_tris = list(tris) 
26            description = 'merged:' + getattr(fid, 'description')         
27            first_file = False
28        else:
29            for tri in tris:
30                verts = [vertex+tri_offset for vertex in tri]
31                out_tris.append(verts)
32
33        tri_offset += fid.dimensions['number_of_points']
34       
35        if verbose:
36            print '  new triangle index offset is ', tri_offset
37           
38        x.extend(list(fid.variables['x'][:]))
39        y.extend(list(fid.variables['y'][:]))
40     
41    points = [[xx, yy] for xx, yy in zip(x, y)]   
42
43
44    # NetCDF file definition
45    fido = NetCDFFile(output, netcdf_mode_w)
46    sww = Write_sww(['elevation'], ['stage', 'xmomentum', 'ymomentum'])
47    sww.store_header(fido, times,
48                             len(out_tris),
49                             len(points),
50                             description=description)
51
52    sww.store_triangulation(fido,
53                                    points,
54                                    num.array(out_tris).astype(num.float32))
55
56
57    # Get names of static quantities
58    static_quantities = {}
59    for name in ['elevation']:
60        static_quantities[name] = fid.variables[name][:]
61   
62    # Store static quantities       
63    self.writer.store_static_quantities(fido, **static_quantities)
64                                       
65    fid.close()
66
67
68
69    domain = Domain(points, out_tris)
70    domain.set_name(output)
71    sww = SWW_file(domain)
72    sww.store_connectivity()
73   
74    for _ in times:
75        sww.store_timestep()
76   
77
78from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular, \
79                                                            rectangular_cross
80from anuga.shallow_water.shallow_water_domain import Domain
81from anuga.file.sww import SWW_file
82
83# Create shallow water domain
84domain = Domain(*rectangular_cross(2, 2))
85domain.set_name('test1')
86sww = SWW_file(domain)
87sww.store_connectivity()
88
89domain = Domain(*rectangular(3, 3))
90domain.set_name('test2')
91sww = SWW_file(domain)
92sww.store_connectivity()
93
94       
95sww_merge(['test1.sww', 'test2.sww'], 'test_out.sww', verbose = True)
Note: See TracBrowser for help on using the repository browser.