source: trunk/anuga_core/source/anuga/file_conversion/sww2pts.py @ 7814

Last change on this file since 7814 was 7814, checked in by hudson, 12 years ago

New filename conventions for file conversion. Filenames must always be passed in with the correct extension.

File size: 6.3 KB
Line 
1
2import numpy as num
3import os
4
5from anuga.coordinate_transforms.geo_reference import Geo_reference
6
7##
8# @brief Convert SWW file to PTS file (at selected points).
9# @param basename_in Stem name of input SWW file.
10# @param basename_out Stem name of output file.
11# @param data_points If given, points where quantity is to be computed.
12# @param quantity Name (or expression) of existing quantity(s) (def: elevation).
13# @param timestep If given, output quantity at that timestep.
14# @param reduction If given, reduce quantity by this factor.
15# @param NODATA_value The NODATA value (default -9999).
16# @param verbose True if this function is to be verbose.
17# @param origin ??
18def sww2pts(name_in, name_out=None,
19            data_points=None,
20            quantity=None,
21            timestep=None,
22            reduction=None,
23            NODATA_value=-9999,
24            verbose=False,
25            origin=None):
26    """Read SWW file and convert to interpolated values at selected points
27
28    The parameter 'quantity' must be the name of an existing quantity or
29    an expression involving existing quantities. The default is 'elevation'.
30
31    if timestep (an index) is given, output quantity at that timestep.
32
33    if reduction is given use that to reduce quantity over all timesteps.
34
35    data_points (Nx2 array) give locations of points where quantity is to
36    be computed.
37    """
38
39    import sys
40    from anuga.geometry.polygon import inside_polygon, outside_polygon
41    from anuga.abstract_2d_finite_volumes.util import \
42             apply_expression_to_dictionary
43    from anuga.geospatial_data.geospatial_data import Geospatial_data
44
45    if quantity is None:
46        quantity = 'elevation'
47
48    if reduction is None:
49        reduction = max
50
51    basename_in, in_ext = os.path.splitext(name_in)
52   
53    if name_out != None:
54        basename_out, out_ext = os.path.splitext(name_out)
55    else:
56        basename_out = basename_in + '_%s' % quantity
57        out_ext = '.pts'
58        name_out = basename_out + out_ext
59
60    if in_ext != '.sww':
61        raise IOError('Input format for %s must be .sww' % name_in)
62
63    if out_ext != '.pts':
64        raise IOError('Output format for %s must be .pts' % name_out)
65
66
67    # Read sww file
68    if verbose: log.critical('Reading from %s' % name_in)
69    from Scientific.IO.NetCDF import NetCDFFile
70    fid = NetCDFFile(name_in)
71
72    # Get extent and reference
73    x = fid.variables['x'][:]
74    y = fid.variables['y'][:]
75    volumes = fid.variables['volumes'][:]
76
77    number_of_timesteps = fid.dimensions['number_of_timesteps']
78    number_of_points = fid.dimensions['number_of_points']
79    if origin is None:
80        # Get geo_reference
81        # sww files don't have to have a geo_ref
82        try:
83            geo_reference = Geo_reference(NetCDFObject=fid)
84        except AttributeError, e:
85            geo_reference = Geo_reference() # Default georef object
86
87        xllcorner = geo_reference.get_xllcorner()
88        yllcorner = geo_reference.get_yllcorner()
89        zone = geo_reference.get_zone()
90    else:
91        zone = origin[0]
92        xllcorner = origin[1]
93        yllcorner = origin[2]
94
95    # FIXME: Refactor using code from file_function.statistics
96    # Something like print swwstats(swwname)
97    if verbose:
98        x = fid.variables['x'][:]
99        y = fid.variables['y'][:]
100        times = fid.variables['time'][:]
101        log.critical('------------------------------------------------')
102        log.critical('Statistics of SWW file:')
103        log.critical('  Name: %s' % swwfile)
104        log.critical('  Reference:')
105        log.critical('    Lower left corner: [%f, %f]' % (xllcorner, yllcorner))
106        log.critical('    Start time: %f' % fid.starttime[0])
107        log.critical('  Extent:')
108        log.critical('    x [m] in [%f, %f], len(x) == %d'
109                     % (num.min(x), num.max(x), len(x.flat)))
110        log.critical('    y [m] in [%f, %f], len(y) == %d'
111                     % (num.min(y), num.max(y), len(y.flat)))
112        log.critical('    t [s] in [%f, %f], len(t) == %d'
113                     % (min(times), max(times), len(times)))
114        log.critical('  Quantities [SI units]:')
115        for name in ['stage', 'xmomentum', 'ymomentum', 'elevation']:
116            q = fid.variables[name][:].flat
117            log.critical('    %s in [%f, %f]' % (name, min(q), max(q)))
118
119    # Get quantity and reduce if applicable
120    if verbose: log.critical('Processing quantity %s' % quantity)
121
122    # Turn NetCDF objects into numeric arrays
123    quantity_dict = {}
124    for name in fid.variables.keys():
125        quantity_dict[name] = fid.variables[name][:]
126
127    # Convert quantity expression to quantities found in sww file
128    q = apply_expression_to_dictionary(quantity, quantity_dict)
129
130    if len(q.shape) == 2:
131        # q has a time component and needs to be reduced along
132        # the temporal dimension
133        if verbose: log.critical('Reducing quantity %s' % quantity)
134
135        q_reduced = num.zeros(number_of_points, num.float)
136        for k in range(number_of_points):
137            q_reduced[k] = reduction(q[:,k])
138        q = q_reduced
139
140    # Post condition: Now q has dimension: number_of_points
141    assert len(q.shape) == 1
142    assert q.shape[0] == number_of_points
143
144    if verbose:
145        log.critical('Processed values for %s are in [%f, %f]'
146                     % (quantity, min(q), max(q)))
147
148    # Create grid and update xll/yll corner and x,y
149    vertex_points = num.concatenate((x[:, num.newaxis], y[:, num.newaxis]), axis=1)
150    assert len(vertex_points.shape) == 2
151
152    # Interpolate
153    from anuga.fit_interpolate.interpolate import Interpolate
154    interp = Interpolate(vertex_points, volumes, verbose=verbose)
155
156    # Interpolate using quantity values
157    if verbose: log.critical('Interpolating')
158    interpolated_values = interp.interpolate(q, data_points).flatten()
159
160    if verbose:
161        log.critical('Interpolated values are in [%f, %f]'
162                     % (num.min(interpolated_values),
163                        num.max(interpolated_values)))
164
165    # Assign NODATA_value to all points outside bounding polygon
166    # (from interpolation mesh)
167    P = interp.mesh.get_boundary_polygon()
168    outside_indices = outside_polygon(data_points, P, closed=True)
169
170    for i in outside_indices:
171        interpolated_values[i] = NODATA_value
172
173    # Store results
174    G = Geospatial_data(data_points=data_points, attributes=interpolated_values)
175
176    G.export_points_file(name_out, absolute = True)
177
178    fid.close()
179
Note: See TracBrowser for help on using the repository browser.