source: inundation/pyvolution/util.py @ 1986

Last change on this file since 1986 was 1986, checked in by ole, 19 years ago

Incorporated caching option into file_function API

File size: 18.1 KB
Line 
1"""This module contains various auxiliary function used by pyvolution.
2
3It is also a clearing house for functions that may later earn a module
4of their own.
5"""
6
7
8#FIXME: Deprecate this shortcut
9from utilities.numerical_tools import ensure_numeric
10
11import utilities.polygon
12from warnings import warn
13
14def point_on_line(*args, **kwargs):
15    """Temporary Interface to new location"""
16
17    msg = 'point_on_line has moved from util.py.  '
18    msg += 'Please use "from utilities.polygon import point_on_line"'
19    warn(msg, DeprecationWarning) 
20
21    return utilities.polygon.point_on_line(*args, **kwargs)     
22   
23def inside_polygon(*args, **kwargs):
24    """Temporary Interface to new location"""
25
26    print 'inside_polygon has moved from util.py.  ',
27    print 'Please use "from utilities.polygon import inside_polygon"'
28
29    return utilities.polygon.inside_polygon(*args, **kwargs)   
30   
31def outside_polygon(*args, **kwargs):
32    """Temporary Interface to new location"""
33
34    print 'outside_polygon has moved from util.py.  ',
35    print 'Please use "from utilities.polygon import outside_polygon"'
36
37    return utilities.polygon.outside_polygon(*args, **kwargs)   
38
39
40def separate_points_by_polygon(*args, **kwargs):
41    """Temporary Interface to new location"""
42
43    print 'separate_points_by_polygon has moved from util.py.  ',
44    print 'Please use "from utilities.polygon import separate_points_by_polygon"'
45
46    return utilities.polygon.separate_points_by_polygon(*args, **kwargs)   
47
48
49from utilities.polygon import Polygon_function #No warning
50
51
52def read_polygon(*args, **kwargs):
53    """Temporary Interface to new location"""
54
55    print 'read_polygon has moved from util.py.  ',
56    print 'Please use "from utilities.polygon import read_polygon"'
57
58    return utilities.polygon.read_polygon(*args, **kwargs)   
59
60
61def populate_polygon(*args, **kwargs):
62    """Temporary Interface to new location"""
63
64    print 'populate_polygon has moved from util.py.  ',
65    print 'Please use "from utilities.polygon import populate_polygon"'
66
67    return utilities.polygon.populate_polygon(*args, **kwargs)   
68
69def read_xya(filename):
70    """Read simple xya file, possibly with a header in the first line, just
71    x y [attributes]
72    separated by whitespace
73
74    Format can be either ASCII or NetCDF
75
76    Return list of points, list of attributes and
77    attribute names if present otherwise None
78    """
79    print "read_xya is obsolete.  Use import_points_file from load_mesh.loadASCII"
80    #FIXME: Probably obsoleted by similar function in load_ASCII
81    #FIXME: Name read_data_points (or see 'load' in pylab)
82
83   
84    from load_mesh.loadASCII import import_points_file
85
86    points_dict = import_points_file(filename)
87    return points_dict['pointlist'], points_dict['attributelist']
88
89   
90
91
92#Normal util stuff
93
94def angle(v):
95    """Compute angle between e1 (the unit vector in the x-direction)
96    and the specified vector
97    """
98
99    from math import acos, pi, sqrt
100    from Numeric import sum, array
101
102    l = sqrt( sum (array(v)**2))
103    v1 = v[0]/l
104    v2 = v[1]/l
105
106
107    theta = acos(v1)
108   
109    #try:
110    #   theta = acos(v1)
111    #except ValueError, e:
112    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
113    #
114    #    #FIXME, hack to avoid acos(1.0) Value error
115    #    # why is it happening?
116    #    # is it catching something we should avoid?
117    #    #FIXME (Ole): When did this happen? We need a unit test to
118    #    #reveal this condition or otherwise remove the hack
119    #   
120    #    s = 1e-6
121    #    if (v1+s >  1.0) and (v1-s < 1.0) :
122    #        theta = 0.0
123    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
124    #        theta = 3.1415926535897931
125    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
126    #          %(v1, theta)
127
128    if v2 < 0:
129        #Quadrant 3 or 4
130        theta = 2*pi-theta
131
132    return theta
133
134
135def anglediff(v0, v1):
136    """Compute difference between angle of vector x0, y0 and x1, y1.
137    This is used for determining the ordering of vertices,
138    e.g. for checking if they are counter clockwise.
139
140    Always return a positive value
141    """
142
143    from math import pi
144
145    a0 = angle(v0)
146    a1 = angle(v1)
147
148    #Ensure that difference will be positive
149    if a0 < a1:
150        a0 += 2*pi
151
152    return a0-a1
153
154
155def mean(x):
156    from Numeric import sum
157    return sum(x)/len(x)
158
159
160
161
162def file_function(filename,
163                  domain = None,
164                  quantities = None,
165                  interpolation_points = None,
166                  verbose = False,
167                  use_cache = False):
168    """Read time history of spatial data from NetCDF file and return
169    a callable object.
170
171    Input variables:
172   
173    filename - Name of sww or tms file
174       
175       If the file has extension 'sww' then it is assumed to be spatio-temporal
176       or temporal and the callable object will have the form f(t,x,y) or f(t)
177       depending on whether the file contains spatial data
178
179       If the file has extension 'tms' then it is assumed to be temporal only
180       and the callable object will have the form f(t)
181
182       Either form will return interpolated values based on the input file
183       using the underlying interpolation_function.
184
185    domain - Associated domain object   
186       If domain is specified, model time (domain.starttime)
187       will be checked and possibly modified.
188   
189       All times are assumed to be in UTC
190       
191       All spatial information is assumed to be in absolute UTM coordinates.
192
193    quantities - the name of the quantity to be interpolated or a
194                 list of quantity names. The resulting function will return
195                 a tuple of values - one for each quantity 
196
197    interpolation_points - list of absolute UTM coordinates for points at
198    which values are sought
199   
200    use_cache: True means that caching of intermediate result of
201               Interpolation_function is attempted
202
203   
204    See Interpolation function for further documentation
205    """
206
207
208    if use_cache is True:
209        try:
210            from caching import cache
211        except:
212            msg = 'Caching was requested, but caching module'+\
213                  'could not be imported'
214            raise msg
215
216
217        f = cache(_file_function,
218                  filename,
219                  {'domain': domain,
220                   'quantities': quantities,
221                   'interpolation_points': interpolation_points,
222                   'verbose': verbose},
223                  dependencies = [filename],
224                  compression = False,
225                  verbose = verbose)
226        #FIXME (Ole): Pass cache arguments, such as compression, in some sort of
227        #structure
228       
229    else:
230        f = _file_function(filename,
231                           domain,
232                           quantities,
233                           interpolation_points,
234                           verbose)           
235
236    return f
237
238
239
240def _file_function(filename,
241                   domain = None,
242                   quantities = None,
243                   interpolation_points = None,
244                   verbose = False):
245    """Internal function
246   
247    See file_function for documentatiton
248    """
249   
250
251    #FIXME (OLE): Should check origin of domain against that of file
252    #In fact, this is where origin should be converted to that of domain
253    #Also, check that file covers domain fully.
254    #If we use the suggested Point_set class for interpolation points
255    #here it would be easier
256
257    #Take into account:
258    #- domain's georef
259    #- sww file's georef
260    #- interpolation points as absolute UTM coordinates
261
262
263    assert type(filename) == type(''),\
264               'First argument to File_function must be a string'
265
266    try:
267        fid = open(filename)
268    except Exception, e:
269        msg = 'File "%s" could not be opened: Error="%s"'\
270                  %(filename, e)
271        raise msg
272
273    line = fid.readline()
274    fid.close()
275
276    if quantities is None: 
277        if domain is not None:
278            quantities = domain.conserved_quantities
279
280
281
282    if line[:3] == 'CDF':
283        return get_netcdf_file_function(filename, domain, quantities,
284                                        interpolation_points,
285                                        verbose = verbose)
286    else:
287        raise 'Must be a NetCDF File'
288
289
290
291def get_netcdf_file_function(filename,
292                             domain=None,
293                             quantity_names=None,
294                             interpolation_points=None,
295                             verbose = False):
296    """Read time history of spatial data from NetCDF sww file and
297    return a callable object f(t,x,y)
298    which will return interpolated values based on the input file.
299
300    If domain is specified, model time (domain.starttime)
301    will be checked and possibly modified
302   
303    All times are assumed to be in UTC
304
305    See Interpolation function for further documetation
306   
307    """
308   
309   
310    #FIXME: Check that model origin is the same as file's origin
311    #(both in UTM coordinates)
312    #If not - modify those from file to match domain
313    #Take this code from e.g. dem2pts in data_manager.py
314    #FIXME: Use geo_reference to read and write xllcorner...
315       
316
317    #FIXME: Maybe move caching out to this level rather than at the
318    #Interpolation_function level (below)
319
320    import time, calendar, types
321    from config import time_format
322    from Scientific.IO.NetCDF import NetCDFFile
323    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
324    from util import ensure_numeric   
325
326    #Open NetCDF file
327    if verbose: print 'Reading', filename
328    fid = NetCDFFile(filename, 'r')
329
330    if type(quantity_names) == types.StringType:
331        quantity_names = [quantity_names]       
332   
333    if quantity_names is None or len(quantity_names) < 1:
334        #If no quantities are specified get quantities from file
335        #x, y, time are assumed as the independent variables so
336        #they are excluded as potentiol quantities
337        quantity_names = []
338        for name in fid.variables.keys():
339            if name not in ['x', 'y', 'time']:
340                quantity_names.append(name)
341
342    if len(quantity_names) < 1:               
343        msg = 'ERROR: At least one independent value must be specified'
344        raise msg
345
346
347    if interpolation_points is not None:
348        interpolation_points = ensure_numeric(interpolation_points, Float)
349
350
351
352    #Now assert that requested quantitites (and the independent ones)
353    #are present in file
354    missing = []
355    for quantity in ['time'] + quantity_names:
356        if not fid.variables.has_key(quantity):
357            missing.append(quantity)
358
359    if len(missing) > 0:
360        msg = 'Quantities %s could not be found in file %s'\
361              %(str(missing), filename)
362        fid.close()
363        raise msg
364
365    #Decide whether this data has a spatial dimension
366    spatial = True
367    for quantity in ['x', 'y']:
368        if not fid.variables.has_key(quantity):
369            spatial = False
370
371    if filename[-3:] == 'tms' and spatial is True:
372        msg = 'Files of type tms must not contain spatial information'
373        raise msg
374
375    if filename[-3:] == 'sww' and spatial is False:
376        msg = 'Files of type sww must contain spatial information'       
377        raise msg       
378
379    #Get first timestep
380    try:
381        starttime = fid.starttime[0]
382    except ValueError:
383        msg = 'Could not read starttime from file %s' %filename
384        raise msg
385
386    #Get variables
387    if verbose: print 'Get variables'   
388    time = fid.variables['time'][:]   
389
390    if spatial:
391        #Get origin
392        xllcorner = fid.xllcorner[0]
393        yllcorner = fid.yllcorner[0]
394        zone = fid.zone[0]       
395
396        x = fid.variables['x'][:]
397        y = fid.variables['y'][:]
398        triangles = fid.variables['volumes'][:]
399
400        x = reshape(x, (len(x),1))
401        y = reshape(y, (len(y),1))
402        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
403
404        if interpolation_points is not None:
405            #Adjust for georef
406            interpolation_points[:,0] -= xllcorner
407            interpolation_points[:,1] -= yllcorner       
408       
409
410
411
412    if domain is not None:
413        #Update domain.startime if it is *earlier* than starttime
414        if starttime > domain.starttime:
415            msg = 'WARNING: Start time as specified in domain (%f)'\
416                  %domain.starttime
417            msg += ' is earlier than the starttime of file %s (%f).'\
418                     %(filename, starttime)
419            msg += ' Modifying domain starttime accordingly.'
420           
421            if verbose: print msg
422            domain.starttime = starttime #Modifying model time
423            if verbose: print 'Domain starttime is now set to %f'\
424               %domain.starttime
425
426
427        #If domain.startime is *later* than starttime,
428        #move time back - relative to domain's time
429        if domain.starttime > starttime:
430            time = time - domain.starttime + starttime
431
432        #FIXME Use method in geo to reconcile
433        #if spatial:
434        #assert domain.geo_reference.xllcorner == xllcorner
435        #assert domain.geo_reference.yllcorner == yllcorner
436        #assert domain.geo_reference.zone == zone       
437       
438    if verbose:
439        print 'File_function data obtained from: %s' %filename
440        print '  References:'
441        #print '    Datum ....' #FIXME
442        if spatial:
443            print '    Lower left corner: [%f, %f]'\
444                  %(xllcorner, yllcorner)
445        print '    Start time:   %f' %starttime               
446       
447   
448    #Produce values for desired data points at
449    #each timestep
450
451    quantities = {}
452    for i, name in enumerate(quantity_names):
453        quantities[name] = fid.variables[name][:]
454    fid.close()
455   
456
457    from least_squares import Interpolation_function
458
459    if not spatial:
460        vertex_coordinates = triangles = interpolation_points = None         
461
462
463    return Interpolation_function(time,
464                                  quantities,
465                                  quantity_names,
466                                  vertex_coordinates,
467                                  triangles,
468                                  interpolation_points,
469                                  verbose = verbose)
470
471
472
473
474def tsh2sww(filename, verbose=True):
475    """
476    to check if a tsh/msh file 'looks' good.
477    """
478
479    #FIXME Move to data_manager
480    from shallow_water import Domain
481    from pmesh2domain import pmesh_to_domain_instance
482    import time, os
483    from data_manager import get_dataobject
484    from os import sep, path
485    #from util import mean
486
487    if verbose == True:print 'Creating domain from', filename
488    domain = pmesh_to_domain_instance(filename, Domain)
489    if verbose == True:print "Number of triangles = ", len(domain)
490
491    domain.smooth = True
492    domain.format = 'sww'   #Native netcdf visualisation format
493    file_path, filename = path.split(filename)
494    filename, ext = path.splitext(filename)
495    domain.filename = filename
496    domain.reduction = mean
497    if verbose == True:print "file_path",file_path
498    if file_path == "":file_path = "."
499    domain.set_datadir(file_path)
500
501    if verbose == True:
502        print "Output written to " + domain.get_datadir() + sep + \
503              domain.filename + "." + domain.format
504    sww = get_dataobject(domain)
505    sww.store_connectivity()
506    sww.store_timestep('stage')
507
508
509
510
511def multiple_replace(text, dictionary):
512    """Multiple replace of words in text
513
514    text:       String to be modified
515    dictionary: Mapping of words that are to be substituted
516
517    Python Cookbook 3.14 page 88 and page 90
518    """
519
520    import re
521   
522    #Create a regular expression from all of the dictionary keys
523    #matching only entire words
524    regex = re.compile(r'\b'+ \
525                       r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \
526                       r'\b' )
527
528    #For each match, lookup the corresponding value in the dictionary
529    return regex.sub(lambda match: dictionary[match.group(0)], text)
530
531
532
533
534def apply_expression_to_dictionary(expression, dictionary):#dictionary):
535    """Apply arbitrary expression to values of dictionary
536
537    Given an expression in terms of the keys, replace key by the
538    corresponding values and evaluate.
539   
540
541    expression: Arbitrary, e.g. arithmetric, expression relating keys
542                from dictionary.
543
544    dictionary: Mapping between symbols in expression and objects that
545                will be evaluated by expression.
546                Values in dictionary must support operators given in
547                expression e.g. by overloading               
548    """
549
550    import types
551    import re
552
553    assert isinstance(expression, basestring)
554    assert type(dictionary) == types.DictType
555
556    #Convert dictionary values to textual representations suitable for eval
557    D = {}
558    for key in dictionary:
559        D[key] = 'dictionary["%s"]' %key
560
561    #Perform substitution of variables   
562    expression = multiple_replace(expression, D)
563
564    #Evaluate and return
565    try:
566        return eval(expression)
567    except NameError, e:
568        msg = 'Expression "%s" could not be evaluated: %s' %(expression, e)
569        raise NameError, msg
570   
571
572
573
574
575
576
577
578####################################################################
579#Python versions of function that are also implemented in util_gateway.c
580#
581
582def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
583    """
584    """
585
586    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
587    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
588    a /= det
589
590    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
591    b /= det
592
593    return a, b
594
595
596def gradient2_python(x0, y0, x1, y1, q0, q1):
597    """Compute radient based on two points and enforce zero gradient
598    in the direction orthogonal to (x1-x0), (y1-y0)
599    """
600
601    #Old code
602    #det = x0*y1 - x1*y0
603    #if det != 0.0:
604    #    a = (y1*q0 - y0*q1)/det
605    #    b = (x0*q1 - x1*q0)/det
606
607    #Correct code (ON)
608    det = (x1-x0)**2 + (y1-y0)**2
609    if det != 0.0:
610        a = (x1-x0)*(q1-q0)/det
611        b = (y1-y0)*(q1-q0)/det
612       
613    return a, b       
614
615
616##############################################
617#Initialise module
618
619from utilities import compile
620if compile.can_use_C_extension('util_ext.c'):
621    from util_ext import gradient, gradient2#, point_on_line
622    #separate_points_by_polygon = separate_points_by_polygon_c
623else:
624    gradient = gradient_python
625    gradient2 = gradient2_python   
626
627
628if __name__ == "__main__":
629    pass
630
631
632
633
Note: See TracBrowser for help on using the repository browser.