Changeset 877

Timestamp:

Feb 14, 2005, 6:09:05 PM (20 years ago)

Author:

ole

Message:

Added conversion from ferret let/lon format to sww (UTM)

File:

• inundation/ga/storm_surge/pyvolution/data_manager.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -996 +996 @@
 
 
+
+
+
+
+def ferret2sww(basefilename, verbose=False,
+               minlat = None, maxlat =None,
+               minlon = None, maxlon =None,
+               mint = None, maxt = None, mean_stage = 0):
+    """Convert 'Ferret' NetCDF format for wave propagation to
+    sww format native to pyvolution.
+
+    Specify only basefilename and read files of the form
+    basefilename_ha.nc, basefilename_ua.nc, basefilename_va.nc containing
+    relative height, x-velocity and y-velocity, respectively.
+
+    Also convert latitude and longitude to UTM. All coordinates are
+    assumed to be given in the GDA94 datum.
+
+
+    min's and max's: If omitted - full extend is used.
+    To include a value min may equal it, while max must exceed it.
+    """ 
+
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, Int, searchsorted, zeros, array
+    precision = Float 
+
+
+    #Get NetCDF data
+    file_h = NetCDFFile(basefilename + '_ha.nc', 'r') #Wave amplitude (cm)
+    file_u = NetCDFFile(basefilename + '_ua.nc', 'r') #Velocity (x) (cm/s)
+    file_v = NetCDFFile(basefilename + '_va.nc', 'r') #Velocity (y) (cm/s)   
+
+    swwname = basefilename + '.sww'
+
+    times = file_h.variables['TIME']
+    latitudes = file_h.variables['LAT']
+    longitudes = file_h.variables['LON']
+
+    if mint == None:
+        jmin = 0
+    else:    
+        jmin = searchsorted(times, mint)
+    
+    if maxt == None:
+        jmax=len(times)
+    else:    
+        jmax = searchsorted(times, maxt)
+        
+    if minlat == None:
+        kmin=0
+    else:
+        kmin = searchsorted(latitudes, minlat)
+
+    if maxlat == None:
+        kmax = len(latitudes)
+    else:
+        kmax = searchsorted(latitudes, maxlat)
+
+    if minlon == None:
+        lmin=0
+    else:    
+        lmin = searchsorted(longitudes, minlon)
+    
+    if maxlon == None:
+        lmax = len(longitudes)
+    else:
+        lmax = searchsorted(longitudes, maxlon)            
+        
+    latitudes = latitudes[kmin:kmax]
+    longitudes = longitudes[lmin:lmax]
+    times = times[jmin:jmax]
+    
+    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]    
+
+    number_of_latitudes = latitudes.shape[0]
+    number_of_longitudes = longitudes.shape[0]
+
+
+    #print times
+    #print latitudes
+    #print longitudes
+
+    #print 'MIN', min(min(min(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()    
+
+    if verbose: print 'Store to SWW file %s' %swwname
+    # NetCDF file definition
+    outfile = NetCDFFile(swwname, 'w')
+        
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'Converted from Ferret files: %s, %s, %s'\
+                          %(basefilename + '_ha.nc',
+                            basefilename + '_ua.nc',
+                            basefilename + '_va.nc')
+
+
+    #For sww compatibility
+    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)
+
+    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))
+
+    # variable definitions
+    outfile.createVariable('x', precision, ('number_of_points',))
+    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',
+                            'number_of_points'))
+
+    outfile.createVariable('xmomentum', precision,
+                           ('number_of_timesteps',
+                            'number_of_points'))
+    
+    outfile.createVariable('ymomentum', precision,
+                           ('number_of_timesteps',
+                            'number_of_points'))
+
+
+    #Store
+    from coordinate_transforms.redfearn import redfearn
+    x = zeros(number_of_points, Float)  #Easting
+    y = zeros(number_of_points, Float)  #Northing
+    #volumes = zeros(number_of_volumes, Int)
+    i = 0
+
+    #Check zone boundaries
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
+
+    vertices = {}
+    for l, lon in enumerate(longitudes):    
+        for k, lat in enumerate(latitudes):
+            vertices[l,k] = i
+            
+            zone, easting, northing = redfearn(lat,lon)
+
+            msg = 'Zone boundary crossed at longitude =', lon
+            assert zone == refzone, msg
+            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
+            x[i] = easting
+            y[i] = northing
+            i += 1
+
+
+    #Construct 2 triangles per 'rectangular' element
+    volumes = []
+    for l in range(number_of_longitudes-1):
+        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]            
+
+            volumes.append([v1,v2,v3]) #Upper element
+            volumes.append([v4,v3,v2]) #Lower element            
+
+    volumes = array(volumes)        
+
+    origin = (min(x), min(y))
+    outfile.minx = origin[0]
+    outfile.miny = origin[1]    
+
+    #print x - origin[0]
+    #print y - origin[1]
+
+    #print len(x)
+    #print number_of_points
+    
+    outfile.variables['x'][:] = x - origin[0]
+    outfile.variables['y'][:] = y - origin[1]
+    outfile.variables['z'][:] = 0.0
+    outfile.variables['elevation'][:] = 0.0    
+    outfile.variables['volumes'][:] = volumes
+    outfile.variables['time'][:] = times    
+
+    #Time stepping
+    stage = outfile.variables['stage']
+    xmomentum = outfile.variables['xmomentum']
+    ymomentum = outfile.variables['ymomentum']        
+
+    for j in range(len(times)):
+        i = 0
+        for l in range(number_of_longitudes):    
+            for k in range(number_of_latitudes):
+                h = 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                
+                i += 1
+        
+    outfile.close()                
+    
+    
+
+
+
 #OBSOLETE STUFF
 #Native checkpoint format.