Index: /issm/trunk-jpl/src/m/plot/applyoptions.py =================================================================== --- /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 21253) @@ -32,16 +32,13 @@ #ax=p.gca() - #font {{{ + # {{{ font fontsize=options.getfieldvalue('fontsize',8) fontweight=options.getfieldvalue('fontweight','normal') fontfamily=options.getfieldvalue('fontfamily','sans-serif') - font={ - 'fontsize' :fontsize, - 'fontweight' :fontweight, - 'family' :fontfamily - } - #}}} - - #title {{{ + font={'fontsize' :fontsize, + 'fontweight' :fontweight, + 'family' :fontfamily} + # }}} + # {{{ title if options.exist('title'): title=options.getfieldvalue('title') @@ -59,7 +56,6 @@ titlefont['weight']=titlefontweight ax.set_title(title,**titlefont) - #}}} - - #xlabel, ylabel, zlabel {{{ + # }}} + # {{{ xlabel, ylabel, zlabel if options.exist('labelfontsize'): labelfontsize=options.getfieldvalue('labelfontsize') @@ -82,7 +78,6 @@ if options.exist('zlabel'): ax.set_zlabel(options.getfieldvalue('zlabel'),**labelfont) - #}}} - - #xticks, yticks, zticks (tick locations) {{{ + # }}} + # {{{ xticks, yticks, zticks (tick locations) if options.exist('xticks'): if options.exist('xticklabels'): @@ -103,7 +98,6 @@ else: ax.set_zticks(options.getfieldvalue('zticks')) - #}}} - - #xticklabels,yticklabels,zticklabels {{{ + # }}} + # {{{ xticklabels,yticklabels,zticklabels if options.getfieldvalue('ticklabels','off')=='off' or options.getfieldvalue('ticklabels',0)==0: options.addfielddefault('xticklabels',[]) @@ -119,11 +113,9 @@ zticklabels=options.getfieldvalue('zticklabels') ax.set_zticklabels(zticklabels) - #}}} - - #ticklabel notation {{{ + # }}} + # {{{ ticklabel notation #ax.ticklabel_format(style='sci',scilimits=(0,0)) - #}}} - - #ticklabelfontsize {{{ + # }}} + # {{{ ticklabelfontsize if options.exist('ticklabelfontsize'): for label in ax.get_xticklabels() + ax.get_yticklabels(): @@ -132,13 +124,12 @@ for label in ax.get_zticklabels(): label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) - #}}} - - #view + # }}} + # {{{ view TOFIX #if int(md.mesh.dimension) == 3 and options.exist('layer'): # #options.getfieldvalue('view') ? # ax=fig.gca(projection='3d') #plt.show() - - #axis {{{ + # }}} + # {{{ axis if options.exist('axis'): if options.getfieldvalue('axis',True)=='off': @@ -147,10 +138,9 @@ p.setp(ax.get_yticklabels(), visible=False) # }}} - - #box + # {{{ box if options.exist('box'): eval(options.getfieldvalue('box')) - - #xlim, ylim, zlim {{{ + # }}} + # {{{ xlim, ylim, zlim if options.exist('xlim'): ax.set_xlim(options.getfieldvalue('xlim')) @@ -159,13 +149,12 @@ if options.exist('zlim'): ax.set_zlim(options.getfieldvalue('zlim')) - #}}} - - #latlon - - #Basinzoom - - #ShowBasins - - #clim {{{ + # }}} + # {{{ latlon TODO + # }}} + # {{{ Basinzoom TODO + # }}} + # {{{ ShowBasins TODO + # }}} + # {{{ clim if options.exist('clim'): lims=options.getfieldvalue('clim') @@ -178,21 +167,25 @@ if len(data)>0: lims=[data.min(),data.max()] else: lims=[0,1] - #}}} - - #shading + # }}} + # {{{ shading TODO #if options.exist('shading'): - - #grid {{{ + # }}} + # {{{ grid if options.exist('grid'): if 'on' in options.getfieldvalue('grid','on'): ax.grid() - #}}} - - #colormap {{{ + # }}} + # {{{ colormap # default sequential colormap defaultmap=truncate_colormap(mpl.cm.gnuplot2,0.1,0.9,128) cmap=options.getfieldvalue('colormap',defaultmap) - norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1]) + if options.exist('log'): + norm = mpl.colors.LogNorm(vmin=lims[0], vmax=lims[1]) + else: + norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1]) options.addfield('colornorm',norm) + if options.exist('cmap_set_bad'): + NaNcolor=options.getfieldvalue('cmap_set_bad','w') + cmap.set_bad(NaNcolor,1.0) cbar_extend=0 if options.exist('cmap_set_over'): @@ -205,113 +198,109 @@ cbar_extend+=2 options.addfield('colormap',cmap) - #}}} - - #contours {{{ + # }}} + # {{{ contours if options.exist('contourlevels'): plot_contour(md,data,options,ax) - #}}} - - #wrapping - - #colorbar {{{ + # }}} + # {{{ wrapping TODO + # }}} + # {{{ colorbar if options.getfieldvalue('colorbar',1)==1: - if cbar_extend==0: - extend='neither' - elif cbar_extend==1: - extend='max' - elif cbar_extend==2: - extend='min' - elif cbar_extend==3: - extend='both' - cb = mpl.colorbar.ColorbarBase(ax.cax, cmap=cmap, norm=norm, extend=extend) - if options.exist('alpha'): - cb.set_alpha(options.getfieldvalue('alpha')) - if options.exist('colorbarnumticks'): - cb.locator=MaxNLocator(nbins=options.getfieldvalue('colorbarnumticks',5)) - else: - cb.locator=MaxNLocator(nbins=5) # default 5 ticks - if options.exist('colorbartickspacing'): - locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing')) - cb.set_ticks(locs) - if options.exist('colorbarlines'): - locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines')) - cb.add_lines(locs,['k' for i in range(len(locs))],npy.ones_like(locs)) - if options.exist('colorbarlineatvalue'): - locs=options.getfieldvalue('colorbarlineatvalue') - colors=options.getfieldvalue('colorbarlineatvaluecolor',['k' for i in range (len(locs))]) - widths=options.getfieldvalue('colorbarlineatvaluewidth',npy.ones_like(locs)) - cb.add_lines(locs,colors,widths) - if options.exist('colorbartitle'): - if options.exist('colorbartitlepad'): - cb.set_label(options.getfieldvalue('colorbartitle'),\ - labelpad=options.getfieldvalue('colorbartitlepad'),fontsize=fontsize) - else: - cb.set_label(options.getfieldvalue('colorbartitle'),fontsize=fontsize) - cb.ax.tick_params(labelsize=fontsize) - cb.solids.set_rasterized(True) - cb.update_ticks() - cb.set_alpha(1) - cb.draw_all() - plt.sca(ax) # return to original axes control - #}}} - - #expdisp {{{ - if options.exist('expdisp'): - filename=options.getfieldvalue('expdisp') - style=options.getfieldvalue('expstyle','k') - linewidth=options.getfieldvalue('explinewidth',1) - for i in xrange(len(filename)): - filenamei=filename[i] - stylei=style[i] - if type(linewidth)==list: - linewidthi=linewidth[i] - else: - linewidthi=linewidth - expdisp(filenamei,ax,linestyle=stylei,linewidth=linewidthi,unitmultiplier=options.getfieldvalue('unit',1)) - #}}} - - #area - - #text {{{ + if cbar_extend==0: + extend='neither' + elif cbar_extend==1: + extend='max' + elif cbar_extend==2: + extend='min' + elif cbar_extend==3: + extend='both' + cb = mpl.colorbar.ColorbarBase(ax.cax, cmap=cmap, norm=norm, extend=extend) + if options.exist('alpha'): + cb.set_alpha(options.getfieldvalue('alpha')) + if options.exist('colorbarnumticks'): + cb.locator=MaxNLocator(nbins=options.getfieldvalue('colorbarnumticks',5)) + else: + cb.locator=MaxNLocator(nbins=5) # default 5 ticks + if options.exist('colorbartickspacing'): + locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing')) + cb.set_ticks(locs) + if options.exist('colorbarlines'): + locs=npy.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines')) + cb.add_lines(locs,['k' for i in range(len(locs))],npy.ones_like(locs)) + if options.exist('colorbarlineatvalue'): + locs=options.getfieldvalue('colorbarlineatvalue') + colors=options.getfieldvalue('colorbarlineatvaluecolor',['k' for i in range (len(locs))]) + widths=options.getfieldvalue('colorbarlineatvaluewidth',npy.ones_like(locs)) + cb.add_lines(locs,colors,widths) + if options.exist('colorbartitle'): + if options.exist('colorbartitlepad'): + cb.set_label(options.getfieldvalue('colorbartitle'), + labelpad=options.getfieldvalue('colorbartitlepad'),fontsize=fontsize) + else: + cb.set_label(options.getfieldvalue('colorbartitle'),fontsize=fontsize) + cb.ax.tick_params(labelsize=fontsize) + cb.solids.set_rasterized(True) + cb.update_ticks() + cb.set_alpha(1) + cb.draw_all() + plt.sca(ax) # return to original axes control + # }}} + # {{{ expdisp + if options.exist('expdisp'): + filename=options.getfieldvalue('expdisp') + style=options.getfieldvalue('expstyle','k') + linewidth=options.getfieldvalue('explinewidth',1) + for i in xrange(len(filename)): + filenamei=filename[i] + stylei=style[i] + if type(linewidth)==list: + linewidthi=linewidth[i] + else: + linewidthi=linewidth + expdisp(filenamei,ax,linestyle=stylei,linewidth=linewidthi,unitmultiplier=options.getfieldvalue('unit',1)) + # }}} + # {{{ area TODO + # }}} + # {{{ text if options.exist('text'): - text=options.getfieldvalue('text') - textx=options.getfieldvalue('textx') - texty=options.getfieldvalue('texty') - textcolor=options.getfieldvalue('textcolor') - textweight=options.getfieldvalue('textweight') - textrotation=options.getfieldvalue('textrotation') - textfontsize=options.getfieldvalue('textfontsize') - for label,x,y,size,color,weight,rotation in zip(text,textx,texty,textfontsize,textcolor,textweight,textrotation): - ax.text(x,y,label,transform=ax.transAxes,fontsize=size,color=color,weight=weight,rotation=rotation) - #}}} - - #north arrow - - #scale ruler - - #streamlines - if options.exist('streamlines'): - plot_streamlines(md,options,ax) - - - #axis positions - - #figure position - - #axes position - - #showregion - - #flat edges of a partition - - #scatter - - #backgroundcolor - - #figurebackgroundcolor - - #lighting - - #point cloud - - #inset + text=options.getfieldvalue('text') + textx=options.getfieldvalue('textx') + texty=options.getfieldvalue('texty') + textcolor=options.getfieldvalue('textcolor') + textweight=options.getfieldvalue('textweight') + textrotation=options.getfieldvalue('textrotation') + textfontsize=options.getfieldvalue('textfontsize') + for label,x,y,size,color,weight,rotation in zip(text,textx,texty,textfontsize,textcolor,textweight,textrotation): + ax.text(x,y,label,transform=ax.transAxes,fontsize=size,color=color,weight=weight,rotation=rotation) + # }}} + # {{{ north arrow TODO + # }}} + # {{{ scale ruler TODO + # }}} + # {{{ streamlines TOFIX + if options.exist('streamlines'): + plot_streamlines(md,options,ax) + # }}} + # {{{ axis positions TODO + # }}} + # {{{ figure position TODO + # }}} + # {{{ axes position TODO + # }}} + # {{{ showregion TODO + # }}} + # {{{ flat edges of a partition TODO + # }}} + # {{{ scatter TODO + # }}} + # {{{ backgroundcolor TODO + # }}} + # {{{ figurebackgroundcolor TODO + # }}} + # {{{ lighting TODO + # }}} + # {{{ point cloud TODO + # }}} + # {{{ inset TODO + # }}} + Index: /issm/trunk-jpl/src/m/plot/checkplotoptions.py =================================================================== --- /issm/trunk-jpl/src/m/plot/checkplotoptions.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/checkplotoptions.py (revision 21253) @@ -1,3 +1,3 @@ -import numpy as npy +import numpy as np def checkplotoptions(md,options): @@ -13,6 +13,5 @@ ''' - - #units + # {{{ units if options.exist('unit'): if 'km' in options.getfieldvalue('unit','km'): @@ -20,26 +19,26 @@ elif '100km' in options.getfieldvalue('unit','100km'): options.changefieldvalue('unit',10**-5) - - #density + # }}} + # {{{ density if options.exist('density'): density=options.getfieldvalue('density') options.changefieldvalue('density',abs(ceil(density))) - - #show section + # }}} + # {{{ show section if options.exist('showsection'): if 'on' in options.getfieldvalue('showsection','on'): options.changefieldvalue('showsection',4) - - #smooth values + # }}} + # {{{ smooth values if options.exist('smooth'): if 'on' in options.getfieldvalue('smooth','on'): options.changefieldvalue('smooth',0) - - #contouronly values + # }}} + # {{{ contouronly values if options.exist('contouronly'): if 'on' in options.getfieldvalue('contouronly','on'): options.changefieldvalue('contouronly',1) - - #colorbar + # }}} + # {{{ colorbar if options.exist('colorbar'): if 'on' in options.getfieldvalue('colorbar','on'): @@ -47,8 +46,7 @@ elif 'off' in options.getfieldvalue('colorbar','off'): options.changefieldvalue('colorbar',0) - - #text + # }}} + # {{{ text if options.exist('text'): - # text values (coerce to list for consistent functionality) textlist=[] @@ -56,10 +54,9 @@ textlist.extend([text] if isinstance(text,str) else text) numtext=len(textlist) - # text position textpos=options.getfieldvalue('textposition',[0.5,0.5]) if not isinstance(textpos,list): raise Exception('textposition should be passed as a list') - if any(isinstance(i,list) for i in textpos): + if any(isinstance(i,list) for i in textpos): textx=[item[0] for item in textpos] texty=[item[1] for item in textpos] @@ -78,5 +75,5 @@ sizelist=[12] if len(sizelist)==1: - sizelist=npy.tile(sizelist,numtext) + sizelist=np.tile(sizelist,numtext) # font color @@ -88,5 +85,5 @@ colorlist=['k'] if len(colorlist)==1: - colorlist=npy.tile(colorlist,numtext) + colorlist=np.tile(colorlist,numtext) # textweight @@ -98,5 +95,5 @@ weightlist=['normal'] if len(weightlist)==1: - weightlist=npy.tile(weightlist,numtext) + weightlist=np.tile(weightlist,numtext) # text rotation @@ -108,5 +105,5 @@ rotationlist=[0] if len(rotationlist)==1: - rotationlist=npy.tile(rotationlist,numtext) + rotationlist=np.tile(rotationlist,numtext) options.changefieldvalue('text',textlist) @@ -117,18 +114,18 @@ options.changefieldvalue('textweight',weightlist) options.changefieldvalue('textrotation',rotationlist) - - #expdisp + # }}} + # {{{ expdisp expdispvaluesarray=[] expstylevaluesarray=[] expstylevalues=[] if options.exist('expstyle'): - expstylevalues=options.getfieldvalue('expstyle') - if type(expstylevalues)==str: - expstylevalues=[expstylevalues] + expstylevalues=options.getfieldvalue('expstyle') + if type(expstylevalues)==str: + expstylevalues=[expstylevalues] if options.exist('expdisp'): expdispvalues=options.getfieldvalue('expdisp') - if type(expdispvalues)==str: - expdispvalues=[expdispvalues] - for i in npy.arange(len(expdispvalues)): + if type(expdispvalues)==str: + expdispvalues=[expdispvalues] + for i in np.arange(len(expdispvalues)): expdispvaluesarray.append(expdispvalues[i]) if len(expstylevalues)>i: @@ -136,14 +133,13 @@ else: expstylevaluesarray.append('-k') - options.changefieldvalue('expstyle',expstylevaluesarray) options.changefieldvalue('expdisp',expdispvaluesarray) - - #latlonnumbering + # }}} + # {{{ latlonnumbering if options.exist('latlonclick'): if 'on' in options.getfieldvalue('latlonclick','on'): options.changefieldvalue('latlonclick',1) - - #northarrow + # }}} + # {{{ northarrow if options.exist('northarrow'): if 'on' in options.getfieldvalue('northarrow','on'): @@ -152,6 +148,6 @@ Ly=max(md.mesh.y)-min(md.mesh.y) options.changefieldvalue('northarrow',[min(md.mesh.x)+1./6.*Lx, min(md.mesh.y)+5./6.*Ly, 1./15.*Ly, 0.25, 1./250.*Ly]) - - #scale ruler + # }}} + # {{{ scale ruler if options.exist('scaleruler'): if 'on' in options.exist('scaleruler','on'): @@ -159,10 +155,8 @@ Ly=max(md.mesh.y)-min(md.mesh.y) options.changefieldvalue('scaleruler',[min(md.mesh.x)+6./8.*Lx, min(md.mesh.y)+1./10.*Ly, 10**(ceil(log10(Lx)))/5, floor(Lx/100), 5]) - - #log scale + # }}} + # {{{ log scale if options.exist('log'): - if options.exist('clim'): - options.changefieldvalue('clim',log(options.getfieldvalue('clim'))/log(options.getfieldvalue('log'))) - options.changefieldvalue('cutoff',log(options.getfieldvalue('cutoff',1.5))/log(options.getfieldvalue('log'))) - + options.changefieldvalue('cutoff',np.log10(options.getfieldvalue('cutoff',1.5))/np.log10(options.getfieldvalue('log'))) + # }}} return options Index: /issm/trunk-jpl/src/m/plot/plot_manager.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_manager.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/plot_manager.py (revision 21253) @@ -39,20 +39,14 @@ #parse options and get a structure of options options=checkplotoptions(md,options) - #get data to be plotted data=options.getfieldvalue('data') #add ticklabel has a default option options.addfielddefault('ticklabels','on') - #initialize plot handle variable - #handle=None - # initialize subplot - #p.subplot(nrows,ncols,i,aspect='equal') - - ##basemap plot + # {{{ basemap plot TOFIX #if options.exist('basemap'): # plot_basemap(md,data,options,nrows,ncols,i) - - #overlay plot + # }}} + # {{{ overlay plot if options.exist('overlay') and overlaysupport: plot_overlay(md,data,options,ax) @@ -60,8 +54,7 @@ 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 + # }}} + # {{{ dealing with special plot (mesh for now) if isinstance(data,(str,unicode)): - # convert string to lower case for a case-insensitive comparison if data.lower()=='mesh': @@ -76,18 +69,28 @@ else: print "WARNING: '%s' is not implemented or is not a valid string for option 'data'" % data - - #elif data in vars(md): - #else: - #print "'data' not a string, plotting model properties yet to be implemented..." - - #Gridded plot - - #Section plot - - #Profile plot + # }}} + # {{{ Gridded plot TODO + # }}} + # {{{ Section plot TODO + # }}} + # {{{ Profile plot TODO + # }}} #process data and model x,y,z,elements,is2d,isplanet=processmesh(md,data,options) data2,datatype=processdata(md,data,options) + #plot unit + plot_unit(x,y,z,elements,data2,is2d,isplanet,datatype,options,ax) + #apply all options + applyoptions(md,data2,options,fig,ax) + + #ground overlay on kml plot_unit + + # Bits and pieces + #initialize plot handle variable + #handle=None + + # initialize subplot + #p.subplot(nrows,ncols,i,aspect='equal') #standard plot @@ -95,9 +98,5 @@ # p.subplot(nrows,ncols,i,aspect='equal') - #plot unit - plot_unit(x,y,z,elements,data2,is2d,isplanet,datatype,options,ax) - - #apply all options - applyoptions(md,data2,options,fig,ax) - - #ground overlay on kml plot_unit + #elif data in vars(md): + #else: + #print "'data' not a string, plotting model properties yet to be implemented..." Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 21253) @@ -1,92 +1,115 @@ from cmaptools import truncate_colormap try: - import pylab as p - import matplotlib as mpl - import matplotlib.pyplot as plt - import numpy as npy + import pylab as p + import matplotlib as mpl + import matplotlib.pyplot as plt + import numpy as np 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) -def plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options,ax): - """ - PLOT_UNIT - unit plot, display data + See also: PLOTMODEL, PLOT_MANAGER + """ + + #edgecolor + edgecolor=options.getfieldvalue('edgecolor','None') + + # colormap + # {{{ give number of colorlevels and transparency + colorlevels=options.getfieldvalue('colorlevels',128) + alpha=options.getfieldvalue('alpha',1) + # }}} + # {{{ define wich colormap to use + try: + defaultmap=plt.cm.viridis + except AttributeError: + print("Viridis can't be found (probably too old Matplotlib) reverting to gnuplot 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) + # }}} + # {{{ Get the colormap limits + if options.exist('clim'): + lims=options.getfieldvalue('clim',[np.amin(data),np.amax(data)]) + elif options.exist('caxis'): + lims=options.getfieldvalue('caxis',[np.amin(data),np.amax(data)]) + else: + if np.amin(data)==np.amax(data): + lims=[np.amin(data)-0.5,np.amax(data)+0.5] + else: + lims=[np.amin(data),np.amax(data)] + # }}} + # {{{ Set the spread of the colormap (default is normal + if options.exist('log'): + norm = mpl.colors.LogNorm(vmin=lims[0], vmax=lims[1]) + else: + norm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1]) + # }}} + + # Plot depending on the datatype + # {{{ data are on elements + if datatype==1: + if is2d: + if options.exist('mask'): + triangles=mpl.tri.Triangulation(x,y,elements,data.mask) + else: + triangles=mpl.tri.Triangulation(x,y,elements) + tri=ax.tripcolor(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,edgecolors=edgecolor) + else: + raise ValueError('plot_unit error: 3D element plot not supported yet') + return + # }}} + # {{{ data are on nodes + elif datatype==2: + if is2d: + if options.exist('mask'): + EltMask=np.asarray([np.any(np.in1d(index,np.where(data.mask))) for index in elements]) + triangles=mpl.tri.Triangulation(x,y,elements,EltMask) + else: + triangles=mpl.tri.Triangulation(x,y,elements) + ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha) +# ax.tricontourf(triangles,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 + # }}} + # {{{ plotting quiver (TODO) + elif datatype==3: + vx=data[:,0] + vy=data[:,1] + #TODO write plot_quiver.py to handle this here + color=np.sqrt(vx**2+vy**2) + scale=options.getfieldvalue('scale',1000) + width=options.getfieldvalue('width',0.005*(np.amax(x)-np.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 + # }}} + # {{{ plotting P1 Patch (TODO) + elif datatype==4: + print 'plot_unit message: P1 patch plot not implemented yet' + return + # }}} + # {{{ plotting P0 Patch (TODO) + elif datatype==5: + print 'plot_unit message: P0 patch plot not implemented yet' + return + # }}} + else: + raise ValueError('datatype=%d not supported' % datatype) - 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/plotmodel.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/plotmodel.py (revision 21253) @@ -13,6 +13,5 @@ def plotmodel(md,*args): - ''' - at command prompt, type 'plotdoc' for additional documentation + ''' at command prompt, type 'plotdoc' for additional documentation ''' Index: /issm/trunk-jpl/src/m/plot/processdata.py =================================================================== --- /issm/trunk-jpl/src/m/plot/processdata.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/processdata.py (revision 21253) @@ -1,108 +1,132 @@ from math import isnan -import numpy as npy +import numpy as np def processdata(md,data,options): - """ - PROCESSDATA - process data to be plotted + """ + 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 + """ + # {{{ Initialisation and grabbing auxiliaries + # 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") + if 'numberofvertices2d' in dir(md.mesh): + numberofvertices2d=md.mesh.numberofvertices2d + numberofelements2d=md.mesh.numberofelements2d + else: + numberofvertices2d=np.nan + numberofelements2d=np.nan + procdata=np.copy(data) + #initialize datatype + datatype=0 + # init patches + # get datasize + if np.ndim(procdata)==1: + datasize=np.array([len(procdata),1]) + else: + datasize=np.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 np.any(np.isnan(procdata)): + lb=np.nanmin(procdata) + ub=np.nanmax(procdata) + if lb==ub: + lb=lb-0.5 + ub=ub+0.5 + nanfill=lb-1 + #procdata[np.isnan(procdata)]=nanfill + procdata=np.ma.array(procdata,mask=np.isnan(procdata)) + 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" + # }}} + # {{{ log + if options.exist('log'): + cutoff=options.getfieldvalue('log',1) + procdata[np.where(procdata1 and datasize[0]!= md.mesh.numberofvertices+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') + # }}} + # {{{ element data + if datasize[0]==md.mesh.numberofelements and datasize[1]==1: + print'ploting elements' + #initialize datatype if non patch + if datatype!=4 and datatype!=5: + datatype=1 + # {{{mask + if options.exist('mask'): + flags=options.getfieldvalue('mask') + hide=np.invert(flags) + if np.size(flags)==md.mesh.numberofvertices: + EltMask=np.asarray([np.any(np.in1d(index,np.where(hide))) for index in md.mesh.elements-1]) + procdata=np.ma.array(procdata,mask=EltMask) + options.addfielddefault('cmap_set_bad','w') + elif np.size(flags)==md.mesh.numberofelements: + procdata=np.ma.array(procdata,mask=hide) + options.addfielddefault('cmap_set_bad','w') + else: + print('plotmodel warning: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements') + # }}} + # }}} + # {{{ node data + if datasize[0]==md.mesh.numberofvertices and datasize[1]==1: + print'ploting nodes' + datatype=2 + # {{{ Mask + if options.exist('mask'): + flags=options.getfieldvalue('mask') + hide=np.invert(flags) + if np.size(flags)==md.mesh.numberofvertices: + procdata=np.ma.array(procdata,mask=hide) + options.addfielddefault('cmap_set_bad','w') + elif np.size(flags)==md.mesh.numberofelements: + NodeMask=np.zeros(np.shape(md.mesh.x),dtype=bool) + HideElt=md.mesh.elements[np.where(hide)[0]]-1 + NodeMask[HideElt]=True + procdata=np.ma.array(procdata,mask=NodeMask) + options.addfielddefault('cmap_set_bad','w') + else: + print('plotmodel warning: mask length not supported yet (supported length are md.mesh.numberofvertices and md.mesh.numberofelements') + # }}} + # }}} + # {{{ 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] - 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 and datasize[0]!= md.mesh.numberofvertices+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? + #mask? #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 + # }}} + # {{{ convert rank-2 array to rank-1 + if np.ndim(procdata)==2 and np.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 + # }}} Index: /issm/trunk-jpl/src/m/plot/processmesh.py =================================================================== --- /issm/trunk-jpl/src/m/plot/processmesh.py (revision 21252) +++ /issm/trunk-jpl/src/m/plot/processmesh.py (revision 21253) @@ -1,4 +1,3 @@ from math import isnan -import MatlabFuncs as m import numpy as npy @@ -13,40 +12,39 @@ """ - #some checks + # {{{ check mesh size parameters if md.mesh.numberofvertices==0: raise ValueError('processmesh error: mesh is empty') if md.mesh.numberofvertices==md.mesh.numberofelements: raise ValueError('processmesh error: the number of elements is the same as the number of nodes') - + # }}} + # {{{ treating non data plots mesh if len(data)==0 or not isinstance(data,dict): - if 'latlon' not in options.getfieldvalue('coord','xy').lower(): #convert to lower case for comparison - x=md.mesh.x - if 'x2d' in dir(md.mesh): x2d=md.mesh.x2d - y=md.mesh.y - if 'y2d' in dir(md.mesh): y2d=md.mesh.x2d + try: + x=md.mesh.x2d + except AttributeError: + x=md.mesh.x + try: + y=md.mesh.x2d + except AttributeError: + y=md.mesh.y else: x=md.mesh.long y=md.mesh.lat - - if 'z' in dir(md.mesh): + try: z=md.mesh.z - else: + except AttributeError: z=npy.zeros_like(md.mesh.x) - if 'elements2d' in dir(md.mesh): - elements2d=md.mesh.elements2d - elements2d=elements2d-1 # subtract one since python indexes from zero - elements=md.mesh.elements - elements=elements-1 + try: + elements2d=md.mesh.elements2d-1 + except AttributeError: + elements=md.mesh.elements-1 #is it a 2D plot? - if md.mesh.dimension()==2: + if md.mesh.dimension()==2 or options.getfieldvalue('layer',0)>=1: is2d=1 else: - if options.getfieldvalue('layer',0)>=1: - is2d=1 - else: - is2d=0 + is2d=0 #layer projection? @@ -54,10 +52,10 @@ if 'latlon' in options.getfieldvalue('coord','xy').lower(): raise ValueError('processmesh error: cannot work with 3D mesh in lat-lon coords') - #we modify the mesh temporarily to a 2D mesh from which the 3D mesh was extruded - x=x2d - y=y2d - z=zeros(size(x2d)) - elements=elements2d - + #we modify the mesh temporarily to a 2D mesh from which the 3D mesh was extruded + #Basile: does not seem necessary as we already picked x as x2d + # x=x2d + # y=y2d + # z=zeros(size(x2d)) + # elements=elements2d else: #Process mesh for plotting