Changeset 1668

Timestamp:

Aug 1, 2005, 5:30:53 PM (20 years ago)

Author:

ole

Message:

Separated Interpolation_function out from File_function and refactored the latter in
terms of the former.
Started unifying the general version with temporal interpolation only.

Location:

inundation/ga/storm_surge/pyvolution

Files:

• test_util.py (modified) (4 diffs)
• util.py (modified) (7 diffs)

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -559 +559 @@
                           interpolation_points = interpolation_points)
         assert allclose(domain.starttime, start+1)
-
         domain.starttime = start
 
 
-
         #Check linear interpolation in time
+        F = file_function(filename + '.sww', domain,
+                          quantities = domain.conserved_quantities,
+                          interpolation_points = interpolation_points)                
         #for id in range(len(interpolation_points)):
         for id in [1]:
@@ -580 +581 @@
 
                 q = F(t, point_id=id)
+                #print i, k, t, q
+                #print ' ', q0
+                #print ' ', q1
+                #print 's', (k*q1 + (6-k)*q0)/6
+                #print
                 assert allclose(q, (k*q1 + (6-k)*q0)/6)
 
@@ -589 +595 @@
 
             t = 90 #Halfway between 60 and 120
-            q = F(t,point_id=id)
+            q = F(t, point_id=id)
             assert allclose( (q120+q60)/2, q )
 
@@ -1083 +1089 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    #suite = unittest.makeSuite(TestCase,'test_inside_polygon_main')
+    #suite = unittest.makeSuite(Test_Util,'test_spatio_temporal_file_function_time')
     suite = unittest.makeSuite(Test_Util,'test')
     runner = unittest.TextTestRunner()

inundation/ga/storm_surge/pyvolution/util.py

@@ -5 +5 @@
 """
 
+#FIXME: Move polygon stuff out
+#FIXME: Finish Interpolation function and get rid of File_function_ASCII
+
+
 def angle(v):
     """Compute angle between e1 (the unit vector in the x-direction)
@@ -17 +21 @@
     v2 = v[1]/l
 
+
     theta = acos(v1)
-    try:
-        theta = acos(v1)
-    except ValueError, e:
-        print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
-
-        #FIXME, hack to avoid acos(1.0) Value error
-        # why is it happening?
-        # is it catching something we should avoid?
-        #FIXME (Ole): When did this happen? We need a unit test to
-        #reveal this condition or otherwise remove the hack
-        
-        s = 1e-6
-        if (v1+s >  1.0) and (v1-s < 1.0) :
-            theta = 0.0
-        elif (v1+s >  -1.0) and (v1-s < -1.0):
-            theta = 3.1415926535897931
-        print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
-              %(v1, theta)
+    
+    #try:
+    #   theta = acos(v1)
+    #except ValueError, e:
+    #    print 'WARNING (util.py): Angle acos(%s) failed: %s' %(str(v1), e)
+    #
+    #    #FIXME, hack to avoid acos(1.0) Value error
+    #    # why is it happening?
+    #    # is it catching something we should avoid?
+    #    #FIXME (Ole): When did this happen? We need a unit test to
+    #    #reveal this condition or otherwise remove the hack
+    #    
+    #    s = 1e-6
+    #    if (v1+s >  1.0) and (v1-s < 1.0) :
+    #        theta = 0.0
+    #    elif (v1+s >  -1.0) and (v1-s < -1.0):
+    #        theta = 3.1415926535897931
+    #    print 'WARNING (util.py): angle v1 is %f, setting acos to %f '\
+    #          %(v1, theta)
 
     if v2 < 0:
@@ -156 +162 @@
         quantities = None
 
-    #Choose format
-    #FIXME: Maybe these can be merged later on
-    #FIXME: We shoulp probably not waste energy on ASCII files -
-    #the different formats are a pain. E.g. the time field.
-    #
+
     if line[:3] == 'CDF':
-        return File_function_NetCDF(filename, domain, quantities,
-                                    interpolation_points, verbose = verbose)
+        return get_netcdf_file_function(filename, domain, quantities,
+                                        interpolation_points, verbose = verbose)
     else:
