Changeset 4295

Timestamp:

Mar 6, 2007, 10:56:39 AM (18 years ago)

Author:

duncan

Message:

Adding mint, maxt to new URS2sww

Location:

anuga_core/source/anuga/shallow_water

Files:

• data_manager.py (modified) (11 diffs)
• test_data_manager.py (modified) (1 diff)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -4587 +4587 @@
     
 def urs_ungridded2sww(basename_in='o', basename_out=None, verbose=False,
-            #mint=None, maxt=None,
+            mint=None, maxt=None,
             mean_stage=0,
             origin = None,
@@ -4594 +4594 @@
     """
     parameters not used!
-    origin
     #mint=None, maxt=None,
     """ 
@@ -4635 +4634 @@
 
     #mesh.export_mesh_file(basename_in + '.tsh')
-    times = []
+    # These are the times of the mux file
+    mux_times = []
     for i in range(a_mux.time_step_count):
-        times.append(a_mux.time_step * i)
+        mux_times.append(a_mux.time_step * i)  
+    mux_times_start_i, mux_times_fin_i = mux2sww_time(mux_times, mint, maxt)
+    times = mux_times[mux_times_start_i:mux_times_fin_i]
+    
+    if mux_times_start_i == mux_times_fin_i:
+        # Close the mux files
+        for quantity, file in map(None, quantities, files_in):
+            mux[quantity].close()
+        msg="Due to mint and maxt there's no time info in the boundary SWW."
+        assert 1==0 , msg #Happy to know a better way of doing this..
         
+    # If this raise is removed there is currently no downstream errors
+           
     points_utm=ensure_numeric(points_utm)
     #print "mesh_dic['generatedpointlist']", mesh_dic['generatedpointlist']
@@ -4654 +4665 @@
 
     outfile = NetCDFFile(swwname, 'w')
-    # For a different way of doing this, check out tsh2sww 
+    # For a different way of doing this, check out tsh2sww
+    # work out sww_times and the index range this covers
     write_sww_header(outfile, times, len(volumes), len(points_utm))
+    outfile.mean_stage = mean_stage
+    outfile.zscale = zscale
     write_sww_triangulation(outfile, points_utm, volumes,
                             elevation, zone,  origin=origin, verbose=verbose)
@@ -4662 +4676 @@
     j = 0
     # Read in a time slice from each mux file and write it to the sww file
-    for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):            
-        write_sww_time_slice(outfile, HA, UA, VA,  elevation, j,
-                              mean_stage=mean_stage, zscale=zscale,
-                              verbose=verbose)
+    for HA, UA, VA in map(None, mux['HA'], mux['UA'], mux['VA']):
+        if j >= mux_times_start_i and j < mux_times_fin_i:
+            write_sww_time_slice(outfile, HA, UA, VA,  elevation,
+                                 j - mux_times_start_i,
+                                 mean_stage=mean_stage, zscale=zscale,
+                                 verbose=verbose)
         j += 1
     outfile.close()
     #
     # Do some conversions while writing the sww file
+def mux2sww_time(mux_times, mint, maxt):
+    """
+    """
+
+    if mint == None:
+        mux_times_start_i = 0
+    else:
+        mux_times_start_i = searchsorted(mux_times, mint)
+       
+    if maxt == None:
+        mux_times_fin_i = len(mux_times)
+    else:
+        maxt += 0.5 # so if you specify a time where there is
+                    # data that time will be included
+        mux_times_fin_i = searchsorted(mux_times, maxt)
+
+    return mux_times_start_i, mux_times_fin_i
 
 def write_sww_header(outfile, times, number_of_volumes, number_of_points ):
@@ -4679 +4712 @@
     number_of_times = len(times)
     
-    #Create new file
     outfile.institution = 'Geoscience Australia'
     outfile.description = 'Converted from XXX'
@@ -4689 +4721 @@
 
     #Start time in seconds since the epoch (midnight 1/1/1970)
-    outfile.starttime = starttime = times[0]
+    if len(times) == 0:
+        outfile.starttime = starttime = 0
+    else:
+        outfile.starttime = starttime = times[0]
 
 
@@ -4775 +4810 @@
 
 
-    # FIXME: Don't store the whole speeds and stage in  memory.
-    # block it?
 def write_sww_time_slices(outfile, has, uas, vas, elevation,
                          mean_stage=0, zscale=1,
@@ -4796 +4829 @@
         j += 1
 
- 
 def write_sww_time_slice(outfile, ha, ua, va, elevation, slice_index,
                          mean_stage=0, zscale=1,
@@ -4804 +4836 @@
     xmomentum = outfile.variables['xmomentum']
     ymomentum = outfile.variables['ymomentum']
-
-    
+   
     w = zscale*ha + mean_stage
     stage[slice_index] = w
@@ -4811 +4842 @@
     xmomentum[slice_index] = ua*h
     ymomentum[slice_index] = va*h
+
         
 class Urs_points:

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -5436 +5436 @@
         os.remove(sww_file)
         
+    def test_urs_ungridded2sww_mint_maxt (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
+                          mint=101, maxt=500)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_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]
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(lat_long[0][0], lat_long[0][1]) 
+        assert allclose([x[0],y[0]], [e,n])
+
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = [200,300,400,500]
+       
+        assert allclose(ensure_numeric(times),
+                        ensure_numeric(times_actual))
+        
+               #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+        assert allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        assert allclose(answer, actual)  #Meters
+
+        # check the stage values, first time step.
+        # These arrays are equal since the Easting values were used as
+        # the stage
+        assert allclose(stage[0], x +tide)  #Meters
+        # check the elevation values.
+        # -ve since urs measures depth, sww meshers height,
+        # these arrays are equal since the northing values were used as
+        # the elevation
+        assert allclose(-elevation, y)  #Meters
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_urs_ungridded2sww_mint_maxtII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        urs_ungridded2sww(base_name, mean_stage=tide, origin =(50,23432,4343),
+                          mint=0, maxt=100000)
+        
+        # now I want to check the sww file ...
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file
+        # Note, the sww info is not gridded.  It is point data.
+        fid = NetCDFFile(sww_file)
+        
+        # Make x and y absolute
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        points = geo_reference.get_absolute(map(None, fid.variables['x'][:],
+                                                fid.variables['y'][:]))
+        points = ensure_numeric(points)
+        x = points[:,0]
+        
+        #Check the time vector
+        times = fid.variables['time'][:]
+        
+        times_actual = [0,100,200,300,400,500]
+        assert allclose(ensure_numeric(times),
+                        ensure_numeric(times_actual))
+        
+        #Check first value
+        stage = fid.variables['stage'][:]
+        assert allclose(stage[0], x +tide)
+        
+        fid.close()
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def test_urs_ungridded2sww_mint_maxtIII (self):
+        
+        #Zone:   50    
+        #Easting:  240992.578  Northing: 7620442.472 
+        #Latitude:   -21  30 ' 0.00000 ''  Longitude: 114  30 ' 0.00000 '' 
+        lat_long = [[-21.5,114.5],[-21,114.5],[-21,115]]
+        time_step_count = 6
+        time_step = 100
+        tide = 9000000
+        base_name, files = self.write_mux(lat_long,
+                                          time_step_count, time_step)
+        try:
+            urs_ungridded2sww(base_name, mean_stage=tide,
+                          origin =(50,23432,4343),
+                          mint=301, maxt=399)
+        except: 
+            pass
+        else:
+            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
+
+        self.delete_mux(files)
+        
     def test_URS_points_needed_and_urs_ungridded2sww(self):
         # This doesn't actually check anything