source: trunk/anuga_core/source/anuga/file_conversion/sww2dem.py @ 7796

Last change on this file since 7796 was 7796, checked in by hudson, 13 years ago

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