Changeset 820 for inundation/ga

Timestamp:

Feb 1, 2005, 5:06:54 AM (20 years ago)

Author:

ole

Message:

Implemented xya2rectangular and tested it

Location:

inundation/ga/storm_surge/pyvolution

Files:

• least_squares.py (modified) (2 diffs)
• mesh_factory.py (modified) (2 diffs)
• test_least_squares.py (modified) (1 diff)
• test_util.py (modified) (1 diff)
• util.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -143 +143 @@
     return vertex_attributes
 
+    
+    
+def xya2rectangular(xya_name, M, N, alpha = DEFAULT_ALPHA):
+    """Fits attributes from xya file to MxN rectangular mesh 
+    
+    Read xya file and create rectangular mesh of resolution MxN such that
+    it covers all points specified in xya file. 
+    
+    FIXME: This may be a temporary function until we decide on 
+    netcdf formats etc
+    """   
+    
+    import util, mesh_factory
+    
+    points, attributes = util.read_xya(xya_name) 
+    
+    #Find extent
+    max_x = min_x = points[0][0]
+    max_y = min_y = points[0][1]
+    for point in points[1:]:
+        x = point[0] 
+        if x > max_x: max_x = x
+        if x < min_x: min_x = x           
+        y = point[1] 
+        if y > max_y: max_y = y
+        if y < min_y: min_y = y           
+    
+    #Create appropriate mesh
+    vertex_coordinates, triangles, boundary =\
+         mesh_factory.rectangular(M, N, max_x-min_x, max_y-min_y, 
+                                (min_x, min_y))
+
+    #Fit attributes to mesh      
+    vertex_attributes = fit_to_mesh(vertex_coordinates,
+                        triangles,
+                        points,
+                        attributes, alpha=alpha)
+
+
+         
+    return vertex_coordinates, triangles, boundary, vertex_attributes
+                         
+    
 
 class Interpolation:
@@ -548 +591 @@
         if n<m and self.alpha == 0.0:
             msg = 'ERROR (least_squares): Too few data points\n'
-            msg += 'There only %d data points. Need at least %d\n' %(n,m)
-            msg += 'Alternatively, increase smoothing parameter alpha' 
+            msg += 'There are only %d data points and alpha == 0. ' %n
+            msg += 'Need at least %d\n' %m
+            msg += 'Alternatively, set smoothing parameter alpha to a small '
+            msg += 'positive value,\ne.g. 1.0e-3.'
             raise msg
 

inundation/ga/storm_surge/pyvolution/mesh_factory.py

@@ -131 +131 @@
     return points, elements, boundary       
 
-#OLD FORMULA FOR SETTING BOUNDARY TAGS - OBSOLETE NOW             
-#     for id, face in M.boundary:
-
-#         e = element_class.instances[id]
-#         x0, y0, x1, y1, x2, y2 = e.get_instance_vertex_coordinates()
-
-#         if face==2:
-#             #Left or right#
-#             if abs(x0-origin[0]) < epsilon:
-#                 M.boundary[(id,face)] = 'left'            
-#             elif abs(origin[0]+lenx-x0) < epsilon:
-#                 M.boundary[(id,face)] = 'right'
-#             else:
-#                 print face, id, id%m, m, n
-#                 raise 'Left or Right Unknown boundary'                            
-#         elif face==1:
-#             #Top or bottom
-#             if x0 > 0.25*lenx and abs(y0-a1-a2*x0-origin[1]) < epsilon or\
-#                x0 <= 0.25*lenx and abs(y0-origin[1]) < epsilon:            
-#                    M.boundary[(id,face)] = 'bottom'
-#             elif abs(origin[1]+leny-y0) < epsilon:
-#                 M.boundary[(id,face)] = 'top'
-#             else:
-#                 print face, id, id%m, m, n
-#                 raise 'Top or Bottom Unknown boundary'  
-#         else:
-#             print face, id, id%m, m, n
-#             raise 'Unknown boundary'          
-#   
-#     return M
-
-
-
 
 def circular(m, n, radius=1.0, center = (0.0, 0.0)):
