source: inundation/pyvolution/util.py @ 1919

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

Implemented arbitrary expressions for sww2dem using code from changeset:1916

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