Changeset 1360

Timestamp:

May 9, 2005, 10:03:30 PM (20 years ago)

Author:

steve

Message:

 

Location:

inundation/ga/storm_surge

Files:

• pyvolution/data_manager.py (modified) (52 diffs)
• pyvolution/domain.py (modified) (1 diff)
• pyvolution/mesh.py (modified) (1 diff)
• pyvolution/netherlands.py (modified) (5 diffs)
• pyvolution/shallow_water.py (modified) (2 diffs)
• pyvolution/test_data_manager.py (modified) (46 diffs)
• zeus/pyvolution.zbi (modified) (previous)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -381 +381 @@
             #other modules (I think).
 
-            #write a filename addon that won't break swollens reader 
+            #write a filename addon that won't break swollens reader
             #(10.sww is bad)
-            filename_ext = '_time_%s'%self.domain.time            
+            filename_ext = '_time_%s'%self.domain.time
             filename_ext = filename_ext.replace('.', '_')
-            #remember the old filename, then give domain a 
+            #remember the old filename, then give domain a
             #name with the extension
             old_domain_filename = self.domain.filename
@@ -539 +539 @@
 
         # Get X, Y and bed elevation Z
-        Q = domain.quantities['elevation']
+        Q = domain.quantities['elevation']
         X,Y,Z,V = Q.get_vertex_values(xy=True,
                                       precision = self.precision)
@@ -559 +559 @@
         """
         from Scientific.IO.NetCDF import NetCDFFile
-        from time import sleep
+        from time import sleep
 
         #Get NetCDF
-        retries = 0
-        file_open = False
-        while not file_open and retries < 10:
-            try:
+        retries = 0
+        file_open = False
+        while not file_open and retries < 10:
+            try:
                 fid = NetCDFFile(self.filename, 'a') #Open existing file
-            except IOError:
-                #This could happen if someone was reading the file.
-                #In that case, wait a while and try again
-                msg = 'Warning (store_timestep): File %s could not be opened'\
-                %self.filename
-                msg += ' - trying again'
-                print msg
-                retries += 1
-                sleep(1)
-            else:
-               file_open = True
-
-        if not file_open:
-            msg = 'File %s could not be opened for append' %self.filename
-            raise msg
+            except IOError:
+                #This could happen if someone was reading the file.
+                #In that case, wait a while and try again
+                msg = 'Warning (store_timestep): File %s could not be opened'\
+                          %self.filename
+                msg += ' - trying again'
+                print msg
+                retries += 1
+                sleep(1)
+            else:
+                file_open = True
+
+            if not file_open:
+                msg = 'File %s could not be opened for append' %self.filename
+            raise msg
 
 
@@ -595 +595 @@
 
         # Get quantity
-        Q = domain.quantities[name]
+        Q = domain.quantities[name]
         A,V = Q.get_vertex_values(xy=False,
                                   precision = self.precision)
@@ -1025 +1025 @@
 
 def sww2asc(basename_in, basename_out = None,
-            quantity = None,                         
+            quantity = None,
             timestep = None,
             reduction = None,
@@ -1061 +1061 @@
     if quantity is given, out values from quantity otherwise default to
     elevation
-    
+
     if timestep (an index) is given, output quantity at that timestep
 
     if reduction is given use that to reduce quantity over all timesteps.
-    
+
     """
     from Numeric import array, Float, concatenate, NewAxis, zeros,\
          sometrue
-    
+
 
     #FIXME: Should be variable
@@ -1076 +1076 @@
     false_northing = 10000000
     NODATA_value = -9999
-    
+
     if quantity is None:
         quantity = 'elevation'
@@ -1085 +1085 @@
     if basename_out is None:
         basename_out = basename_in + '_%s' %quantity
-        
+
     swwfile = basename_in + '.sww'
     ascfile = basename_out + '.asc'
-    prjfile = basename_out + '.prj'    
+    prjfile = basename_out + '.prj'
 
 
@@ -1103 +1103 @@
     ymin = min(y); ymax = max(y)
     xmin = min(x); xmax = max(x)
-    
+
     number_of_timesteps = fid.dimensions['number_of_timesteps']
     number_of_points = fid.dimensions['number_of_points']
     if origin is None:
