source: anuga_core/source/anuga/file_conversion/sww2dem.py @ 7744

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

Split off sww2dem from data_manager.
File size: 13.9 KB
Line 
1# external modules
2import numpy as num
3
4# ANUGA modules
5from anuga.abstract_2d_finite_volumes.util import remove_lone_verts     
6from anuga.coordinate_transforms.geo_reference import Geo_reference
7from anuga.utilities.system_tools import get_vars_in_expression
8import anuga.utilities.log as log
9
10
11######
12# formula mappings
13######
14
15quantity_formula = {'momentum':'(xmomentum**2 + ymomentum**2)**0.5',
16                    'depth':'stage-elevation',
17                    'speed': \
18 '(xmomentum**2 + ymomentum**2)**0.5/(stage-elevation+1.e-6/(stage-elevation))'}
19
20
21
22# Default block size for sww2dem()
23DEFAULT_BLOCK_SIZE = 10000
24
25def sww2dem(basename_in, basename_out=None,
26            quantity=None, # defaults to elevation
27            reduction=None,
28            cellsize=10,
29            number_of_decimal_places=None,
30            NODATA_value=-9999,
31            easting_min=None,
32            easting_max=None,
33            northing_min=None,
34            northing_max=None,
35            verbose=False,
36            origin=None,
37            datum='WGS84',
38            format='ers',
39            block_size=None):
40    """Read SWW file and convert to Digitial Elevation model format
41    (.asc or .ers)
42
43    Example (ASC):
44    ncols         3121
45    nrows         1800
46    xllcorner     722000
47    yllcorner     5893000
48    cellsize      25
49    NODATA_value  -9999
50    138.3698 137.4194 136.5062 135.5558 ..........
51
52    The number of decimal places can be specified by the user to save
53    on disk space requirements by specifying in the call to sww2dem.
54
55    Also write accompanying file with same basename_in but extension .prj
56    used to fix the UTM zone, datum, false northings and eastings.
57
58    The prj format is assumed to be as
59
60    Projection    UTM
61    Zone          56
62    Datum         WGS84
63    Zunits        NO
64    Units         METERS
65    Spheroid      WGS84
66    Xshift        0.0000000000
67    Yshift        10000000.0000000000
68    Parameters
69
70
71    The parameter quantity must be the name of an existing quantity or
72    an expression involving existing quantities. The default is
73    'elevation'. Quantity is not a list of quantities.
74
75    if timestep (an index) is given, output quantity at that timestep
76
77    if reduction is given and its an index, output quantity at that timestep. If reduction is given
78    and is a built in function, use that to reduce quantity over all timesteps.
79
80    datum
81
82    format can be either 'asc' or 'ers'
83    block_size - sets the number of slices along the non-time axis to
84                 process in one block.
85    """
86
87    import sys
88    import types
89
90    from anuga.geometry.polygon import inside_polygon, outside_polygon
91    from anuga.abstract_2d_finite_volumes.util import \
92         apply_expression_to_dictionary
93
94    msg = 'Format must be either asc or ers'
95    assert format.lower() in ['asc', 'ers'], msg
96
97    false_easting = 500000
98    false_northing = 10000000
99
100    if quantity is None:
101        quantity = 'elevation'
102   
103    if reduction is None:
104        reduction = max
105
106    if basename_out is None:
107        basename_out = basename_in + '_%s' % quantity
108
109    if quantity_formula.has_key(quantity):
110        quantity = quantity_formula[quantity]
111
112    if number_of_decimal_places is None:
113        number_of_decimal_places = 3
114
115    if block_size is None:
116        block_size = DEFAULT_BLOCK_SIZE
117
118    # Read SWW file
119    swwfile = basename_in + '.sww'
120    demfile = basename_out + '.' + format
121
122    # Read sww file
123    if verbose:
124        log.critical('Reading from %s' % swwfile)
125        log.critical('Output directory is %s' % basename_out)
126
127    from Scientific.IO.NetCDF import NetCDFFile
128    fid = NetCDFFile(swwfile)
129
130    #Get extent and reference
131    x = fid.variables['x'][:]
132    y = fid.variables['y'][:]
133    volumes = fid.variables['volumes'][:]
134    if type(reduction) is not types.BuiltinFunctionType:
135        times = fid.variables['time'][reduction]
136    else:
137        times = fid.variables['time'][:]
138
139    number_of_timesteps = fid.dimensions['number_of_timesteps']
140    number_of_points = fid.dimensions['number_of_points']
141
142    if origin is None:
143        # Get geo_reference
144        # sww files don't have to have a geo_ref
145        try:
146            geo_reference = Geo_reference(NetCDFObject=fid)
147        except AttributeError, e:
148            geo_reference = Geo_reference() # Default georef object
149
150        xllcorner = geo_reference.get_xllcorner()
151        yllcorner = geo_reference.get_yllcorner()
152        zone = geo_reference.get_zone()
153    else:
154        zone = origin[0]
155        xllcorner = origin[1]
156        yllcorner = origin[2]
157
158    # FIXME: Refactor using code from Interpolation_function.statistics
159    # (in interpolate.py)
160    # Something like print swwstats(swwname)
161    if verbose:
162        log.critical('------------------------------------------------')
163        log.critical('Statistics of SWW file:')
164        log.critical('  Name: %s' % swwfile)
165        log.critical('  Reference:')
166        log.critical('    Lower left corner: [%f, %f]' % (xllcorner, yllcorner))
167        if type(reduction) is not types.BuiltinFunctionType:
168            log.critical('    Time: %f' % times)
169        else:
170            log.critical('    Start time: %f' % fid.starttime[0])
171        log.critical('  Extent:')
172        log.critical('    x [m] in [%f, %f], len(x) == %d'
173                     %(num.min(x), num.max(x), len(x.flat)))
174        log.critical('    y [m] in [%f, %f], len(y) == %d'
175                     % (num.min(y), num.max(y), len(y.flat)))
176        if type(reduction) is not types.BuiltinFunctionType:
177            log.critical('    t [s] = %f, len(t) == %d' % (times, 1))
178        else:
179            log.critical('    t [s] in [%f, %f], len(t) == %d'
180                         % (min(times), max(times), len(times)))
181        log.critical('  Quantities [SI units]:')
182       
183        # Comment out for reduced memory consumption
184        for name in ['stage', 'xmomentum', 'ymomentum']:
185            q = fid.variables[name][:].flatten()
186            if type(reduction) is not types.BuiltinFunctionType:
187                q = q[reduction*len(x):(reduction+1)*len(x)]
188            if verbose: log.critical('    %s in [%f, %f]'
189                                     % (name, min(q), max(q)))
190        for name in ['elevation']:
191            q = fid.variables[name][:].flatten()
192            if verbose: log.critical('    %s in [%f, %f]'
193                                     % (name, min(q), max(q)))
194
195    # Get the variables in the supplied expression.
196    # This may throw a SyntaxError exception.
197    var_list = get_vars_in_expression(quantity)
198
199    # Check that we have the required variables in the SWW file.
200    missing_vars = []
201    for name in var_list:
202        try:
203            _ = fid.variables[name]
204        except:
205            missing_vars.append(name)
206    if missing_vars:
207        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
208               % (quantity, swwfile))
209        raise Exception, msg
210
211    # Create result array and start filling, block by block.
212    result = num.zeros(number_of_points, num.float)
213
214    for start_slice in xrange(0, number_of_points, block_size):
215        # Limit slice size to array end if at last block
216        end_slice = min(start_slice + block_size, number_of_points)
217       
218        # Get slices of all required variables
219        q_dict = {}
220        for name in var_list:
221            # check if variable has time axis
222            if len(fid.variables[name].shape) == 2:
223                q_dict[name] = fid.variables[name][:,start_slice:end_slice]
224            else:       # no time axis
225                q_dict[name] = fid.variables[name][start_slice:end_slice]
226
227        # Evaluate expression with quantities found in SWW file
228        res = apply_expression_to_dictionary(quantity, q_dict)
229
230        if len(res.shape) == 2:
231            new_res = num.zeros(res.shape[1], num.float)
232            for k in xrange(res.shape[1]):
233                if type(reduction) is not types.BuiltinFunctionType:
234                    new_res[k] = res[reduction,k]
235                else:
236                    new_res[k] = reduction(res[:,k])
237            res = new_res
238
239        result[start_slice:end_slice] = res
240                                   
241    # Post condition: Now q has dimension: number_of_points
242    assert len(result.shape) == 1
243    assert result.shape[0] == number_of_points
244
245    if verbose:
246        log.critical('Processed values for %s are in [%f, %f]'
247                     % (quantity, min(result), max(result)))
248
249    # Create grid and update xll/yll corner and x,y
250    # Relative extent
251    if easting_min is None:
252        xmin = min(x)
253    else:
254        xmin = easting_min - xllcorner
255
256    if easting_max is None:
257        xmax = max(x)
258    else:
259        xmax = easting_max - xllcorner
260
261    if northing_min is None:
262        ymin = min(y)
263    else:
264        ymin = northing_min - yllcorner
265
266    if northing_max is None:
267        ymax = max(y)
268    else:
269        ymax = northing_max - yllcorner
270
271    msg = 'xmax must be greater than or equal to xmin.\n'
272    msg += 'I got xmin = %f, xmax = %f' %(xmin, xmax)
273    assert xmax >= xmin, msg
274
275    msg = 'ymax must be greater than or equal to xmin.\n'
276    msg += 'I got ymin = %f, ymax = %f' %(ymin, ymax)
277    assert ymax >= ymin, msg
278
279    if verbose: log.critical('Creating grid')
280    ncols = int((xmax-xmin)/cellsize) + 1
281    nrows = int((ymax-ymin)/cellsize) + 1
282
283    # New absolute reference and coordinates
284    newxllcorner = xmin + xllcorner
285    newyllcorner = ymin + yllcorner
286
287    x = x + xllcorner - newxllcorner
288    y = y + yllcorner - newyllcorner
289
290    vertex_points = num.concatenate ((x[:,num.newaxis], y[:,num.newaxis]), axis=1)
291    assert len(vertex_points.shape) == 2
292
293    grid_points = num.zeros ((ncols*nrows, 2), num.float)
294
295    for i in xrange(nrows):
296        if format.lower() == 'asc':
297            yg = i * cellsize
298        else:
299            # this will flip the order of the y values for ers
300            yg = (nrows-i) * cellsize
301
302        for j in xrange(ncols):
303            xg = j * cellsize
304            k = i*ncols + j
305
306            grid_points[k, 0] = xg
307            grid_points[k, 1] = yg
308
309    # Interpolate
310    from anuga.fit_interpolate.interpolate import Interpolate
311
312    # Remove loners from vertex_points, volumes here
313    vertex_points, volumes = remove_lone_verts(vertex_points, volumes)
314    # export_mesh_file('monkey.tsh',{'vertices':vertex_points, 'triangles':volumes})
315    interp = Interpolate(vertex_points, volumes, verbose = verbose)
316
317    # Interpolate using quantity values
318    if verbose: log.critical('Interpolating')
319    grid_values = interp.interpolate(result, grid_points).flatten()
320
321    if verbose:
322        log.critical('Interpolated values are in [%f, %f]'
323                     % (num.min(grid_values), num.max(grid_values)))
324
325    # Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
326    P = interp.mesh.get_boundary_polygon()
327    outside_indices = outside_polygon(grid_points, P, closed=True)
328
329    for i in outside_indices:
330        grid_values[i] = NODATA_value
331
332    if format.lower() == 'ers':
333        # setup ERS header information
334        grid_values = num.reshape(grid_values, (nrows, ncols))
335        header = {}
336        header['datum'] = '"' + datum + '"'
337        # FIXME The use of hardwired UTM and zone number needs to be made optional
338        # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
339        header['projection'] = '"UTM-' + str(zone) + '"'
340        header['coordinatetype'] = 'EN'
341        if header['coordinatetype'] == 'LL':
342            header['longitude'] = str(newxllcorner)
343            header['latitude'] = str(newyllcorner)
344        elif header['coordinatetype'] == 'EN':
345            header['eastings'] = str(newxllcorner)
346            header['northings'] = str(newyllcorner)
347        header['nullcellvalue'] = str(NODATA_value)
348        header['xdimension'] = str(cellsize)
349        header['ydimension'] = str(cellsize)
350        header['value'] = '"' + quantity + '"'
351        #header['celltype'] = 'IEEE8ByteReal'  #FIXME: Breaks unit test
352
353        #Write
354        if verbose: log.critical('Writing %s' % demfile)
355
356        import ermapper_grids
357
358        ermapper_grids.write_ermapper_grid(demfile, grid_values, header)
359
360        fid.close()
361    else:
362        #Write to Ascii format
363        #Write prj file
364        prjfile = basename_out + '.prj'
365
366        if verbose: log.critical('Writing %s' % prjfile)
367        prjid = open(prjfile, 'w')
368        prjid.write('Projection    %s\n' %'UTM')
369        prjid.write('Zone          %d\n' %zone)
370        prjid.write('Datum         %s\n' %datum)
371        prjid.write('Zunits        NO\n')
372        prjid.write('Units         METERS\n')
373        prjid.write('Spheroid      %s\n' %datum)
374        prjid.write('Xshift        %d\n' %false_easting)
375        prjid.write('Yshift        %d\n' %false_northing)
376        prjid.write('Parameters\n')
377        prjid.close()
378
379        if verbose: log.critical('Writing %s' % demfile)
380
381        ascid = open(demfile, 'w')
382
383        ascid.write('ncols         %d\n' %ncols)
384        ascid.write('nrows         %d\n' %nrows)
385        ascid.write('xllcorner     %d\n' %newxllcorner)
386        ascid.write('yllcorner     %d\n' %newyllcorner)
387        ascid.write('cellsize      %f\n' %cellsize)
388        ascid.write('NODATA_value  %d\n' %NODATA_value)
389
390        #Get bounding polygon from mesh
391        #P = interp.mesh.get_boundary_polygon()
392        #inside_indices = inside_polygon(grid_points, P)
393
394        # change printoptions so that a long string of zeros in not
395        # summarized as [0.0, 0.0, 0.0, ... 0.0, 0.0, 0.0]
396        #printoptions = num.get_printoptions()
397        #num.set_printoptions(threshold=sys.maxint)
398
399        format = '%.'+'%g' % number_of_decimal_places +'e'
400        for i in range(nrows):
401            if verbose and i % ((nrows+10)/10) == 0:
402                log.critical('Doing row %d of %d' % (i, nrows))
403
404            base_index = (nrows-i-1)*ncols
405
406            slice = grid_values[base_index:base_index+ncols]
407
408            num.savetxt(ascid, slice.reshape(1,ncols), format, ' ' )
409           
410       
411        #Close
412        ascid.close()
413        fid.close()
414
415        return basename_out
Note: See TracBrowser for help on using the repository browser.