1 | import numpy as num |
---|
2 | |
---|
3 | from Scientific.IO.NetCDF import NetCDFFile |
---|
4 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
5 | from anuga.file.sww import SWW_file, Write_sww |
---|
6 | |
---|
7 | def 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 | |
---|
78 | from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular, \ |
---|
79 | rectangular_cross |
---|
80 | from anuga.shallow_water.shallow_water_domain import Domain |
---|
81 | from anuga.file.sww import SWW_file |
---|
82 | |
---|
83 | # Create shallow water domain |
---|
84 | domain = Domain(*rectangular_cross(2, 2)) |
---|
85 | domain.set_name('test1') |
---|
86 | sww = SWW_file(domain) |
---|
87 | sww.store_connectivity() |
---|
88 | |
---|
89 | domain = Domain(*rectangular(3, 3)) |
---|
90 | domain.set_name('test2') |
---|
91 | sww = SWW_file(domain) |
---|
92 | sww.store_connectivity() |
---|
93 | |
---|
94 | |
---|
95 | sww_merge(['test1.sww', 'test2.sww'], 'test_out.sww', verbose = True) |
---|