Changeset 851 for inundation/ga/storm_surge/pyvolution/util.py

Timestamp:

Feb 8, 2005, 7:07:56 PM (20 years ago)

Author:

ole

Message:

First cut at a spatio-temporal boundary.
Need to improve the test, though.

File:

• inundation/ga/storm_surge/pyvolution/util.py (modified) (7 diffs)

inundation/ga/storm_surge/pyvolution/util.py

@@ -123 +123 @@
     quantities = ['stage', 'xmomentum', 'ymomentum']
     
-    interpolation_points decides at which points interpolated quantities are to be computed
+    interpolation_points decides at which points interpolated 
+    quantities are to be computed whenever object is called.
+    
+    
+    
     If None, return average value
     """
@@ -133 +137 @@
         will be checked and possibly modified
         
-
-
         All times are assumed to be in UTC
         """
@@ -213 +215 @@
         msg = 'File %s must list time as a monotonuosly ' %self.filename
         msg += 'increasing sequence'
-        assert alltrue( time[1:] - time[:-1] > 0 ), msg 
+        assert alltrue(time[1:] - time[:-1] > 0 ), msg  
         
 
@@ -227 +229 @@
             
             ####self.Q = zeros( (len(time),
-             
-            for i, t in enumerate(time[:]):
+            self.T = time[:] 
+            self.index = 0 #Initial index       
+
+            self.values = {}                
+            for name in self.quantities:
+                self.values[name] = zeros( (len(self.T), 
+                                            len(interpolation_points)), 
+                                            Float)                  
+            
+            for i, t in enumerate(self.T):
     
                 #Interpolate quantities at this timestep 
@@ -238 +248 @@
                                          verbose = False)
 
-                for name in self.quantities:                            
-                    print t, name, interpol.interpolate(fid.variables[name][i,:])                    
+                for name in self.quantities:
+                    self.values[name][i, :] =\
+                    interpol.interpolate(fid.variables[name][i,:])   
 
         else:
             pass 
 
-        self.T = time[:]
-        
-        return
-        
-        fid = open(self.filename)        
-        lines = fid.readlines()
-        fid.close()
-        
-        N = len(lines)
-        d = self.number_of_values
-
-        T = zeros(N, Float)       #Time
-        Q = zeros((N, d), Float)  #Values
-        
-        for i, line in enumerate(lines):
-            fields = line.split(',')
-            realtime = calendar.timegm(time.strptime(fields[0], time_format))
-
-            T[i] = realtime - self.starttime
-
-            for j, value in enumerate(fields[1].split()):
-                Q[i, j] = float(value)
-
-
-
-        self.T = T     #Time
-        self.Q = Q     #Values
-        self.index = 0 #Initial index
-
-        
+    def __repr__(self):
+        return 'File function'
+
+    def __call__(self, t, x=None, y=None, point_id = None):
+        """Evaluate f(t, point_id)
+        
+        Inputs:
+          t: time - Model time (tau = domain.starttime-self.starttime+t)  
+                    must lie within existing timesteps   
+          point_id: index of one of the preprocessed points. 
+                    If point_id is None all preprocessed points are computed
+
+          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? 
+        """
+
+        from math import pi, cos, sin, sqrt
+        from Numeric import zeros, Float
+
+        #Find time tau such that
+        #
+        # domain.starttime + t == self.starttime + tau
+        
+        if self.domain is not None:
+            tau = self.domain.starttime-self.starttime+t
+        else:
+            tau = t
+                
+        msg = 'Time interval derived from file %s (%s:%s) does not match model time: %s'\
+              %(self.filename, self.T[0], self.T[1], tau)
+        if tau < self.T[0]: raise msg
+        if tau > self.T[-1]: raise msg        
+
+        oldindex = self.index #Time index
+
+        #Find time slot
+        while tau > self.T[self.index]: self.index += 1
+        while tau < self.T[self.index]: self.index -= 1
+
+        #t is now between index and index+1
+        ratio = (tau - self.T[self.index])/\
+                (self.T[self.index+1] - self.T[self.index])
+
+                
+        #Compute interpolated values
+        q = zeros( len(self.quantities), Float)
+        
+        for i, name in enumerate(self.quantities):
+            Q = self.values[name]    
+
+            q[i] = Q[self.index, point_id] +\
+                 ratio*(Q[self.index+1, point_id] - Q[self.index, point_id]) 
+
+            
+        #Return vector of interpolated values
+        return q
+                
              
 class File_function_ASCII:
@@ -430 +470 @@
         return 'File function'
 
-        
-
-    def __call__(self, t, x=None, y=None):
+    def __call__(self, t, x=None, y=None, point_id=None):
         """Evaluate f(t,x,y)
 
@@ -438 +476 @@
         result is a vector of same length as x and y
         FIXME: Naaaa
+        
+        FIXME: Rethink semantics when x,y are vectors.
         """