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

Last change on this file since 8143 was 8143, checked in by wilsonr, 13 years ago

Removed '@brief' comments.

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