Changeset 4700

Timestamp:

Sep 4, 2007, 5:32:22 PM (18 years ago)

Author:

sexton

Message:

working on converting matlab script to matplotlib

File:

• anuga_validation/convergence_study/animatesww2d_alt.py (modified) (2 diffs)

anuga_validation/convergence_study/animatesww2d_alt.py

@@ -17 +17 @@
     ymomentum = fid.variables['ymomentum']
     
-    
+
+    # find indices for where y = 0
+
+    # 
     # Set up plot
 
@@ -64 +67 @@
 #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')