Changeset 21440

Timestamp:

12/14/16 03:36:13 (8 years ago)

Author:

bdef

Message:

NEW: adding 3D treatment for mesh and clean way to remove the colorbar

Location:

issm/trunk-jpl/src/m/plot

Files:

• plot_manager.py (modified) (3 diffs)
• plot_mesh.py (modified) (2 diffs)
• plotmodel.py (modified) (4 diffs)
• processmesh.py (modified) (1 diff)

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

@@ -24 +24 @@
         from plot_overlay import plot_overlay
 
-def plot_manager(md,options,fig,ax):
+def plot_manager(md,options,fig,axgrid,gridindex):
         '''
         PLOT_MANAGER - distribute the plots called by plotmodel
@@ -46 +46 @@
         options.addfielddefault('ticklabels','on')
 
+        ax=axgrid[gridindex]
         # {{{ basemap plot TOFIX
         #if options.exist('basemap'):
@@ -60 +61 @@
         if isinstance(data,(str,unicode)):
                 if data=='mesh': 
-                        plot_mesh(md,options,fig,ax)
+                        plot_mesh(md,options,fig,axgrid,gridindex)
+
                         #fig.delaxes(fig.axes[1]) # hack to remove colorbar after the fact
                         return

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

@@ -3 +3 @@
 except ImportError:
         print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled"
-
+import numpy as np
+import matplotlib as mpl
 from processmesh import processmesh
 from applyoptions import applyoptions
-
-def plot_mesh(md,options,fig,ax):
+from matplotlib.patches import Polygon
+from mpl_toolkits.mplot3d.art3d import Line3DCollection
+def plot_mesh(md,options,fig,axgrid,gridindex):
         '''
         PLOT_MESH - plot model mesh
@@ -16 +18 @@
                 See also: PLOTMODEL
         '''
+        x,y,z,elements,is2d,isplanet=processmesh(md,'mesh',options)
 
-        x,y,z,elements,is2d,isplanet=processmesh(md,[],options)
-
+        ax=axgrid[gridindex]
+        fig.delaxes(axgrid.cbar_axes[gridindex])
+        
         if is2d:
                 ax.triplot(x,y,elements)
         else:
-                print 'WARNING: only 2D mesh plot is currently implemented'
-        
+                fig.delaxes(ax)
+                geom=int(str(axgrid.get_geometry()[0])+str(axgrid.get_geometry()[1])+str(gridindex+1))
+                ax = fig.add_subplot(geom,projection='3d')
+
+                AB=elements[:,0:2]
+                BC=elements[:,1:3]
+                CA=np.vstack((elements[:,2],elements[:,0])).T
+                DE=elements[:,3:5]
+                EF=elements[:,4:]
+                FD=np.vstack((elements[:,5],elements[:,3])).T
+                AD=np.vstack((elements[:,0],elements[:,3])).T
+                BE=np.vstack((elements[:,1],elements[:,4])).T
+                CF=np.vstack((elements[:,2],elements[:,5])).T
+
+                tmpa=np.vstack((AB,BC,CA,DE,EF,FD,AD,BE,CF))
+                #deleting segments that appear multiple times
+                tmpb = np.ascontiguousarray(tmpa).view(np.dtype((np.void, tmpa.dtype.itemsize * tmpa.shape[1])))
+                _, idx = np.unique(tmpb, return_index=True)
+                triel= tmpa[idx]
+                
+                for triangle in triel:
+                        tri=zip(x[triangle],y[triangle],z[triangle])
+                        pl3=Line3DCollection([tri],edgecolor='r')
+                        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)])
         #apply options
         options.addfielddefault('title','Mesh')

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

@@ -2 +2 @@
 from plotoptions import plotoptions
 from plotdoc import plotdoc
+from plot_manager import plot_manager
+from math import ceil, sqrt
 
 try:
         import pylab as p
         import matplotlib.pyplot as plt
-        from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid
+        from mpl_toolkits.axes_grid1 import ImageGrid, AxesGrid
+        from mpl_toolkits.mplot3d import Axes3D
 except ImportError:
         print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled"
-
-from plot_manager import plot_manager
-from math import ceil, sqrt
 
 def plotmodel(md,*args):
@@ -83 +83 @@
                 cbar_pad=options.list[0].getfieldvalue('colorbarpad','2.5%') # None or %
 
-                axgrid=ImageGrid(fig, 111,
+                axgrid=ImageGrid(fig,111,
                                 nrows_ncols=(nrows,ncols),
                                 ngrids=plotnum,
@@ -94 +94 @@
                                 cbar_location=cbar_location,
                                 cbar_size=cbar_size,
-                                cbar_pad=cbar_pad
-                                )
+                                cbar_pad=cbar_pad)
 
                 if cbar_mode=='None':
@@ -101 +100 @@
                                 fig._axstack.remove(ax)
 
-                for i in xrange(numberofplots):
-                        plot_manager(options.list[i].getfieldvalue('model',md),options.list[i],fig,axgrid[i])
-
+                for i,ax in enumerate(axgrid.axes_all):
+                        plot_manager(options.list[i].getfieldvalue('model',md),options.list[i],fig,axgrid,i)
                 fig.show()
         else:

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

@@ -15 +15 @@
         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')
+        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
-                        try:
-                                x=md.mesh.x2d
-                        except AttributeError:
-                                x=md.mesh.x
-                        try:
-                                y=md.mesh.y2d
-                        except AttributeError:                          
-                                y=md.mesh.y
-                else:
-                        x=md.mesh.long
-                        y=md.mesh.lat
-                try:
-                        z=md.mesh.z
-                except AttributeError:
-                        z=np.zeros_like(md.mesh.x)
+  # {{{ treating coordinates
+
+        try:
+                z=md.mesh.z
+        except AttributeError:
+                z=np.zeros(np.shape(md.mesh.x))
+        elements=md.mesh.elements-1
+        
+        if options.getfieldvalue('layer',0)>=1:
+                x=md.mesh.x2d
+                y=md.mesh.y2d
+                z=np.zeros(np.shape(x))
+                elements=md.mesh.elements2d-1
+        elif 'latlon' in options.getfieldvalue('coord','xy'):
+                x=md.mesh.long
+                y=md.mesh.lat
+        else:
+                x=md.mesh.x
+                y=md.mesh.y
+
+        #is it a 2D plot?
+        if md.mesh.dimension()==2 or options.getfieldvalue('layer',0)>=1:
+                is2d=1
+        else:
+                is2d=0
                 
-                try:
-                        elements=md.mesh.elements2d-1
-                except AttributeError:
-                        elements=md.mesh.elements-1
-
-                #is it a 2D plot?
-                if md.mesh.dimension()==2 or options.getfieldvalue('layer',0)>=1:
-                        is2d=1
-                else:
-                        is2d=0
-
-                #layer projection?
-                if options.getfieldvalue('layer',0)>=1:
-                         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
-                         z=np.zeros(np.size(x))
-                         elements=elements
-        else:
-                #Process mesh for plotting 
-                if md.mesh.dimension()==2:
-                        is2d=1
-                else:
-                        # process polycollection here for 3D plot
-                        is2d=0
         #units
         if options.exist('unit'):