<p><b>mhoffman@lanl.gov</b> 2012-06-14 23:13:30 -0600 (Thu, 14 Jun 2012)</p><p>BRANCH COMMIT -- land ice<br>
A python script that can be used to plot most MPAS fields on planar meshes.<br>
Useful for quickly visualizing output on remote machines (provided they have matplotlib, numpy, and either of Scientific or netCDF4 modules installed).<br>
Could be used for output or input from other cores besides land ice with little to no modification.<br>
</p><hr noshade><pre><font color="gray">Added: branches/land_ice_projects/grid_tools/plot_mpas_field.py
===================================================================
--- branches/land_ice_projects/grid_tools/plot_mpas_field.py (rev 0)
+++ branches/land_ice_projects/grid_tools/plot_mpas_field.py 2012-06-15 05:13:30 UTC (rev 1983)
@@ -0,0 +1,147 @@
+#!/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
+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')
+
+
+parser = OptionParser()
+parser.add_option("-f", "--file", dest="filename", help="file to visualize", metavar="FILE")
+parser.add_option("-v", "--var", dest="variable", help="variable to visualize; use 'uMag' for velocity magnitude", metavar="VAR")
+parser.add_option("-t", "--time", dest="time", help="time step to visualize (0 based)", metavar="TIME")
+parser.add_option("-l", "--level", dest="vertlevel", help="vertical level to visualize (0 based)", metavar="LEVEL")
+parser.add_option("-s", "--save", action="store_true", dest="saveimages", help="include this flag to save plots as files")
+parser.add_option("-n", "--nodisp", action="store_true", dest="hidefigs", help="include this flag to not display plots (usually used with -s)")
+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")
+
+options, args = parser.parse_args()
+
+print "** Gathering information. (Invoke with 'python plot_mpas_field.py --help' for more details.)"
+if not options.filename:
+ print "No filename provided. Using output.nc."
+ options.filename = "output.nc"
+
+if not options.time:
+ print "No time provided. Using time 0, if time is needed."
+ time_slice = 0
+else:
+ time_slice = int(options.time)
+ print "Using time level " + options.time
+
+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']
+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[:]
+
+# Get the needed slice and determine the plot title. Make some assumptions about how dimensions are arranged:
+plottitle = varname
+if 'Time' in dims:
+ plottitle = plottitle + ', at time ' + str(time_slice)
+ if 'nVertLevels' in dims:
+ plottitle = plottitle + ', for layer ' + str(vert_level)
+ var_slice = var[time_slice,:,vert_level]
+ else:
+ var_slice = var[time_slice,:]
+else:
+ if 'nVertLevels' in dims:
+ plottitle = plottitle + ', for layer ' + str(vert_level)
+ var_slice = var[:,vert_level]
+ else:
+ var_slice = var[:]
+
+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.'
+fig = plt.figure(1, facecolor='w')
+ax = fig.add_subplot(111, aspect='equal')
+plt.scatter(xCell[:], yCell[:], 3, var_slice, marker='o', edgecolors='none', vmin=minval, vmax=maxval)
+plt.colorbar()
+plt.title( plottitle )
+plt.draw()
+print "Plotted " + varname + "."
+
+if options.saveimages:
+ plotname = varname + '.' + options.filename + '.png'
+ plt.savefig(plotname)
+ print 'Saved plot as ' + plotname
+
+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/plot_mpas_field.py
___________________________________________________________________
Added: svn:executable
+ *
</font>
</pre>