<p><b>mhoffman@lanl.gov</b> 2012-06-15 15:55:13 -0600 (Fri, 15 Jun 2012)</p><p>BRANCH COMMIT -- land ice<br>
Created python script to animate fields in an MPAS output file.<br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice_projects/grid_tools/animate_mpas_field.py
===================================================================
--- branches/land_ice_projects/grid_tools/animate_mpas_field.py (rev 0)
+++ branches/land_ice_projects/grid_tools/animate_mpas_field.py 2012-06-15 21:55:13 UTC (rev 1995)
@@ -0,0 +1,166 @@
+#!/usr/bin/python
+# Script capable of plotting most MPAS fields. Invoke with 'python plot_mpas_field.py --help for details about how to use.
+# Matt Hoffman, June 14, 2012
+
+import numpy, time
+from optparse import OptionParser
+import matplotlib.pyplot as plt
+try:
+ from Scientific.IO.NetCDF import NetCDFFile
+ netCDF_module = 'Scientific.IO.NetCDF'
+except ImportError:
+ try:
+ from netCDF4 import Dataset as NetCDFFile
+ netCDF_module = 'netCDF4'
+ except ImportError:
+ print 'Unable to import any of the following python modules:'
+ print ' Scientific.IO.NetCDF </font>
<font color="blue"> netcdf4 '
+ print 'One of them must be installed.'
+ raise ImportError('No netCDF module found')
+
+
+print "** Gathering information. (Invoke with 'python plot_mpas_field.py --help' for more details. All arguments are optional)"
+parser = OptionParser()
+parser.add_option("-f", "--file", dest="filename", help="file to visualize; default: output.nc", metavar="FILE")
+parser.add_option("-v", "--var", dest="variable", help="variable to visualize; use 'uMag' for velocity magnitude; default: thickness", metavar="VAR")
+parser.add_option("-l", "--level", dest="vertlevel", help="vertical level to visualize (0 based); default: 0", metavar="LEVEL")
+parser.add_option("-s", "--save", action="store_true", dest="saveimages", help="include this flag to save animation as series of .png files")
+parser.add_option("-n", "--nodisp", action="store_true", dest="hidefigs", help="include this flag to not display plot after animation is complete")
+parser.add_option("-x", "--max", dest="maximum", help="maximum for color bar", metavar="MAX")
+parser.add_option("-m", "--min", dest="minimum", help="minimum for color bar", metavar="MIN")
+parser.add_option("-g", "--log", action="store_true", dest="log", help="include this flag to plot on log scale", metavar="LOG")
+parser.add_option("-k", "--mask", action="store_true", dest="mask", help="include this flag to mask non-ice cells with gray", metavar="MASK")
+
+options, args = parser.parse_args()
+
+if not options.filename:
+ print "No filename provided. Using output.nc."
+ options.filename = "output.nc"
+
+if not options.vertlevel:
+ print "No vertical level provided. Using level 0 (surface), if level is needed."
+ vert_level = 0
+else:
+ vert_level = int(options.vertlevel)
+ print "Using vertical level " + options.vertlevel
+
+if not options.variable:
+ print "No variable provided. Using 'thickness'."
+ varname = 'thickness'
+else:
+ varname = options.variable
+ print "Using variable " + varname
+
+
+f = NetCDFFile(options.filename,'r')
+
+# Get grid stuff
+#times = f.variables['xtime'] # Not needed unless trying to print actual time stamp
+xCell = f.variables['xCell']
+yCell = f.variables['yCell']
+xEdge = f.variables['xEdge']
+yEdge = f.variables['yEdge']
+angleEdge = f.variables['angleEdge']
+
+# get the requested variable
+if varname == 'uMag':
+ uReconstructX = f.variables['uReconstructX']
+ uReconstructY = f.variables['uReconstructY']
+ var = numpy.zeros(uReconstructX.shape)
+ var[:] = (uReconstructX[:] ** 2 + uReconstructY[:] **2)**0.5
+ dims = uReconstructX.dimensions
+else:
+ var = f.variables[varname]
+ dims = var.dimensions
+
+# Determine what the appropriate X & Y values are:
+if 'nCells' in dims:
+ x = xCell[:]
+ y = yCell[:]
+elif 'nEdges' in dims:
+ x = xEdge[:]
+ y = yEdge[:]
+elif 'nVertices' in dims:
+ x = xVertex[:]
+ y = yVertex[:]
+
+if options.mask:
+ thickness = f.variables['thickness']
+
+# Get the needed slice and determine the plot title. Make some assumptions about how dimensions are arranged:
+plottitle = varname
+if 'Time' in dims:
+ time_length = var.shape[0] # Assume time is the first dimension
+ if 'nVertLevels' in dims:
+ plottitle = plottitle + ', for layer ' + str(vert_level)
+ var_slice = var[:,:,vert_level]
+ else:
+ var_slice = var[:,:]
+else:
+ print "Time is not a dimension of this variable. Unable to animate it!"
+ sys.exit()
+
+if options.log:
+ var_slice = numpy.log10(var_slice + 1.0e-5)
+ plottitle = plottitle + ', log scale'
+
+# Determine color axis max & min values to use.
+if not options.maximum:
+ maxval = var_slice[:].max()
+else:
+ maxval = float(options.maximum)
+
+if not options.minimum:
+ minval = var_slice[:].min()
+else:
+ minval = float(options.minimum)
+
+
+# MAKE THE PLOT
+print '** Beginning to create plot.'
+plt.ion()
+fig = plt.figure(1, facecolor='w', figsize=(8, 6), dpi=200)
+ax = fig.add_subplot(111, aspect='equal')
+
+markersize = 70
+markershape = 'h'
+
+plt.scatter(xCell[:], yCell[:], markersize, var_slice[0,:], marker=markershape, edgecolors='none', vmin=minval, vmax=maxval)
+plt.colorbar()
+if options.mask:
+ maskindices = numpy.nonzero(thickness[0,:] <= 0.0)[:]
+ plt.scatter(numpy.take(xCell,maskindices), numpy.take(yCell,maskindices), markersize+2, '0.9', marker=markershape, edgecolors='none')
+plt.title( plottitle + ' at time 0' )
+plt.draw()
+if options.saveimages:
+ plotname = varname + '.000.' + options.filename + '.png'
+ plt.savefig(plotname)
+ print 'Saved plot as ' + plotname
+
+
+for t in range(1, time_length):
+ ax.cla()
+ plt.scatter(xCell[:], yCell[:], markersize, var_slice[t,:], marker=markershape, edgecolors='none', vmin=minval, vmax=maxval)
+ if options.mask:
+ maskindices = numpy.nonzero(thickness[t,:] <= 0.0)[:]
+ plt.scatter(numpy.take(xCell,maskindices), numpy.take(yCell,maskindices), markersize+2, '0.9', marker=markershape, edgecolors='none')
+ plt.title( plottitle + ' at time ' + str(t) )
+ plt.draw()
+ time.sleep(0.05)
+ if options.saveimages:
+ plotname = varname + '.' + '{0:04d}'.format(t) + '.' + options.filename + '.png'
+ plt.savefig(plotname)
+ print 'Saved plot as ' + plotname
+
+plt.ioff()
+print "Plotted " + varname + "."
+
+
+if options.hidefigs:
+ print "Plot display disabled with -n argument."
+else:
+ print 'Showing plot... Close plot window to exit.'
+ plt.show()
+
+f.close()
+
Property changes on: branches/land_ice_projects/grid_tools/animate_mpas_field.py
___________________________________________________________________
Added: svn:executable
+ *
</font>
</pre>