Index: /issm/trunk-jpl/src/m/classes/plotoptions.py =================================================================== --- /issm/trunk-jpl/src/m/classes/plotoptions.py (revision 23562) +++ /issm/trunk-jpl/src/m/classes/plotoptions.py (revision 23563) @@ -52,10 +52,10 @@ print "WARNING: option number %d is not a string and will be ignored." % (i+1) - #get figure number + #get figure number self.figurenumber=rawoptions.getfieldvalue('figure',1) rawoptions.removefield('figure',0) - #get number of subplots - numberofplots=Counter(x for sublist in rawlist for x in sublist if isinstance(x,(str,unicode)))['data'] + #get number of subplots + numberofplots=len([1 for sublist in rawlist for x in sublist if str(x)=='data']) self.numberofplots=numberofplots @@ -76,5 +76,5 @@ #if alloptions flag is on, apply to all plots if (allflag and 'data' not in rawlist[i][0] and '#' not in rawlist[i][0]): - + for j in xrange(numberofplots): self.list[j].addfield(rawlist[i][0],rawlist[i][1]) @@ -106,5 +106,5 @@ raise ValueError('error: in option i-j both i and j must be integers') for j in xrange(int(nums[0])-1,int(nums[1])): - self.list[j].addfield(field,rawlist[i][1]) + self.list[j].addfield(field,rawlist[i][1]) # Deal with #i @@ -115,5 +115,5 @@ self.list[int(plotnum)-1].addfield(field,rawlist[i][1]) else: - + #go through all subplots and assign key-value pairs j=0 Index: /issm/trunk-jpl/src/m/plot/applyoptions.py =================================================================== --- /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 23562) +++ /issm/trunk-jpl/src/m/plot/applyoptions.py (revision 23563) @@ -120,5 +120,5 @@ for label in ax.get_xticklabels() + ax.get_yticklabels(): label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) - if int(md.mesh.dimension)==3: + if int(md.mesh.dimension)==3: for label in ax.get_zticklabels(): label.set_fontsize(options.getfieldvalue('ticklabelfontsize')) @@ -142,5 +142,5 @@ # }}} # {{{ xlim, ylim, zlim - if options.exist('xlim'): + if options.exist('xlim'): ax.set_xlim(options.getfieldvalue('xlim')) if options.exist('ylim'): @@ -164,7 +164,7 @@ options.addfielddefault('clim',lims) else: - if len(data)>0: + if len(data)>0: lims=[data.min(),data.max()] - else: + else: lims=[0,1] # }}} @@ -231,17 +231,16 @@ cb.solids.set_rasterized(True) cb.update_ticks() - cb.set_alpha(1) 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) + 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 + # {{{ expdisp if options.exist('expdisp'): expdisp(ax,options) @@ -291,3 +290,2 @@ # {{{ inset TODO # }}} - Index: /issm/trunk-jpl/src/m/plot/plot_manager.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_manager.py (revision 23562) +++ /issm/trunk-jpl/src/m/plot/plot_manager.py (revision 23563) @@ -96,5 +96,5 @@ #apply all options applyoptions(md,data2,options,fig,axgrid,gridindex) - + #ground overlay on kml plot_unit Index: /issm/trunk-jpl/src/m/plot/plot_overlay.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_overlay.py (revision 23562) +++ /issm/trunk-jpl/src/m/plot/plot_overlay.py (revision 23563) @@ -5,13 +5,11 @@ import matplotlib.pyplot as plt import matplotlib as mpl +import os try: - from mpl_toolkits.basemap import Basemap + from mpl_toolkits.basemap import Basemap except ImportError: - print 'Basemap toolkit not installed' - -import os - + print 'Basemap toolkit not installed' try: - from osgeo import gdal + from osgeo import gdal except ImportError: print 'osgeo/gdal for python not installed, plot_overlay is disabled' @@ -25,12 +23,11 @@ x,y,z,elements,is2d,isplanet=processmesh(md,[],options) - - if data=='none' or data==None: + try: + data,datatype=processdata(md,data,options) + imageonly=0 + except (TypeError,ValueError):#that should catch None and 'none' but may also catch unwanted errors imageonly=1 data=np.float('nan')*np.ones((md.mesh.numberofvertices,)) datatype=1 - else: - imageonly=0 - data,datatype=processdata(md,data,options) if not is2d: @@ -50,5 +47,4 @@ ymin=trans[3]+gtif.RasterYSize*trans[5] ymax=trans[3] - # allow supplied image to have limits smaller than basemap or model limits x0=max(min(xlim),xmin) @@ -58,9 +54,9 @@ inputname='temp.tif' os.system('gdal_translate -quiet -projwin ' + str(x0) + ' ' + str(y1) + ' ' + str(x1) + ' ' + str(y0) + ' ' + image+ ' ' + inputname) - + gtif=gdal.Open(inputname) arr=gtif.ReadAsArray() #os.system('rm -rf ./temp.tif') - + if gtif.RasterCount>=3: # RGB array r=gtif.GetRasterBand(1).ReadAsArray() @@ -84,7 +80,7 @@ plt.hist(arr.flatten(),bins=256,range=(0.,1.)) plt.title('histogram of overlay image, use for setting overlaylims') - plt.show() + plt.show() plt.sca(ax) # return to original axes/figure - + # get parameters from cropped geotiff trans=gtif.GetGeoTransform() @@ -94,6 +90,6 @@ ymax=trans[3] dx=trans[1] - dy=trans[5] - + dy=trans[5] + xarr=np.arange(xmin,xmax,dx) yarr=np.arange(ymin,ymax,-dy) # -dy since origin='upper' (not sure how robust this is) @@ -102,31 +98,35 @@ norm=mpl.colors.Normalize(vmin=overlaylims[0],vmax=overlaylims[1]) + pc=ax.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm) + if options.exist('basemap'): # create coordinate grid in map projection units (for plotting) - lat,lon=xy2ll(xlim,ylim,-1,0,71) - #plt.sca(ax) - width=xmax-xmin - height=ymax-ymin - lat_0,lon_0=xy2ll(xmin+width/2.,ymin+height/2.,-1,0,71) - m=Basemap(projection='spstere', - llcrnrlon=lon[0],llcrnrlat=lat[0],urcrnrlon=lon[1],urcrnrlat=lat[1], - epsg=3031, - resolution='c') - #width=width,height=height,lon_0=lon_0,lat_0=lat_0, - #lat_0=-90,lon_0=0,lat_ts=-71, - #llcrnrx=x0,llcrnry=y0,urcrnrx=x1,urcrnry=y1) - #test - #m.ax=ax - meridians=np.arange(-180.,181.,1.) - parallels=np.arange(-80.,80.,1.) - m.drawparallels(parallels,labels=[0,0,1,1]) # labels=[left,right,top,bottom] - m.drawmeridians(meridians,labels=[1,1,0,0]) - m.drawcoastlines() - pc=m.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm, ax=ax) + if md.mesh.epsg==3413: + hemisphere=1 + st_lat=70 + lon_0=45 + elif md.mesh.epsg==3031: + hemisphere=-1 + st_lat=71 + lon_0=0 + else: + hemisphere=raw_input('epsg code {} is not supported chose your hemisphere (1 for North, -1 for south)'.format(mesh.epsg)) - else: - pc=ax.pcolormesh(xg, yg, np.flipud(arr), cmap=mpl.cm.Greys, norm=norm) - - #rasterization? + lat,lon=xy2ll(xlim,ylim,hemisphere,lon_0,st_lat) + extent=[np.diff(xlim)[0],np.diff(ylim)[0]] + center=[lon[0]+np.diff(lon)[0]*0.5,lat[0]+np.diff(lat)[0]*0.5] + m=Basemap(llcrnrlon=lon[0],llcrnrlat=lat[0],urcrnrlon=lon[1],urcrnrlat=lat[1], + lon_0=center[0],lat_0=center[1],#width=extent[0],height=extent[1],# + epsg=md.mesh.epsg,anchor='NW', + resolution='i',ax=ax) + + meridians=np.arange(np.floor(lon[0]),np.ceil(lon[1]),1.) + parallels=np.arange(np.floor(lat[0]),np.ceil(lat[1]),1.) + m.drawparallels(parallels,labels=[1,0,0,0],ax=ax) # labels=[left,right,top,bottom] + m.drawmeridians(meridians,labels=[0,0,1,0],ax=ax) + m.drawcoastlines(ax=ax) + m.drawmapboundary(ax=ax) + + #rasterization? if options.getfieldvalue('rasterized',0): pc.set_rasterized(True) Index: /issm/trunk-jpl/src/m/plot/plot_streamlines.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_streamlines.py (revision 23562) +++ /issm/trunk-jpl/src/m/plot/plot_streamlines.py (revision 23563) @@ -12,5 +12,5 @@ def plot_streamlines(md,options,ax): ''' - plot streamlines on a figure, using by default vx and vy components in md.initialization. + plot streamlines on a figure, using by default vx and vy components in md.initialization. Usage: @@ -25,5 +25,5 @@ streamlineswidthscale: scaling multiplier for linewidth scaled by velocity streamlinesarrowsize: size of arrows on lines (default 1) - + ''' @@ -36,5 +36,5 @@ arrowsize=options.getfieldvalue('streamlinesarrowsize',1) - #process mesh and data + #process mesh and data x,y,z,elements,is2d,isplanet=processmesh(md,vx,options) u,datatype=processdata(md,vx,options) Index: /issm/trunk-jpl/src/m/plot/plot_unit.py =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 23562) +++ /issm/trunk-jpl/src/m/plot/plot_unit.py (revision 23563) @@ -11,9 +11,9 @@ except ImportError: 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) @@ -36,6 +36,10 @@ colorlevels=options.getfieldvalue('colorlevels',128) alpha=options.getfieldvalue('alpha',1) - # }}} - # {{{ define wich colormap to use + if alpha<1: + antialiased=True + else: + antialiased=False + # }}} + # {{{ define wich colormap to use try: defaultmap=plt.cm.get_cmap('viridis',colorlevels) @@ -51,5 +55,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: @@ -83,7 +87,8 @@ options.addfield('colornorm',norm) # }}} - + # Plot depending on the datatype # {{{ data are on elements + if datatype==1: if is2d: @@ -141,6 +146,8 @@ #raise ValueError('plot_unit error: 3D element plot not supported yet') return + # }}} # {{{ data are on nodes + elif datatype==2: if is2d: @@ -152,7 +159,11 @@ #tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha) if options.exist('log'): - tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha) + 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: - tri=ax.tricontourf(triangles,data,colorlevels,cmap=cmap,norm=norm,alpha=alpha,extend='both') + 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) @@ -162,5 +173,5 @@ 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 @@ -177,5 +188,5 @@ pl3.set_color(color) ax.add_collection3d(pl3) - + #deal with prism faces triface=np.vstack((elements[:,0:3],elements[:,3:6])) @@ -191,5 +202,5 @@ pl3.set_color(color) ax.add_collection3d(pl3) - + ax.set_xlim([min(x),max(x)]) ax.set_ylim([min(y),max(y)]) @@ -197,4 +208,5 @@ #raise ValueError('plot_unit error: 3D element plot not supported yet') return + # }}} # {{{ plotting quiver @@ -223,3 +235,2 @@ else: raise ValueError('datatype=%d not supported' % datatype) -