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

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

Added extra file conversion classes to file_conversion.

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