Wednesday, November 14, 2007

Sample movie for GETM

Here is a sample script for generating a series of figures (that can be turned into a movie with ffmpeg, see [below](http://pong.tamu.edu/~rob/?p=100)). The case used is the box_curvilinear test case.



'''
Make an animation of surface temperature and velocity from the
box_curvilinear GETM test case
'''

# Copyright R. Hetland on 2007-11-15.
# All rights reserved.
# Version 1.0.0 on 2007-11-15

from numpy import *
import pylab as pl
import netCDF4

nc3d = netCDF4.Dataset('box_curvilinear.3d.nc')
nctopo = netCDF4.Dataset('topo.nc')

# load in cell verticies
xx = nctopo.variables['xx'][:]
yx = nctopo.variables['yx'][:]

# load in cell centers
xc = nc3d.variables['xc'][:]
yc = nc3d.variables['yc'][:]

# load in model time
t = nc3d.variables['time'][:]

kidx = -5 # vertical slice index (-1 = uppermost slice)
Nvec = 3 # decimation of velocity field for vector plot

# Set figure font properties
mono_font = pl.matplotlib.font_manager.FontProperties()
mono_font.set_family('monospace')
mono_font.set_size(16)

title_font = pl.matplotlib.font_manager.FontProperties()
title_font.set_style('oblique')
title_font.set_size(12)

# Create figure and axes
fig = pl.figure(figsize=(15, 5))
pl.ioff()

for tidx in range(len(t)):
temp = nc3d.variables['temp'][tidx, kidx, :, :]
temp = ma.masked_where(temp==-9999, temp)

uu = nc3d.variables['uu'][tidx, kidx, :, :]
uu[uu==-9999.0] = nan
vv = nc3d.variables['vv'][tidx, kidx, :, :]
vv[vv==-9999.0] = nan

fig.clf()
ax = fig.add_axes([0, 0, 1, 1])
cax = fig.add_axes([0.4, 0.2, 0.2, 0.03])
ax.set_axis_off()
pcm = ax.pcolor(xx, yx, temp, vmin=0, vmax=20)
q = ax.quiver(xc[::Nvec, ::Nvec], yc[::Nvec, ::Nvec],
uu[::Nvec, ::Nvec], vv[::Nvec, ::Nvec],
units='width', scale=10.0, width=0.002,
pivot='middle')

pl.quiverkey(q, 0.5, 0.3, 0.3, r'0.3 m s$^{-1}$')

cb = pl.colorbar(pcm, cax=cax, orientation='horizontal')
cb.set_label(r'Temperature [$^\circ$C]')

ax.text(0.1, 0.9, '%5.3f days' % (t[tidx]/86400.0),
ha='left', va='bottom',
transform = ax.transAxes,
fontproperties=mono_font)

ax.text(0.95, 0.95, r'box_curvilinear GETM test case',
ha='right', va='bottom',
transform = ax.transAxes,
fontproperties=title_font)
ax.set_aspect(1.0)

pl.savefig('frame_%03d.png' % tidx, dpi=80)

No comments:

Post a Comment