[7817] | 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) |
---|