Changeset 1083 for inundation/ga/storm_surge/pyvolution/data_manager.py

Timestamp:

Mar 16, 2005, 1:01:22 PM (20 years ago)

Author:

prow

Message:

added sww2domain. pretty sure it works ok.

File:

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

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -1533 +1533 @@
 
 
+def sww2domain(filename,t=None):
+    """Read sww Net CDF file containing Shallow Water Wave simulation
+
+    Quantities stage, elevation, xmomentum and ymomentum.
+
+    The momentum is not always stored.
+
+    """
+    from Scientific.IO.NetCDF import NetCDFFile
+    from domain import Domain
+    from Numeric import asarray, transpose
+    print 'Reading from ', filename
+    fid = NetCDFFile(filename, 'r')    #Open existing file for read
+    time = fid.variables['time']       #Timesteps
+    if t is None:
+        t = time[-1]
+    time_interp = get_time_interp(time,t)
+################################
+########################################    
+    # Get the variables as Numeric arrays
+    x = fid.variables['x'][:]             #x-coordinates of vertices
+    y = fid.variables['y'][:]             #y-coordinates of vertices
+    elevation = fid.variables['elevation']     #Elevation
+    stage = fid.variables['stage']     #Water level
+    xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
+    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction    
+#################################
+    xllcorner = fid.xllcorner[0]
+    yllcorner = fid.yllcorner[0]
+    starttime = fid.starttime[0]
+    zone = fid.zone
+    volumes = fid.variables['volumes'][:] #Connectivity
+    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
+#
+    conserved_quantities = []
+    interpolated_quantities = {}
+    other_quantities = []
+#
+    print '    interpolating quantities'
+    for quantity in fid.variables.keys():
+        dimensions = fid.variables[quantity].dimensions
+        if 'number_of_timesteps' in dimensions:
+            conserved_quantities.append(quantity)
+            interpolated_quantities[quantity]=\
+                  interpolated_quantity(fid.variables[quantity][:],time_interp)
+        else: other_quantities.append(quantity)
+#
+    other_quantities.remove('x')
+    other_quantities.remove('y')
+    other_quantities.remove('z')
+    other_quantities.remove('volumes')
+#
+    conserved_quantities.remove('time')
+#
+    print other_quantities
+    print conserved_quantities
+    print '    building domain'
+    domain = Domain(coordinates, volumes,\
+                    conserved_quantities = conserved_quantities,\
+                    other_quantities = other_quantities,zone=zone,\
+                    xllcorner=xllcorner, yllcorner=yllcorner)
+    domain.starttime=starttime
+    domain.time=t
+    return domain
+
+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  
+    if ratio > 0:
+        q = (1-ratio)*Q[index]+ ratio*Q[index+1]
+    else:
+         q = Q[index]
+    #Return vector of interpolated values
+    return q
+
+def get_time_interp(time,t=None):
+    #Finds the ratio and index for time interpolation.
+    #It is borrowed from previous pyvolution code.
+    if t is None:
+        t=time[-1]
+        index = -1
+        ratio = 0.
+    else:
+        T = time
+        tau = t
+        index=0
+        msg = 'Time interval derived from file %s [%s:%s]'\
+            %('FIXMEfilename', T[0], T[-1])
+        msg += ' does not match model time: %s' %tau
+        if tau < time[0]: raise msg
+        if tau > time[-1]: raise msg
+        while tau > time[index]: index += 1
+        while tau < time[index]: index -= 1
+        if tau == time[index]:
+            #Protect against case where tau == time[-1] (last time)
+            # - also works in general when tau == time[i]
+            ratio = 0
+        else:
+            #t is now between index and index+1            
+            ratio = (tau - time[index])/(time[index+1] - time[index])
+    return (index,ratio)
+
+
+
 #OBSOLETE STUFF
 #Native checkpoint format.