-        
+
         # get geo_reference
         #sww files don't have to have a geo_ref
@@ -1114 +1114 @@
         except AttributeError, e:
             geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
-        
+
         xllcorner = geo_reference.get_xllcorner()
         yllcorner = geo_reference.get_yllcorner()
@@ -1122 +1122 @@
         xllcorner = origin[1]
         yllcorner = origin[2]
-        
+
 
     #Get quantity and reduce if applicable
@@ -1129 +1129 @@
     if quantity.lower() == 'depth':
         q = fid.variables['stage'][:] - fid.variables['elevation'][:]
-    else:  
+    else:
         q = fid.variables[quantity][:]
 
-    
+
     if len(q.shape) == 2:
-        if verbose: print 'Reducing quantity %s' %quantity            
-        q_reduced = zeros( number_of_points, Float ) 
+        if verbose: print 'Reducing quantity %s' %quantity
+        q_reduced = zeros( number_of_points, Float )
 
         for k in range(number_of_points):
@@ -1144 +1144 @@
     #Now q has dimension: number_of_points
 
-    
+
     #Write prj file
     if verbose: print 'Writing %s' %prjfile
     prjid = open(prjfile, 'w')
     prjid.write('Projection    %s\n' %'UTM')
-    prjid.write('Zone          %d\n' %zone)                
+    prjid.write('Zone          %d\n' %zone)
     prjid.write('Datum         %s\n' %datum)
     prjid.write('Zunits        NO\n')
@@ -1166 +1166 @@
     newxllcorner = xmin+xllcorner
     newyllcorner = ymin+yllcorner
-    
+
     x = x+xllcorner-newxllcorner
     y = y+yllcorner-newyllcorner
@@ -1172 +1172 @@
     vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
     assert len(vertex_points.shape) == 2
-    
-
-    from Numeric import zeros, Float 
+
+
+    from Numeric import zeros, Float
     grid_points = zeros ( (ncols*nrows, 2), Float )
 
 
     for i in xrange(nrows):
-        yg = i*cellsize        
+        yg = i*cellsize
         for j in xrange(ncols):
             xg = j*cellsize
@@ -1187 +1187 @@
             grid_points[k,0] = xg
             grid_points[k,1] = yg
-            
+
 
     #Interpolate
-    
+
     from least_squares import Interpolation
     from util import inside_polygon
@@ -1197 +1197 @@
     #take forever. With expand_search  set to False, some grid points might
     #miss out....
-    
+
     interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
                            precrop = False, expand_search = False,
@@ -1207 +1207 @@
 
     #Write
-    if verbose: print 'Writing %s' %ascfile    
+    if verbose: print 'Writing %s' %ascfile
     ascid = open(ascfile, 'w')
 
@@ -1215 +1215 @@
     ascid.write('yllcorner     %d\n' %newyllcorner)
     ascid.write('cellsize      %f\n' %cellsize)
-    ascid.write('NODATA_value  %d\n' %NODATA_value)                    
+    ascid.write('NODATA_value  %d\n' %NODATA_value)
 
 
@@ -1221 +1221 @@
     P = interp.mesh.get_boundary_polygon()
     inside_indices = inside_polygon(grid_points, P)
-    
+
     for i in range(nrows):
         if verbose and i%((nrows+10)/10)==0:
             print 'Doing row %d of %d' %(i, nrows)
-            
+
         for j in range(ncols):
             index = (nrows-i-1)*ncols+j
-            
+
             if sometrue(inside_indices == index):
-                ascid.write('%f ' %grid_values[index])                
+                ascid.write('%f ' %grid_values[index])
             else:
                 ascid.write('%d ' %NODATA_value)
-                
+
         ascid.write('\n')
-        
+
     #Close
     ascid.close()
     fid.close()
-    
+
 
 def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,
@@ -1407 +1407 @@
                fail_on_NaN = True,
                NaN_filler = 0,
-               elevation = None,
+               elevation = None,
                inverted_bathymetry = False
                ): #FIXME: Bathymetry should be obtained
                                   #from MOST somehow.
                                   #Alternatively from elsewhere
