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

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

Improved docs for sww_merge.

File size: 4.1 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        The sww files to be merged must have exactly the same timesteps.
22       
23        swwfiles is a list of .sww files to merge.
24        output is the output filename, including .sww extension.
25        verbose True to log output information
26    """
27   
28    static_quantities = ['elevation']
29    dynamic_quantities = ['stage', 'xmomentum', 'ymomentum']
30   
31    first_file = True
32    tri_offset = 0
33    for filename in swwfiles:
34        if verbose:
35            print 'Reading file ', filename, ':'   
36   
37        fid = NetCDFFile(filename, netcdf_mode_r)
38        tris = fid.variables['volumes'][:]       
39         
40        if first_file:
41            times = fid.variables['time'][:]
42            x = []
43            y = []
44            out_tris = list(tris) 
45            out_s_quantities = {}
46            out_d_quantities = {}
47           
48            for quantity in static_quantities:
49                out_s_quantities[quantity] = []
50
51            # Quantities are stored as a 2D array of timesteps x data.
52            for quantity in dynamic_quantities:
53                out_d_quantities[quantity] = [ [] for _ in range(len(times))]
54                 
55            description = 'merged:' + getattr(fid, 'description')         
56            first_file = False
57        else:
58            for tri in tris:
59                # Advance new tri indices to point at newly appended points.
60                verts = [vertex+tri_offset for vertex in tri]
61                out_tris.append(verts)
62
63        num_pts = fid.dimensions['number_of_points']
64        tri_offset += num_pts
65       
66        if verbose:
67            print '  new triangle index offset is ', tri_offset
68           
69        x.extend(list(fid.variables['x'][:]))
70        y.extend(list(fid.variables['y'][:]))
71       
72        # Grow the list of static quantities associated with the x,y points
73        for quantity in static_quantities:
74            out_s_quantities[quantity].extend(fid.variables[quantity][:])
75           
76        #Collate all dynamic quantities according to their timestep
77        for quantity in dynamic_quantities:
78            time_chunks = fid.variables[quantity][:]
79            for i, time_chunk in enumerate(time_chunks):
80                out_d_quantities[quantity][i].extend(time_chunk)           
81   
82    # Mash all points into a single big list   
83    points = [[xx, yy] for xx, yy in zip(x, y)]
84    fid.close()
85   
86    # Write out the SWW file
87    fido = NetCDFFile(output, netcdf_mode_w)
88    sww = Write_sww(static_quantities, dynamic_quantities)
89    sww.store_header(fido, times,
90                             len(out_tris),
91                             len(points),
92                             description=description,
93                             sww_precision=netcdf_float)
94
95
96    sww.store_triangulation(fido, points, out_tris)
97       
98    sww.store_static_quantities(fido, verbose=verbose, **out_s_quantities)
99
100    # Write out all the dynamic quantities for each timestep
101    for q in dynamic_quantities:
102        q_values = out_d_quantities[q]
103        for i, time_slice in enumerate(q_values):
104            fido.variables[q][i] = num.array(time_slice, netcdf_float)
105       
106        # This updates the _range values
107        q_range = fido.variables[q + Write_sww.RANGE][:]
108        q_values_min = num.min(q_values)
109        if q_values_min < q_range[0]:
110            fido.variables[q + Write_sww.RANGE][0] = q_values_min
111        q_values_max = num.max(q_values)
112        if q_values_max > q_range[1]:
113            fido.variables[q + Write_sww.RANGE][1] = q_values_max       
114
115                                       
116    fido.close()
117   
118
Note: See TracBrowser for help on using the repository browser.