source: development/steve/sww_viewer_vtk.py @ 2386

Last change on this file since 2386 was 2386, checked in by steve, 18 years ago
File size: 6.4 KB
Line 
1"""
2Created: 11/01/2006
3Author: John Jakeman
4
5Program Visualises .sww files
6.sww: Netcdf format for storing model output f(t,x,y)
7"""
8import time as systime
9def update_viewer(time_interp, yeildstep, t):
10    stage = interpolated_quantity(fid.variables['stage'][:],time_interp)
11    domain.set_quantity('stage',stage)
12   
13    ratio = time_interp[1]
14    domain.set_time(t+ratio*yieldstep)
15   
16    #domain.visualiser.update_quantity('stage')
17    #domain.visualiser.update_timer()
18    if domain.visualise is True:
19        domain.visualiser.redraw_ready.set()
20        domain.visualiser.idle.wait()
21        domain.visualiser.idle.clear()
22        domain.visualiser.unpaused.wait()   
23   
24    print 'time = ',domain.write_time()
25    #systime.sleep(0.5)
26
27import project
28from pyvolution.data_manager import get_time_interp, interpolated_quantity
29from caching import cache
30filename = project.filename[:-3] + '.sww' #filename is of form filename.py
31#filename = 'run_new_meribula.sww'
32#filename = 'results29_01_2006.sww'
33filename = 'cairns250m.sww'
34verbose = True
35very_verbose = False
36boundary = None
37t = 0.0
38fail_if_NaN=False
39NaN_filler=0
40
41"""
42The following is borrowed from data_manager.sww2domain
43"""
44NaN=9.969209968386869e+036
45#initialise NaN.
46
47from Scientific.IO.NetCDF import NetCDFFile
48from pyvolution.shallow_water_vtk import Domain
49from Numeric import asarray, transpose, resize
50
51if verbose: print 'Reading from ', filename
52fid = NetCDFFile(filename, 'r')    #Open existing file for read
53time = fid.variables['time']       #Timesteps
54
55#if t is None:
56#    t = time[-1]
57time_interp = 0, 0.0
58#time_interp = get_time_interp(time,t)
59
60# Get the variables as Numeric arrays
61x = fid.variables['x'][:]             #x-coordinates of vertices
62y = fid.variables['y'][:]             #y-coordinates of vertices
63elevation = fid.variables['elevation']     #Elevation
64stage = fid.variables['stage']     #Water level
65xmomentum = fid.variables['xmomentum']   #Momentum in the x-direction
66ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
67
68starttime = fid.starttime[0]
69volumes = fid.variables['volumes'][:] #Connectivity
70coordinates=transpose(asarray([x.tolist(),y.tolist()]))
71
72conserved_quantities = []
73interpolated_quantities = {}
74other_quantities = []
75
76# get geo_reference
77#sww files don't have to have a geo_ref
78try:
79    geo_reference = Geo_reference(NetCDFObject=fid)
80except: #AttributeError, e:
81    geo_reference = None
82
83if verbose: print '    getting quantities'
84for quantity in fid.variables.keys():
85    dimensions = fid.variables[quantity].dimensions
86    if 'number_of_timesteps' in dimensions:
87        conserved_quantities.append(quantity)
88        interpolated_quantities[quantity]=\
89              interpolated_quantity(fid.variables[quantity][:],time_interp)
90    else:
91        other_quantities.append(quantity)
92
93other_quantities.remove('x')
94other_quantities.remove('y')
95other_quantities.remove('z')
96other_quantities.remove('volumes')
97
98conserved_quantities.remove('time')
99
100if verbose: print '    building domain'
101#    From domain.Domain:
102#    domain = Domain(coordinates, volumes,\
103#                    conserved_quantities = conserved_quantities,\
104#                    other_quantities = other_quantities,zone=zone,\
105#                    xllcorner=xllcorner, yllcorner=yllcorner)
106
107#   From shallow_water.Domain:
108coordinates=coordinates.tolist()
109volumes=volumes.tolist()
110#FIXME:should this be in mesh?(peter row)
111if fid.smoothing == 'Yes': unique = False
112else: unique = True
113if unique:
114    coordinates,volumes,boundary=weed(coordinates,volumes,boundary)
115
116
117domain = cache(Domain, (coordinates, volumes, boundary))
118
119if not boundary is None:
120    domain.boundary = boundary
121
122domain.geo_reference = geo_reference
123
124domain.starttime=float(starttime)+float(t)
125domain.time=0.0
126
127for quantity in other_quantities:
128    try:
129        NaN = fid.variables[quantity].missing_value
130    except:
131        pass #quantity has no missing_value number
132    X = fid.variables[quantity][:]
133    if very_verbose:
134        print '       ',quantity
135        print '        NaN =',NaN
136        print '        max(X)'
137        print '       ',max(X)
138        print '        max(X)==NaN'
139        print '       ',max(X)==NaN
140        print ''
141    if (max(X)==NaN) or (min(X)==NaN):
142        if fail_if_NaN:
143            msg = 'quantity "%s" contains no_data entry'%quantity
144            raise msg
145        else:
146            data = (X<>NaN)
147            X = (X*data)+(data==0)*NaN_filler
148    if unique:
149        X = resize(X,(len(X)/3,3))
150    domain.set_quantity(quantity,X)
151#
152for quantity in conserved_quantities:
153    try:
154        NaN = fid.variables[quantity].missing_value
155    except:
156        pass #quantity has no missing_value number
157    X = interpolated_quantities[quantity]
158    if very_verbose:
159        print '       ',quantity
160        print '        NaN =',NaN
161        print '        max(X)'
162        print '       ',max(X)
163        print '        max(X)==NaN'
164        print '       ',max(X)==NaN
165        print ''
166    if (max(X)==NaN) or (min(X)==NaN):
167        if fail_if_NaN:
168            msg = 'quantity "%s" contains no_data entry'%quantity
169            raise msg
170        else:
171            data = (X<>NaN)
172            X = (X*data)+(data==0)*NaN_filler
173    if unique:
174        X = resize(X,(X.shape[0]/3,3))
175    domain.set_quantity(quantity,X)
176
177
178"""
179The following initialises visualiser and uses .sww file to update stage
180at each timestep
181"""
182
183domain.initialise_visualiser()
184#domain.visualiser.setup_all()
185domain.visualiser.scale_z['stage'] = 20.0
186domain.visualiser.coloring['stage'] = True
187domain.visualiser.scale_z['elevation'] = 20.0
188domain.visualiser.coloring['elevation'] = False
189domain.visualiser.updating['elevation'] = True
190#domain.visualiser.setup['stage'] = False
191#domain.visualiser.updating['stage'] = False
192
193
194domain.visualiser.start()
195# Things go haywire if we start evolving before the vis is ready
196domain.visualiser.idle.wait()
197domain.visualiser.idle.clear()
198
199index = 0 
200ratio = 0.0
201num_subt = 1 # use 1 if only want to display at yieldstep timeinterval
202
203if len(time) > 1:
204    yieldstep = time[1]-time[0]
205else:
206    yieldstep = 0.0
207
208for t in time[:]:
209    if t < time[-1]:
210        if index%10 == 0:
211            for subt in range(0,num_subt):
212                ratio = 1-float(num_subt-subt)/float(num_subt)
213                #print ratio
214                time_interp = index, ratio
215                update_viewer(time_interp, yieldstep, t)
216    else:
217        time_interp = index, 0.0
218        update_viewer(time_interp, yieldstep, t)
219           
220    index = index + 1
221   
222print 'Finished'   
223fid.close()
224
Note: See TracBrowser for help on using the repository browser.