Changeset 1884

Timestamp:

Oct 10, 2005, 11:41:26 AM (19 years ago)

Author:

ole

Message:

Made interpolation_points in file_function absolute UTM coordinates and wrote test

Location:

inundation/pyvolution

Files:

• least_squares.py (modified) (1 diff)
• test_util.py (modified) (1 diff)
• util.py (modified) (5 diffs)

inundation/pyvolution/least_squares.py

@@ -939 +939 @@
         from Numeric import array, zeros, Float, alltrue, concatenate,\
              reshape, ArrayType
+
 
         from util import mean, ensure_numeric

inundation/pyvolution/test_util.py

@@ -472 +472 @@
         os.remove(filename)
 
-        
+
+
+    def test_spatio_temporal_file_function_different_origin(self):
+        """Test that spatio temporal file function performs the correct
+        interpolations in both time and space where space is offset by
+        xllcorner and yllcorner
+        NetCDF version (x,y,t dependency)        
+        """
+        import time
+
+        #Create sww file of simple propagation from left to right
+        #through rectangular domain
+        from shallow_water import Domain, Dirichlet_boundary
+        from mesh_factory import rectangular
+        from Numeric import take, concatenate, reshape
+
+
+        from coordinate_transforms.geo_reference import Geo_reference
+        xllcorner = 2048
+        yllcorner = 11000
+        zone = 2
+
+        #Create basic mesh and shallow water domain
+        points, vertices, boundary = rectangular(3, 3)
+        domain1 = Domain(points, vertices, boundary,
+                         geo_reference = Geo_reference(xllcorner = xllcorner,
+                                                       yllcorner = yllcorner))
+        
+
+        from util import mean
+        domain1.reduction = mean
+        domain1.smooth = True #NOTE: Mimic sww output where each vertex has
+                              # only one value.
+
+        domain1.default_order = 2
+        domain1.store = True
+        domain1.set_datadir('.')
+        domain1.set_name('spatio_temporal_boundary_source_%d' %(id(self)))
+        domain1.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+        #Bed-slope, friction and IC at vertices (and interpolated elsewhere)
+        domain1.set_quantity('elevation', 0)
+        domain1.set_quantity('friction', 0)
+        domain1.set_quantity('stage', 0)
+
+        # Boundary conditions
+        B0 = Dirichlet_boundary([0,0,0])
+        B6 = Dirichlet_boundary([0.6,0,0])
+        domain1.set_boundary({'left': B6, 'top': B6, 'right': B0, 'bottom': B0})
+        domain1.check_integrity()
+
+        finaltime = 8
+        #Evolution
+        for t in domain1.evolve(yieldstep = 0.1, finaltime = finaltime):
+            pass
+            #domain1.write_time()
+
+
+        #Now read data from sww and check
+        from Scientific.IO.NetCDF import NetCDFFile
+        filename = domain1.get_name() + '.' + domain1.format
+        fid = NetCDFFile(filename)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        time = fid.variables['time'][:]
+
+        #Take stage vertex values at last timestep on diagonal
+        #Diagonal is identified by vertices: 0, 5, 10, 15
+
+        timestep = len(time)-1 #Last timestep
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0 = (D[0] + D[1])/2
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
+
+        #And the midpoints are found now
+        Dx = take(reshape(x, (16,1)), [0,5,10,15])
+        Dy = take(reshape(y, (16,1)), [0,5,10,15])
+
+        diag = concatenate( (Dx, Dy), axis=1)
+        d_midpoints = (diag[1:] + diag[:-1])/2
+
+
+        #Adjust for georef - make interpolation points absolute
+        d_midpoints[:,0] += xllcorner
+        d_midpoints[:,1] += yllcorner                
+
+        #Let us see if the file function can find the correct
+        #values at the midpoints at the last timestep:
+        f = file_function(filename, domain1,
+                          interpolation_points = d_midpoints)
+
+        q = f(timestep/10., point_id=0); assert allclose(r0, q)
+        q = f(timestep/10., point_id=1); assert allclose(r1, q)
+        q = f(timestep/10., point_id=2); assert allclose(r2, q)
+
+
+        ##################
+        #Now do the same for the first timestep
+
+        timestep = 0 #First timestep
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0 = (D[0] + D[1])/2
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
+
+        #Let us see if the file function can find the correct
+        #values
+        q = f(0, point_id=0); assert allclose(r0, q)
+        q = f(0, point_id=1); assert allclose(r1, q)
+        q = f(0, point_id=2); assert allclose(r2, q)
+
+
+        ##################
+        #Now do it again for a timestep in the middle
+
+        timestep = 33
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0 = (D[0] + D[1])/2
+        r1 = (D[1] + D[2])/2
+        r2 = (D[2] + D[3])/2
+
+        q = f(timestep/10., point_id=0); assert allclose(r0, q)
+        q = f(timestep/10., point_id=1); assert allclose(r1, q)
+        q = f(timestep/10., point_id=2); assert allclose(r2, q)
+
+
+        ##################
+        #Now check temporal interpolation
+        #Halfway between timestep 15 and 16
+
+        timestep = 15
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0_0 = (D[0] + D[1])/2
+        r1_0 = (D[1] + D[2])/2
+        r2_0 = (D[2] + D[3])/2
+
+        #
+        timestep = 16
+        d_stage = reshape(take(stage[timestep, :], [0,5,10,15]), (4,1))
+        d_uh = reshape(take(xmomentum[timestep, :], [0,5,10,15]), (4,1))
+        d_vh = reshape(take(ymomentum[timestep, :], [0,5,10,15]), (4,1))
+        D = concatenate( (d_stage, d_uh, d_vh), axis=1)
+
+        #Reference interpolated values at midpoints on diagonal at
+        #this timestep are
+        r0_1 = (D[0] + D[1])/2
+        r1_1 = (D[1] + D[2])/2
+        r2_1 = (D[2] + D[3])/2
+
+        # The reference values are
+        r0 = (r0_0 + r0_1)/2
+        r1 = (r1_0 + r1_1)/2
+        r2 = (r2_0 + r2_1)/2
+
+        q = f((timestep - 0.5)/10., point_id=0); assert allclose(r0, q)
+        q = f((timestep - 0.5)/10., point_id=1); assert allclose(r1, q)
+        q = f((timestep - 0.5)/10., point_id=2); assert allclose(r2, q)
+
+        ##################
+        #Finally check interpolation 2 thirds of the way
+        #between timestep 15 and 16
+
+        # The reference values are
+        r0 = (r0_0 + 2*r0_1)/3
+        r1 = (r1_0 + 2*r1_1)/3
+        r2 = (r2_0 + 2*r2_1)/3
+
+        #And the file function gives
+        q = f((timestep - 1.0/3)/10., point_id=0); assert allclose(r0, q)
+        q = f((timestep - 1.0/3)/10., point_id=1); assert allclose(r1, q)
+        q = f((timestep - 1.0/3)/10., point_id=2); assert allclose(r2, q)
+
+        fid.close()
+        import os
+        os.remove(filename)
+
+        
+
 
     def test_spatio_temporal_file_function_time(self):