@@ -394 +361 @@
     return points, elements                             
 
-
-
-# def contracting_channel_mesh(m, n, x1 = 0.0, x2 = 1./3., x3 = 2./3., x4 = 1.0,
-#                             y1 =0.0, y4 = -1./4., y8 = 1.0,
-#                              origin = (0.0, 0.0), point_class=Point,
-#                              element_class=Triangle, mesh_class=Mesh):
-#     """Setup a oblique grid of triangles
-#     with m segments in the x-direction and n segments in the y-direction
-
-#     Triangle refers to the actual class or subclass to be instantiated:
-#     e.g. if Volume is a subclass of Triangle,
-#     this function can be invoked with the keywords
-#       oblique_mesh(...,Triangle=Volume, Mesh=Domain)
-#     """
-
-#     from Numeric import array
-#     from visual import rate#
-
-#     import math
-
-#     from config import epsilon
-
-#     E = m*n*2        #Number of triangular elements
-#     P = (m+1)*(n+1)  #Number of initial vertices
-
-#     initialise_consecutive_datastructure(P+E, E,
-#                                          point_class,
-#                                          element_class,
-#                                          mesh_class)
-
-#     deltax = (x4 - x1)/float(m)
-#     deltay = (y8 - y1)/float(n)
-#     a = y4 - y1
-    
-#     if y8 - a <= y1 + a:
-#         print a,y1,y4
-#         raise 'Constriction is too large reduce a'
-            
-#     y2 = y1
-#     y3 = y4
-#     x5 = x4
-#     y5 = y8 - a
-#     x6 = x3
-#     y6 = y5
-#     x7 = x2
-#     y7 = y8
-#     x8 = x1
-    
-#     a2 = a/(x3 - x2)
-#     a1 = a - a*x3/(x3 - x2)
-#     a4 = (-a + a2*(x7 - x6))/(x6 - x7)/y7
-#     a3 = (y7 - a1 - x7*a2 - a4*x7*y7)/y7
-    
-#     # Dictionary of vertex objects
-#     vertices = {}
-
-#     for i in range(m+1):
-#         x = deltax*i
-#         for j in range(n+1):     
-#             y = deltay*j
-#             if x > x2:
-#                 if x < x3:
-#                     y = a1 + a2*x + a3*y + a4*x*y
-#                 else:
-#                     y = a + y*(y5 - y4)/(y8 - y1)   
-                
-#             vertices[i,j] = Point(x + origin[0],y + origin[1])
-
-#     # Construct 2 elements per element        
-#     elements = []
-#     for i in range(m):
-#         for j in range(n):
-#             v1 = vertices[i,j+1]
-#             v2 = vertices[i,j]            
-#             v3 = vertices[i+1,j+1]            
-#             v4 = vertices[i+1,j]
- 
-#             elements.append(element_class(v4,v3,v2)) #Lower                
-#             elements.append(element_class(v1,v2,v3)) #Upper    
-
-#     M = mesh_class(elements)
-
-#     #Set a default tagging
-
-#     for id, face in M.boundary:
-
-#         e = element_class.instances[id]
-#         x_0, y_0, x_1, y_1, x_2, y_2 = e.get_instance_vertex_coordinates()
-#         lenx = x4 - x1
-#         if face==2:
-#             #Left or right#
-#             if abs(x_0-origin[0]) < epsilon:
-#                 M.boundary[(id,face)] = 'left'            
-#             elif abs(origin[0]+lenx-x_0) < epsilon:
-#                 M.boundary[(id,face)] = 'right'
-#             else:
-#                 print face, id, id%m, m, n
-#                 raise 'Left or Right Unknown boundary'                            
-#         elif face==1:
-#             #Top or bottom
-#             if x_0 <= x2 and abs(y_0-y1-origin[1]) < epsilon or\
-#                x_0 > x3 and abs(y_0-y3-origin[1]) < epsilon or\
-#                x_0 >  x2 and x_0 <= x3 and abs(y_0-(y2+(y3-y2)*(x_0-x2)/(x3-x2)+origin[1])) < epsilon:            
-#                 M.boundary[(id,face)] = 'bottom'
-#             elif x_0 <= x7 and abs(y_0-y8-origin[1]) < epsilon or\
-#                x_0 > x6 and abs(y_0-y6-origin[1]) < epsilon or\
-#                x_0 > x7 and x_0 <= x6 and abs(y_0-(y7+(y6-y7)*(x_0-x7)/(x6-x7)+origin[1])) < epsilon:            
-#                 M.boundary[(id,face)] = 'top'
-#             else:
-#                 print face, id, id%m, m, n 
-#                 raise 'Top or Bottom Unknown boundary'  
-#         else:
-#             print face, id, id%m, m, n
-#             raise 'Unknown boundary'          
-
-            
-#       #  print id, face, M.boundary[(id,face)],x_0,y_0,x_1,y_1,x_2,y_2
-    
-#     return M
-
-
-
 # #Map from edge number to indices of associated vertices
 # edge_map = ((1,2), (0,2), (0,1))

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -488 +488 @@
 
 
+        
+    def test_xya2rectangular(self):
+
+        import time, os
+        FN = 'xyatest' + str(time.time()) + '.xya'
+        fid = open(FN, 'w')
+        fid.write('%f %f %f %f\n' %(1,1,2,4) )
+        fid.write('%f %f %f %f\n' %(1,3,4,8) )
+        fid.write('%f %f %f %f\n' %(3,1,4,8) )          
+        fid.close()
+        
+        points, triangles, boundary, attributes = xya2rectangular(FN, 4, 4)
+        
+        
+        
+
+        z1 = [2, 4]
+        z2 = [4, 8]
+        z3 = [4, 8]
+        data_coords = [ [1,1], [1,3], [3,1] ]
+        z = [z1, z2, z3]
+
+        ref = fit_to_mesh(points, triangles, data_coords, z)            
+        
+        assert allclose(attributes, ref)
+
+        
+        
 
     #Tests of smoothing matrix    

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -336 +336 @@
         assert allclose( 2*sin(120*pi/600)/3 + sin(60*pi/600)/3, q[2] )
         