-                                  #or, as a last resort, 
-                                  #specified here. 
-                                  #The value of -100 will work 
+                                  #or, as a last resort,
+                                  #specified here.
+                                  #The value of -100 will work
                                   #for the Wollongong tsunami
                                   #scenario but is very hacky
@@ -1440 +1440 @@
     ferret2sww uses grid points as vertices in a triangular grid
     counting vertices from lower left corner upwards, then right
-    """ 
+    """
 
     import os
     from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import Float, Int, Int32, searchsorted, zeros, array
-    precision = Float 
+    precision = Float
 
 
@@ -1452 +1452 @@
     file_h = NetCDFFile(basename_in + '_ha.nc', 'r') #Wave amplitude (cm)
     file_u = NetCDFFile(basename_in + '_ua.nc', 'r') #Velocity (x) (cm/s)
-    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)   
-    file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)   
+    file_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)
+    file_e = NetCDFFile(basename_in + '_e.nc', 'r') #Elevation (z) (m)
 
     if basename_out is None:
         swwname = basename_in + '.sww'
     else:
-        swwname = basename_out + '.sww'        
+        swwname = basename_out + '.sww'
 
     times = file_h.variables['TIME']
@@ -1466 +1466 @@
     if mint == None:
         jmin = 0
-    else:    
+    else:
         jmin = searchsorted(times, mint)
-    
+
     if maxt == None:
         jmax=len(times)
-    else:    
+    else:
         jmax = searchsorted(times, maxt)
-        
+
     if minlat == None:
         kmin=0
@@ -1486 +1486 @@
     if minlon == None:
         lmin=0
-    else:    
+    else:
         lmin = searchsorted(longitudes, minlon)
-    
+
     if maxlon == None:
         lmax = len(longitudes)
     else:
-        lmax = searchsorted(longitudes, maxlon)            
-
-
-
-    times = times[jmin:jmax]        
+        lmax = searchsorted(longitudes, maxlon)
+
+
+
+    times = times[jmin:jmax]
     latitudes = latitudes[kmin:kmax]
     longitudes = longitudes[lmin:lmax]
@@ -1531 +1531 @@
     else:
         nan_e = None
-    
+
     #Cleanup
     from Numeric import sometrue
-    
+
     missing = (amplitudes == nan_ha)
     if sometrue (missing):
@@ -1574 +1574 @@
     #######
 
-                
+
 
     number_of_times = times.shape[0]
@@ -1582 +1582 @@
     assert amplitudes.shape[0] == number_of_times
     assert amplitudes.shape[1] == number_of_latitudes
-    assert amplitudes.shape[2] == number_of_longitudes          
+    assert amplitudes.shape[2] == number_of_longitudes
 
     if verbose:
@@ -1596 +1596 @@
         print '    t in [%f, %f], len(t) == %d'\
               %(min(times.flat), max(times.flat), len(times.flat))
-        
-        q = amplitudes.flat        
+
+        q = amplitudes.flat
         name = 'Amplitudes (ha) [cm]'
         print '  %s in [%f, %f]' %(name, min(q), max(q))
-        
+
         q = uspeed.flat
-        name = 'Speeds (ua) [cm/s]'        
+        name = 'Speeds (ua) [cm/s]'
         print '  %s in [%f, %f]' %(name, min(q), max(q))
-        
+
         q = vspeed.flat
-        name = 'Speeds (va) [cm/s]'                
-        print '  %s in [%f, %f]' %(name, min(q), max(q))                
+        name = 'Speeds (va) [cm/s]'
+        print '  %s in [%f, %f]' %(name, min(q), max(q))
 
         q = elevations.flat
-        name = 'Elevations (e) [m]'                
-        print '  %s in [%f, %f]' %(name, min(q), max(q))                
+        name = 'Elevations (e) [m]'
+        print '  %s in [%f, %f]' %(name, min(q), max(q))
 
 
@@ -1617 +1617 @@
     number_of_points = number_of_latitudes*number_of_longitudes
     number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
-    
+
 
     file_h.close()
     file_u.close()
-    file_v.close()    
-    file_e.close()    
+    file_v.close()
+    file_e.close()
 
 
     # NetCDF file definition
     outfile = NetCDFFile(swwname, 'w')
-        
+
     #Create new file
     outfile.institution = 'Geoscience Australia'
@@ -1638 +1638 @@
 
     #For sww compatibility
