source: trunk/anuga_core/source/anuga/file/sts.py @ 7871

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

Cleaned up some code.

File size: 12.1 KB
Line 
1import numpy as num
2import anuga.utilities.log as log
3from Scientific.IO.NetCDF import NetCDFFile
4
5from anuga.config import max_float
6from anuga.config import netcdf_float, netcdf_float32, netcdf_int
7from anuga.utilities.numerical_tools import ensure_numeric
8from anuga.coordinate_transforms.geo_reference import Geo_reference, \
9     ensure_geo_reference
10from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
11
12class Write_sts:
13    """ A class to write STS files.
14    """
15    sts_quantities = ['stage', 'xmomentum', 'ymomentum']
16    RANGE = '_range'
17    EXTREMA = ':extrema'
18
19    ##
20    # @brief Instantiate this STS writer class.
21    def __init__(self):
22        pass
23
24    ##
25    # @brief Write a header to the underlying data file.
26    # @param outfile Handle to open file to write.
27    # @param times A list of the time slice times *or* a start time.
28    # @param number_of_points The number of URS gauge sites.
29    # @param description Description string to write into the STS file.
30    # @param sts_precision Format of data to write (netcdf constant ONLY).
31    # @param verbose True if this function is to be verbose.
32    # @note If 'times' is a list, the info will be made relative.
33    def store_header(self,
34                     outfile,
35                     times,
36                     number_of_points,
37                     description='Converted from URS mux2 format',
38                     sts_precision=netcdf_float32,
39                     verbose=False):
40        """
41        outfile - the name of the file that will be written
42        times - A list of the time slice times OR a start time
43        Note, if a list is given the info will be made relative.
44        number_of_points - the number of urs gauges sites
45        """
46
47        outfile.institution = 'Geoscience Australia'
48        outfile.description = description
49
50        try:
51            revision_number = get_revision_number()
52        except:
53            revision_number = None
54
55        # Allow None to be stored as a string
56        outfile.revision_number = str(revision_number)
57
58        # Start time in seconds since the epoch (midnight 1/1/1970)
59        # This is being used to seperate one number from a list.
60        # what it is actually doing is sorting lists from numeric arrays.
61        if isinstance(times, (list, num.ndarray)):
62            number_of_times = len(times)
63            times = ensure_numeric(times)
64            if number_of_times == 0:
65                starttime = 0
66            else:
67                starttime = times[0]
68                times = times - starttime  #Store relative times
69        else:
70            number_of_times = 0
71            starttime = times
72
73        outfile.starttime = starttime
74
75        # Dimension definitions
76        outfile.createDimension('number_of_points', number_of_points)
77        outfile.createDimension('number_of_timesteps', number_of_times)
78        outfile.createDimension('numbers_in_range', 2)
79
80        # Variable definitions
81        outfile.createVariable('permutation', netcdf_int, ('number_of_points',))
82        outfile.createVariable('x', sts_precision, ('number_of_points',))
83        outfile.createVariable('y', sts_precision, ('number_of_points',))
84        outfile.createVariable('elevation', sts_precision, \
85                                    ('number_of_points',))
86
87        q = 'elevation'
88        outfile.createVariable(q + Write_sts.RANGE, sts_precision,
89                               ('numbers_in_range',))
90
91        # Initialise ranges with small and large sentinels.
92        # If this was in pure Python we could have used None sensibly
93        outfile.variables[q + Write_sts.RANGE][0] = max_float  # Min
94        outfile.variables[q + Write_sts.RANGE][1] = -max_float # Max
95
96        self.write_dynamic_quantities(outfile, Write_sts.sts_quantities, times)
97
98    ##
99    # @brief
100    # @param outfile
101    # @param points_utm
102    # @param elevation
103    # @param zone
104    # @param new_origin
105    # @param points_georeference
106    # @param verbose True if this function is to be verbose.
107    def store_points(self,
108                     outfile,
109                     points_utm,
110                     elevation, zone=None, new_origin=None,
111                     points_georeference=None, verbose=False):
112
113        """
114        points_utm - currently a list or array of the points in UTM.
115        points_georeference - the georeference of the points_utm
116
117        How about passing new_origin and current_origin.
118        If you get both, do a convertion from the old to the new.
119
120        If you only get new_origin, the points are absolute,
121        convert to relative
122
123        if you only get the current_origin the points are relative, store
124        as relative.
125
126        if you get no georefs create a new georef based on the minimums of
127        points_utm.  (Another option would be to default to absolute)
128
129        Yes, and this is done in another part of the code.
130        Probably geospatial.
131
132        If you don't supply either geo_refs, then supply a zone. If not
133        the default zone will be used.
134
135        precondition:
136             header has been called.
137        """
138
139        number_of_points = len(points_utm)
140        points_utm = num.array(points_utm)
141
142        # given the two geo_refs and the points, do the stuff
143        # described in the method header
144        points_georeference = ensure_geo_reference(points_georeference)
145        new_origin = ensure_geo_reference(new_origin)
146
147        if new_origin is None and points_georeference is not None:
148            points = points_utm
149            geo_ref = points_georeference
150        else:
151            if new_origin is None:
152                new_origin = Geo_reference(zone, min(points_utm[:,0]),
153                                                 min(points_utm[:,1]))
154            points = new_origin.change_points_geo_ref(points_utm,
155                                                      points_georeference)
156            geo_ref = new_origin
157
158        # At this stage I need a georef and points
159        # the points are relative to the georef
160        geo_ref.write_NetCDF(outfile)
161
162        x = points[:,0]
163        y = points[:,1]
164        z = outfile.variables['elevation'][:]
165
166        if verbose:
167            log.critical('------------------------------------------------')
168            log.critical('More Statistics:')
169            log.critical('  Extent (/lon):')
170            log.critical('    x in [%f, %f], len(lat) == %d'
171                         % (min(x), max(x), len(x)))
172            log.critical('    y in [%f, %f], len(lon) == %d'
173                         % (min(y), max(y), len(y)))
174            log.critical('    z in [%f, %f], len(z) == %d'
175                         % (min(elevation), max(elevation), len(elevation)))
176            log.critical('geo_ref: %s' % str(geo_ref))
177            log.critical('------------------------------------------------')
178
179        z = resize(bath_grid,outfile.variables['elevation'][:].shape)
180        outfile.variables['x'][:] = points[:,0] #- geo_ref.get_xllcorner()
181        outfile.variables['y'][:] = points[:,1] #- geo_ref.get_yllcorner()
182        #outfile.variables['z'][:] = elevation
183        outfile.variables['elevation'][:] = elevation  #FIXME HACK4
184
185        # This updates the _range values
186        q = 'elevation'
187        outfile.variables[q + Write_sts.RANGE][0] = min(elevation)
188        outfile.variables[q + Write_sts.RANGE][1] = max(elevation)
189
190    ##
191    # @brief Store quantity data in the underlying file.
192    # @param outfile
193    # @param sts_precision
194    # @param slice_index
195    # @param time
196    # @param verboseTrue if this function is to be verbose.
197    # @param **quant Extra keyword args.
198    def store_quantities(self, outfile, sts_precision=num.float32,
199                         slice_index=None, time=None,
200                         verbose=False, **quant):
201        """Write the quantity info.
202
203        **quant is extra keyword arguments passed in. These must be
204          the sts quantities, currently; stage.
205
206        if the time array is already been built, use the slice_index
207        to specify the index.
208
209        Otherwise, use time to increase the time dimension
210
211        Maybe make this general, but the viewer assumes these quantities,
212        so maybe we don't want it general - unless the viewer is general
213
214        precondition:
215            triangulation and header have been called.
216        """
217
218        if time is not None:
219            file_time = outfile.variables['time']
220            slice_index = len(file_time)
221            file_time[slice_index] = time
222
223        # Write the conserved quantities from Domain.
224        # Typically stage, xmomentum, ymomentum
225        # other quantities will be ignored, silently.
226        # Also write the ranges: stage_range
227        for q in Write_sts.sts_quantities:
228            if not quant.has_key(q):
229                msg = 'STS file can not write quantity %s' % q
230                raise NewQuantity, msg
231            else:
232                q_values = quant[q]
233                outfile.variables[q][slice_index] = \
234                                q_values.astype(sts_precision)
235
236                # This updates the _range values
237                q_range = outfile.variables[q + Write_sts.RANGE][:]
238                q_values_min = num.min(q_values)
239                if q_values_min < q_range[0]:
240                    outfile.variables[q + Write_sts.RANGE][0] = q_values_min
241                q_values_max = num.max(q_values)
242                if q_values_max > q_range[1]:
243                    outfile.variables[q + Write_sts.RANGE][1] = q_values_max
244
245
246
247    def write_dynamic_quantities(self, outfile, quantities,
248                    times, precis = netcdf_float32, verbose = False):   
249        """
250            Write out given quantities to file.
251        """
252        for q in quantities:
253            outfile.createVariable(q, precis, ('number_of_timesteps',
254                                                      'number_of_points'))
255            outfile.createVariable(q + Write_sts.RANGE, precis,
256                                   ('numbers_in_range',))
257
258            # Initialise ranges with small and large sentinels.
259            # If this was in pure Python we could have used None sensibly
260            outfile.variables[q+Write_sts.RANGE][0] = max_float  # Min
261            outfile.variables[q+Write_sts.RANGE][1] = -max_float # Max
262
263        # Doing sts_precision instead of Float gives cast errors.
264        outfile.createVariable('time', netcdf_float, ('number_of_timesteps',))
265
266        if isinstance(times, (list, num.ndarray)):
267            outfile.variables['time'][:] = times    # Store time relative
268
269        if verbose:
270            log.critical('------------------------------------------------')
271            log.critical('Statistics:')
272            log.critical('    t in [%f, %f], len(t) == %d'
273                         % (num.min(times), num.max(times), len(times.flat)))
274
275
276
277##
278# @brief Create a list of points defining a boundary from an STS file.
279# @param stsname Stem of path to the STS file to read.
280# @return A list of boundary points.
281def create_sts_boundary(sts_filename):
282    """Create a list of points defining a boundary from an STS file.
283
284    Create boundary segments from .sts file. Points can be stored in
285    arbitrary order within the .sts file. The order in which the .sts points
286    make up the boundary are given in order.txt file
287
288    FIXME:
289    Point coordinates are stored in relative eastings and northings.
290    But boundary is produced in absolute coordinates
291    """
292
293    if sts_filename.endswith('.sts'):
294        stsname_postfixed = sts_filename
295    else:
296        stsname_postfixed = sts_filename + '.sts'
297
298    try:
299        fid = NetCDFFile(stsname_postfixed, netcdf_mode_r)
300    except IOError:
301        msg = 'Cannot open %s' % stsname_postfixed
302        raise IOError(msg)
303
304    xllcorner = fid.xllcorner[0]
305    yllcorner = fid.yllcorner[0]
306
307    #Points stored in sts file are normalised to [xllcorner,yllcorner] but
308    #we cannot assume that boundary polygon will be. At least the
309    #additional points specified by the user after this function is called
310    x = fid.variables['x'][:] + xllcorner
311    y = fid.variables['y'][:] + yllcorner
312
313    x = num.reshape(x, (len(x), 1))
314    y = num.reshape(y, (len(y), 1))
315    sts_points = num.concatenate((x,y), axis=1)
316
317    return sts_points.tolist()
318
319
320
Note: See TracBrowser for help on using the repository browser.