+
+    def test_xya(self):
+        import time, os
+        FN = 'xyatest' + str(time.time()) + '.xya'
+        fid = open(FN, 'w')
+        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
+        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
+        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )            
+        fid.close()
+        
+        points, attributes = read_xya(FN)
+        
+        assert allclose(points, [ [0,1], [1,0], [1,1] ])
+        assert allclose(attributes, [ [10,20,30], [30,20,10], 
+                [40.2,40.3,40.4] ])
+        
+        os.remove(FN)
         
         

inundation/ga/storm_surge/pyvolution/util.py

@@ -80 +80 @@
       return False  
     
+             
+class File_function:
+    """Read time series from file and return a callable object:
+    f(t,x,y)
+    which will return interpolated values based on the input file.
+
+    The file format is assumed to be either two fields separated by a comma:
+
+        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
+
+    or four comma separated fields
+
+        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
+
+    In either case, the callable object will return a tuple of
+    interpolated values, one each value stated in the file.
+
+    
+    E.g
+
+      31/08/04 00:00:00, 1.328223 0 0
+      31/08/04 00:15:00, 1.292912 0 0
+
+    will provide a time dependent function f(t,x=None,y=None),
+    ignoring x and y, which returns three values per call.
+
+    
+    NOTE: At this stage the function is assumed to depend on
+    time only, i.e  no spatial dependency!!!!!
+    When that is need we can use the least_squares interpolation.
+    """
+
+    
+    def __init__(self, filename, domain=None):
+        """Initialise callable object from a file.
+        See docstring for class File_function
+
+        If domain is specified, model time (domain,starttime)
+        will be checked and possibly modified
+
+        ALl times are assumed to be in UTC
+        """
+
+        import time, calendar
+        from Numeric import array
+        from config import time_format
+
+        assert type(filename) == type(''),\
+               'First argument to File_function must be a string'
+
+
+        try:
+            fid = open(filename)
+        except Exception, e:
+            msg = 'File "%s" could not be opened: Error="%s"'\
+                  %(filename, e)
+            raise msg
+
+
+        line = fid.readline()
+        fid.close()
+        fields = line.split(',')
+        msg = 'File %s must have the format date, value0 value1 value2 ...'
+        msg += ' or date, x, y, value0 value1 value2 ...'
+        mode = len(fields)
+        assert mode in [2,4], msg
+
+        try:
+            starttime = calendar.timegm(time.strptime(fields[0], time_format))
+        except ValueError:    
+            msg = 'First field in file %s must be' %filename
+            msg += ' date-time with format %s.\n' %time_format
+            msg += 'I got %s instead.' %fields[0] 
+            raise msg
+
+
+        #Split values
+        values = []
+        for value in fields[mode-1].split():
+            values.append(float(value))
+
+        q = array(values)
+        
+        msg = 'ERROR: File must contain at least one independent value'
+        assert len(q.shape) == 1, msg
+            
+        self.number_of_values = len(q)
+
+        self.filename = filename
+        self.starttime = starttime
+        self.domain = domain
+
+        if domain is not None:
+            if domain.starttime is None:
+                domain.starttime = self.starttime
+            else:
+                msg = 'WARNING: Start time as specified in domain (%s)'\
+                      %domain.starttime
+                msg += ' is earlier than the starttime of file %s: %s.'\
+                         %(self.filename, self.starttime)
+                msg += 'Modifying starttime accordingly.'
+                if self.starttime > domain.starttime:
+                    #FIXME: Print depending on some verbosity setting
+                    #print msg
+                    domain.starttime = self.starttime #Modifying model time
+
+        if mode == 2:        
+            self.read_times() #Now read all times in.
+        else:
+            raise 'x,y dependecy not yet implemented'
+
+        
+    def read_times(self):
+        """Read time and values 
+        """
+        from Numeric import zeros, Float, alltrue
+        from config import time_format
+        import time, calendar
+        
+        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)
+
+        msg = 'File %s must list time as a monotonuosly ' %self.filename
+        msg += 'increasing sequence'
+        assert alltrue( T[1:] - T[:-1] > 0 ), msg
+
+        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):
+        """Evaluate f(t,x,y)
+
+        FIXME: x, y dependency not yet implemented except that
+        result is a vector of same length as x and y
+        """
+
+        from math import pi, cos, sin, sqrt
+        
+
+        #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
+
+        #Find 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 = self.Q[self.index,:] +\
+            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
+
+        #Return vector of interpolated values
+        if x == None and y == None:
+            return q
+        else:
+            from Numeric import ones, Float
+            #Create one constant column for each value
+            N = len(x)
+            assert len(y) == N, 'x and y must have same length'
+
+            res = []
+            for col in q:
+                res.append(col*ones(N, Float))
+            
+            return res
+
+def read_xya(file_name):
+    """Read simple xya file, no header, just
+    x y [attributes]
+    separated by whitespace
+    
+    Return list of points and list of attributes
+    """
+    
+    fid = open(file_name)
+    lines = fid.readlines()
+    
+    points = [] 
+    attributes = []
+    for line in lines:
+        fields = line.strip().split()
+        
+        points.append( (float(fields[0]),  float(fields[1])) )
+        attributes.append( [float(z) for z in fields[2:] ] )
+        
+    return points, attributes
+    
+
+#####################################
+#POLYFON STUFF
+#
+#FIXME: ALl these should be put into new module polygon.py
 
 def inside_polygon(point, polygon, closed = True):
