Changeset 21444

Timestamp:

12/15/16 03:14:02 (8 years ago)

Author:

bdef

Message:

NEW: adding 3D ability for standard data plot (may need to shift to mayavi on a later stage)

Location:

issm/trunk-jpl/src/m/plot

Files:

• applyoptions.py (modified) (2 diffs)
• plot_manager.py (modified) (1 diff)
• plot_mesh.py (modified) (5 diffs)
• plot_unit.py (modified) (5 diffs)
• plotdoc.py (modified) (3 diffs)
• plotmodel.py (modified) (1 diff)

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

@@ -10 +10 @@
         from mpl_toolkits.mplot3d import Axes3D
         import matplotlib as mpl
-        import pylab as p
         import matplotlib.pyplot as plt
 except ImportError:
         print "could not import pylab, matplotlib has not been installed, no plotting capabilities enabled"
 
-def applyoptions(md,data,options,fig,ax):
+def applyoptions(md,data,options,fig,axgrid,gridindex):
         '''
         APPLYOPTIONS - apply options to current plot
@@ -30 +29 @@
         # get handle to current figure and axes instance
         #fig = p.gcf()
-        #ax=p.gca()
+        ax=     axgrid[gridindex]
 
         # {{{ font

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

@@ -93 +93 @@
         data2,datatype=processdata(md,data,options)
         #plot unit
-        plot_unit(x,y,z,elements,data2,is2d,isplanet,datatype,options,ax)
+        plot_unit(x,y,z,elements,data2,is2d,isplanet,datatype,options,fig,axgrid,gridindex)
         #apply all options
-        applyoptions(md,data2,options,fig,ax)
+        applyoptions(md,data2,options,fig,axgrid,gridindex)
         
         #ground overlay on kml plot_unit

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

@@ -9 +9 @@
 from matplotlib.patches import Polygon
 from mpl_toolkits.mplot3d.art3d import Line3DCollection
+from mpl_toolkits.axes_grid1 import inset_locator
+from mpl_toolkits.mplot3d import Axes3D
 def plot_mesh(md,options,fig,axgrid,gridindex):
         '''
@@ -25 +27 @@
         if is2d:
                 ax.triplot(x,y,elements)
+
         else:
-                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')
-
+                # 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')
+                ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D)
+                
+                
                 AB=elements[:,0:2]
                 BC=elements[:,1:3]
@@ -39 +44 @@
                 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
@@ -50 +55 @@
                         pl3=Line3DCollection([tri],edgecolor='r')
                         ax.add_collection3d(pl3)
-
+                        
                 ax.set_xlim([min(x),max(x)])
                 ax.set_ylim([min(y),max(y)])
@@ -58 +63 @@
         options.addfielddefault('colorbar','off')
         options.addfielddefault('ticklabels','on')
-        applyoptions(md,[],options,fig,ax)
+        applyoptions(md,[],options,fig,axgrid,gridindex)

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

@@ -1 +1 @@
 from cmaptools import truncate_colormap
 from plot_quiver import plot_quiver
+import numpy as  np
 try:
-        import pylab as p
         import matplotlib as mpl
         import matplotlib.pyplot as plt
-        import numpy as  np
+        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"
         
-def plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options,ax):
+def plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options,fig,axgrid,gridindex):
         """
         PLOT_UNIT - unit plot, display data
@@ -18 +20 @@
         See also: PLOTMODEL, PLOT_MANAGER
         """
-        
+        #if we are plotting 3d replace the current axis
+        if not is2d:
+                ax=inset_locator.inset_axes(axgrid[gridindex],width='100%',height='100%',loc=3,borderpad=0,axes_class=Axes3D)
+        else:
+                ax=axgrid[gridindex]
+
         #edgecolor
         edgecolor=options.getfieldvalue('edgecolor','None')
@@ -42 +49 @@
         options.addfield('colormap',cmap)
         # }}}   
-        # {{{ if plotting only one of several layers reduce dataset
+        # {{{ if plotting only one of several layers reduce dataset, same for surface
         if options.getfieldvalue('layer',0)>=1:
