source: anuga_core/source/anuga/abstract_2d_finite_volumes/file_function.py @ 7690

Last change on this file since 7690 was 7690, checked in by hudson, 15 years ago

Check for complex polygons raises an exception.
File size: 19.0 KB
Line 
1#!/usr/bin/env python
2"""File function
3Takes a file as input, and returns it as a mathematical function.
4For example, you can load an arbitrary 2D heightfield mesh, and treat it as a function as so:
5
6F = file_function('my_mesh.sww', ...)
7evaluated_point = F(x, y)
8
9Values will be interpolated across the surface of the mesh. Holes in the mesh have an undefined value.
10
11"""
12
13import numpy as num
14
15from anuga.geospatial_data.geospatial_data import ensure_absolute
16from Scientific.IO.NetCDF import NetCDFFile
17from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a
18from anuga.utilities.numerical_tools import ensure_numeric
19
20import anuga.utilities.log as log
21
22
23
24##
25# @brief Read time history of data from NetCDF file, return callable object.
26# @param filename  Name of .sww or .tms file.
27# @param domain Associated domain object.
28# @param quantities Name of quantity to be interpolated or a list of names.
29# @param interpolation_points List of absolute UTM coordinates for points
30#                             (N x 2) or geospatial object or
31#                             points file name at which values are sought.
32# @param time_thinning
33# @param verbose True if this function is to be verbose.
34# @param use_cache True means that caching of intermediate result is attempted.
35# @param boundary_polygon
36# @param output_centroids if True, data for the centroid of the triangle will be output
37# @return A callable object.
38def file_function(filename,
39                  domain=None,
40                  quantities=None,
41                  interpolation_points=None,
42                  time_thinning=1,
43                  time_limit=None,
44                  verbose=False,
45                  use_cache=False,
46                  boundary_polygon=None,
47                  output_centroids=False):
48    """Read time history of spatial data from NetCDF file and return
49    a callable object.
50
51    Input variables:
52   
53    filename - Name of sww, tms or sts file
54       
55       If the file has extension 'sww' then it is assumed to be spatio-temporal
56       or temporal and the callable object will have the form f(t,x,y) or f(t)
57       depending on whether the file contains spatial data
58
59       If the file has extension 'tms' then it is assumed to be temporal only
60       and the callable object will have the form f(t)
61
62       Either form will return interpolated values based on the input file
63       using the underlying interpolation_function.
64
65    domain - Associated domain object   
66       If domain is specified, model time (domain.starttime)
67       will be checked and possibly modified.
68   
69       All times are assumed to be in UTC
70       
71       All spatial information is assumed to be in absolute UTM coordinates.
72
73    quantities - the name of the quantity to be interpolated or a
74                 list of quantity names. The resulting function will return
75                 a tuple of values - one for each quantity
76                 If quantities are None, the default quantities are
77                 ['stage', 'xmomentum', 'ymomentum']
78                 
79
80    interpolation_points - list of absolute UTM coordinates for points (N x 2)
81    or geospatial object or points file name at which values are sought
82
83    time_thinning -
84
85    verbose -
86
87    use_cache: True means that caching of intermediate result of
88               Interpolation_function is attempted
89
90    boundary_polygon -
91
92   
93    See Interpolation function in anuga.fit_interpolate.interpolation for
94    further documentation
95    """
96
97    # FIXME (OLE): Should check origin of domain against that of file
98    # In fact, this is where origin should be converted to that of domain
99    # Also, check that file covers domain fully.
100
101    # Take into account:
102    # - domain's georef
103    # - sww file's georef
104    # - interpolation points as absolute UTM coordinates
105
106    if quantities is None:
107        if verbose:
108            msg = 'Quantities specified in file_function are None,'
109            msg += ' so I will use stage, xmomentum, and ymomentum in that order'
110            log.critical(msg)
111        quantities = ['stage', 'xmomentum', 'ymomentum']
112
113    # Use domain's startime if available
114    if domain is not None:   
115        domain_starttime = domain.get_starttime()
116    else:
117        domain_starttime = None
118
119    # Build arguments and keyword arguments for use with caching or apply.
120    args = (filename,)
121
122    # FIXME (Ole): Caching this function will not work well
123    # if domain is passed in as instances change hash code.
124    # Instead we pass in those attributes that are needed (and return them
125    # if modified)
126    kwargs = {'quantities': quantities,
127              'interpolation_points': interpolation_points,
128              'domain_starttime': domain_starttime,
129              'time_thinning': time_thinning,     
130              'time_limit': time_limit,                                 
131              'verbose': verbose,
132              'boundary_polygon': boundary_polygon,
133              'output_centroids': output_centroids}
134
135    # Call underlying engine with or without caching
136    if use_cache is True:
137        try:
138            from caching import cache
139        except:
140            msg = 'Caching was requested, but caching module'+\
141                  'could not be imported'
142            raise msg
143
144        f, starttime = cache(_file_function,
145                             args, kwargs,
146                             dependencies=[filename],
147                             compression=False,                 
148                             verbose=verbose)
149    else:
150        f, starttime = apply(_file_function,
151                             args, kwargs)
152
153    #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
154    #structure
155
156    f.starttime = starttime
157    f.filename = filename
158   
159    if domain is not None:
160        #Update domain.startime if it is *earlier* than starttime from file
161        if starttime > domain.starttime:
162            msg = 'WARNING: Start time as specified in domain (%f)'\
163                  %domain.starttime
164            msg += ' is earlier than the starttime of file %s (%f).'\
165                     %(filename, starttime)
166            msg += ' Modifying domain starttime accordingly.'
167           
168            if verbose: log.critical(msg)
169
170            domain.set_starttime(starttime) #Modifying model time
171
172            if verbose: log.critical('Domain starttime is now set to %f'
173                                     % domain.starttime)
174    return f
175
176
177##
178# @brief ??
179# @param filename  Name of .sww or .tms file.
180# @param domain Associated domain object.
181# @param quantities Name of quantity to be interpolated or a list of names.
182# @param interpolation_points List of absolute UTM coordinates for points
183#                             (N x 2) or geospatial object or
184#                             points file name at which values are sought.
185# @param time_thinning
186# @param verbose True if this function is to be verbose.
187# @param use_cache True means that caching of intermediate result is attempted.
188# @param boundary_polygon
189def _file_function(filename,
190                   quantities=None,
191                   interpolation_points=None,
192                   domain_starttime=None,
193                   time_thinning=1,
194                   time_limit=None,
195                   verbose=False,
196                   boundary_polygon=None,
197                   output_centroids=False):
198    """Internal function
199   
200    See file_function for documentatiton
201    """
202
203    assert type(filename) == type(''),\
204               'First argument to File_function must be a string'
205
206    try:
207        fid = open(filename)
208    except Exception, e:
209        msg = 'File "%s" could not be opened: Error="%s"' % (filename, e)
210        raise msg
211
212    # read first line of file, guess file type
213    line = fid.readline()
214    fid.close()
215
216    if line[:3] == 'CDF':
217        return get_netcdf_file_function(filename,
218                                        quantities,
219                                        interpolation_points,
220                                        domain_starttime,
221                                        time_thinning=time_thinning,
222                                        time_limit=time_limit,
223                                        verbose=verbose,
224                                        boundary_polygon=boundary_polygon,
225                                        output_centroids=output_centroids)
226    else:
227        # FIXME (Ole): Could add csv file here to address Ted Rigby's
228        # suggestion about reading hydrographs.
229        # This may also deal with the gist of ticket:289
230        raise 'Must be a NetCDF File'
231
232
233##
234# @brief ??
235# @param filename  Name of .sww or .tms file.
236# @param quantity_names Name of quantity to be interpolated or a list of names.
237# @param interpolation_points List of absolute UTM coordinates for points
238#                             (N x 2) or geospatial object or
239#                             points file name at which values are sought.
240# @param domain_starttime Start time from domain object.
241# @param time_thinning ??
242# @param verbose True if this function is to be verbose.
243# @param boundary_polygon ??
244# @return A callable object.
245def get_netcdf_file_function(filename,
246                             quantity_names=None,
247                             interpolation_points=None,
248                             domain_starttime=None,                           
249                             time_thinning=1,                 
250                             time_limit=None,           
251                             verbose=False,
252                             boundary_polygon=None,
253                             output_centroids=False):
254    """Read time history of spatial data from NetCDF sww file and
255    return a callable object f(t,x,y)
256    which will return interpolated values based on the input file.
257
258    Model time (domain_starttime)
259    will be checked, possibly modified and returned
260   
261    All times are assumed to be in UTC
262
263    See Interpolation function for further documetation
264    """
265
266    # FIXME: Check that model origin is the same as file's origin
267    # (both in UTM coordinates)
268    # If not - modify those from file to match domain
269    # (origin should be passed in)
270    # Take this code from e.g. dem2pts in data_manager.py
271    # FIXME: Use geo_reference to read and write xllcorner...
272
273    import time, calendar, types
274    from anuga.config import time_format
275
276    # Open NetCDF file
277    if verbose: log.critical('Reading %s' % filename)
278
279    fid = NetCDFFile(filename, netcdf_mode_r)
280
281    if type(quantity_names) == types.StringType:
282        quantity_names = [quantity_names]       
283
284    if quantity_names is None or len(quantity_names) < 1:
285        msg = 'No quantities are specified in file_function'
286        raise Exception, msg
287 
288    if interpolation_points is not None:
289        interpolation_points = ensure_absolute(interpolation_points)
290        msg = 'Points must by N x 2. I got %d' % interpolation_points.shape[1]
291        assert interpolation_points.shape[1] == 2, msg
292
293    # Now assert that requested quantitites (and the independent ones)
294    # are present in file
295    missing = []
296    for quantity in ['time'] + quantity_names:
297        if not fid.variables.has_key(quantity):
298            missing.append(quantity)
299
300    if len(missing) > 0:
301        msg = 'Quantities %s could not be found in file %s'\
302              % (str(missing), filename)
303        fid.close()
304        raise Exception, msg
305
306    # Decide whether this data has a spatial dimension
307    spatial = True
308    for quantity in ['x', 'y']:
309        if not fid.variables.has_key(quantity):
310            spatial = False
311
312    if filename[-3:] == 'tms' and spatial is True:
313        msg = 'Files of type tms must not contain spatial  information'
314        raise msg
315
316    if filename[-3:] == 'sww' and spatial is False:
317        msg = 'Files of type sww must contain spatial information'       
318        raise msg
319
320    if filename[-3:] == 'sts' and spatial is False:
321        #What if mux file only contains one point
322        msg = 'Files of type sts must contain spatial information'       
323        raise msg
324
325    if filename[-3:] == 'sts' and boundary_polygon is None:
326        #What if mux file only contains one point
327        msg = 'Files of type sts require boundary polygon'       
328        raise msg
329
330    # Get first timestep
331    try:
332        starttime = fid.starttime[0]
333    except ValueError:
334        msg = 'Could not read starttime from file %s' % filename
335        raise msg
336
337    # Get variables
338    # if verbose: log.critical('Get variables'    )
339    time = fid.variables['time'][:]   
340    # FIXME(Ole): Is time monotoneous?
341
342    # Apply time limit if requested
343    upper_time_index = len(time)   
344    msg = 'Time vector obtained from file %s has length 0' % filename
345    assert upper_time_index > 0, msg
346   
347    if time_limit is not None:
348        # Adjust given time limit to given start time
349        time_limit = time_limit - starttime
350
351
352        # Find limit point
353        for i, t in enumerate(time):
354            if t > time_limit:
355                upper_time_index = i
356                break
357               
358        msg = 'Time vector is zero. Requested time limit is %f' % time_limit
359        assert upper_time_index > 0, msg
360
361        if time_limit < time[-1] and verbose is True:
362            log.critical('Limited time vector from %.2fs to %.2fs'
363                         % (time[-1], time_limit))
364
365    time = time[:upper_time_index]
366
367
368   
369   
370    # Get time independent stuff
371    if spatial:
372        # Get origin
373        xllcorner = fid.xllcorner[0]
374        yllcorner = fid.yllcorner[0]
375        zone = fid.zone[0]       
376
377        x = fid.variables['x'][:]
378        y = fid.variables['y'][:]
379        if filename.endswith('sww'):
380            triangles = fid.variables['volumes'][:]
381
382        x = num.reshape(x, (len(x),1))
383        y = num.reshape(y, (len(y),1))
384        vertex_coordinates = num.concatenate((x,y), axis=1) #m x 2 array
385
386        if boundary_polygon is not None:
387            # Remove sts points that do not lie on boundary
388            # FIXME(Ole): Why don't we just remove such points from the list of points and associated data?
389            # I am actually convinced we can get rid of neighbour_gauge_id altogether as the sts file is produced using the ordering file.
390            # All sts points are therefore always present in the boundary. In fact, they *define* parts of the boundary.
391            boundary_polygon=ensure_numeric(boundary_polygon)
392            boundary_polygon[:,0] -= xllcorner
393            boundary_polygon[:,1] -= yllcorner
394            temp=[]
395            boundary_id=[]
396            gauge_id=[]
397            for i in range(len(boundary_polygon)):
398                for j in range(len(x)):
399                    if num.allclose(vertex_coordinates[j],
400                                    boundary_polygon[i], 1e-4):
401                        #FIXME:
402                        #currently gauges lat and long is stored as float and
403                        #then cast to double. This cuases slight repositioning
404                        #of vertex_coordinates.
405                        temp.append(boundary_polygon[i])
406                        gauge_id.append(j)
407                        boundary_id.append(i)
408                        break
409            gauge_neighbour_id=[]
410            for i in range(len(boundary_id)-1):
411                if boundary_id[i]+1==boundary_id[i+1]:
412                    gauge_neighbour_id.append(i+1)
413                else:
414                    gauge_neighbour_id.append(-1)
415            if boundary_id[len(boundary_id)-1]==len(boundary_polygon)-1 \
416               and boundary_id[0]==0:
417                gauge_neighbour_id.append(0)
418            else:
419                gauge_neighbour_id.append(-1)
420            gauge_neighbour_id=ensure_numeric(gauge_neighbour_id)
421
422            if len(num.compress(gauge_neighbour_id>=0,gauge_neighbour_id)) \
423               != len(temp)-1:
424                msg='incorrect number of segments'
425                raise msg
426            vertex_coordinates=ensure_numeric(temp)
427            if len(vertex_coordinates)==0:
428                msg = 'None of the sts gauges fall on the boundary'
429                raise msg
430        else:
431            gauge_neighbour_id=None
432
433        if interpolation_points is not None:
434            # Adjust for georef
435            interpolation_points[:,0] -= xllcorner
436            interpolation_points[:,1] -= yllcorner       
437    else:
438        gauge_neighbour_id=None
439       
440    if domain_starttime is not None:
441        # If domain_startime is *later* than starttime,
442        # move time back - relative to domain's time
443        if domain_starttime > starttime:
444            time = time - domain_starttime + starttime
445
446        # FIXME Use method in geo to reconcile
447        # if spatial:
448        # assert domain.geo_reference.xllcorner == xllcorner
449        # assert domain.geo_reference.yllcorner == yllcorner
450        # assert domain.geo_reference.zone == zone       
451       
452    if verbose:
453        log.critical('File_function data obtained from: %s' % filename)
454        log.critical('  References:')
455        if spatial:
456            log.critical('    Lower left corner: [%f, %f]'
457                         % (xllcorner, yllcorner))
458        log.critical('    Start time:   %f' % starttime)
459       
460   
461    # Produce values for desired data points at
462    # each timestep for each quantity
463    quantities = {}
464    for i, name in enumerate(quantity_names):
465        quantities[name] = fid.variables[name][:]
466        if boundary_polygon is not None:
467            #removes sts points that do not lie on boundary
468            quantities[name] = num.take(quantities[name], gauge_id, axis=1)
469           
470    # Close sww, tms or sts netcdf file         
471    fid.close()
472
473    from anuga.fit_interpolate.interpolate import Interpolation_function
474
475    if not spatial:
476        vertex_coordinates = triangles = interpolation_points = None
477    if filename[-3:] == 'sts':#added
478        triangles = None
479        #vertex coordinates is position of urs gauges
480
481    if verbose:
482        log.critical('Calling interpolation function')
483       
484    # Return Interpolation_function instance as well as
485    # starttime for use to possible modify that of domain
486    return (Interpolation_function(time,
487                                   quantities,
488                                   quantity_names,
489                                   vertex_coordinates,
490                                   triangles,
491                                   interpolation_points,
492                                   time_thinning=time_thinning,
493                                   verbose=verbose,
494                                   gauge_neighbour_id=gauge_neighbour_id,
495                                   output_centroids=output_centroids),
496            starttime)
497
498    # NOTE (Ole): Caching Interpolation function is too slow as
499    # the very long parameters need to be hashed.
Note: See TracBrowser for help on using the repository browser.