Index: development/steve/sww_viewer_vtk.py
===================================================================
--- development/steve/sww_viewer_vtk.py	(revision 2386)
+++ development/steve/sww_viewer_vtk.py	(revision 2386)
@@ -0,0 +1,224 @@
+"""
+Created: 11/01/2006
+Author: John Jakeman
+
+Program Visualises .sww files
+.sww: Netcdf format for storing model output f(t,x,y)
+"""
+import time as systime
+def update_viewer(time_interp, yeildstep, t):
+    stage = interpolated_quantity(fid.variables['stage'][:],time_interp)
+    domain.set_quantity('stage',stage)
+    
+    ratio = time_interp[1]
+    domain.set_time(t+ratio*yieldstep)
+    
+    #domain.visualiser.update_quantity('stage')
+    #domain.visualiser.update_timer()
+    if domain.visualise is True:
+        domain.visualiser.redraw_ready.set()
+        domain.visualiser.idle.wait()
+        domain.visualiser.idle.clear()
+        domain.visualiser.unpaused.wait()    
+    
+    print 'time = ',domain.write_time()
+    #systime.sleep(0.5)
+
+import project
+from pyvolution.data_manager import get_time_interp, interpolated_quantity
+from caching import cache
+filename = project.filename[:-3] + '.sww' #filename is of form filename.py
+#filename = 'run_new_meribula.sww'
+#filename = 'results29_01_2006.sww'
+filename = 'cairns250m.sww'
+verbose = True
+very_verbose = False
+boundary = None
+t = 0.0
+fail_if_NaN=False
+NaN_filler=0
+
+"""
+The following is borrowed from data_manager.sww2domain
+"""
+NaN=9.969209968386869e+036
+#initialise NaN.
+
+from Scientific.IO.NetCDF import NetCDFFile
+from pyvolution.shallow_water_vtk import Domain
+from Numeric import asarray, transpose, resize
+
+if verbose: 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 = 0, 0.0
+#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
+
+starttime = fid.starttime[0]
+volumes = fid.variables['volumes'][:] #Connectivity
+coordinates=transpose(asarray([x.tolist(),y.tolist()]))
+
+conserved_quantities = []
+interpolated_quantities = {}
+other_quantities = []
+
+# get geo_reference
+#sww files don't have to have a geo_ref
+try:
+    geo_reference = Geo_reference(NetCDFObject=fid)
+except: #AttributeError, e:
+    geo_reference = None
+
+if verbose: print '    getting 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')
+
+if verbose: print '    building domain'
+#    From domain.Domain:
+#    domain = Domain(coordinates, volumes,\
+#                    conserved_quantities = conserved_quantities,\
+#                    other_quantities = other_quantities,zone=zone,\
+#                    xllcorner=xllcorner, yllcorner=yllcorner)
+
+#   From shallow_water.Domain:
+coordinates=coordinates.tolist()
+volumes=volumes.tolist()
+#FIXME:should this be in mesh?(peter row)
+if fid.smoothing == 'Yes': unique = False
+else: unique = True
+if unique:
+    coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
+
+
+domain = cache(Domain, (coordinates, volumes, boundary))
+
+if not boundary is None:
+    domain.boundary = boundary
+
+domain.geo_reference = geo_reference
+
+domain.starttime=float(starttime)+float(t)
+domain.time=0.0
+
+for quantity in other_quantities:
+    try:
+        NaN = fid.variables[quantity].missing_value
+    except:
+        pass #quantity has no missing_value number
+    X = fid.variables[quantity][:]
+    if very_verbose:
+        print '       ',quantity
+        print '        NaN =',NaN
+        print '        max(X)'
+        print '       ',max(X)
+        print '        max(X)==NaN'
+        print '       ',max(X)==NaN
+        print ''
+    if (max(X)==NaN) or (min(X)==NaN):
+        if fail_if_NaN:
+            msg = 'quantity "%s" contains no_data entry'%quantity
+            raise msg
+        else:
+            data = (X<>NaN)
+            X = (X*data)+(data==0)*NaN_filler
+    if unique:
+        X = resize(X,(len(X)/3,3))
+    domain.set_quantity(quantity,X)
+#
+for quantity in conserved_quantities:
+    try:
+        NaN = fid.variables[quantity].missing_value
+    except:
+        pass #quantity has no missing_value number
+    X = interpolated_quantities[quantity]
+    if very_verbose:
+        print '       ',quantity
+        print '        NaN =',NaN
+        print '        max(X)'
+        print '       ',max(X)
+        print '        max(X)==NaN'
+        print '       ',max(X)==NaN
+        print ''
+    if (max(X)==NaN) or (min(X)==NaN):
+        if fail_if_NaN:
+            msg = 'quantity "%s" contains no_data entry'%quantity
+            raise msg
+        else:
+            data = (X<>NaN)
+            X = (X*data)+(data==0)*NaN_filler
+    if unique:
+        X = resize(X,(X.shape[0]/3,3))
+    domain.set_quantity(quantity,X)
+
+
+"""
+The following initialises visualiser and uses .sww file to update stage
+at each timestep
+"""
+
+domain.initialise_visualiser()
+#domain.visualiser.setup_all()
+domain.visualiser.scale_z['stage'] = 20.0
+domain.visualiser.coloring['stage'] = True
+domain.visualiser.scale_z['elevation'] = 20.0
+domain.visualiser.coloring['elevation'] = False
+domain.visualiser.updating['elevation'] = True
+#domain.visualiser.setup['stage'] = False
+#domain.visualiser.updating['stage'] = False
+
+
+domain.visualiser.start()
+# Things go haywire if we start evolving before the vis is ready
+domain.visualiser.idle.wait()
+domain.visualiser.idle.clear()
+
+index = 0 
+ratio = 0.0
+num_subt = 1 # use 1 if only want to display at yieldstep timeinterval
+
+if len(time) > 1:
+    yieldstep = time[1]-time[0]
+else:
+    yieldstep = 0.0
+
+for t in time[:]:
+    if t < time[-1]:
+        if index%10 == 0:
+            for subt in range(0,num_subt):
+                ratio = 1-float(num_subt-subt)/float(num_subt)
+                #print ratio
+                time_interp = index, ratio
+                update_viewer(time_interp, yieldstep, t)
+    else:
+        time_interp = index, 0.0
+        update_viewer(time_interp, yieldstep, t)
+            
+    index = index + 1
+    
+print 'Finished'   
+fid.close()
+
