Changeset 7866 for trunk/anuga_core/source/anuga/file_conversion

Timestamp:

Jun 22, 2010, 12:04:32 PM (15 years ago)

Author:

James Hudson

Message:

More swb tests passing. Cleaned up some pylint errors.

Location:

trunk/anuga_core/source/anuga/file_conversion

Files:

• file_conversion.py (modified) (5 diffs)
• test_sww2dem.py (modified) (1 diff)
• test_urs2sts.py (modified) (1 diff)

trunk/anuga_core/source/anuga/file_conversion/file_conversion.py

@@ -91 +91 @@
 # @param quantity_names 
 # @param time_as_seconds 
-def timefile2netcdf(file_text, quantity_names=None, time_as_seconds=False):
+def timefile2netcdf(file_text, file_out = None, quantity_names=None, \
+                                time_as_seconds=False):
     """Template for converting typical text files with time series to
     NetCDF tms file.
@@ -111 +112 @@
     will provide a time dependent function f(t) with three attributes
 
-    filename is assumed to be the rootname with extenisons .txt and .sww
+    filename is assumed to be the rootname with extensions .txt/.tms and .sww
     """
 
@@ -121 +122 @@
         raise IOError('Input file %s should be of type .txt.' % file_text)
 
-    filename = file_text[:-4]
+    if file_out == None:
+        file_out = file_text[:-4] + '.tms'
+
     fid = open(file_text)
     line = fid.readline()
@@ -181 +184 @@
             Q[i, j] = float(value)
 
-    msg = 'File %s must list time as a monotonuosly ' % filename
+    msg = 'File %s must list time as a monotonuosly ' % file_text
     msg += 'increasing sequence'
     assert num.alltrue(T[1:] - T[:-1] > 0), msg
@@ -188 +191 @@
     from Scientific.IO.NetCDF import NetCDFFile
 
-    fid = NetCDFFile(filename + '.tms', netcdf_mode_w)
+    fid = NetCDFFile(file_out, netcdf_mode_w)
 
     fid.institution = 'Geoscience Australia'

trunk/anuga_core/source/anuga/file_conversion/test_sww2dem.py

@@ -1520 +1520 @@
         os.remove(ascfile)
         os.remove(swwfile2)
