Changeset 23563

Timestamp:

12/20/18 08:43:50 (6 years ago)

Author:

bdef

Message:

CHG: updating overlay for plotmodel

Location:

issm/trunk-jpl/src/m

Files:

• classes/plotoptions.py (modified) (4 diffs)
• plot/applyoptions.py (modified) (5 diffs)
• plot/plot_manager.py (modified) (1 diff)
• plot/plot_overlay.py (modified) (7 diffs)
• plot/plot_streamlines.py (modified) (3 diffs)
• plot/plot_unit.py (modified) (11 diffs)

issm/trunk-jpl/src/m/classes/plotoptions.py

@@ -52 +52 @@
                                 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 +76 @@
                         #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 +106 @@
                                                         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 +115 @@
                                                 self.list[int(plotnum)-1].addfield(field,rawlist[i][1])
                         else:
-                                
+
                                 #go through all subplots and assign key-value pairs
                                 j=0

issm/trunk-jpl/src/m/plot/applyoptions.py

@@ -120 +120 @@
                 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 +142 @@
         # }}}
         # {{{ xlim, ylim, zlim
-        if options.exist('xlim'):
+        if options.exist('xlim'):
                 ax.set_xlim(options.getfieldvalue('xlim'))
         if options.exist('ylim'):
@@ -164 +164 @@
                 options.addfielddefault('clim',lims)
         else:
-                if len(data)>0: 
+                if len(data)>0:
                         lims=[data.min(),data.max()]
-                else: 
+                else:
                         lims=[0,1]
         # }}}
@@ -231 +231 @@
                 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 +290 @@
         # {{{ inset TODO
         # }}}
-        

issm/trunk-jpl/src/m/plot/plot_manager.py

@@ -96 +96 @@
         #apply all options
         applyoptions(md,data2,options,fig,axgrid,gridindex)
-        
+
         #ground overlay on kml plot_unit
 

issm/trunk-jpl/src/m/plot/plot_overlay.py

@@ -5 +5 @@
 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 +23 @@
 
         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 +47 @@
         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 +54 @@
         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 +80 @@
                 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 +90 @@
         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 +98 @@
         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)

issm/trunk-jpl/src/m/plot/plot_streamlines.py

@@ -12 +12 @@
 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 +25 @@
         streamlineswidthscale: scaling multiplier for linewidth scaled by velocity
         streamlinesarrowsize: size of arrows on lines (default 1)
-        
+
     '''
 
@@ -36 +36 @@
     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)

issm/trunk-jpl/src/m/plot/plot_unit.py

@@ -11 +11 @@
 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 +36 @@
         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 +55 @@
                 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 +87 @@
         options.addfield('colornorm',norm)
         # }}}
-        
+
         # Plot depending on the datatype
         # {{{ data are on elements
+
         if datatype==1:
                 if is2d:
@@ -141 +146 @@
                         #raise ValueError('plot_unit error: 3D element plot not supported yet')
                 return 
+
         # }}}
         # {{{ data are on nodes
+
         elif datatype==2:
                 if is2d:
@@ -152 +159 @@
                                 #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 +173 @@
                         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 +188 @@
                                 pl3.set_color(color)
                                 ax.add_collection3d(pl3)
-                                
+
                         #deal with prism faces
                         triface=np.vstack((elements[:,0:3],elements[:,3:6]))
@@ -191 +202 @@
                                 pl3.set_color(color)
                                 ax.add_collection3d(pl3)
-                                
+
                         ax.set_xlim([min(x),max(x)])
                         ax.set_ylim([min(y),max(y)])
@@ -197 +208 @@
                         #raise ValueError('plot_unit error: 3D element plot not supported yet')
                 return
+
         # }}}
         # {{{ plotting quiver
@@ -223 +235 @@
         else:
                 raise ValueError('datatype=%d not supported' % datatype)
-