Changeset 1044

Timestamp:

Mar 8, 2005, 10:19:53 PM (20 years ago)

Author:

ole

Message:

Added temporary fix for momentum computation in ferret2sww

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (18 diffs)
• test_data_manager.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -1129 +1129 @@
                minlon = None, maxlon = None,
                mint = None, maxt = None, mean_stage = 0,
-               origin = None, zscale = 1):
+               origin = None, zscale = 1, 
+               mean_bathymetry = -100): #FIXME: Bathymetry should be obtained
+                                        #from MOST somehow.
+                                        #Alternatively from elsewhere
+                                        #or, as a last resort, 
+                                        #specified here. 
+                                        #The value of -100 will work 
+                                        #for the Wollongong tsunami
+                                        #scenario but is very hacky
     """Convert 'Ferret' NetCDF format for wave propagation to
     sww format native to pyvolution.
@@ -1153 +1161 @@
     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 
 
 
@@ -1165 +1173 @@
     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_v = NetCDFFile(basename_in + '_va.nc', 'r') #Velocity (y) (cm/s)   
 
     if basename_out is None:
         swwname = basename_in + '.sww'
     else:
-        swwname = basename_out + '.sww'
+        swwname = basename_out + '.sww'        
 
     times = file_h.variables['TIME']
@@ -1178 +1186 @@
     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
@@ -1198 +1206 @@
     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]
@@ -1214 +1222 @@
     amplitudes = file_h.variables['HA'][jmin:jmax, kmin:kmax, lmin:lmax]
     xspeed = file_u.variables['UA'][jmin:jmax, kmin:kmax, lmin:lmax]
-    yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]
+    yspeed = file_v.variables['VA'][jmin:jmax, kmin:kmax, lmin:lmax]    
 
     number_of_times = times.shape[0]
@@ -1222 +1230 @@
     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          
 
     #print times
@@ -1229 +1237 @@
 
     #print 'MIN', min(min(min(amplitudes)))
-    #print 'MAX', max(max(max(amplitudes)))
+    #print 'MAX', max(max(max(amplitudes)))    
 
     #print number_of_latitudes, number_of_longitudes
     number_of_points = number_of_latitudes*number_of_longitudes
     number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
-
+    
     #print file_h.dimensions.keys()
-    #print file_h.variables.keys()
+    #print file_h.variables.keys()    
 
     file_h.close()
     file_u.close()
-    file_v.close()
+    file_v.close()    
 
 
@@ -1246 +1254 @@
     # NetCDF file definition
     outfile = NetCDFFile(swwname, 'w')
-
+        
     #Create new file
     outfile.institution = 'Geoscience Australia'
@@ -1256 +1264 @@
 
     #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 = times[0]
-
+    
     # dimension definitions
     outfile.createDimension('number_of_volumes', number_of_volumes)
@@ -1267 +1275 @@
     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))
@@ -1276 +1284 @@
     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',
@@ -1294 +1302 @@
                            ('number_of_timesteps',
                             'number_of_points'))
-
+    
     outfile.createVariable('ymomentum', precision,
                            ('number_of_timesteps',
@@ -1308 +1316 @@
 
     #Check zone boundaries
-    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0])    
 
     vertices = {}
-    for k, lat in enumerate(latitudes):
-        for l, lon in enumerate(longitudes):
+    for k, lat in enumerate(latitudes):    
+        for l, lon in enumerate(longitudes):    
 
             vertices[l,k] = i
 
-            zone, easting, northing = redfearn(lat,lon)
+            zone, easting, northing = redfearn(lat,lon)                
 
             msg = 'Zone boundary crossed at longitude =', lon
@@ -1331 +1339 @@
         for k in range(number_of_latitudes-1):
             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)
@@ -1347 +1355 @@
         zone = origin[0]
         xllcorner = origin[1]
-        yllcorner = origin[2]
-
-
+        yllcorner = origin[2]                
+
+    
     outfile.xllcorner = xllcorner
     outfile.yllcorner = yllcorner
@@ -1357 +1365 @@
     outfile.variables['y'][:] = y - yllcorner
     outfile.variables['z'][:] = 0.0
-    outfile.variables['elevation'][:] = 0.0
-    outfile.variables['time'][:] = times
+    outfile.variables['elevation'][:] = 0.0 #Grrrrrrr
+    outfile.variables['time'][:] = times        
     outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
-
+    
 
 
@@ -1366 +1374 @@
     stage = outfile.variables['stage']
     xmomentum = outfile.variables['xmomentum']
-    ymomentum = outfile.variables['ymomentum']
-
+    ymomentum = outfile.variables['ymomentum']        
+
+    if mean_bathymetry is not None:
+        z = mean_bathymetry
+    else:
+        pass
+        #FIXME: z should be obtained from MOST and passed in here       
+    
     for j in range(len(times)):
         i = 0
         for k in range(number_of_latitudes):
-            for l in range(number_of_longitudes):
-                h = zscale*amplitudes[j,k,l]/100 + mean_stage
-                stage[j,i] = h
-                xmomentum[j,i] = xspeed[j,k,l]/100*h
-                ymomentum[j,i] = yspeed[j,k,l]/100*h
+            for l in range(number_of_longitudes):                    
+                w = zscale*amplitudes[j,k,l]/100 + mean_stage
+                stage[j,i] = w
+
+                h = w-z  
+    
+                xmomentum[j,i] = xspeed[j,k,l]/100*h
+                ymomentum[j,i] = yspeed[j,k,l]/100*h                
                 i += 1
-
-    outfile.close()
+        
+    outfile.close()                
+    
 
 

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -602 +602 @@
 
         fid = NetCDFFile('small_ha.nc')
-
         first_value = fid.variables['HA'][:][0,0,0]
         fourth_value = fid.variables['HA'][:][0,0,3]
+
 
         #Call conversion (with zero origin)
@@ -627 +627 @@
         #Check first value
         stage = fid.variables['stage'][:]
-
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]       
+
+        #print ymomentum
+                
         assert allclose(stage[0,0], first_value/100)  #Meters