-                plotlayer=options.getfieldvalue('layer',1)
+                plotlayer=options.getfieldvalue('layer',0)
                 if datatype==1:
                         slicesize=np.shape(elements)[0]
@@ -84 +91 @@
                         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')
+                        #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(recface).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
+                        _, idx = np.unique(tmp, return_index=True)
+                        recel= recface[idx]
+                        recindex=eltind[idx]
+                        for i,rectangle in enumerate(recel):
+                                rec=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 = np.unique(tmp, return_index=True)
+                        triel= triface[idx]
+                        triindex=eltind[idx]
+                        for i,triangle in enumerate(triel):
+                                tri=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 
         # }}}
@@ -99 +150 @@
                                 ax.triplot(x,y,elements,color=edgecolor)
                 else:
-                        raise ValueError('plot_unit error: 3D node plot not supported yet')
+                        #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(recface).view(np.dtype((np.void, recface.dtype.itemsize * recface.shape[1])))
+                        _, idx = np.unique(tmp, return_index=True)
+                        recel= recface[idx]                                     
+                        for rectangle in recel:
+                                rec=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 = np.unique(tmp, return_index=True)
+                        triel= triface[idx]
+                        for triangle in triel:
+                                tri=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
         # }}}

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

@@ -8 +8 @@
                                         'mesh':' draw mesh using trisurf',
                                         'BC':' this will draw all the boundary conditions (Dirichlet and Neumann).',
-                                        'elementnumbering':' numbering of elements (matlab indices)'}
+                                        'elementnumbering':' numbering of elements (matlab indices)',
+                                        '3D disclaimer':'3D is implemented with plot3d for now this is not optimal and may change to mayavi at some point. The impelementation is on the development side for now so expect some issue and question your plotting before you results.'}
         TODOdata={'basal_drag':' plot the basal drag on the bed (in kPa) based on the velocity in md.initialization',
                                                 'basal_dragx or basal_dragy' :' plot a component of the basal drag on the bed (in kPa)',
@@ -18 +19 @@
                                                 'driving_stress':' plot the driving stress (in kPa)',
                                                 'elements_type':' model used for each element',
-                                                
                                                 'highlightvertices':' to highlight vertices (use highlight option to enter the vertex list',
                                                 'referential':' stressbalance referential',
@@ -47 +47 @@
         
         pyoptions={'axis':" show ('on') or hide ('off') axes",
-                                                 'caxis':" modify  colorbar range. (array of type [a b] where b>=a)",
+                                                 'caxis':" modify  colorbar range. (array of type [a b] where b>=a)",
                                                  'colorlevels':" N, number of levels to use",
                                                  'colorbar':" add colorbar (string 'on','off' or 'one')",
                                                  'axes_pad':" spacing between axes (default is 0.25)",
                                                  'colorbartitle':" colorbar title (string)",
-                                                 'colorbarticks':" set colorbar ticks manually (list)",
-                                                 'colorbarfontsize':" specify colorbar fontsize",
+                                                 'colorbarticks':" set colorbar ticks manually (list)",
+                                                 'colorbarfontsize':" specify colorbar fontsize",
                                                  'colormap':" change the default colormap ('viridis' is the default)",
-                                                 'contourlevels':" N or [value1,...] add the contours of the specified values or N contours",
+                                                 'contourlevels':" N or [value1,...] add the contours of the specified values or N contours",
                                                  'streamlines':" TOFIX argument does nothing",
                                                  'edgecolor':" color of mesh edges. RGB tuple or standard string",

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

@@ -58 +58 @@
                 if options.list[0].exist('figsize'):
                         figsize=options.list[0].getfieldvalue('figsize')
-                        fig=plt.figure(figurenumber,figsize=(figsize[0],figsize[1]),tight_layout=True)
+                        fig=plt.figure(figurenumber,figsize=(figsize[0],figsize[1]))#,tight_layout=True)
                 else:
-                        fig=plt.figure(figurenumber,tight_layout=True)
+                        fig=plt.figure(figurenumber)#,tight_layout=True)
                 fig.clf()