+        #FIXME Should be phased out
         return File_function_ASCII(filename, domain, quantities,
                                    interpolation_points)
-
-
+    
+    
+
+
+
+
+
+
+
+
+def get_netcdf_file_function(filename, domain=None, quantity_names=None,
+                             interpolation_points=None, verbose = False):
+    """Read time history of spatial data from NetCDF sww file and
+    return a callable object f(t,x,y)
+    which will return interpolated values based on the input file.
+
+    If domain is specified, model time (domain.starttime)
+    will be checked and possibly modified
+    
+    All times are assumed to be in UTC
+
+    See Interpolation function for further documetation
+    
+
+    """
+    
+    #verbose = True
+    
+    #FIXME: Check that model origin is the same as file's origin
+    #(both in UTM coordinates)
+    #If not - modify those from file to match domain
+    #Take this code from e.g. dem2pts in data_manager.py
+    #FIXME: Use geo_reference to read and write xllcorner...
+        
+
+    import time, calendar
+    from config import time_format
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import array, zeros, Float, alltrue, concatenate, reshape        
+
+    #Open NetCDF file
+    if verbose: print 'Reading', filename
+    fid = NetCDFFile(filename, 'r')
+
+    if quantity_names is None or len(quantity_names) < 1:
+        msg = 'ERROR: At least one independent value must be specified'
+        raise msg
+
+    missing = []
+    for quantity in ['x', 'y', 'time'] + quantity_names:
+        if not fid.variables.has_key(quantity):
+            missing.append(quantity)
+
+    if len(missing) > 0:
+        msg = 'Quantities %s could not be found in file %s'\
+              %(str(missing), filename)
+        raise msg
+
+
+    #Get first timestep
+    try:
+        starttime = fid.starttime[0]
+    except ValueError:
+        msg = 'Could not read starttime from file %s' %filename
+        raise msg
+
+    #Get origin
+    xllcorner = fid.xllcorner[0]
+    yllcorner = fid.yllcorner[0]
+
+
+    #Get variables
+    if verbose: print 'Get variables'
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    triangles = fid.variables['volumes'][:]
+    time = fid.variables['time'][:]
+
+    x = reshape(x, (len(x),1))
+    y = reshape(y, (len(y),1))
+    vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+
+    if domain is not None:
+        #Update domain.startime if it is *earlier* than starttime
+        if starttime > domain.starttime:
+            msg = 'WARNING: Start time as specified in domain (%f)'\
+                  %domain.starttime
+            msg += ' is earlier than the starttime of file %s (%f).'\
+                     %(filename, starttime)
+            msg += ' Modifying domain starttime accordingly.'
+            
+            if verbose: print msg
+            domain.starttime = starttime #Modifying model time
+            if verbose: print 'Domain starttime is now set to %f'\
+               %domain.starttime
+
+
+        #If domain.startime is *later* than starttime,
+        #move time back - relative to domain's time
+        if domain.starttime > starttime:
+            time = time - domain.starttime + starttime
+
+
+    if verbose:
+        print 'File_function data obtained from: %s' %filename
+        print '  References:'
+        #print '    Datum ....' #FIXME
+        print '    Lower left corner: [%f, %f]'\
+              %(xllcorner, yllcorner)
+        print '    Start time:   %f' %starttime                
+        
+    
+    #Produce values for desired data points at
+    #each timestep
+
+    quantities = {}
+    for i, name in enumerate(quantity_names):
+        quantities[name] = fid.variables[name][:]
+    fid.close()
+    
+
+    return Interpolation_function(vertex_coordinates,
+                                  triangles,
+                                  time,
+                                  quantities,
+                                  quantity_names,  
+                                  interpolation_points,
+                                  alpha = 0,
+                                  verbose = verbose)
+
+
+
+
+
+class Interpolation_function:
+    """Interpolation_function - creates callable object f(t, id) or f(t,x,y)
+    which is interpolated from time series defined at vertices of
+    triangular mesh 
+
+    Input
+        vertex_coordinates: mx2 array of coordinates (Float)
+        triangles:          nx3 array of indices into vertex_coordinates (Int)
+        time:               px1 array of times (Float)               
+        quantities:         Dictionary of pxm values
+        quantity_names:
+        interpolation_points:
+        alpha:
+        verbose:
+    
+
+    
+    The quantities returned by the callable object are specified by
+    the list quantities which must contain the names of the
+    quantities to be returned and also reflect the order, e.g. for
+    the shallow water wave equation, on would have
+    quantities = ['stage', 'xmomentum', 'ymomentum']
+
+    The parameter interpolation_points decides at which points interpolated
+    quantities are to be computed whenever object is called.
+    If None, return average value
+    """
+
+    
+
+    
+    def __init__(self,
+                 vertex_coordinates,
+                 triangles,
+                 time,
+                 quantities,
+                 quantity_names = None,  
+                 interpolation_points = None,
+                 alpha = 0,
+                 verbose = False):
+        """
+        """
+
+        from Numeric import array, zeros, Float, alltrue, concatenate,\
+             reshape, ArrayType
+
+        from config import time_format
+        from least_squares import Interpolation
+
+        #Check
+        msg = 'Time must be a monotonuosly '
+        msg += 'increasing sequence %s' %time
+        assert alltrue(time[1:] - time[:-1] > 0 ), msg
+
+
+        #Check if quantities is a single array only
+        if type(quantities) == ArrayType:
+            quantity_names = 'Attribute'
+            quantities = {quantity_names: quantities}            
+
+        #Use keys if no names specified
+        if quantity_names is not None:
+            self.quantity_names = quantity_names
+        else:
+            self.quantity_names = quantities.keys()
+
+
+        self.interpolation_points = interpolation_points
+        self.T = time[:]  #Time assumed to be relative to starttime
+        self.index = 0    #Initial time index
+        self.values = {}
+            
+
+        if interpolation_points is not None:
+            #Spatial interpolation
+
+            try:
+                interpolation_points = ensure_numeric(interpolation_points)
+            except:
+                msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
+                msg += 'I got: %s.' %( str(interpolation_points)[:60] + '...')
+                raise msg
+
+
+            for name in quantity_names:
+                self.values[name] = zeros( (len(self.T),
+                                            len(interpolation_points)),
+                                            Float)
+
+            #Build interpolator
+            interpol = Interpolation(vertex_coordinates,
+                                     triangles,
+                                     point_coordinates = interpolation_points,
+                                     alpha = 0,
+                                     verbose = verbose)
+
+
+            if verbose: print 'Interpolate'
+            for i, t in enumerate(self.T):
+                #Interpolate quantities at this timestep
+                if verbose: print ' time step %d of %d' %(i, len(self.T))
+                for name in quantity_names:
+                    self.values[name][i, :] =\
+                    interpol.interpolate(quantities[name][i,:])
+
+            #Report
+            if verbose:
+                x = vertex_coordinates[:,0]
+                y = vertex_coordinates[:,1]                
+                
+                print '------------------------------------------------'
+                print 'Interpolation_function statistics:'
+                print '  Extent:'
+                print '    x in [%f, %f], len(x) == %d'\
+                      %(min(x.flat), max(x.flat), len(x.flat))
+                print '    y in [%f, %f], len(y) == %d'\
+                      %(min(y.flat), max(y.flat), len(y.flat))
+                print '    t in [%f, %f], len(t) == %d'\
+                      %(min(self.T), max(self.T), len(self.T))
+                print '  Quantities:'
+                for name in quantity_names:
+                    q = quantities[name][:].flat
+                    print '    %s in [%f, %f]' %(name, min(q), max(q))
+
+
+                print '  Interpolation points (xi, eta):'\
+                      'number of points == %d ' %interpolation_points.shape[0]
+                print '    xi in [%f, %f]' %(min(interpolation_points[:,0]),
+                                             max(interpolation_points[:,0]))
+                print '    eta in [%f, %f]' %(min(interpolation_points[:,1]),
+                                              max(interpolation_points[:,1]))
+                print '  Interpolated quantities (over all timesteps):'
+                for name in quantity_names:
+                    q = self.values[name][:].flat
+                    print '    %s at interpolation points in [%f, %f]'\
+                          %(name, min(q), max(q))
+                print '------------------------------------------------'
+        else:
+            #Return an average, making this a time series
+            #FIXME: Needs a test
+            for name in quantity_names:
+                self.values[name] = avg(quantities[name][i,:])
+
+
+    def __repr__(self):
+        return 'Interpolation function'
+
+    def __call__(self, t, point_id = None, x = None, y = None):
+        """Evaluate f(t), f(t, point_id) or f(t, x, y)
+
+        Inputs:
+          t: time - Model time. Must lie within existing timesteps
+          point_id: index of one of the preprocessed points.
+                    If point_id is None all preprocessed points are computed
+                    
+          If no x,y info is present, point_id and x,y arguments are ignored
+          making f a function of time only.
+          
+          FIXME: point_id could also be a slice
+          FIXME: One could allow arbitrary x, y coordinates
+          (requires computation of a new interpolator)
+          FIXME: What if x and y are vectors?
+          FIXME: Does point_id or x,y take precedence?
+          FIXME: What about f(x,y) without t
+          
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import zeros, Float
+
+
+        msg = 'Time interval [%s:%s]' %(self.T[0], self.T[1])
+        msg += ' does not match model time: %s\n' %t
+        if t < self.T[0]: raise msg
+        if t > self.T[-1]: raise msg
+
+        oldindex = self.index #Time index
+
+        #Find current time slot
+        while t > self.T[self.index]: self.index += 1
+        while t < self.T[self.index]: self.index -= 1
+
+        if t == self.T[self.index]:
+            #Protect against case where t == T[-1] (last time)
+            # - also works in general when t == T[i]
+            ratio = 0
+        else:
+            #t is now between index and index+1
+            ratio = (t - self.T[self.index])/\
+                    (self.T[self.index+1] - self.T[self.index])
+
+        #Compute interpolated values
+        q = zeros(len(self.quantity_names), Float)
+
+        for i, name in enumerate(self.quantity_names):
+            Q = self.values[name]
+
+            if point_id is None:
+                Q0 = Q[self.index]
+
+                try:
+                    Q1 = Q[self.index+1]
+                except:
+                    pass                
+            else:    
+                Q0 = Q[self.index, point_id]
+
+                try:
+                    Q1 = Q[self.index+1, point_id]
+                except:
+                    pass
+                
+
+            #Linear temporal interpolation    
+            if ratio > 0:
+                q[i] = Q0 + ratio*(Q1 - Q0)
+            else:
+                q[i] = Q0
+            
+
+            #if ratio > 0:
+            #    q[i] = Q[self.index, point_id] +\
+            #           ratio*(Q[self.index+1, point_id] -\
+            #                  Q[self.index, point_id])
+            #else:
+            #    q[i] = Q[self.index, point_id]
+
+        #Return vector of interpolated values
+        return q
+
+
+
+
+
+
+
+
+
+
+#FIXME Obsolete
 class File_function_NetCDF:
     """Read time history of spatial data from NetCDF sww file and
@@ -203 +578 @@
         #If not - modify those from file to match domain
         #Take this code from e.g. dem2pts in data_manager.py
+        #FIXME: Use geo_reference to read and write xllcorner...
         
 
@@ -281 +657 @@
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
-        z = fid.variables['elevation'][:]
+        #z = fid.variables['elevation'][:]
         triangles = fid.variables['volumes'][:]
         time = fid.variables['time'][:]
@@ -410 +786 @@
             tau = self.domain.starttime-self.starttime+t
         else:
-            #print 'DOMAIN IS NONE!!!!!!!!!'
             tau = t
 
-        #print 'D start', self.domain.starttime
-        #print 'S start', self.starttime
-        #print 't', t
-        #print 'tau', tau
 
         msg = 'Time interval derived from file %s [%s:%s]'\
@@ -464 +835 @@
         return q
 
-
+#FIXME Obsolete
 class File_function_ASCII:
     """Read time series from file and return a callable object: