Index: /issm/trunk-jpl/src/m/plot/plot_overlay.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_overlay.py (revision 19431) +++ /issm/trunk-jpl/src/m/plot/plot_overlay.py (revision 19432) @@ -36,6 +36,4 @@ 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'): Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 19431) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 19432) @@ -1,82 +1,92 @@ from cmaptools import truncate_colormap try: - import pylab as p - import matplotlib as mpl - import matplotlib.pyplot as plt + import pylab as p + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as npy except ImportError: - print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled" + print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled" def plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options,ax): - """ - PLOT_UNIT - unit plot, display data - - Usage: - plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options) - - See also: PLOTMODEL, PLOT_MANAGER - """ - - #edgecolor - edgecolor=options.getfieldvalue('edgecolor','None') - - #number of colorlevels for plots - colorlevels=options.getfieldvalue('colorlevels',128) - - alpha=options.getfieldvalue('alpha',1) - - #colormap - # default sequential colormap - defaultmap=truncate_colormap(mpl.cm.gnuplot2,0.1,0.9,128) - cmap=options.getfieldvalue('colormap',defaultmap) - if options.exist('cmap_set_over'): - over=options.getfieldvalue('cmap_set_over','0.5') - cmap.set_over(over) - if options.exist('cmap_set_under'): - under=options.getfieldvalue('cmap_set_under','0.5') - cmap.set_under(under) - - #normalize colormap if clim/caxis specified - if options.exist('clim'): - lims=options.getfieldvalue('clim',[min(data),max(data)]) - elif options.exist('caxis'): - lims=options.getfieldvalue('caxis',[min(data),max(data)]) - else: - if min(data)==max(data): - lims=[min(data)-0.5,max(data)+0.5] - else: - lims=[min(data),max(data)] - norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1]) - if datatype==1: - #element plot - if is2d: - tri=ax.tripcolor(x,y,elements,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,edgecolors=edgecolor) - else: - raise ValueError('plot_unit error: 3D element plot not supported yet') - return - - elif datatype==2: - #node plot - if is2d: - tri=ax.tricontourf(x,y,elements,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,extend='both') - if edgecolor != 'None': - ax.triplot(x,y,elements,color=edgecolor) - else: - raise ValueError('plot_unit error: 3D node plot not supported yet') - return - - elif datatype==3: - print 'plot_unit message: quiver plot not implemented yet' - return - - elif datatype==4: - #P1 patch plot - print 'plot_unit message: P1 patch plot not implemented yet' - return - - elif datatype==5: - print 'plot_unit message: P0 patch plot not implemented yet' - return - - else: - raise ValueError('datatype=%d not supported' % datatype) - + """ + PLOT_UNIT - unit plot, display data + + Usage: + plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options) + + See also: PLOTMODEL, PLOT_MANAGER + """ + + #edgecolor + edgecolor=options.getfieldvalue('edgecolor','None') + + #number of colorlevels for plots + colorlevels=options.getfieldvalue('colorlevels',128) + + alpha=options.getfieldvalue('alpha',1) + + #colormap + # default sequential colormap + defaultmap=truncate_colormap(mpl.cm.gnuplot2,0.1,0.9,128) + cmap=options.getfieldvalue('colormap',defaultmap) + if options.exist('cmap_set_over'): + over=options.getfieldvalue('cmap_set_over','0.5') + cmap.set_over(over) + if options.exist('cmap_set_under'): + under=options.getfieldvalue('cmap_set_under','0.5') + cmap.set_under(under) + + #normalize colormap if clim/caxis specified + if options.exist('clim'): + lims=options.getfieldvalue('clim',[npy.amin(data),npy.amax(data)]) + elif options.exist('caxis'): + lims=options.getfieldvalue('caxis',[npy.amin(data),npy.amax(data)]) + else: + if npy.amin(data)==npy.amax(data): + lims=[npy.amin(data)-0.5,npy.amax(data)+0.5] + else: + lims=[npy.amin(data),npy.amax(data)] + norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1]) + if datatype==1: + #element plot + if is2d: + tri=ax.tripcolor(x,y,elements,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,edgecolors=edgecolor) + else: + raise ValueError('plot_unit error: 3D element plot not supported yet') + return + + elif datatype==2: + #node plot + if is2d: + tri=ax.tricontourf(x,y,elements,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,extend='both') + if edgecolor != 'None': + ax.triplot(x,y,elements,color=edgecolor) + else: + raise ValueError('plot_unit error: 3D node plot not supported yet') + return + + elif datatype==3: + vx=data[:,0] + vy=data[:,1] + #TODO write plot_quiver.py to handle this here + color=npy.sqrt(vx**2+vy**2) + scale=options.getfieldvalue('scale',1000) + width=options.getfieldvalue('width',0.005*(npy.amax(x)-npy.amin(y))) + headwidth=options.getfieldvalue('headwidth',3) + headlength=options.getfieldvalue('headlength',5) + Q=ax.quiver(x,y,vx,vy,color,cmap=cmap,norm=norm,scale=scale, + width=width,headwidth=headwidth,headlength=headlength) + return + + elif datatype==4: + #P1 patch plot + print 'plot_unit message: P1 patch plot not implemented yet' + return + + elif datatype==5: + print 'plot_unit message: P0 patch plot not implemented yet' + return + + else: + raise ValueError('datatype=%d not supported' % datatype) + Index: /issm/trunk-jpl/src/m/plot/processdata.py =================================================================== --- /issm/trunk-jpl/src/m/plot/processdata.py (revision 19431) +++ /issm/trunk-jpl/src/m/plot/processdata.py (revision 19432) @@ -3,97 +3,106 @@ def processdata(md,data,options): - """ - PROCESSDATA - process data to be plotted - - datatype = 1 -> elements - datatype = 2 -> nodes - datatype = 3 -> node quivers - datatype = 4 -> patch - - Usage: - data,datatype=processdata(md,data,options); - - See also: PLOTMODEL, PROCESSMESH - """ - - #check format - if (len(data)==0 or (len(data)==1 and not isinstance(data,dict) and isnan(data).all())): - raise ValueError("processdata error message: 'data' provided is empty") - - #needed later on - if 'numberofvertices2d' in dir(md.mesh): - numberofvertices2d=md.mesh.numberofvertices2d - numberofelements2d=md.mesh.numberofelements2d - else: - numberofvertices2d=npy.nan - numberofelements2d=npy.nan - - procdata=npy.copy(data) - - #process patch - - #initialize datatype - datatype=0 - - #get datasize - if npy.ndim(procdata)==1: - datasize=npy.array([len(procdata),1]) - else: - datasize=npy.shape(procdata) - - #process NaN's if any - nanfill=options.getfieldvalue('nan',-9999) - if npy.any(npy.isnan(procdata)): - lb=npy.min(data[~npy.isnan(data)]) - ub=npy.max(data[~npy.isnan(data)]) - if lb==ub: - lb=lb-0.5 - ub=ub+0.5 - nanfill=lb-1 - procdata[npy.isnan(procdata)]=nanfill - options.addfielddefault('clim',[lb,ub]) - options.addfielddefault('cmap_set_under','1') - print "WARNING: nan's treated as", nanfill, "by default. Change using pairoption 'nan',nan_fill_value in plotmodel call" - - #non-patch processing - - #element data - if datasize[0]==md.mesh.numberofelements and datasize[1]==1: - - #initialize datatype if non patch - if datatype!=4 and datatype!=5: - datatype=1 - - #mask? - - #log? - - #node data - if datasize[0]==md.mesh.numberofvertices and datasize[1]==1: - datatype=2 - - #spc time series? - if datasize[0]==md.mesh.numberofvertices+1: - datatype=2 - spccol=options.getfieldvalue('spccol',0) - print 'multiple-column spc field; specify column to plot using option "spccol"' - print 'column ', spccol, ' plotted for time: ', procdata[-1,spccol] - procdata=procdata[0:-1,spccol] - - #mask? - - #log? - - #layer projection? - - #control arrow density if quiver plot - - #convert rank-2 array to rank-1 - if npy.ndim(procdata)==2 and npy.shape(procdata)[1]==1: - procdata=procdata.reshape(-1,) - - #if datatype is still zero, error out - if datatype==0: - raise ValueError("processdata error: data provided not recognized or not supported") - else: - return procdata, datatype + """ + PROCESSDATA - process data to be plotted + + datatype = 1 -> elements + datatype = 2 -> nodes + datatype = 3 -> node quivers + datatype = 4 -> patch + + Usage: + data,datatype=processdata(md,data,options); + + See also: PLOTMODEL, PROCESSMESH + """ + + #check format + if (len(data)==0 or (len(data)==1 and not isinstance(data,dict) and isnan(data).all())): + raise ValueError("processdata error message: 'data' provided is empty") + + #needed later on + if 'numberofvertices2d' in dir(md.mesh): + numberofvertices2d=md.mesh.numberofvertices2d + numberofelements2d=md.mesh.numberofelements2d + else: + numberofvertices2d=npy.nan + numberofelements2d=npy.nan + + procdata=npy.copy(data) + + #process patch + + #initialize datatype + datatype=0 + + #get datasize + if npy.ndim(procdata)==1: + datasize=npy.array([len(procdata),1]) + else: + datasize=npy.shape(procdata) + if len(datasize)>2: + raise ValueError('data passed to plotmodel has more than 2 dimensions; check that column vectors are rank-1') + + #process NaN's if any + nanfill=options.getfieldvalue('nan',-9999) + if npy.any(npy.isnan(procdata)): + lb=npy.min(data[~npy.isnan(data)]) + ub=npy.max(data[~npy.isnan(data)]) + if lb==ub: + lb=lb-0.5 + ub=ub+0.5 + nanfill=lb-1 + procdata[npy.isnan(procdata)]=nanfill + options.addfielddefault('clim',[lb,ub]) + options.addfielddefault('cmap_set_under','1') + print "WARNING: nan's treated as", nanfill, "by default. Change using pairoption 'nan',nan_fill_value in plotmodel call" + + #quiver plot + if datasize[1]>1: + if datasize[0]==md.mesh.numberofvertices and datasize[1]==2: + datatype=3 + else: + raise ValueError('plotmodel error message: data should have two columns of length md.mesh.numberofvertices for a quiver plot') + + #non-patch processing + + #element data + if datasize[0]==md.mesh.numberofelements and datasize[1]==1: + + #initialize datatype if non patch + if datatype!=4 and datatype!=5: + datatype=1 + + #mask? + + #log? + + #node data + if datasize[0]==md.mesh.numberofvertices and datasize[1]==1: + datatype=2 + + #spc time series? + if datasize[0]==md.mesh.numberofvertices+1: + datatype=2 + spccol=options.getfieldvalue('spccol',0) + print 'multiple-column spc field; specify column to plot using option "spccol"' + print 'column ', spccol, ' plotted for time: ', procdata[-1,spccol] + procdata=procdata[0:-1,spccol] + + #mask? + + #log? + + #layer projection? + + #control arrow density if quiver plot + + #convert rank-2 array to rank-1 + if npy.ndim(procdata)==2 and npy.shape(procdata)[1]==1: + procdata=procdata.reshape(-1,) + + #if datatype is still zero, error out + if datatype==0: + raise ValueError("processdata error: data provided not recognized or not supported") + else: + return procdata, datatype