source: inundation/pyvolution/util.py @ 1961

Last change on this file since 1961 was 1961, checked in by ole, 18 years ago

Ensure interpolation_points are converted to Float in file_function

File size: 16.6 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
160def file_function(filename,
161                  domain = None,
162                  quantities = None,
163                  interpolation_points = None, verbose = False):
164    """Read time history of spatial data from NetCDF file and return
165    a callable object.
166
167    Input variables:
168   
169    filename - Name of sww or tms file
170       
171       If the file has extension 'sww' then it is assumed to be spatio-temporal
172       or temporal and the callable object will have the form f(t,x,y) or f(t)
173       depending on whether the file contains spatial data
174
175       If the file has extension 'tms' then it is assumed to be temporal only
176       and the callable object will have the form f(t)
177
178       Either form will return interpolated values based on the input file
179       using the underlying interpolation_function.
180
181    domain - Associated domain object   
182       If domain is specified, model time (domain.starttime)
183       will be checked and possibly modified.
184   
185       All times are assumed to be in UTC
186       
187       All spatial information is assumed to be in absolute UTM coordinates.
188
189    quantities - the name of the quantity to be interpolated or a
190                 list of quantity names. The resulting function will return
191                 a tuple of values - one for each quantity 
192
193    interpolation_points - list of absolute UTM coordinates for points at
194    which values are sought
195   
196
197   
198    See Interpolation function for further documentation
199    """
200   
201
202    #FIXME (OLE): Should check origin of domain against that of file
203    #In fact, this is where origin should be converted to that of domain
204    #Also, check that file covers domain fully.
205    #If we use the suggested Point_set class for interpolation points
206    #here it would be easier
207
208    #Take into account:
209    #- domain's georef
210    #- sww file's georef
211    #- interpolation points as absolute UTM coordinates
212
213
214    assert type(filename) == type(''),\
215               'First argument to File_function must be a string'
216
217    try:
218        fid = open(filename)
219    except Exception, e:
220        msg = 'File "%s" could not be opened: Error="%s"'\
221                  %(filename, e)
222        raise msg
223
224    line = fid.readline()
225    fid.close()
226
227    if quantities is None: 
228        if domain is not None:
229            quantities = domain.conserved_quantities
230
231
232
233    if line[:3] == 'CDF':
234        return get_netcdf_file_function(filename, domain, quantities,
235                                        interpolation_points, verbose = verbose)
236    else:
237        raise 'Must be a NetCDF File'
238
239
240
241def get_netcdf_file_function(filename, domain=None, quantity_names=None,
242                             interpolation_points=None, verbose = False):
243    """Read time history of spatial data from NetCDF sww file and
244    return a callable object f(t,x,y)
245    which will return interpolated values based on the input file.
246
247    If domain is specified, model time (domain.starttime)
248    will be checked and possibly modified
249   
250    All times are assumed to be in UTC
251
252    See Interpolation function for further documetation
253   
254
255    """
256   
257    #verbose = True
258   
259    #FIXME: Check that model origin is the same as file's origin
260    #(both in UTM coordinates)
261    #If not - modify those from file to match domain
262    #Take this code from e.g. dem2pts in data_manager.py
263    #FIXME: Use geo_reference to read and write xllcorner...
264       
265
266    import time, calendar, types
267    from config import time_format
268    from Scientific.IO.NetCDF import NetCDFFile
269    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
270    from util import ensure_numeric   
271
272    #Open NetCDF file
273    if verbose: print 'Reading', filename
274    fid = NetCDFFile(filename, 'r')
275
276    if type(quantity_names) == types.StringType:
277        quantity_names = [quantity_names]       
278   
279    if quantity_names is None or len(quantity_names) < 1:
280        #If no quantities are specified get quantities from file
281        #x, y, time are assumed as the independent variables so
282        #they are excluded as potentiol quantities
283        quantity_names = []
284        for name in fid.variables.keys():
285            if name not in ['x', 'y', 'time']:
286                quantity_names.append(name)
287
288    if len(quantity_names) < 1:               
289        msg = 'ERROR: At least one independent value must be specified'
290        raise msg
291
292
293    if interpolation_points is not None:
294        interpolation_points = ensure_numeric(interpolation_points, Float)
295
296
297
298    #Now assert that requested quantitites (and the independent ones)
299    #are present in file
300    missing = []
301    for quantity in ['time'] + quantity_names:
302        if not fid.variables.has_key(quantity):
303            missing.append(quantity)
304
305    if len(missing) > 0:
306        msg = 'Quantities %s could not be found in file %s'\
307              %(str(missing), filename)
308        fid.close()
309        raise msg
310
311    #Decide whether this data has a spatial dimension
312    spatial = True
313    for quantity in ['x', 'y']:
314        if not fid.variables.has_key(quantity):
315            spatial = False
316
317    if filename[-3:] == 'tms' and spatial is True:
318        msg = 'Files of type tms must not contain spatial information'
319        raise msg
320
321    if filename[-3:] == 'sww' and spatial is False:
322        msg = 'Files of type sww must contain spatial information'       
323        raise msg       
324
325    #Get first timestep
326    try:
327        starttime = fid.starttime[0]
328    except ValueError:
329        msg = 'Could not read starttime from file %s' %filename
330        raise msg
331
332    #Get variables
333    if verbose: print 'Get variables'   
334    time = fid.variables['time'][:]   
335
336    if spatial:
337        #Get origin
338        xllcorner = fid.xllcorner[0]
339        yllcorner = fid.yllcorner[0]
340        zone = fid.zone[0]       
341
342        x = fid.variables['x'][:]
343        y = fid.variables['y'][:]
344        triangles = fid.variables['volumes'][:]
345
346        x = reshape(x, (len(x),1))
347        y = reshape(y, (len(y),1))
348        vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
349
350        if interpolation_points is not None:
351            #Adjust for georef
352            interpolation_points[:,0] -= xllcorner
353            interpolation_points[:,1] -= yllcorner       
354       
355
356
357
358    if domain is not None:
359        #Update domain.startime if it is *earlier* than starttime
360        if starttime > domain.starttime:
361            msg = 'WARNING: Start time as specified in domain (%f)'\
362                  %domain.starttime
363            msg += ' is earlier than the starttime of file %s (%f).'\
364                     %(filename, starttime)
365            msg += ' Modifying domain starttime accordingly.'
366           
367            if verbose: print msg
368            domain.starttime = starttime #Modifying model time
369            if verbose: print 'Domain starttime is now set to %f'\
370               %domain.starttime
371
372
373        #If domain.startime is *later* than starttime,
374        #move time back - relative to domain's time
375        if domain.starttime > starttime:
376            time = time - domain.starttime + starttime
377
378        #FIXME Use method in geo to reconcile
379        #if spatial:
380        #assert domain.geo_reference.xllcorner == xllcorner
381        #assert domain.geo_reference.yllcorner == yllcorner
382        #assert domain.geo_reference.zone == zone       
383       
384    if verbose:
385        print 'File_function data obtained from: %s' %filename
386        print '  References:'
387        #print '    Datum ....' #FIXME
388        if spatial:
389            print '    Lower left corner: [%f, %f]'\
390                  %(xllcorner, yllcorner)
391        print '    Start time:   %f' %starttime               
392       
393   
394    #Produce values for desired data points at
395    #each timestep
396
397    quantities = {}
398    for i, name in enumerate(quantity_names):
399        quantities[name] = fid.variables[name][:]
400    fid.close()
401   
402
403    from least_squares import Interpolation_function
404
405    if not spatial:
406        vertex_coordinates = triangles = interpolation_points = None         
407       
408    return Interpolation_function(time,
409                                  quantities,
410                                  quantity_names,
411                                  vertex_coordinates,
412                                  triangles,
413                                  interpolation_points,
414                                  verbose = verbose)
415   
416
417
418
419
420
421 
422
423
424def tsh2sww(filename, verbose=True):
425    """
426    to check if a tsh/msh file 'looks' good.
427    """
428
429    #FIXME Move to data_manager
430    from shallow_water import Domain
431    from pmesh2domain import pmesh_to_domain_instance
432    import time, os
433    from data_manager import get_dataobject
434    from os import sep, path
435    #from util import mean
436
437    if verbose == True:print 'Creating domain from', filename
438    domain = pmesh_to_domain_instance(filename, Domain)
439    if verbose == True:print "Number of triangles = ", len(domain)
440
441    domain.smooth = True
442    domain.format = 'sww'   #Native netcdf visualisation format
443    file_path, filename = path.split(filename)
444    filename, ext = path.splitext(filename)
445    domain.filename = filename
446    domain.reduction = mean
447    if verbose == True:print "file_path",file_path
448    if file_path == "":file_path = "."
449    domain.set_datadir(file_path)
450
451    if verbose == True:
452        print "Output written to " + domain.get_datadir() + sep + \
453              domain.filename + "." + domain.format
454    sww = get_dataobject(domain)
455    sww.store_connectivity()
456    sww.store_timestep('stage')
457
458
459
460
461def multiple_replace(text, dictionary):
462    """Multiple replace of words in text
463
464    text:       String to be modified
465    dictionary: Mapping of words that are to be substituted
466
467    Python Cookbook 3.14 page 88 and page 90
468    """
469
470    import re
471   
472    #Create a regular expression from all of the dictionary keys
473    #matching only entire words
474    regex = re.compile(r'\b'+ \
475                       r'\b|\b'.join(map(re.escape, dictionary.keys()))+ \
476                       r'\b' )
477
478    #For each match, lookup the corresponding value in the dictionary
479    return regex.sub(lambda match: dictionary[match.group(0)], text)
480
481
482
483
484def apply_expression_to_dictionary(expression, dictionary):#dictionary):
485    """Apply arbitrary expression to values of dictionary
486
487    Given an expression in terms of the keys, replace key by the
488    corresponding values and evaluate.
489   
490
491    expression: Arbitrary, e.g. arithmetric, expression relating keys
492                from dictionary.
493
494    dictionary: Mapping between symbols in expression and objects that
495                will be evaluated by expression.
496                Values in dictionary must support operators given in
497                expression e.g. by overloading               
498    """
499
500    import types
501    import re
502
503    assert isinstance(expression, basestring)
504    assert type(dictionary) == types.DictType
505
506    #Convert dictionary values to textual representations suitable for eval
507    D = {}
508    for key in dictionary:
509        D[key] = 'dictionary["%s"]' %key
510
511    #Perform substitution of variables   
512    expression = multiple_replace(expression, D)
513
514    #Evaluate and return
515    try:
516        return eval(expression)
517    except NameError, e:
518        msg = 'Expression "%s" could not be evaluated: %s' %(expression, e)
519        raise NameError, msg
520   
521
522
523
524
525
526
527
528####################################################################
529#Python versions of function that are also implemented in util_gateway.c
530#
531
532def gradient_python(x0, y0, x1, y1, x2, y2, q0, q1, q2):
533    """
534    """
535
536    det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)
537    a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
538    a /= det
539
540    b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
541    b /= det
542
543    return a, b
544
545
546def gradient2_python(x0, y0, x1, y1, q0, q1):
547    """Compute radient based on two points and enforce zero gradient
548    in the direction orthogonal to (x1-x0), (y1-y0)
549    """
550
551    #Old code
552    #det = x0*y1 - x1*y0
553    #if det != 0.0:
554    #    a = (y1*q0 - y0*q1)/det
555    #    b = (x0*q1 - x1*q0)/det
556
557    #Correct code (ON)
558    det = (x1-x0)**2 + (y1-y0)**2
559    if det != 0.0:
560        a = (x1-x0)*(q1-q0)/det
561        b = (y1-y0)*(q1-q0)/det
562       
563    return a, b       
564
565
566##############################################
567#Initialise module
568
569from utilities import compile
570if compile.can_use_C_extension('util_ext.c'):
571    from util_ext import gradient, gradient2#, point_on_line
572    #separate_points_by_polygon = separate_points_by_polygon_c
573else:
574    gradient = gradient_python
575    gradient2 = gradient2_python   
576
577
578if __name__ == "__main__":
579    pass
580
581
582
583
Note: See TracBrowser for help on using the repository browser.