@@ -286 +518 @@
 
     return polygon      
-             
-             
-class File_function:
-    """Read time series from file and return a callable object:
-    f(t,x,y)
-    which will return interpolated values based on the input file.
-
-    The file format is assumed to be either two fields separated by a comma:
-
-        time [DD/MM/YY hh:mm:ss], value0 value1 value2 ...
-
-    or four comma separated fields
-
-        time [DD/MM/YY hh:mm:ss], x, y, value0 value1 value2 ...
-
-    In either case, the callable object will return a tuple of
-    interpolated values, one each value stated in the file.
-
-    
-    E.g
-
-      31/08/04 00:00:00, 1.328223 0 0
-      31/08/04 00:15:00, 1.292912 0 0
-
-    will provide a time dependent function f(t,x=None,y=None),
-    ignoring x and y, which returns three values per call.
-
-    
-    NOTE: At this stage the function is assumed to depend on
-    time only, i.e  no spatial dependency!!!!!
-    When that is need we can use the least_squares interpolation.
-    """
-
-    
-    def __init__(self, filename, domain=None):
-        """Initialise callable object from a file.
-        See docstring for class File_function
-
-        If domain is specified, model time (domain,starttime)
-        will be checked and possibly modified
-
-        ALl times are assumed to be in UTC
-        """
-
-        import time, calendar
-        from Numeric import array
-        from config import time_format
-
-        assert type(filename) == type(''),\
-               'First argument to File_function must be a string'
-
-
-        try:
-            fid = open(filename)
-        except Exception, e:
-            msg = 'File "%s" could not be opened: Error="%s"'\
-                  %(filename, e)
-            raise msg
-
-
-        line = fid.readline()
-        fid.close()
-        fields = line.split(',')
-        msg = 'File %s must have the format date, value0 value1 value2 ...'
-        msg += ' or date, x, y, value0 value1 value2 ...'
-        mode = len(fields)
-        assert mode in [2,4], msg
-
-        try:
-            starttime = calendar.timegm(time.strptime(fields[0], time_format))
-        except ValueError:    
-            msg = 'First field in file %s must be' %filename
-            msg += ' date-time with format %s.\n' %time_format
-            msg += 'I got %s instead.' %fields[0] 
-            raise msg
-
-
-        #Split values
-        values = []
-        for value in fields[mode-1].split():
-            values.append(float(value))
-
-        q = array(values)
-        
-        msg = 'ERROR: File must contain at least one independent value'
-        assert len(q.shape) == 1, msg
-            
-        self.number_of_values = len(q)
-
-        self.filename = filename
-        self.starttime = starttime
-        self.domain = domain
-
-        if domain is not None:
-            if domain.starttime is None:
-                domain.starttime = self.starttime
-            else:
-                msg = 'WARNING: Start time as specified in domain (%s)'\
-                      %domain.starttime
-                msg += ' is earlier than the starttime of file %s: %s.'\
-                         %(self.filename, self.starttime)
-                msg += 'Modifying starttime accordingly.'
-                if self.starttime > domain.starttime:
-                    #FIXME: Print depending on some verbosity setting
-                    #print msg
-                    domain.starttime = self.starttime #Modifying model time
-
-        if mode == 2:        
-            self.read_times() #Now read all times in.
-        else:
-            raise 'x,y dependecy not yet implemented'
-
-        
-    def read_times(self):
-        """Read time and values 
-        """
-        from Numeric import zeros, Float, alltrue
-        from config import time_format
-        import time, calendar
-        
-        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)
-
-        msg = 'File %s must list time as a monotonuosly ' %self.filename
-        msg += 'increasing sequence'
-        assert alltrue( T[1:] - T[:-1] > 0 ), msg
-
-        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):
-        """Evaluate f(t,x,y)
-
-        FIXME: x, y dependency not yet implemented except that
-        result is a vector of same length as x and y
-        """
-
-        from math import pi, cos, sin, sqrt
-        
-
-        #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
-
-        #Find 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 = self.Q[self.index,:] +\
-            ratio*(self.Q[self.index+1,:] - self.Q[self.index,:]) 
-
-        #Return vector of interpolated values
-        if x == None and y == None:
-            return q
-        else:
-            from Numeric import ones, Float
-            #Create one constant column for each value
-            N = len(x)
-            assert len(y) == N, 'x and y must have same length'
-
-            res = []
-            for col in q:
-                res.append(col*ones(N, Float))
-            
-            return res
-    
-
+    
 def populate_polygon(polygon, number_of_points, seed = None):
     """Populate given polygon with uniformly distributed points.