inundation/pyvolution/util.py

@@ -155 +155 @@
     All times are assumed to be in UTC
 
-    All spatial information is assumed to be in UTM coordinates.
+    All spatial information is assumed to be in absolute UTM coordinates.
 
     See Interpolation function for further documentation
@@ -170 +170 @@
     #- domain's georef
     #- sww file's georef
-    #- interpolation points georef
+    #- interpolation points as absolute UTM coordinates
 
 
@@ -228 +228 @@
     from config import time_format
     from Scientific.IO.NetCDF import NetCDFFile
-    from Numeric import array, zeros, Float, alltrue, concatenate, reshape        
+    from Numeric import array, zeros, Float, alltrue, concatenate, reshape
+    from util import ensure_numeric    
 
     #Open NetCDF file
@@ -251 +252 @@
 
 
+    if interpolation_points is not None:
+        interpolation_points = ensure_numeric(interpolation_points)
+
+
+
     #Now assert that requested quantitites (and the independent ones)
     #are present in file 
@@ -302 +308 @@
         y = reshape(y, (len(y),1))
         vertex_coordinates = concatenate((x,y), axis=1) #m x 2 array
+
+        if interpolation_points is not None:
+            #Adjust for georef
+            interpolation_points[:,0] -= xllcorner
+            interpolation_points[:,1] -= yllcorner        
+        
+