Changeset 5462 for anuga_core/source/anuga/shallow_water/test_data_manager.py

Timestamp:

Jul 4, 2008, 11:48:37 AM (15 years ago)

Author:

ole

Message:

John and Ole implemented urs2sts using an ordering file creating smaller and pre-ordered sts files for use with boundary conditions.

File:

• anuga_core/source/anuga/shallow_water/test_data_manager.py (modified) (4 diffs)

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -6319 +6319 @@
                                       va=va)
 
-        urs2sts(base_name,mean_stage=tide,verbose=False)
+        urs2sts(base_name,
+                basename_out=base_name, 
+                mean_stage=tide,verbose=False)
 
         # now I want to check the sts file ...
@@ -6441 +6443 @@
 
         # Call urs2sts with multiple mux files
-        urs2sts([base_nameI, base_nameII], weights=[1.0,1.0],
+        urs2sts([base_nameI, base_nameII], 
+                basename_out=base_nameI, 
+                weights=[1.0,1.0],
                 mean_stage=tide,
                 verbose=False)
@@ -6534 +6538 @@
         self.delete_mux(filesII)
         os.remove(sts_file)
+        
+        
+    def test_urs2sts_ordering(self):
+        """Test multiple sources with ordering file
+        """
+        
+        tide=0
+        time_step_count = 3
+        time_step = 2
+        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
+        n=len(lat_long_points)
+        first_tstep=ones(n,Int)
+        first_tstep[0]+=1
+        first_tstep[2]+=1
+        last_tstep=(time_step_count)*ones(n,Int)
+        last_tstep[0]-=1
+
+        gauge_depth=20*ones(n,Float)
+        ha=2*ones((n,time_step_count),Float)
+        ha[0]=arange(0,time_step_count)
+        ha[1]=arange(time_step_count,2*time_step_count)
+        ha[2]=arange(2*time_step_count,3*time_step_count)
+        ha[3]=arange(3*time_step_count,4*time_step_count)
+        ua=5*ones((n,time_step_count),Float)
+        va=-10*ones((n,time_step_count),Float)
+
+        # Create two identical mux files to be combined by urs2sts
+        base_nameI, filesI = self.write_mux2(lat_long_points,
+                                             time_step_count, time_step,
+                                             first_tstep, last_tstep,
+                                             depth=gauge_depth,
+                                             ha=ha,
+                                             ua=ua,
+                                             va=va)
+
+        base_nameII, filesII = self.write_mux2(lat_long_points,
+                                               time_step_count, time_step,
+                                               first_tstep, last_tstep,
+                                               depth=gauge_depth,
+                                               ha=ha,
+                                               ua=ua,
+                                               va=va)
+
+                                               
+        # Create ordering file
+        permutation = [3,0,2]
+
+        # Do it wrongly at first and check that exception is being raised
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][0], 
+                                             lat_long_points[index][1]))
+        order_fid.close()
+        
+        try:
+            urs2sts([base_nameI, base_nameII], 
+                    basename_out=base_nameI, 
+                    ordering_filename=ordering_filename,
+                    weights=[1.0,1.0],
+                    mean_stage=tide,
+                    verbose=False)  
+            os.remove(ordering_filename)            
+        except:
+            pass
+        else:
+            msg = 'Should have caught wrong lat longs'
+            raise Exception, msg
+
+
+        
+            
+        # Then do it right
+        _, ordering_filename = tempfile.mkstemp('')
+        order_fid = open(ordering_filename, 'w')  
+        order_fid.write('index, longitude, latitude\n')
+        for index in permutation:
+            order_fid.write('%d, %f, %f\n' %(index, 
+                                             lat_long_points[index][1], 
+                                             lat_long_points[index][0]))
+        order_fid.close()
+        
+            
+
+                                               
+        # Call urs2sts with multiple mux files
+        urs2sts([base_nameI, base_nameII], 
+                basename_out=base_nameI, 
+                ordering_filename=ordering_filename,
+                weights=[1.0,1.0],
+                mean_stage=tide,
+                verbose=False)
+
+        # now I want to check the sts file ...
+        sts_file = base_nameI + '.sts'
+
+        #Let's interrogate the sts file
+        # Note, the sts info is not gridded.  It is point data.
+        fid = NetCDFFile(sts_file)
+
+        # Make x and y absolute
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, x, y))
+        points = ensure_numeric(points)
+
+        x = points[:,0]
+        y = points[:,1]
+
+        for i, index in enumerate(permutation):
+            # Check that STS points are stored in the correct order
+            
+            # Work out the UTM coordinates sts point i
+            zone, e, n = redfearn(lat_long_points[index][0], 
+                                  lat_long_points[index][1])             
+            
+            assert allclose([x[i],y[i]], [e,n])
+            
+                        
+        # Check the time vector
+        times = fid.variables['time'][:]
+
+        times_actual = []
+        for i in range(time_step_count):
+            times_actual.append(time_step * i)
+
+        assert allclose(ensure_numeric(times),
+                        ensure_numeric(times_actual))
+                        
+
+        # Check sts values
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        # Set original data used to write mux file to be zero when gauges are
+        # not recdoring
+        
+        ha[0][0]=0.0
+        ha[0][time_step_count-1]=0.0
+        ha[2][0]=0.0
+        ua[0][0]=0.0
+        ua[0][time_step_count-1]=0.0
+        ua[2][0]=0.0
+        va[0][0]=0.0
+        va[0][time_step_count-1]=0.0
+        va[2][0]=0.0;
+
+        
+        # The stage stored in the .sts file should be the sum of the stage
+        # in the two mux2 files because both have weights = 1. In this case
+        # the mux2 files are the same so stage == 2.0 * ha
+        #print 2.0*transpose(ha) - stage 
+        
+        ha_permutation = take(ha, permutation) 
+        ua_permutation = take(ua, permutation)         
+        va_permutation = take(va, permutation)                 
+        gauge_depth_permutation = take(gauge_depth, permutation)                         
+
+        
+        assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_ua *(stage+depth)
+
+        depth=zeros((len(lat_long_points),time_step_count),Float)
+        for i in range(len(lat_long_points)):
+            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
+            #2.0*ha necessary because using two files with weights=1 are used
+            
+        depth_permutation = take(depth, permutation)                     
+        
+
+        # The xmomentum stored in the .sts file should be the sum of the ua
+        # in the two mux2 files multiplied by the depth.
+        assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum) 
+
+        #Check the momentums - va
+        #momentum = velocity*(stage-elevation)
+        # elevation = - depth
+        #momentum = velocity_va *(stage+depth)
+
+        # The ymomentum stored in the .sts file should be the sum of the va
+        # in the two mux2 files multiplied by the depth.
+        assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
+
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        assert allclose(-gauge_depth_permutation, elevation)  #Meters
+
+        fid.close()
+        self.delete_mux(filesI)
+        self.delete_mux(filesII)
+        os.remove(sts_file)
+        
 
     def test_file_boundary_stsI(self):
@@ -6563 +6769 @@
 
         sts_file=base_name
-        urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
+        urs2sts(base_name,
+                sts_file,
+                mean_stage=tide,verbose=False)
         self.delete_mux(files)