source: anuga_validation/convergence_study/animatesww2d_alt.py @ 4700

def animatesww2d_alt(sww_filename, movie_filename, range):
    """Plot cross section of model output
    """

    # Read SWW file     
    from Scientific.IO.NetCDF import NetCDFFile
    fid = NetCDFFile(sww_filename, 'r')
    
    x = fid.variables['x']
    y = fid.variables['y']
    volumes = fid.variables['volumes']
    
    elevation = fid.variables['elevation']
    time = fid.variables['time']
    stage = fid.variables['stage']
    xmomentum = fid.variables['xmomentum']
    ymomentum = fid.variables['ymomentum']
    

    # find indices for where y = 0

    # 
    # Set up plot

    spacing=10 # grid interval in metres
    hmax=5     # ??

    
#figure(1);

#if nargin > 1
#    figure(1);set(gcf,'Position',[200 400 800 300]);
#    mov=avifile(movie_filename,'FPS',5)
#end

#i=find(V_y==0);

#[xi,yi]=meshgrid(min(V_x):spacing:max(V_x),min(V_y):spacing:max(V_y));

#elev=griddata(V_x,V_y,V_elevation,xi,yi);

#if nargin < 3
#    range=[min(V_x) max(V_x) min(min(V_stage))-1 max(max(V_stage))+1];
#end

#m=V_stage(1,i);

#for t=1:length(V_time)
    
#    m=max(m,V_stage(t,i));
#    plot(V_x(i),V_stage(t,i),V_x(i),V_elevation(i),V_x(i),m,'.');

#    axis(range)
    

#    title(num2str(V_time(t)));drawnow;

#    if nargin>1
#        F=getframe(gcf);
#        mov=addframe(mov,F);   
#    end
   
#end

#if nargin > 1
#    mov=close(mov);
#end

'''
if __name__ == '__main__':

    import sys
    animatesww2d_alt(sys.argv[1], 'cross_section_plot', 5)

'''

'''
How do I make a movie with matplotlib?
If you want to take an animated plot and turn it into a movie, the best approach is to save a series of image files (eg PNG) and use an external tool to convert them to a movie. There is a matplotlib tutorial on this subject here. You can use mencoder, which is part of the mplayer suite for this 
#fps (frames per second) controls the play speed
mencoder 'mf://*.png' -mf type=png:fps=10 -ovc \
   lavc -lavcopts vcodec=wmv2 -oac copy -o animation.avi

The swiss army knife of image tools, ImageMagick's convert, works for this as well.
Here is a simple example script that saves some PNGs, makes them into a movie, and then cleans up. 

import os, sys
from pylab import *

files = []
figure(figsize=(5,5))
ax = subplot(111)
for i in range(50):  # 50 frames
    cla()
    imshow(rand(5,5), interpolation='nearest')
    fname = '_tmp%03d.png'%i
    print 'Saving frame', fname
    savefig(fname)
    files.append(fname)

print 'Making movie animation.mpg - this make take a while'
os.system("mencoder 'mf://_tmp*.png' -mf type=png:fps=10 \
  -ovc lavc -lavcopts vcodec=wmv2 -oac copy -o animation.mpg")

# cleanup
for fname in files: os.remove(fname)

'''

def animatesww2d(swwfile):
    """Read in sww file and plot cross section of model output
    """

    sww_quantity = ['stage', 'elevation', 'xmomentum', 'ymomentum']

    try:
        fid = open(swwfile)
    except Exception, e:
        msg = 'File "%s" could not be opened: Error="%s"'\
              %(swwfile, e)
        raise msg

    print 'swwfile', swwfile

    # interpolation points are for y = 0
    m = 100 # number of increments in x
    max_x = 100000.
    x = 0
    points = []
    x_points = []
    for i in range(m):
        x += max_x/m
        points.append([x,0])
        x_points.append(x)
            
    time_thinning = 1
    use_cache = True
    verbose = True

    from anuga.abstract_2d_finite_volumes.util import file_function
    
    f = file_function(swwfile,
                      quantities = sww_quantity,
                      interpolation_points = points,
                      time_thinning = time_thinning,
                      verbose = verbose,
                      use_cache = use_cache)

    n = len(f.get_time()) # number of time steps

    from Numeric import zeros, Float
    from math import sqrt
    stage = zeros((n,m), Float)
    elevation = zeros((n,m), Float)
    xmom = zeros((n,m), Float)
    ymom = zeros((n,m), Float)
    momenta = zeros((n,m), Float)
    speed = zeros((n,m), Float)
    max_stages = []
    max_stage = None

    for k, location in enumerate(x_points):
        for i, t in enumerate(f.get_time()):
            w = f(t, point_id = k)[0]
            z = f(t, point_id = k)[1]
            uh = f(t, point_id = k)[2]
            vh = f(t, point_id = k)[3]
            depth = w-z      
            m = sqrt(uh*uh + vh*vh)
            if depth < 0.001:
                vel = 0.0
            else:
                vel = m / (depth + 1.e-6/depth)

            stage[i,k] = w
            elevation[i,k] = z
            xmom[i,k] = uh
            ymom[i,k] = vh
            momenta[i,k] = m
            speed[i,k] = vel

            if w > max_stage: max_stage = w
            max_stages.append(max_stage)

    from pylab import plot, xlabel, ylabel, savefig, close, hold
    hold(False)
    for i in range(n):
        plot(x_points,stage[i,:])
        xlabel('x')
        ylabel('stage (m)')
        #title('')
        name = 'time_%i' %i
        savefig(name)

    close('all')
    
animatesww2d('myexample2.sww')