Index: /issm/trunk-jpl/src/m/plot/applyoptions.py =================================================================== --- /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 23919) +++ /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 23920) @@ -1,4 +1,4 @@ import numpy as np -from cmaptools import truncate_colormap +from cmaptools import getcolormap from plot_contour import plot_contour from plot_streamlines import plot_streamlines @@ -6,286 +6,287 @@ try: - from matplotlib.ticker import MaxNLocator - from mpl_toolkits.axes_grid1 import make_axes_locatable - from mpl_toolkits.mplot3d import Axes3D - import matplotlib as mpl - import matplotlib.pyplot as plt + from matplotlib.ticker import MaxNLocator + from mpl_toolkits.axes_grid1 import make_axes_locatable + from mpl_toolkits.mplot3d import Axes3D + import matplotlib as mpl + import matplotlib.pyplot as plt 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 applyoptions(md,data,options,fig,axgrid,gridindex): - ''' - APPLYOPTIONS - apply options to current plot - - 'plotobj' is the object returned by the specific plot call used to - render the data. This object is used for adding a colorbar. - - Usage: - applyoptions(md,data,options) - - See also: PLOTMODEL, PARSE_OPTIONS - ''' - - # get handle to current figure and axes instance - #fig = p.gcf() - ax= axgrid[gridindex] - - # {{{ 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 - if options.exist('title'): - title=options.getfieldvalue('title') - if options.exist('titlefontsize'): - titlefontsize=options.getfieldvalue('titlefontsize') - else: - titlefontsize=fontsize - if options.exist('titlefontweight'): - titlefontweight=options.getfieldvalue('titlefontweight') - else: - titlefontweight=fontweight - #title font - titlefont=font.copy() - titlefont['size']=titlefontsize - titlefont['weight']=titlefontweight - ax.set_title(title,**titlefont) - # }}} - # {{{ xlabel, ylabel, zlabel - if options.exist('labelfontsize'): - labelfontsize=options.getfieldvalue('labelfontsize') - else: - labelfontsize=fontsize - if options.exist('labelfontweight'): - labelfontweight=options.getfieldvalue('labelfontweight') - else: - labelfontweight=fontweight - - #font dict for labels - labelfont=font.copy() - labelfont['fontsize']=labelfontsize - labelfont['fontweight']=labelfontweight - - if options.exist('xlabel'): - ax.set_xlabel(options.getfieldvalue('xlabel'),**labelfont) - if options.exist('ylabel'): - ax.set_ylabel(options.getfieldvalue('ylabel'),**labelfont) - if options.exist('zlabel'): - ax.set_zlabel(options.getfieldvalue('zlabel'),**labelfont) - # }}} - # {{{ xticks, yticks, zticks (tick locations) - if options.exist('xticks'): - if options.exist('xticklabels'): - xticklabels=options.getfieldvalue('xticklabels') - ax.set_xticks(options.getfieldvalue('xticks'),xticklabels) - else: - ax.set_xticks(options.getfieldvalue('xticks')) - if options.exist('yticks'): - if options.exist('yticklabels'): - yticklabels=options.getfieldvalue('yticklabels') - ax.set_yticks(options.getfieldvalue('yticks'),yticklabels) - else: - ax.set_yticks(options.getfieldvalue('yticks')) - if options.exist('zticks'): - if options.exist('zticklabels'): - zticklabels=options.getfieldvalue('zticklabels') - ax.set_zticks(options.getfieldvalue('zticks'),zticklabels) - else: - ax.set_zticks(options.getfieldvalue('zticks')) - # }}} - # {{{ xticklabels,yticklabels,zticklabels - if options.getfieldvalue('ticklabels','off')=='off' or options.getfieldvalue('ticklabels',0)==0: - options.addfielddefault('xticklabels',[]) - options.addfielddefault('yticklabels',[]) - # TODO check if ax has a z-axis (e.g. is 3D) - if options.exist('xticklabels'): - xticklabels=options.getfieldvalue('xticklabels') - ax.set_xticklabels(xticklabels) - if options.exist('yticklabels'): - yticklabels=options.getfieldvalue('yticklabels') - ax.set_yticklabels(yticklabels) - if options.exist('zticklabels'): - zticklabels=options.getfieldvalue('zticklabels') - ax.set_zticklabels(zticklabels) - # }}} - # {{{ ticklabel notation - #ax.ticklabel_format(style='sci',scilimits=(0,0)) - # }}} - # {{{ ticklabelfontsize - if options.exist('ticklabelfontsize'): - for label in ax.get_xticklabels() + ax.get_yticklabels(): - label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) - if int(md.mesh.dimension)==3: - for label in ax.get_zticklabels(): - label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) - # }}} - # {{{ view TOFIX - #if int(md.mesh.dimension) == 3 and options.exist('layer'): - # #options.getfieldvalue('view') ? - # ax=fig.gca(projection='3d') - #plt.show() - # }}} - # {{{ axis - if options.exist('axis'): - if options.getfieldvalue('axis',True)=='off': - ax.ticklabel_format(style='plain') - p.setp(ax.get_xticklabels(), visible=False) - p.setp(ax.get_yticklabels(), visible=False) - # }}} - # {{{ box - if options.exist('box'): - eval(options.getfieldvalue('box')) - # }}} - # {{{ xlim, ylim, zlim - if options.exist('xlim'): - ax.set_xlim(options.getfieldvalue('xlim')) - if options.exist('ylim'): - ax.set_ylim(options.getfieldvalue('ylim')) - if options.exist('zlim'): - ax.set_zlim(options.getfieldvalue('zlim')) - # }}} - # {{{ latlon TODO - # }}} - # {{{ Basinzoom TODO - # }}} - # {{{ ShowBasins TODO - # }}} - # {{{ clim - if options.exist('clim'): - lims=options.getfieldvalue('clim') - assert len(lims)==2, 'error, clim should be passed as a list of length 2' - elif options.exist('caxis'): - lims=options.getfieldvalue('caxis') - assert len(lims)==2, 'error, caxis should be passed as a list of length 2' - options.addfielddefault('clim',lims) - else: - if len(data)>0: - lims=[data.min(),data.max()] - else: - lims=[0,1] - # }}} - # {{{ shading TODO - #if options.exist('shading'): - # }}} - # {{{ grid - if options.exist('grid'): - if 'on' in options.getfieldvalue('grid','on'): - ax.grid() - # }}} - # {{{ colormap - if options.exist('colornorm'): - norm=options.getfieldvalue('colornorm') - if options.exist('colormap'): - cmap=options.getfieldvalue('colormap') - cbar_extend=0 - if options.exist('cmap_set_over'): - cbar_extend+=1 - if options.exist('cmap_set_under'): - cbar_extend+=2 - # }}} - # {{{ contours - if options.exist('contourlevels'): - plot_contour(md,data,options,ax) - # }}} - # {{{ 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=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing')) - cb.set_ticks(locs) - if options.exist('colorbarlines'): - locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines')) - cb.add_lines(locs,['k' for i in range(len(locs))],np.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',np.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.draw_all() - if options.exist('colorbarfontsize'): - colorbarfontsize=options.getfieldvalue('colorbarfontsize') - cb.ax.tick_params(labelsize=colorbarfontsize) - # cb.set_ticks([0,-10]) - # cb.set_ticklabels([-10,0,10]) - if options.exist('colorbarticks'): - colorbarticks=options.getfieldvalue('colorbarticks') - cb.set_ticks(colorbarticks) - plt.sca(ax) # return to original axes control - # }}} - # {{{ expdisp - if options.exist('expdisp'): - expdisp(ax,options) - # }}} - # {{{ 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 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 - # }}} + ''' + APPLYOPTIONS - apply options to current plot + + 'plotobj' is the object returned by the specific plot call used to + render the data. This object is used for adding a colorbar. + + Usage: + applyoptions(md,data,options) + + See also: PLOTMODEL, PARSE_OPTIONS + ''' + + # get handle to current figure and axes instance + #fig = p.gcf() + ax= axgrid[gridindex] + + # {{{ 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 + if options.exist('title'): + title=options.getfieldvalue('title') + if options.exist('titlefontsize'): + titlefontsize=options.getfieldvalue('titlefontsize') + else: + titlefontsize=fontsize + if options.exist('titlefontweight'): + titlefontweight=options.getfieldvalue('titlefontweight') + else: + titlefontweight=fontweight + #title font + titlefont=font.copy() + titlefont['size']=titlefontsize + titlefont['weight']=titlefontweight + ax.set_title(title,**titlefont) + # }}} + # {{{ xlabel, ylabel, zlabel + if options.exist('labelfontsize'): + labelfontsize=options.getfieldvalue('labelfontsize') + else: + labelfontsize=fontsize + if options.exist('labelfontweight'): + labelfontweight=options.getfieldvalue('labelfontweight') + else: + labelfontweight=fontweight + + #font dict for labels + labelfont=font.copy() + labelfont['fontsize']=labelfontsize + labelfont['fontweight']=labelfontweight + + if options.exist('xlabel'): + ax.set_xlabel(options.getfieldvalue('xlabel'),**labelfont) + if options.exist('ylabel'): + ax.set_ylabel(options.getfieldvalue('ylabel'),**labelfont) + if options.exist('zlabel'): + ax.set_zlabel(options.getfieldvalue('zlabel'),**labelfont) + # }}} + # {{{ xticks, yticks, zticks (tick locations) + if options.exist('xticks'): + if options.exist('xticklabels'): + xticklabels=options.getfieldvalue('xticklabels') + ax.set_xticks(options.getfieldvalue('xticks'),xticklabels) + else: + ax.set_xticks(options.getfieldvalue('xticks')) + if options.exist('yticks'): + if options.exist('yticklabels'): + yticklabels=options.getfieldvalue('yticklabels') + ax.set_yticks(options.getfieldvalue('yticks'),yticklabels) + else: + ax.set_yticks(options.getfieldvalue('yticks')) + if options.exist('zticks'): + if options.exist('zticklabels'): + zticklabels=options.getfieldvalue('zticklabels') + ax.set_zticks(options.getfieldvalue('zticks'),zticklabels) + else: + ax.set_zticks(options.getfieldvalue('zticks')) + # }}} + # {{{ xticklabels,yticklabels,zticklabels + if options.getfieldvalue('ticklabels','off')=='off' or options.getfieldvalue('ticklabels',0)==0: + options.addfielddefault('xticklabels',[]) + options.addfielddefault('yticklabels',[]) + # TODO check if ax has a z-axis (e.g. is 3D) + if options.exist('xticklabels'): + xticklabels=options.getfieldvalue('xticklabels') + ax.set_xticklabels(xticklabels) + if options.exist('yticklabels'): + yticklabels=options.getfieldvalue('yticklabels') + ax.set_yticklabels(yticklabels) + if options.exist('zticklabels'): + zticklabels=options.getfieldvalue('zticklabels') + ax.set_zticklabels(zticklabels) + # }}} + # {{{ ticklabel notation + #ax.ticklabel_format(style='sci',scilimits=(0,0)) + # }}} + # {{{ ticklabelfontsize + if options.exist('ticklabelfontsize'): + for label in ax.get_xticklabels() + ax.get_yticklabels(): + label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) + if int(md.mesh.dimension)==3: + for label in ax.get_zticklabels(): + label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) + # }}} + # {{{ view TOFIX + #if int(md.mesh.dimension) == 3 and options.exist('layer'): + # #options.getfieldvalue('view') ? + # ax=fig.gca(projection='3d') + #plt.show() + # }}} + # {{{ axis + if options.exist('axis'): + if options.getfieldvalue('axis',True)=='off': + ax.ticklabel_format(style='plain') + p.setp(ax.get_xticklabels(), visible=False) + p.setp(ax.get_yticklabels(), visible=False) + # }}} + # {{{ box + if options.exist('box'): + eval(options.getfieldvalue('box')) + # }}} + # {{{ xlim, ylim, zlim + if options.exist('xlim'): + ax.set_xlim(options.getfieldvalue('xlim')) + if options.exist('ylim'): + ax.set_ylim(options.getfieldvalue('ylim')) + if options.exist('zlim'): + ax.set_zlim(options.getfieldvalue('zlim')) + # }}} + # {{{ latlon TODO + # }}} + # {{{ Basinzoom TODO + # }}} + # {{{ ShowBasins TODO + # }}} + # {{{ clim + if options.exist('clim'): + lims=options.getfieldvalue('clim') + assert len(lims)==2, 'error, clim should be passed as a list of length 2' + elif options.exist('caxis'): + lims=options.getfieldvalue('caxis') + assert len(lims)==2, 'error, caxis should be passed as a list of length 2' + options.addfielddefault('clim',lims) + else: + if len(data)>0: + lims=[data.min(),data.max()] + else: + lims=[0,1] + # }}} + # {{{ shading TODO + #if options.exist('shading'): + # }}} + # {{{ grid + if options.exist('grid'): + if 'on' in options.getfieldvalue('grid','on'): + ax.grid() + # }}} + # {{{ colormap + if options.exist('colornorm'): + norm=options.getfieldvalue('colornorm') + if options.exist('colormap'): + cmap=getcolormap(options) + cbar_extend=0 + if options.exist('cmap_set_over'): + cbar_extend+=1 + if options.exist('cmap_set_under'): + cbar_extend+=2 + # }}} + # {{{ contours + if options.exist('contourlevels'): + plot_contour(md,data,options,ax) + # }}} + # {{{ 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=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbartickspacing')) + cb.set_ticks(locs) + if options.exist('colorbarlines'): + locs=np.arange(lims[0],lims[1]+1,options.getfieldvalue('colorbarlines')) + cb.add_lines(locs,['k' for i in range(len(locs))],np.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',np.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.draw_all() + if options.exist('colorbarfontsize'): + colorbarfontsize=options.getfieldvalue('colorbarfontsize') + cb.ax.tick_params(labelsize=colorbarfontsize) + # cb.set_ticks([0,-10]) + # cb.set_ticklabels([-10,0,10]) + if options.exist('colorbarticks'): + colorbarticks=options.getfieldvalue('colorbarticks') + cb.set_ticks(colorbarticks) + plt.sca(ax) # return to original axes control + # }}} + # {{{ expdisp + if options.exist('expdisp'): + expdisp(ax,options) + # }}} + # {{{ 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 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/colormaps/cmaptools.py =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/cmaptools.py (revision 23919) +++ /issm/trunk-jpl/src/m/plot/colormaps/cmaptools.py (revision 23920) @@ -1,8 +1,59 @@ import numpy as np +from demmap import * try: import matplotlib as mpl + import matplotlib.pyplot as plt + from matplotlib.colors import ListedColormap, rgb_to_hsv, hsv_to_rgb except ImportError: print('cannot import matplotlib, no plotting capabilities enabled') + +def getcolormap(options): + ''' + get colormap from options and apply + + default: viridis + supported: + matplotlib defaults (see: pyplot.colormaps()) + Rignot + demmap(50,-300,1200) + demmap(50,-300,1200,'ibcao') + + Usage: + cmap = getcolormap(options) + ''' + + map_name = options.getfieldvalue('colormap') + cmap = 'viridis' + + # already a valid colormap, the name of a valid colormap, or empty (use default) + if type(map_name) == mpl.colors.ListedColormap: + return map_name + elif map_name in plt.colormaps(): + return map_name + elif map_name == '': + return cmap + + # if we don't have a matching colormap, build one + if map_name == 'Rignot': + alpha=options.getfieldvalue('alpha',1) + cmap = np.array((np.linspace(0,1,128,False),np.ones(128,),np.ones(128,))).T + cmap[:,1] =np.maximum(np.minimum((0.1+cmap[:,0]**(1/alpha)),1),0) + cmap = hsv_to_rgb(cmap) + # construct a colormap object from an array of shape (n,3/4) + cmap = ListedColormap(cmap) + + #elif map_name == 'Ala': + + else: + # map is a library or executable function that constructs a colormap, + # function must be imported above + try: + cmap = ListedColormap(eval(map_name)) + except: + raise RuntimeError("getcolormap: Error: provided colormap must be supported map or map-constructing function with syntax: 'jet' or 'function(args)'") + + return cmap + def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100): Index: /issm/trunk-jpl/src/m/plot/colormaps/demmap.py =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/demmap.py (revision 23920) +++ /issm/trunk-jpl/src/m/plot/colormaps/demmap.py (revision 23920) @@ -0,0 +1,75 @@ +import numpy as np +from MatlabFuncs import strcmpi +from landcolor import landcolor +from seacolor import seacolor +from ibcao import ibcao + +def demmap(ncolors,minZ,maxZ,colorscheme='dem'): + '''DEMMAP - concatenate sea and land color deping on zmin and zmax + + Usage: + cmap = demmap(n,zmin,zmax,colorscheme) + + Example: + cmap = demmap(50,-300,1200) + cmap = demmap(50,-300,1200,'dem') + cmap = demmap(50,-300,1200,'ibcao') + ''' + + if type(colorscheme) != str: + raise RuntimeError('demmap: Error: optional argument "colorscheme" should be a string') + + # determine appropriate number of sea and land colors + if minZ == maxZ: + maxZ = minZ+1 + + cmn = minZ + cmx = maxZ + + # determine appropriate number of sea and land colors + if minZ >= 0: + nsea = 0 + nland = ncolors + elif maxZ <= 0: + nland = 0 + nsea = ncolors + else: + # find optimal ratio of land to sea colors + maxminratio = maxZ/abs(minZ) + n1 = np.floor(ncolors/2) + n2 = np.ceil(ncolors/2) + if maxminratio>1: + sea = np.arange(1,n1+1) + land = np.arange(ncolors-1,n2-1,-1) + else: + land = np.arange(1,n1+1) + sea = np.arange(ncolors-1,n2-1,-1) + + ratio = land/sea + errors = abs(ratio - maxminratio) / maxminratio + indx = np.where(errors == min(errors)) + nsea = sea[indx] + nland = land[indx] + + # determine color limits + seaint = abs(minZ)/nsea + landint = maxZ/nland + if seaint >= landint: + interval = seaint + else: + interval = landint + + cmn = -nsea*interval*(1 + 1e-9) # zero values treated as land + cmx = nland*interval + + + clim = [cmn,cmx] + + if strcmpi(colorscheme,'dem'): + # concatenate and transpose to match matplotlib's colormap format + cmap = np.concatenate((seacolor(nsea),landcolor(nland)**1.3),axis=1).T + elif strcmpi(colorscheme,'ibcao'): + cmap = ibcao(nsea,nland) + + return cmap + Index: /issm/trunk-jpl/src/m/plot/colormaps/ibcao.py =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/ibcao.py (revision 23920) +++ /issm/trunk-jpl/src/m/plot/colormaps/ibcao.py (revision 23920) @@ -0,0 +1,82 @@ +import numpy as np + +def ibcao(nsea,nland): + '''IBCAO - IBCAO color map + + Usage: + map = ibcap(nsea,nland) + ''' + + Jsea = [ + 0.18039,0.29020,0.57255, + 0.18039,0.29020,0.57255, + 0.05882,0.44314,0.65490, + 0.05882,0.44314,0.65490, + 0.02745,0.49804,0.73725, + 0.02745,0.49804,0.73725, + 0.01176,0.54510,0.78824, + 0.01176,0.54510,0.78824, + 0.00784,0.63529,0.83922, + 0.00784,0.63529,0.83922, + 0.06667,0.71765,0.86667, + 0.06667,0.71765,0.86667, + 0.17647,0.75294,1.00000, + 0.17647,0.75294,1.00000, + 0.23529,0.76471,0.85882, + 0.23529,0.76471,0.85882, + 0.24314,0.76471,0.83922, + 0.24314,0.76471,0.83922, + 0.25882,0.76078,0.81176, + 0.25882,0.76078,0.81176, + 0.27451,0.76078,0.76078, + 0.27451,0.76078,0.76078, + 0.41961,0.78431,0.74902, + 0.41961,0.78431,0.74902, + 0.60000,0.83137,0.74902, + 0.60000,0.83137,0.74902 + ] + + Jland = [ + 0.85098,0.84314,0.30588, + 0.85098,0.84314,0.30588, + 0.93333,0.89020,0.41961, + 0.93333,0.89020,0.41961, + 0.93725,0.80784,0.35686, + 0.93725,0.80784,0.35686, + 0.89804,0.74510,0.31765, + 0.89804,0.74510,0.31765, + 0.85098,0.63922,0.21961, + 0.85098,0.63922,0.21961, + 0.75686,0.55294,0.22353, + 0.75686,0.55294,0.22353, + 0.71765,0.50980,0.22353, + 0.71765,0.50980,0.22353, + 0.68627,0.48235,0.21961, + 0.68627,0.48235,0.21961, + 0.65490,0.45882,0.21569, + 0.65490,0.45882,0.21569, + 0.58824,0.39608,0.20392, + 0.58824,0.39608,0.20392, + 1.00000,1.00000,1.00000 + ] + + # Jsea and Jland are each a series of r,g,b triples, reshape them as such + + lsea = int(len(Jsea)/3) + Jsea = np.array(Jsea).reshape(lsea,3) + a = np.linspace(1,lsea,nsea) + b = np.arange(1,lsea+1) + # interpolate color on each channel r,g,b + ysea = np.array([np.interp(a, b, Jsea[:,i]) for i in range(3)]) + + lland = int(len(Jland)/3) + Jland = np.array(Jland).reshape(lland,3) + a = np.linspace(1,lland,nland) + b = np.arange(1,lland+1) + # interpolate color on each channel r,g,b + yland = np.array([np.interp(a, b, Jland[:,i]) for i in range(3)]) + + # concatenate and transpose to match matplotlib's colormap format + map = np.concatenate((ysea,yland),axis=1).T + + return map Index: /issm/trunk-jpl/src/m/plot/colormaps/landcolor.py =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/landcolor.py (revision 23920) +++ /issm/trunk-jpl/src/m/plot/colormaps/landcolor.py (revision 23920) @@ -0,0 +1,81 @@ +import numpy as np + +def landcolor(n=256): + '''LANDCOLOR Land colormap''' + + J = [ + 0.095678,0.53427,0.21682, + 0.15785,0.5979,0.23274, + 0.21286,0.64673,0.2514, + 0.26411,0.68789,0.27268, + 0.32959,0.72416,0.31308, + 0.39794,0.75695,0.36038, + 0.46153,0.7871,0.40624, + 0.52108,0.81516,0.45135, + 0.57702,0.84152,0.49547, + 0.62973,0.86645,0.53891, + 0.67946,0.89016,0.58187, + 0.72647,0.91282,0.62427, + 0.77095,0.93455,0.66619, + 0.81306,0.95546,0.70772, + 0.85292,0.97563,0.7489, + 0.89066,0.99514,0.78976, + 0.88379,0.98595,0.77038, + 0.86389,0.96758,0.73236, + 0.84615,0.94972,0.69623, + 0.8303,0.93233,0.66186, + 0.81612,0.91536,0.6291, + 0.80341,0.8988,0.59784, + 0.79201,0.8826,0.56795, + 0.78191,0.86676,0.53946, + 0.7729,0.85123,0.51224, + 0.76479,0.83602,0.48615, + 0.75747,0.8211,0.46111, + 0.75084,0.80645,0.43704, + 0.74506,0.79206,0.41414, + 0.73981,0.77792,0.39211, + 0.73501,0.76401,0.37089, + 0.73068,0.75033,0.35052, + 0.72683,0.73685,0.33106, + 0.72042,0.72074,0.31228, + 0.71032,0.70085,0.29417, + 0.69761,0.67821,0.27694, + 0.68489,0.65558,0.26026, + 0.67235,0.63313,0.24418, + 0.65997,0.61082,0.22889, + 0.64775,0.58874,0.21406, + 0.63568,0.56689,0.19983, + 0.62376,0.54527,0.18622, + 0.61197,0.52391,0.17299, + 0.60033,0.50283,0.16046, + 0.58881,0.48203,0.14832, + 0.57742,0.46151,0.13667, + 0.56616,0.44133,0.12555, + 0.55502,0.4214,0.11472, + 0.54398,0.4019,0.10456, + 0.53306,0.38266,0.094633, + 0.52226,0.36382,0.085242, + 0.51155,0.3453,0.076179, + 0.50095,0.32714,0.067515, + 0.49045,0.30938,0.059259, + 0.48005,0.29193,0.051294, + 0.46973,0.27495,0.043796, + 0.45951,0.25823,0.0365, + 0.44938,0.24206,0.029715, + 0.43934,0.22609,0.023063, + 0.42938,0.21074,0.016949, + 0.41951,0.19556,0.010917, + 0.40971,0.18105,0.0054326, + 0.4,0.16667,0 + ] + + # J is a series of r,g,b triples, reshape it as such + l = int(len(J)/3) + J = np.array(J).reshape(l,3) + a = np.linspace(1,l,n) + b = np.arange(1,l+1) + + # interpolate color on each channel r,g,b + y = np.array([np.interp(a, b, J[:,i]) for i in range(3)]) + + return y Index: /issm/trunk-jpl/src/m/plot/colormaps/seacolor.py =================================================================== --- /issm/trunk-jpl/src/m/plot/colormaps/seacolor.py (revision 23920) +++ /issm/trunk-jpl/src/m/plot/colormaps/seacolor.py (revision 23920) @@ -0,0 +1,62 @@ +import numpy as np + +def seacolor(n=256): + '''SEACOLOR Sea colormap''' + + J = [ + 0.0392, 0, 0.4745, + 0.1020, 0, 0.5373, + 0.1020, 0, 0.5373, + 0.1490, 0, 0.5961, + 0.1490, 0, 0.5961, + 0.1059, 0.0118, 0.6510, + 0.1059, 0.0118, 0.6510, + 0.0627, 0.0235, 0.7059, + 0.0627, 0.0235, 0.7059, + 0.0196, 0.0353, 0.7569, + 0.0196, 0.0353, 0.7569, + 0, 0.0549, 0.7961, + 0, 0.0549, 0.7961, + 0, 0.0863, 0.8235, + 0, 0.0863, 0.8235, + 0, 0.1176, 0.8471, + 0, 0.1176, 0.8471, + 0, 0.1529, 0.8745, + 0, 0.1529, 0.8745, + 0.0471, 0.2667, 0.9059, + 0.0471, 0.2667, 0.9059, + 0.1020, 0.4000, 0.9412, + 0.1020, 0.4000, 0.9412, + 0.0745, 0.4588, 0.9569, + 0.0745, 0.4588, 0.9569, + 0.0549, 0.5216, 0.9765, + 0.0549, 0.5216, 0.9765, + 0.0824, 0.6196, 0.9882, + 0.0824, 0.6196, 0.9882, + 0.1176, 0.6980, 1.0000, + 0.1176, 0.6980, 1.0000, + 0.1686, 0.7294, 1.0000, + 0.1686, 0.7294, 1.0000, + 0.2157, 0.7569, 1.0000, + 0.2157, 0.7569, 1.0000, + 0.2549, 0.7843, 1.0000, + 0.2549, 0.7843, 1.0000, + 0.3098, 0.8235, 1.0000, + 0.3098, 0.8235, 1.0000, + 0.3686, 0.8745, 1.0000, + 0.3686, 0.8745, 1.0000, + 0.5412, 0.8902, 1.0000, + 0.5412, 0.8902, 1.0000, + 0.7373, 0.9020, 1.0000 + ] + + # J is a series of r,g,b triples, reshape it as such + l = int(len(J)/3) + J = np.array(J).reshape(l,3) + a = np.linspace(1,l,n) + b = np.arange(1,l+1) + + # interpolate color on each channel r,g,b + y = np.array([np.interp(a, b, J[:,i]) for i in range(3)]) + + return y Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 23919) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 23920) @@ -1,236 +1,239 @@ -from cmaptools import truncate_colormap +from cmaptools import getcolormap,truncate_colormap from plot_quiver import plot_quiver from scipy.interpolate import griddata import numpy as np try: - import matplotlib as mpl - import matplotlib.pyplot as plt - from mpl_toolkits.axes_grid1 import inset_locator - from mpl_toolkits.mplot3d import Axes3D - from mpl_toolkits.mplot3d.art3d import Poly3DCollection + import matplotlib as mpl + import matplotlib.pyplot as plt + from mpl_toolkits.axes_grid1 import inset_locator + from mpl_toolkits.mplot3d import Axes3D + from mpl_toolkits.mplot3d.art3d import Poly3DCollection 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,fig,axgrid,gridindex): - """ - PLOT_UNIT - unit plot, display data - - Usage: - plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options) - - See also: PLOTMODEL, PLOT_MANAGER - """ - #if we are plotting 3d replace the current axis - if not is2d: - axgrid[gridindex].axis('off') - ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D) - ax.set_axis_bgcolor((0.7,0.7,0.7)) - else: - ax=axgrid[gridindex] - - #edgecolor - edgecolor=options.getfieldvalue('edgecolor','None') - - # colormap - # {{{ give number of colorlevels and transparency - colorlevels=options.getfieldvalue('colorlevels',128) - alpha=options.getfieldvalue('alpha',1) - if alpha<1: - antialiased=True - else: - antialiased=False - # }}} - # {{{ define wich colormap to use - try: - defaultmap=plt.cm.get_cmap('viridis',colorlevels) - except AttributeError: - print("Viridis can't be found (probably too old Matplotlib) reverting to gnuplot colormap") - defaultmap=truncate_colormap('gnuplot2',0.1,0.9,colorlevels) - 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) - options.addfield('colormap',cmap) - # }}} - # {{{ if plotting only one of several layers reduce dataset, same for surface - if options.getfieldvalue('layer',0)>=1: - plotlayer=options.getfieldvalue('layer',0) - if datatype==1: - slicesize=np.shape(elements)[0] - elif datatype in [2,3]: - slicesize=len(x) - data=data[(plotlayer-1)*slicesize:plotlayer*slicesize] - # }}} - # {{{ 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]) - 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) - # }}} - - # 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: - #first deal with colormap - loccmap = plt.cm.ScalarMappable(cmap=cmap) - loccmap.set_array([min(data),max(data)]) - loccmap.set_clim(vmin=min(data),vmax=max(data)) - - #dealing with prism sides - recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T - eltind=np.arange(0,np.shape(elements)[0]) - recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T)) - eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0]))) - recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T)) - eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0]))) - tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1]))) - _, idx, recur = np.unique(tmp, return_index=True, return_counts=True) - recel= recface[idx[np.where(recur==1)]] - recindex=eltind[idx[np.where(recur==1)]] - for i,rectangle in enumerate(recel): - rec=list(zip(x[rectangle],y[rectangle],z[rectangle])) - pl3=Poly3DCollection([rec]) - color=loccmap.to_rgba(data[recindex[i]]) - pl3.set_edgecolor(color) - pl3.set_color(color) - ax.add_collection3d(pl3) - - #dealing with prism bases - triface=np.vstack((elements[:,0:3],elements[:,3:6])) - eltind=np.arange(0,np.shape(elements)[0]) - eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0]))) - tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1]))) - _, idx,recur = np.unique(tmp, return_index=True,return_counts=True) - #we keep only top and bottom elements - triel= triface[idx[np.where(recur==1)]] - triindex=eltind[idx[np.where(recur==1)]] - for i,triangle in enumerate(triel): - tri=list(zip(x[triangle],y[triangle],z[triangle])) - pl3=Poly3DCollection([tri]) - color=loccmap.to_rgba(data[triindex[i]]) - pl3.set_edgecolor(color) - pl3.set_color(color) - ax.add_collection3d(pl3) - - ax.set_xlim([min(x),max(x)]) - ax.set_ylim([min(y),max(y)]) - ax.set_zlim([min(z),max(z)]) - - #raise ValueError('plot_unit error: 3D element plot not supported yet') - return - - # }}} - # {{{ data are on nodes - - elif datatype==2: - if is2d: - if np.ma.is_masked(data): - 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) - #tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha) - if options.exist('log'): - if alpha<1:#help with antialiasing - tri=ax.tricontour(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=0.1,antialiased=antialiased) - tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,antialiased=antialiased) - else: - if alpha<1:#help with antialiasing - tri=ax.tricontour(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=0.1,antialiased=antialiased) - tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,extend='both',antialiased=antialiased) - if edgecolor != 'None': - ax.triplot(x,y,elements,color=edgecolor) - else: - #first deal with the colormap - loccmap = plt.cm.ScalarMappable(cmap=cmap) - loccmap.set_array([min(data),max(data)]) - loccmap.set_clim(vmin=min(data),vmax=max(data)) - - #deal with prism sides - recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T - recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T)) - recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T)) - tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1]))) - _, idx, recur = np.unique(tmp, return_index=True, return_counts=True) - recel= recface[idx[np.where(recur==1)]] - for rectangle in recel: - rec=list(zip(x[rectangle],y[rectangle],z[rectangle])) - pl3=Poly3DCollection([rec]) - color=loccmap.to_rgba(np.mean(data[rectangle])) - pl3.set_edgecolor(color) - pl3.set_color(color) - ax.add_collection3d(pl3) - - #deal with prism faces - triface=np.vstack((elements[:,0:3],elements[:,3:6])) - tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1]))) - _, idx,recur = np.unique(tmp, return_index=True,return_counts=True) - #we keep only top and bottom elements - triel= triface[idx[np.where(recur==1)]] - for triangle in triel: - tri=list(zip(x[triangle],y[triangle],z[triangle])) - pl3=Poly3DCollection([tri]) - color=loccmap.to_rgba(np.mean(data[triangle])) - pl3.set_edgecolor(color) - pl3.set_color(color) - ax.add_collection3d(pl3) - - ax.set_xlim([min(x),max(x)]) - ax.set_ylim([min(y),max(y)]) - ax.set_zlim([min(z),max(z)]) - #raise ValueError('plot_unit error: 3D element plot not supported yet') - return - - # }}} - # {{{ plotting quiver - elif datatype==3: - if is2d: - Q=plot_quiver(x,y,data,options,ax) - else: - raise ValueError('plot_unit error: 3D node plot not supported yet') - 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) + """ + PLOT_UNIT - unit plot, display data + + Usage: + plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options) + + See also: PLOTMODEL, PLOT_MANAGER + """ + #if we are plotting 3d replace the current axis + if not is2d: + axgrid[gridindex].axis('off') + ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D) + ax.set_axis_bgcolor((0.7,0.7,0.7)) + else: + ax=axgrid[gridindex] + + #edgecolor + edgecolor=options.getfieldvalue('edgecolor','None') + + # colormap + # {{{ give number of colorlevels and transparency + colorlevels=options.getfieldvalue('colorlevels',128) + alpha=options.getfieldvalue('alpha',1) + if alpha<1: + antialiased=True + else: + antialiased=False + # }}} + # {{{ define wich colormap to use + try: + defaultmap=plt.cm.get_cmap('viridis',colorlevels) + except AttributeError: + print("Viridis can't be found (probably too old Matplotlib) reverting to gnuplot colormap") + defaultmap=truncate_colormap('gnuplot2',0.1,0.9,colorlevels) + if not options.exist('colormap'): + cmap=defaultmap + else: + cmap=getcolormap(options) + 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) + options.addfield('colormap',cmap) + # }}} + # {{{ if plotting only one of several layers reduce dataset, same for surface + if options.getfieldvalue('layer',0)>=1: + plotlayer=options.getfieldvalue('layer',0) + if datatype==1: + slicesize=np.shape(elements)[0] + elif datatype in [2,3]: + slicesize=len(x) + data=data[(plotlayer-1)*slicesize:plotlayer*slicesize] + # }}} + # {{{ 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]) + 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) + # }}} + + # 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: + #first deal with colormap + loccmap = plt.cm.ScalarMappable(cmap=cmap) + loccmap.set_array([min(data),max(data)]) + loccmap.set_clim(vmin=min(data),vmax=max(data)) + + #dealing with prism sides + recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T + eltind=np.arange(0,np.shape(elements)[0]) + recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T)) + eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0]))) + recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T)) + eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0]))) + tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1]))) + _, idx, recur = np.unique(tmp, return_index=True, return_counts=True) + recel= recface[idx[np.where(recur==1)]] + recindex=eltind[idx[np.where(recur==1)]] + for i,rectangle in enumerate(recel): + rec=list(zip(x[rectangle],y[rectangle],z[rectangle])) + pl3=Poly3DCollection([rec]) + color=loccmap.to_rgba(data[recindex[i]]) + pl3.set_edgecolor(color) + pl3.set_color(color) + ax.add_collection3d(pl3) + + #dealing with prism bases + triface=np.vstack((elements[:,0:3],elements[:,3:6])) + eltind=np.arange(0,np.shape(elements)[0]) + eltind=np.hstack((eltind,np.arange(0,np.shape(elements)[0]))) + tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1]))) + _, idx,recur = np.unique(tmp, return_index=True,return_counts=True) + #we keep only top and bottom elements + triel= triface[idx[np.where(recur==1)]] + triindex=eltind[idx[np.where(recur==1)]] + for i,triangle in enumerate(triel): + tri=list(zip(x[triangle],y[triangle],z[triangle])) + pl3=Poly3DCollection([tri]) + color=loccmap.to_rgba(data[triindex[i]]) + pl3.set_edgecolor(color) + pl3.set_color(color) + ax.add_collection3d(pl3) + + ax.set_xlim([min(x),max(x)]) + ax.set_ylim([min(y),max(y)]) + ax.set_zlim([min(z),max(z)]) + + #raise ValueError('plot_unit error: 3D element plot not supported yet') + return + + # }}} + # {{{ data are on nodes + + elif datatype==2: + if is2d: + if np.ma.is_masked(data): + 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) + #tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha) + if options.exist('log'): + if alpha<1:#help with antialiasing + tri=ax.tricontour(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=0.1,antialiased=antialiased) + tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,antialiased=antialiased) + else: + if alpha<1:#help with antialiasing + tri=ax.tricontour(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=0.1,antialiased=antialiased) + tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,extend='both',antialiased=antialiased) + if edgecolor != 'None': + ax.triplot(x,y,elements,color=edgecolor) + else: + #first deal with the colormap + loccmap = plt.cm.ScalarMappable(cmap=cmap) + loccmap.set_array([min(data),max(data)]) + loccmap.set_clim(vmin=min(data),vmax=max(data)) + + #deal with prism sides + recface=np.vstack((elements[:,0],elements[:,1],elements[:,4],elements[:,3])).T + recface=np.vstack((recface,np.vstack((elements[:,1],elements[:,2],elements[:,5],elements[:,4])).T)) + recface=np.vstack((recface,np.vstack((elements[:,2],elements[:,0],elements[:,3],elements[:,5])).T)) + tmp = np.ascontiguousarray(np.sort(recface)).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1]))) + _, idx, recur = np.unique(tmp, return_index=True, return_counts=True) + recel= recface[idx[np.where(recur==1)]] + for rectangle in recel: + rec=list(zip(x[rectangle],y[rectangle],z[rectangle])) + pl3=Poly3DCollection([rec]) + color=loccmap.to_rgba(np.mean(data[rectangle])) + pl3.set_edgecolor(color) + pl3.set_color(color) + ax.add_collection3d(pl3) + + #deal with prism faces + triface=np.vstack((elements[:,0:3],elements[:,3:6])) + tmp = np.ascontiguousarray(triface).view(np.dtype((np.void, triface.dtype.itemsize * triface.shape[1]))) + _, idx,recur = np.unique(tmp, return_index=True,return_counts=True) + #we keep only top and bottom elements + triel= triface[idx[np.where(recur==1)]] + for triangle in triel: + tri=list(zip(x[triangle],y[triangle],z[triangle])) + pl3=Poly3DCollection([tri]) + color=loccmap.to_rgba(np.mean(data[triangle])) + pl3.set_edgecolor(color) + pl3.set_color(color) + ax.add_collection3d(pl3) + + ax.set_xlim([min(x),max(x)]) + ax.set_ylim([min(y),max(y)]) + ax.set_zlim([min(z),max(z)]) + #raise ValueError('plot_unit error: 3D element plot not supported yet') + return + + # }}} + # {{{ plotting quiver + elif datatype==3: + if is2d: + Q=plot_quiver(x,y,data,options,ax) + else: + raise ValueError('plot_unit error: 3D node plot not supported yet') + 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)