+
+
+    def test_export_grid(self):
+        """
+        test_export_grid(self):
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        try:
+            os.remove('teg*.sww')
+        except:
+            pass
+
+        #Setup
+        self.domain.set_name('teg')
+
+        prjfile = self.domain.get_name() + '_elevation.prj'
+        ascfile = self.domain.get_name() + '_elevation.asc'
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        fid.close()
+
+        #Export to ascii/prj files
+        sww2dem_batch(self.domain.get_name(),
+                quantities = 'elevation',
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def test_export_gridII(self):
+        """
+        test_export_gridII(self):
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        try:
+            os.remove('teg*.sww')
+        except:
+            pass
+
+        #Setup
+        self.domain.set_name('tegII')
+
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep()
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep()
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]        
+
+        #print 'stage', stage
+        #print 'xmom', xmomentum
+        #print 'ymom', ymomentum        
+
+        fid.close()
+
+        #Export to ascii/prj files
+        if True:
+            sww2dem_batch(self.domain.get_name(),
+                        quantities = ['elevation', 'depth'],
+                        cellsize = cellsize,
+                        verbose = self.verbose,
+                        format = 'asc')
+
+        else:
+            sww2dem_batch(self.domain.get_name(),
+                quantities = ['depth'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+
+            export_grid(self.domain.get_name(),
+                quantities = ['elevation'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+        prjfile = self.domain.get_name() + '_elevation.prj'
+        ascfile = self.domain.get_name() + '_elevation.asc'
+        
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #print "ascfile", ascfile
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        
+        #Check asc file
+        ascfile = self.domain.get_name() + '_depth.asc'
+        prjfile = self.domain.get_name() + '_depth.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])                
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+
+    def test_export_gridIII(self):
+        """
+        test_export_gridIII
+        Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        try:
+            os.remove('teg*.sww')
+        except:
+            pass
+
+        #Setup
+        
+        self.domain.set_name('tegIII')
+
+        swwfile = self.domain.get_name() + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        self.domain.set_quantity('stage', 1.0)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+        
+        sww = SWW_file(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep() #'stage')
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep() #'stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, netcdf_mode_r)
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+        fid.close()
+
+        #Export to ascii/prj files
+        extra_name_out = 'yeah'
+        if True:
+            sww2dem_batch(self.domain.get_name(),
+                        quantities = ['elevation', 'depth'],
+                        extra_name_out = extra_name_out,
+                        cellsize = cellsize,
+                        verbose = self.verbose,
+                        format = 'asc')
+
+        else:
+            sww2dem_batch(self.domain.get_name(),
+                quantities = ['depth'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+
+            sww2dem_batch(self.domain.get_name(),
+                quantities = ['elevation'],
+                cellsize = cellsize,
+                verbose = self.verbose,
+                format = 'asc')
+
+        prjfile = self.domain.get_name() + '_elevation_yeah.prj'
+        ascfile = self.domain.get_name() + '_elevation_yeah.asc'
+        
+        #Check asc file
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #print "ascfile", ascfile
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                #print " -i*cellsize - y",  -i*cellsize - y
+                #print "float(L[i])", float(L[i])
+                assert num.allclose(float(L[i]), -i*cellsize - y)
+                
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        
+        #Check asc file
+        ascfile = self.domain.get_name() + '_depth_yeah.asc'
+        prjfile = self.domain.get_name() + '_depth_yeah.prj'
+        ascid = open(ascfile)
+        lines = ascid.readlines()
+        ascid.close()
+
+        L = lines[2].strip().split()
+        assert L[0].strip().lower() == 'xllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 308500)
+
+        L = lines[3].strip().split()
+        assert L[0].strip().lower() == 'yllcorner'
+        assert num.allclose(float(L[1].strip().lower()), 6189000)
+
+        #Check grid values
+        for j in range(5):
+            L = lines[6+j].strip().split()
+            y = (4-j) * cellsize
+            for i in range(5):
+                assert num.allclose(float(L[i]), 1 - (-i*cellsize - y))
+
+        #Cleanup
+        os.remove(prjfile)
+        os.remove(ascfile)
+        os.remove(swwfile)
+
+    def test_export_grid_bad(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        try:
+            sww2dem_batch('a_small_round-egg',
+                        quantities = ['elevation', 'depth'],
+                        cellsize = 99,
+                        verbose = self.verbose,
+                        format = 'asc')
+        except IOError:
+            pass
+        else:
+            self.failUnless(0 ==1,  'Bad input did not throw exception error!')
         
 

trunk/anuga_core/source/anuga/file_conversion/test_urs2sts.py

@@ -1958 +1958 @@
         os.remove(meshname)
         
-
+        
+        
+    def test_file_boundary_sts_time_limit(self):
+        """test_file_boundary_stsIV_sinewave_ordering(self):
+        Read correct points from ordering file and apply sts to boundary
+        This one uses a sine wave and compares to time boundary
+        
+        This one tests that times used can be limited by upper limit
+        """
+        
+        lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]]
+        bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]]
+        tide = 0.35
+        time_step_count = 50
+        time_step = 0.1
+        times_ref = num.arange(0, time_step_count*time_step, time_step)
+        
+        n=len(lat_long_points)
+        first_tstep=num.ones(n,num.int)
+        last_tstep=(time_step_count)*num.ones(n,num.int)
+        
+        gauge_depth=20*num.ones(n,num.float)
+        
+        ha1=num.ones((n,time_step_count),num.float)
+        ua1=3.*num.ones((n,time_step_count),num.float)
+        va1=2.*num.ones((n,time_step_count),num.float)
+        for i in range(n):
+            ha1[i]=num.sin(times_ref)
+        
+        
+        base_name, files = self.write_mux2(lat_long_points,
+                                           time_step_count, time_step,
+                                           first_tstep, last_tstep,
+                                           depth=gauge_depth,
+                                           ha=ha1,
+                                           ua=ua1,
+                                           va=va1)
+
+        # Write order file
+        file_handle, order_base_name = tempfile.mkstemp("")
+        os.close(file_handle)
+        os.remove(order_base_name)
+        d=","
+        order_file=order_base_name+'order.txt'
+        fid=open(order_file,'w')
+        
+        # Write Header
+        header='index, longitude, latitude\n'
+        fid.write(header)
+        indices=[3,0,1]
+        for i in indices:
+            line=str(i)+d+str(lat_long_points[i][1])+d+\
+                str(lat_long_points[i][0])+"\n"
+            fid.write(line)
+        fid.close()
+
+        sts_file=base_name
+        urs2sts(base_name, basename_out=sts_file,
+                ordering_filename=order_file,
+                mean_stage=tide,
+                verbose=False)
+        self.delete_mux(files)
+        
+        
+        
+        # Now read the sts file and check that values have been stored correctly.
+        fid = NetCDFFile(sts_file + '.sts')
+
+        # Check the time vector
+        times = fid.variables['time'][:]
+        starttime = fid.starttime[0]
+        #print times
+        #print starttime
+
+        # Check sts quantities
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+        elevation = fid.variables['elevation'][:]
+
+        
+
+        # Create beginnings of boundary polygon based on sts_boundary
+        boundary_polygon = create_sts_boundary(base_name)
+        
+        os.remove(order_file)
+
+        # Append the remaining part of the boundary polygon to be defined by
+        # the user
+        bounding_polygon_utm=[]
+        for point in bounding_polygon:
+            zone,easting,northing=redfearn(point[0],point[1])
+            bounding_polygon_utm.append([easting,northing])
+
+        boundary_polygon.append(bounding_polygon_utm[3])
+        boundary_polygon.append(bounding_polygon_utm[4])
+
+        #print 'boundary_polygon', boundary_polygon
+        
+
+        assert num.allclose(bounding_polygon_utm,boundary_polygon)
+
+
+        extent_res=1000000
+        meshname = 'urs_test_mesh' + '.tsh'
+        interior_regions=None
+        boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]}
+        
+        # have to change boundary tags from last example because now bounding
+        # polygon starts in different place.
+        create_mesh_from_regions(boundary_polygon,
+                                 boundary_tags=boundary_tags,
+                                 maximum_triangle_area=extent_res,
+                                 filename=meshname,
+                                 interior_regions=interior_regions,
+                                 verbose=False)
+        
+        domain_fbound = Domain(meshname)
+        domain_fbound.set_quantity('stage', tide)
+        
+        
+        Bf = File_boundary(sts_file+'.sts', 
+                           domain_fbound, 
+                           boundary_polygon=boundary_polygon)
+        time_vec = Bf.F.get_time()
+        assert num.allclose(Bf.F.starttime, starttime)
+        assert num.allclose(time_vec, times_ref)                                   
+        
+        for time_limit in [0.1, 0.2, 0.5, 1.0, 2.2, 3.0, 4.3, 6.0, 10.0]:
+            Bf = File_boundary(sts_file+'.sts', 
+                               domain_fbound, 
+                               time_limit=time_limit+starttime,
+                               boundary_polygon=boundary_polygon)
+        
+            time_vec = Bf.F.get_time()
+            assert num.allclose(Bf.F.starttime, starttime)            
+            assert num.alltrue(time_vec < time_limit)
+            
+            
+        try:    
+            Bf = File_boundary(sts_file+'.sts', 
+                               domain_fbound, 
+                               time_limit=-1+starttime,
+                               boundary_polygon=boundary_polygon)            
+            time_vec = Bf.F.get_time()    
+            print time_vec    
+        except AssertionError:
+            pass
+        else:
+            raise Exception, 'Should have raised Exception here'
 
 #-------------------------------------------------------------