source: anuga_work/development/steve/sww_viewer_vtk.py @ 7818

Last change on this file since 7818 was 3514, checked in by duncan, 18 years ago

Hi all,
I'm doing a change in the anuga structure, moving the code to

\anuga_core\source\anuga

After you have done an svn update, the PYTHONPATH has to be changed to;
PYTHONPATH = anuga_core/source/

This is part of changes required to make installation of anuga quicker and reducing the size of our sandpits.

If any imports are broken, try fixing them. With adding anuga. to them for example. If this seems to have really broken things, email/phone me.

Cheers
Duncan

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 anuga.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 anuga.pyvolution.data_manager.sww2domain
43"""
44NaN=9.969209968386869e+036
45#initialise NaN.
46
47from Scientific.IO.NetCDF import NetCDFFile
48from anuga.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.