source: trunk/anuga_core/source/anuga/file_conversion/sww2array.py @ 8879

Last change on this file since 8879 was 8847, checked in by steve, 12 years ago

Small changes to sww2array looking for memory leak
File size: 10.5 KB
Line 
1"""
2    Module to convert SWW to DEM files.
3"""
4
5# external modules
6import os
7import numpy as num
8
9# ANUGA modules
10from anuga.abstract_2d_finite_volumes.util import remove_lone_verts     
11from anuga.coordinate_transforms.geo_reference import Geo_reference
12from anuga.utilities.system_tools import get_vars_in_expression
13import anuga.utilities.log as log
14from anuga.utilities.file_utils import get_all_swwfiles
15
16
17######
18# formula mappings
19######
20
21quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
22                    'depth':'stage-elevation',
23                    'speed': \
24 '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
25
26
27
28# Default block size for sww2dem()
29DEFAULT_BLOCK_SIZE = 100000
30
31def sww2array(name_in,
32            quantity=None, # defaults to elevation
33            reduction=None,
34            cellsize=10,
35            number_of_decimal_places=None,
36            NODATA_value=-9999.0,
37            easting_min=None,
38            easting_max=None,
39            northing_min=None,
40            northing_max=None,
41            verbose=False,
42            origin=None,
43            datum='WGS84',
44            block_size=None):
45    """Read SWW file and convert to a numpy array (can be stored to a png file later)
46
47
48    The parameter quantity must be the name of an existing quantity or
49    an expression involving existing quantities. The default is
50    'elevation'. Quantity is not a list of quantities.
51
52    If reduction is given and it's an index, sww2array will output the quantity at that time-step.
53    If reduction is given and it's a built in function (eg max, min, mean), then that
54    function is used to reduce the quantity over all time-steps. If reduction is not given,
55    reduction is set to "max" by default.
56
57    datum
58
59
60    block_size - sets the number of slices along the non-time axis to
61                 process in one block.
62    """
63
64    import sys
65    import types
66
67    from anuga.geometry.polygon import inside_polygon, outside_polygon
68    from anuga.abstract_2d_finite_volumes.util import \
69         apply_expression_to_dictionary
70
71    basename_in, in_ext = os.path.splitext(name_in)
72
73    if in_ext != '.sww':
74        raise IOError('Input format for %s must be .sww' % name_in)
75
76
77    false_easting = 500000
78    false_northing = 10000000
79
80    if quantity is None:
81        quantity = 'elevation'
82   
83    if reduction is None:
84        reduction = max
85
86    if quantity_formula.has_key(quantity):
87        quantity = quantity_formula[quantity]
88
89    if number_of_decimal_places is None:
90        number_of_decimal_places = 3
91
92    if block_size is None:
93        block_size = DEFAULT_BLOCK_SIZE
94
95    assert(isinstance(block_size, (int, long, float)))
96
97    # Read sww file
98    if verbose:
99        log.critical('Reading from %s' % name_in)
100
101
102    from anuga.file.netcdf import NetCDFFile
103    fid = NetCDFFile(name_in)
104
105    #Get extent and reference
106    x = num.array(fid.variables['x'], num.float)
107    y = num.array(fid.variables['y'], num.float)
108    volumes = num.array(fid.variables['volumes'], num.int)
109    if type(reduction) is not types.BuiltinFunctionType:
110        times = fid.variables['time'][reduction]
111    else:
112        times = fid.variables['time'][:]
113
114    number_of_timesteps = fid.dimensions['number_of_timesteps']
115    number_of_points = fid.dimensions['number_of_points']
116
117    if origin is None:
118        # Get geo_reference
119        # sww files don't have to have a geo_ref
120        try:
121            geo_reference = Geo_reference(NetCDFObject=fid)
122        except AttributeError, e:
123            geo_reference = Geo_reference() # Default georef object
124
125        xllcorner = geo_reference.get_xllcorner()
126        yllcorner = geo_reference.get_yllcorner()
127        zone = geo_reference.get_zone()
128    else:
129        zone = origin[0]
130        xllcorner = origin[1]
131        yllcorner = origin[2]
132
133    # FIXME: Refactor using code from Interpolation_function.statistics
134    # (in interpolate.py)
135    # Something like print swwstats(swwname)
136    if verbose:
137        log.critical('------------------------------------------------')
138        log.critical('Statistics of SWW file:')
139        log.critical('  Name: %s' % name_in)
140        log.critical('  Reference:')
141        log.critical('    Lower left corner: [%f, %f]' % (xllcorner, yllcorner))
142        if type(reduction) is not types.BuiltinFunctionType:
143            log.critical('    Time: %f' % times)
144        else:
145            log.critical('    Start time: %f' % fid.starttime[0])
146        log.critical('  Extent:')
147        log.critical('    x [m] in [%f, %f], len(x) == %d'
148                     %(num.min(x), num.max(x), len(x.flat)))
149        log.critical('    y [m] in [%f, %f], len(y) == %d'
150                     % (num.min(y), num.max(y), len(y.flat)))
151        if type(reduction) is not types.BuiltinFunctionType:
152            log.critical('    t [s] = %f, len(t) == %d' % (times, 1))
153        else:
154            log.critical('    t [s] in [%f, %f], len(t) == %d'
155                         % (min(times), max(times), len(times)))
156        log.critical('  Quantities [SI units]:')
157       
158        # Comment out for reduced memory consumption
159        for name in ['stage', 'xmomentum', 'ymomentum']:
160            q = fid.variables[name][:].flatten()
161            if type(reduction) is not types.BuiltinFunctionType:
162                q = q[reduction*len(x):(reduction+1)*len(x)]
163            if verbose: log.critical('    %s in [%f, %f]'
164                                     % (name, min(q), max(q)))
165        for name in ['elevation']:
166            q = fid.variables[name][:].flatten()
167            if verbose: log.critical('    %s in [%f, %f]'
168                                     % (name, min(q), max(q)))
169
170    # Get the variables in the supplied expression.
171    # This may throw a SyntaxError exception.
172    var_list = get_vars_in_expression(quantity)
173
174    # Check that we have the required variables in the SWW file.
175    missing_vars = []
176    for name in var_list:
177        try:
178            _ = fid.variables[name]
179        except KeyError:
180            missing_vars.append(name)
181    if missing_vars:
182        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
183               % (quantity, str(missing_vars), name_in))
184        raise Exception, msg
185
186    # Create result array and start filling, block by block.
187    result = num.zeros(number_of_points, num.float)
188
189    if verbose:
190        msg = 'Slicing sww file, num points: ' + str(number_of_points)
191        msg += ', block size: ' + str(block_size)
192        log.critical(msg)
193
194    for start_slice in xrange(0, number_of_points, block_size):
195        # Limit slice size to array end if at last block
196        end_slice = min(start_slice + block_size, number_of_points)
197       
198        # Get slices of all required variables
199        if type(reduction) is not types.BuiltinFunctionType:
200            q_dict = {}
201            for name in var_list:
202                # check if variable has time axis
203                if len(fid.variables[name].shape) == 2:
204                    print 'avoiding large array'
205                    q_dict[name] = fid.variables[name][reduction,start_slice:end_slice]
206                else:       # no time axis
207                    q_dict[name] = fid.variables[name][start_slice:end_slice]
208
209            # Evaluate expression with quantities found in SWW file
210            res = apply_expression_to_dictionary(quantity, q_dict)
211
212#            if len(res.shape) == 2:
213#                new_res = num.zeros(res.shape[1], num.float)
214#                for k in xrange(res.shape[1]):
215#                    if type(reduction) is not types.BuiltinFunctionType:
216#                        new_res[k] = res[k]
217#                    else:
218#                        new_res[k] = reduction(res[:,k])
219#                res = new_res
220        else:
221            q_dict = {}
222            for name in var_list:
223                # check if variable has time axis
224                if len(fid.variables[name].shape) == 2:
225                    q_dict[name] = fid.variables[name][:,start_slice:end_slice]
226                else:       # no time axis
227                    q_dict[name] = fid.variables[name][start_slice:end_slice]
228
229            # Evaluate expression with quantities found in SWW file
230            res = apply_expression_to_dictionary(quantity, q_dict)
231
232            if len(res.shape) == 2:
233                new_res = num.zeros(res.shape[1], num.float)
234                for k in xrange(res.shape[1]):
235                    if type(reduction) is not types.BuiltinFunctionType:
236                        new_res[k] = res[reduction,k]
237                    else:
238                        new_res[k] = reduction(res[:,k])
239                res = new_res
240
241        result[start_slice:end_slice] = res
242                                   
243    # Post condition: Now q has dimension: number_of_points
244    assert len(result.shape) == 1
245    assert result.shape[0] == number_of_points
246
247    if verbose:
248        log.critical('Processed values for %s are in [%f, %f]'
249                     % (quantity, min(result), max(result)))
250
251    # Create grid and update xll/yll corner and x,y
252    # Relative extent
253    if easting_min is None:
254        xmin = min(x)
255    else:
256        xmin = easting_min - xllcorner
257
258    if easting_max is None:
259        xmax = max(x)
260    else:
261        xmax = easting_max - xllcorner
262
263    if northing_min is None:
264        ymin = min(y)
265    else:
266        ymin = northing_min - yllcorner
267
268    if northing_max is None:
269        ymax = max(y)
270    else:
271        ymax = northing_max - yllcorner
272
273    msg = 'xmax must be greater than or equal to xmin.\n'
274    msg += 'I got xmin = %f, xmax = %f' %(xmin, xmax)
275    assert xmax >= xmin, msg
276
277    msg = 'ymax must be greater than or equal to xmin.\n'
278    msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
279    assert ymax >= ymin, msg
280
281    if verbose: log.critical('Creating grid')
282    ncols = int((xmax-xmin)/cellsize) + 1
283    nrows = int((ymax-ymin)/cellsize) + 1
284
285    # New absolute reference and coordinates
286    newxllcorner = xmin + xllcorner
287    newyllcorner = ymin + yllcorner
288
289    x = x + xllcorner - newxllcorner
290    y = y + yllcorner - newyllcorner
291
292
293    grid_values = num.zeros( (nrows*ncols, ), num.float)
294    #print '---',grid_values.shape
295
296    num_tri =  len(volumes)
297    norms = num.zeros(6*num_tri, num.float)
298
299    #Use fasr method to calc grid values
300    from calc_grid_values_ext import calc_grid_values
301
302    calc_grid_values(nrows, ncols, cellsize, NODATA_value,
303                     x,y, norms, volumes, result, grid_values)
304
305
306    fid.close()
307   
308    #print outside_indices
309
310    if verbose:
311        log.critical('Interpolated values are in [%f, %f]'
312                     % (num.min(grid_values), num.max(grid_values)))
313
314
315
316    return grid_values.reshape(nrows,ncols)[::-1,:]
317
318
319
320
321
Note: See TracBrowser for help on using the repository browser.