-    outfile.smoothing = 'Yes' 
+    outfile.smoothing = 'Yes'
     outfile.order = 1
-            
+
     #Start time in seconds since the epoch (midnight 1/1/1970)
     outfile.starttime = starttime = times[0]
     times = times - starttime  #Store relative times
-    
+
     # dimension definitions
     outfile.createDimension('number_of_volumes', number_of_volumes)
@@ -1650 +1650 @@
     outfile.createDimension('number_of_vertices', 3)
     outfile.createDimension('number_of_points', number_of_points)
-                            
-                
+
+
     #outfile.createDimension('number_of_timesteps', len(times))
     outfile.createDimension('number_of_timesteps', len(times))
@@ -1659 +1659 @@
     outfile.createVariable('y', precision, ('number_of_points',))
     outfile.createVariable('elevation', precision, ('number_of_points',))
-            
+
     #FIXME: Backwards compatibility
     outfile.createVariable('z', precision, ('number_of_points',))
     #################################
-        
+
     outfile.createVariable('volumes', Int, ('number_of_volumes',
                                             'number_of_vertices'))
-    
+
     outfile.createVariable('time', precision,
                            ('number_of_timesteps',))
-        
+
     outfile.createVariable('stage', precision,
                            ('number_of_timesteps',
@@ -1677 +1677 @@
                            ('number_of_timesteps',
                             'number_of_points'))
-    
+
     outfile.createVariable('ymomentum', precision,
                            ('number_of_timesteps',
@@ -1691 +1691 @@
     if verbose: print 'Making triangular grid'
     #Check zone boundaries
-    refzone, _, _ = redfearn(latitudes[0],longitudes[0])    
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
 
     vertices = {}
-    i = 0    
+    i = 0
     for k, lat in enumerate(latitudes):       #Y direction
         for l, lon in enumerate(longitudes):  #X direction
@@ -1700 +1700 @@
             vertices[l,k] = i
 
-            zone, easting, northing = redfearn(lat,lon)                
+            zone, easting, northing = redfearn(lat,lon)
 
             msg = 'Zone boundary crossed at longitude =', lon
@@ -1715 +1715 @@
         for k in range(number_of_latitudes-1): #Y direction
             v1 = vertices[l,k+1]
-            v2 = vertices[l,k]            
-            v3 = vertices[l+1,k+1]            
-            v4 = vertices[l+1,k]            
+            v2 = vertices[l,k]
+            v3 = vertices[l+1,k+1]
+            v4 = vertices[l+1,k]
 
             volumes.append([v1,v2,v3]) #Upper element
-            volumes.append([v4,v3,v2]) #Lower element            
-
-    volumes = array(volumes)        
+            volumes.append([v4,v3,v2]) #Lower element
+
+    volumes = array(volumes)
 
     if origin == None:
-        zone = refzone        
+        zone = refzone
         xllcorner = min(x)
         yllcorner = min(y)
@@ -1731 +1731 @@
         zone = origin[0]
         xllcorner = origin[1]
-        yllcorner = origin[2]                
-
-    
+        yllcorner = origin[2]
+
+
     outfile.xllcorner = xllcorner
     outfile.yllcorner = yllcorner
@@ -1747 +1747 @@
             z = elevations
         #FIXME: z should be obtained from MOST and passed in here
-  
+
     from Numeric import resize
     z = resize(z,outfile.variables['z'][:].shape)
@@ -1756 +1756 @@
     outfile.variables['time'][:] = times   #Store time relative
     outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
-    
+
 
 
@@ -1762 +1762 @@
     stage = outfile.variables['stage']
     xmomentum = outfile.variables['xmomentum']
-    ymomentum = outfile.variables['ymomentum'] 
+    ymomentum = outfile.variables['ymomentum']
 
     if verbose: print 'Converting quantities'
@@ -1781 +1781 @@
     if verbose:
         x = outfile.variables['x'][:]
-        y = outfile.variables['y'][:]        
+        y = outfile.variables['y'][:]
         print '------------------------------------------------'
         print 'Statistics of output file:'
@@ -1800 +1800 @@
             q = outfile.variables[name][:].flferret2swwat
             print '    %s in [%f, %f]' %(name, min(q), max(q))
-            
-        
-
-        
-    outfile.close()                
-    
+
+
+
+
+    outfile.close()
+
 
 
@@ -1850 +1850 @@
     Usage: domain = sww2domain('file.sww',t=time (default = last time in file))
 
-    Boundary is not recommended if domian.smooth is not selected, as it       
-    uses unique coordinates, but not unique boundaries. This means that 
+    Boundary is not recommended if domian.smooth is not selected, as it
+    uses unique coordinates, but not unique boundaries. This means that
     the boundary file will not be compatable with the coordiantes, and will
     give a different final boundary, or crash.
@@ -1891 +1891 @@
     except: #AttributeError, e:
         geo_reference = None
-        
+
     if verbose: print '    getting quantities'
     for quantity in fid.variables.keys():
@@ -1929 +1929 @@
     if not boundary is None:
         domain.boundary = boundary
-    
+
     domain.geo_reference = geo_reference
 
@@ -1988 +1988 @@
 
 def interpolated_quantity(saved_quantity,time_interp):
-        
+
     #given an index and ratio, interpolate quantity with respect to time.
     index,ratio = time_interp
-    Q = saved_quantity  
+    Q = saved_quantity
     if ratio > 0:
         q = (1-ratio)*Q[index]+ ratio*Q[index+1]
@@ -2022 +2022 @@
             ratio = 0
         else:
-            #t is now between index and index+1            
+            #t is now between index and index+1
             ratio = (tau - time[index])/(time[index+1] - time[index])
     return (index,ratio)
@@ -2042 +2042 @@
             same_point[i]=point
             #to change all point i references to point j
-        else: 
+        else:
             point_dict[point]=i
             same_point[i]=point
@@ -2079 +2079 @@
             else: new_boundary[(point0,point1)]=label
 
-        boundary = new_boundary               
+        boundary = new_boundary
 
     return coordinates,volumes,boundary
-        
-        
+
+
 
 #OBSOLETE STUFF
@@ -2230 +2230 @@
 
     The momentum is not always stored.
-    
+
     """
     from Scientific.IO.NetCDF import NetCDFFile
@@ -2243 +2243 @@
     stage = fid.variables['stage']     #Water level
     #xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
-    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction    
+    #ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
 
     volumes = fid.variables['volumes'] #Connectivity

inundation/ga/storm_surge/pyvolution/domain.py

@@ -16 +16 @@
     def __init__(self, coordinates, vertices, boundary = None,
                  conserved_quantities = None, other_quantities = None,
-                 tagged_elements = None, geo_reference = None):
+                 tagged_elements = None, geo_reference = None, use_inscribed_circle=False):
 
         Mesh.__init__(self, coordinates, vertices, boundary,
-                      tagged_elements, geo_reference)
+                      tagged_elements, geo_reference, use_inscribed_circle)
 
         from Numeric import zeros, Float

inundation/ga/storm_surge/pyvolution/mesh.py

@@ -130 +130 @@
                 c = sqrt((x2-x0)**2+(y2-y0)**2)
 
-                self.radii[i]=self.areas[i]/(2*(a+b+c))
+                self.radii[i]=2.0*self.areas[i]/(a+b+c)
 
 

inundation/ga/storm_surge/pyvolution/netherlands.py

@@ -63 +63 @@
                 z[i] = -x[i]/20+0.4
 
-            if (x[i] - 0.72)**2 + (y[i] - 0.4)**2 < 0.05**2:# or\
-                   #x[i] > 0.75 and x[i] < 0.9 and y[i] > 0.35 and y[i] < 0.45:
-                z[i] = -x[i]/20+0.4
 
             #Wall
@@ -83 +80 @@
 N = 130 #size = 33800
 N = 600 #Size = 720000
-N = 40
+N = 60
 
 #N = 15
@@ -92 +89 @@
 
 #Create shallow water domain
-domain = Domain(points, vertices, boundary)
+domain = Domain(points, vertices, boundary, use_inscribed_circle=True)
 
 domain.check_integrity()
@@ -105 +102 @@
 
 #Set bed-slope and friction
-inflow_stage = 0.4
+inflow_stage = 0.1
 manning = 0.02
 Z = Weir(inflow_stage)
@@ -120 +117 @@
 else:
     domain.visualise = True
-    domain.visualise_color_stage = True
+    domain.visualise_color_stage = False
     domain.visualise_timer = True
     domain.checkpoint = False
     domain.store = False
+
 
 

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -63 +63 @@
 
     def __init__(self, coordinates, vertices, boundary = None,
-                 tagged_elements = None, geo_reference = None):
+                 tagged_elements = None, geo_reference = None,
+                 use_inscribed_circle=False):
 
         conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
@@ -69 +70 @@
         Generic_domain.__init__(self, coordinates, vertices, boundary,
                                 conserved_quantities, other_quantities,
-                                tagged_elements, geo_reference)
+                                tagged_elements, geo_reference, use_inscribed_circle)
 
         from config import minimum_allowed_height, g

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -263 +263 @@
 
 
-        Q = self.domain.quantities['stage']
+        Q = self.domain.quantities['stage']
         Q0 = Q.vertex_values[:,0]
         Q1 = Q.vertex_values[:,1]
@@ -321 +321 @@
 
         #Check values
-        Q = self.domain.quantities['stage']
+        Q = self.domain.quantities['stage']
         Q0 = Q.vertex_values[:,0]
         Q1 = Q.vertex_values[:,1]
@@ -335 +335 @@
         assert allclose(A[4], (Q1[0] + Q2[1] + Q0[2] +\
                                  Q0[5] + Q2[6] + Q1[7])/6)
-
-
-
 
 
@@ -380 +377 @@
 
         #Check values
-        Q = self.domain.quantities['stage']
+        Q = self.domain.quantities['stage']
         Q0 = Q.vertex_values[:,0]
         Q1 = Q.vertex_values[:,1]
@@ -432 +429 @@
 
             fid.close()
+
         os.remove(self.domain.writer.filename)
 
@@ -596 +594 @@
         #Setup
         self.domain.filename = 'datatest'
-        
+
         prjfile = self.domain.filename + '_elevation.prj'
-        ascfile = self.domain.filename + '_elevation.asc'        
+        ascfile = self.domain.filename + '_elevation.asc'
         swwfile = self.domain.filename + '.sww'
-        
+
         self.domain.set_datadir('.')
         self.domain.format = 'sww'
@@ -607 +605 @@
 
         self.domain.geo_reference = Geo_reference(56,308500,6189000)
-        
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
@@ -630 +628 @@
 
         #Export to ascii/prj files
-        sww2asc(self.domain.filename, 
-                quantity = 'elevation',                         
+        sww2asc(self.domain.filename,
+                quantity = 'elevation',
                 cellsize = cellsize,
                 verbose = False)
@@ -655 +653 @@
         L = lines[3].strip().split()
         assert L[0].strip().lower() == 'zunits'
-        assert L[1].strip().lower() == 'no'                        
+        assert L[1].strip().lower() == 'no'
 
         L = lines[4].strip().split()
         assert L[0].strip().lower() == 'units'
-        assert L[1].strip().lower() == 'meters'                
+        assert L[1].strip().lower() == 'meters'
 
         L = lines[5].strip().split()
@@ -671 +669 @@
         L = lines[7].strip().split()
         assert L[0].strip().lower() == 'yshift'
-        assert L[1].strip().lower() == '10000000'        
+        assert L[1].strip().lower() == '10000000'
 
         L = lines[8].strip().split()
         assert L[0].strip().lower() == 'parameters'
-        
+
 
         #Check asc file
         ascid = open(ascfile)
         lines = ascid.readlines()
-        ascid.close()        
+        ascid.close()
 
         L = lines[0].strip().split()
@@ -688 +686 @@
         L = lines[1].strip().split()
         assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'        
+        assert L[1].strip().lower() == '5'
 
         L = lines[2].strip().split()
@@ -704 +702 @@
         L = lines[5].strip().split()
         assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'        
+        assert L[1].strip().lower() == '-9999'
 
         #Check grid values
         for j in range(5):
-            L = lines[6+j].strip().split()            
+            L = lines[6+j].strip().split()
             y = (4-j) * cellsize
             for i in range(5):
                 assert allclose(float(L[i]), -i*cellsize - y)
-        
+
 
         fid.close()
@@ -718 +716 @@
         #Cleanup
         os.remove(prjfile)
-        os.remove(ascfile)        
+        os.remove(ascfile)
         os.remove(swwfile)
 
@@ -735 +733 @@
         #Setup
         self.domain.filename = 'datatest'
-        
+
         prjfile = self.domain.filename + '_stage.prj'
-        ascfile = self.domain.filename + '_stage.asc'                
+        ascfile = self.domain.filename + '_stage.asc'
         swwfile = self.domain.filename + '.sww'
-        
+
         self.domain.set_datadir('.')
         self.domain.format = 'sww'
@@ -746 +744 @@
 
         self.domain.geo_reference = Geo_reference(56,308500,6189000)
-        
-        
+
+
         sww = get_dataobject(self.domain)
         sww.store_connectivity()
@@ -770 +768 @@
 
         #Export to ascii/prj files
-        sww2asc(self.domain.filename, 
-                quantity = 'stage',                         
+        sww2asc(self.domain.filename,
+                quantity = 'stage',
                 cellsize = cellsize,
                 reduction = min)
@@ -779 +777 @@
         ascid = open(ascfile)
         lines = ascid.readlines()
-        ascid.close()        
+        ascid.close()
 
         L = lines[0].strip().split()
@@ -787 +785 @@
         L = lines[1].strip().split()
         assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'        
+        assert L[1].strip().lower() == '5'
 
         L = lines[2].strip().split()
@@ -803 +801 @@
         L = lines[5].strip().split()
         assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'        
-
-        
+        assert L[1].strip().lower() == '-9999'
+
+
         #Check grid values (where applicable)
         for j in range(5):
             if j%2 == 0:
-                L = lines[6+j].strip().split()            
+                L = lines[6+j].strip().split()
                 jj = 4-j
                 for i in range(5):
@@ -825 +823 @@
         #Cleanup
         os.remove(prjfile)
-        os.remove(ascfile)        
+        os.remove(ascfile)
         #os.remove(swwfile)
 
@@ -837 +835 @@
         This test includes the writing of missing values
         """
-        
+
         import time, os
         from Numeric import array, zeros, allclose, Float, concatenate
@@ -848 +846 @@
         points = [                        [1.0, 1.0],
                               [0.5, 0.5], [1.0, 0.5],
-                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]            
+                  [0.0, 0.0], [0.5, 0.0], [1.0, 0.0]]
 
         vertices = [ [4,1,3], [5,2,4], [1,4,2], [2,0,1]]
-        
+
         #Create shallow water domain
         domain = Domain(points, vertices)
@@ -858 +856 @@
 
         #Set some field values
-        domain.set_quantity('elevation', lambda x,y: -x-y)        
+        domain.set_quantity('elevation', lambda x,y: -x-y)
         domain.set_quantity('friction', 0.03)
 
@@ -885 +883 @@
 
         domain.filename = 'datatest'
-        
+
         prjfile = domain.filename + '_elevation.prj'
-        ascfile = domain.filename + '_elevation.asc'                
+        ascfile = domain.filename + '_elevation.asc'
         swwfile = domain.filename + '.sww'
-        
+
         domain.set_datadir('.')
         domain.format = 'sww'
@@ -895 +893 @@
 
         domain.geo_reference = Geo_reference(56,308500,6189000)
-        
+
         sww = get_dataobject(domain)
         sww.store_connectivity()
@@ -916 +914 @@
         except AttributeError, e:
             geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
-        
+
         #Export to ascii/prj files
-        sww2asc(domain.filename, 
-                quantity = 'elevation',                         
+        sww2asc(domain.filename,
+                quantity = 'elevation',
                 cellsize = cellsize,
                 verbose = False)
@@ -927 +925 @@
         ascid = open(ascfile)
         lines = ascid.readlines()
-        ascid.close()        
+        ascid.close()
 
         L = lines[0].strip().split()
@@ -935 +933 @@
         L = lines[1].strip().split()
         assert L[0].strip().lower() == 'nrows'
-        assert L[1].strip().lower() == '5'        
+        assert L[1].strip().lower() == '5'
 
         L = lines[2].strip().split()
@@ -951 +949 @@
         L = lines[5].strip().split()
         assert L[0].strip() == 'NODATA_value'
-        assert L[1].strip().lower() == '-9999'        
+        assert L[1].strip().lower() == '-9999'
 
 
         #Check grid values
         for j in range(5):
-            L = lines[6+j].strip().split()            
+            L = lines[6+j].strip().split()
             y = (4-j) * cellsize
             for i in range(5):
@@ -971 +969 @@
         #Cleanup
         os.remove(prjfile)
-        os.remove(ascfile)        
+        os.remove(ascfile)
         os.remove(swwfile)
 
@@ -1019 +1017 @@
         stage = fid.variables['stage'][:]
         xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]       
+        ymomentum = fid.variables['ymomentum'][:]
 
         #print ymomentum
-                
+
         assert allclose(stage[0,0], first_value/100)  #Meters
 
@@ -1109 +1107 @@
         # LAT = -34.5, -34.33333, -34.16667, -34 ;
         # ELEVATION = [-1 -2 -3 -4
-        #             -5 -6 -7 -8 
+        #             -5 -6 -7 -8
         #              ...
         #              ...      -16]
@@ -1146 +1144 @@
             fid.variables[lat_name].units='degrees_north'
             fid.variables[lat_name].assignValue(h2_list)
-        
+
             fid.createDimension('TIME',2)
             fid.createVariable('TIME','d',('TIME',))
@@ -1153 +1151 @@
             fid.variables['TIME'].assignValue([0.,1.])
             if fid == fid3: break
-        
+
 
         for fid in [fid4]:
@@ -1173 +1171 @@
         name[fid3]='VA'
         name[fid4]='ELEVATION'
-        
+
         units = {}
         units[fid1]='cm'
@@ -1185 +1183 @@
         values[fid3]=[[[13., 12.,11.], [10.,9.,8.]],[[7.,6.,5.],[4.,3.,2.]]]
         values[fid4]=[[-3000,-3100,-3200],[-4000,-5000,-6000]]
-        
+
         for fid in [fid1,fid2,fid3]:
           fid.createVariable(name[fid],'d',('TIME',lat_name,long_name))
@@ -1201 +1199 @@
             fid.variables[name[fid]].missing_value = -99999999.
 
-        
+
         fid1.sync(); fid1.close()
         fid2.sync(); fid2.close()
@@ -1210 +1208 @@
         fid2 = NetCDFFile('test_e.nc','r')
         fid3 = NetCDFFile('test_va.nc','r')
-        
+
 
         first_amp = fid1.variables['HA'][:][0,0,0]
@@ -1240 +1238 @@
         stage = fid.variables['stage'][:]
         xmomentum = fid.variables['xmomentum'][:]
-        ymomentum = fid.variables['ymomentum'][:]       
+        ymomentum = fid.variables['ymomentum'][:]
 
         #print ymomentum
@@ -1343 +1341 @@
         yiel=0.01
         points, vertices, boundary = rectangular(10,10)
-        
+
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
@@ -1465 +1463 @@
             #print bit
             assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
-        
+
 
     def test_sww2domain2(self):
@@ -1507 +1505 @@
 
         domain.check_integrity()
-        
+
         for t in domain.evolve(yieldstep = 1, finaltime = 2.0):
             pass
             #domain.write_time()
 
-        
+
 
         ##################################
@@ -1563 +1561 @@
         coordinates2,volumes2,boundary2=weed(coordinates1,volumes1,boundary1)
 
-        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}        
+        points2 = {(0.,0.):None,(1.,0.):None,(1.,1.):None,(2.,0.):None}
 
         assert len(points2)==len(coordinates2)
@@ -1577 +1575 @@
 
 
-     #FIXME This fails - smooth makes the comparism too hard for allclose 
+     #FIXME This fails - smooth makes the comparism too hard for allclose
     def ztest_sww2domain3(self):
         ################################################
@@ -1590 +1588 @@
         yiel=0.01
         points, vertices, boundary = rectangular(10,10)
-        
+
         #Create shallow water domain
         domain = Domain(points, vertices, boundary)
@@ -1714 +1712 @@
             print bit
             assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
-        
+
 
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Data_Manager,'test')    
+    suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2domain')
     runner = unittest.TextTestRunner()