Index: /issm/trunk-jpl/src/m/plot/plot_manager.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 17653)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 17654)
@@ -1,3 +1,2 @@
-
 try:
 	import pylab as p
@@ -11,4 +10,5 @@
 from plot_unit import plot_unit
 from applyoptions import applyoptions
+from plot_overlay import plot_overlay
 
 def plot_manager(md,options,subplotwidth,nlines,ncols,i):
@@ -24,8 +24,21 @@
 	#parse options and get a structure of options
 	options=checkplotoptions(md,options)
-	#print 'options:', options
 
 	#get data to be plotted
 	data=options.getfieldvalue('data');
+
+	#initialize plot handle variable
+	handle=None
+
+	##basemap plot
+	#if options.exist('basemap'):
+	#	plot_basemap(md,data,options,nlines,ncols,i)
+
+	#overlay plot
+	if options.exist('overlay'):
+		handle=plot_overlay(md,data,options,nlines,ncols,i)
+		options.addfielddefault('alpha',0.5)
+		options.addfielddefault('xlim',[min(md.mesh.x),max(md.mesh.x)])
+		options.addfielddefault('ylim',[min(md.mesh.y),max(md.mesh.y)])
 
 	#figure out if this is a special plot
@@ -36,4 +49,7 @@
 			plot_mesh(md,options,nlines,ncols,i)
 			return
+		elif data.lower()=='none':
+			print 'no data provided to plot (TODO: write plot_none.py)'
+			return
 		else:
 			print "WARNING: '%s' is not implemented or is not a valid string for option 'data'" % data
@@ -42,6 +58,4 @@
 	#else:
 		#print "'data' not a string, plotting model properties yet to be implemented..."
-
-	#Overlay plot
 
 	#Gridded plot
@@ -56,5 +70,6 @@
 
 	#standard plot
-	p.subplot(nlines,ncols,i,aspect='equal')
+	if not handle:
+		p.subplot(nlines,ncols,i,aspect='equal')
 
 	#plot unit
Index: /issm/trunk-jpl/src/m/plot/plot_overlay.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_overlay.py	(revision 17654)
+++ /issm/trunk-jpl/src/m/plot/plot_overlay.py	(revision 17654)
@@ -0,0 +1,85 @@
+import numpy as npy
+from processmesh import processmesh
+from processdata import processdata
+from osgeo import gdal
+import matplotlib.pyplot as plt
+import matplotlib as mpl
+import os
+
+def plot_overlay(md,data,options,rows,cols,i):
+	'''
+	Function for plotting a georeferenced image.  This function is called
+	from within the plotmodel code.
+	'''
+
+	x,y,z,elements,is2d,isplanet=processmesh(md,[],options)
+
+	if data=='none' or data==None:
+		imageonly=1
+		data=npy.float('nan')*npy.ones((md.mesh.numberofvertices,))
+		datatype=1
+	else:
+		imageonly=0
+		data,datatype=processdata(md,data,options)
+
+	if not is2d:
+		raise StandardError('overlay plot not supported for 3D meshes, project on a 2D layer first')
+	if datatype==3:
+		raise StandardError('overlay not yet supported for quiver plots')	
+
+	if not options.exist('geotiff_name'):
+		raise StandardError('overlay error: provide geotiff_name with path to geotiff file')
+	geotiff=options.getfieldvalue('geotiff_name')
+
+	xlim=options.getfieldvalue('xlim',[min(md.mesh.x),max(md.mesh.x)])
+	ylim=options.getfieldvalue('ylim',[min(md.mesh.y),max(md.mesh.y)])
+
+	gtif=gdal.Open(geotiff)
+	trans=gtif.GetGeoTransform()
+	xmin=trans[0]
+	xmax=trans[0]+gtif.RasterXSize*trans[1]
+	ymin=trans[3]+gtif.RasterYSize*trans[5]
+	ymax=trans[3]
+	
+	# allow supplied geotiff to have limits smaller than basemap or model limits
+	x0=max(min(xlim),xmin)
+	x1=min(max(xlim),xmax)
+	y0=max(min(ylim),ymin)
+	y1=min(max(ylim),ymax)
+	inputname='temp.tif'
+	os.system('gdal_translate -quiet -projwin ' + str(x0) + ' ' + str(y1) + ' ' + str(x1) + ' ' + str(y0) + ' ' + geotiff + ' ' + inputname)
+	
+	gtif=gdal.Open(inputname)
+	arr=gtif.ReadAsArray()
+	#os.system('rm -rf ./temp.tif')
+	
+	# get parameters from cropped geotiff
+	trans=gtif.GetGeoTransform()
+	xmin=trans[0]
+	xmax=trans[0]+gtif.RasterXSize*trans[1]
+	ymin=trans[3]+gtif.RasterYSize*trans[5]
+	ymax=trans[3]
+	dx=trans[1]
+	dy=trans[5]	
+	
+	xarr=npy.arange(xmin,xmax,dx)
+	yarr=npy.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is)
+	xg,yg=npy.meshgrid(xarr,yarr)
+	if options.exist('basemap'):
+		# get handle to basemap instance
+		# handle= functiontogethandle()
+		# create coordinate grid in map projection units (for plotting)
+		lats,lons=xy2ll(xg,yg,-1,0,71)
+		xgmap,ygmap=m(lons,lats) # map projection units returned by basemap instance
+	else:
+		xgmap=xg
+		ygmap=yg
+		handle=plt.gca()
+	
+	# for single band geotiffs
+	overlaylims=options.getfieldvalue('overlaylims',[min(arr.ravel()),max(arr.ravel())])
+
+	norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1])
+	handle.pcolormesh(xgmap, ygmap, npy.flipud(arr), cmap=mpl.cm.Greys, norm=norm)
+	plt.axes().set_aspect('equal','box')
+	return handle
Index: /issm/trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 17653)
+++ /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 17654)
@@ -21,5 +21,7 @@
    #number of colorlevels for plots
    colorlevels=options.getfieldvalue('colorlevels',256)
-   
+
+   alpha=options.getfieldvalue('alpha',1)
+
    #colormap
    cmap=options.getfieldvalue('colormap',mpl.cm.gnuplot2)
@@ -50,5 +52,5 @@
    	#node plot
    	if is2d:
-   		p.tricontourf(x,y,elements,data,colorlevels,cmap=cmap,norm=norm)
+   		p.tricontourf(x,y,elements,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha)
    		if edgecolor != 'None':
    			p.triplot(x,y,elements,color=edgecolor)