Index: /issm/trunk-jpl/src/m/plot/applyoptions.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 27197)
+++ /issm/trunk-jpl/src/m/plot/applyoptions.py	(revision 27198)
@@ -34,5 +34,5 @@
     fontsize = options.getfieldvalue('fontsize', 8)
     fontweight = options.getfieldvalue('fontweight', 'normal')
-    fontfamily = options.getfieldvalue('fontfamily', 'sans - serif')
+    fontfamily = options.getfieldvalue('fontfamily', 'sans-serif')
     font = {
         'fontsize': fontsize,
@@ -186,18 +186,10 @@
     if options.exist('cmap_set_under'):
         cbar_extend += 2
-    # }}}
-    # {{{ contours
-    if options.exist('contourlevels'):
-        plot_contour(md, data, options, ax)
-    # }}}
-    # {{{ edgeoverlay
-    if options.exist('edgeoverlay'):
-        edgedata = options.getfieldvalue('edgeoverlay')
-        plot_edgeoverlay(md, edgedata, options, ax)
-    # }}}
-    # {{{ wrapping TODO
-    # }}}
-    # {{{ colorbar
-    if options.getfieldvalue('colorbar', 1) == 1:
+
+    # }}}
+    # {{{ colorbar extension
+    if options.exist('cbar_extend'):
+        extend = options.getfieldvalue('cbar_extend', 'neither')
+    else:
         if cbar_extend == 0:
             extend = 'neither'
@@ -208,5 +200,19 @@
         elif cbar_extend == 3:
             extend = 'both'
-
+        options.addfielddefault('cbar_extend', extend)
+    # }}}
+    # {{{ contours
+    if options.exist('contourlevels'):
+        plot_contour(md, data, options, ax)
+    # }}}
+    # {{{ edgeoverlay
+    if options.exist('edgeoverlay'):
+        edgedata = options.getfieldvalue('edgeoverlay')
+        plot_edgeoverlay(md, edgedata, options, ax)
+    # }}}
+    # {{{ wrapping TODO
+    # }}}
+    # {{{ colorbar
+    if options.getfieldvalue('colorbar', 1) == 1:
         cb = mpl.colorbar.ColorbarBase(ax.cax, cmap=cmap, norm=norm, extend=extend)
         if options.exist('alpha'):
Index: /issm/trunk-jpl/src/m/plot/plot_edgeoverlay.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_edgeoverlay.py	(revision 27197)
+++ /issm/trunk-jpl/src/m/plot/plot_edgeoverlay.py	(revision 27198)
@@ -3,4 +3,6 @@
 import matplotlib.pyplot as plt
 from processmesh import processmesh
+from matplotlib import collections
+from scipy.stats import percentileofscore
 
 
@@ -13,32 +15,119 @@
     See also: PLOTMODEL'''
 
+    # if md.mesh.numberofedges not in np.shape(datain):
+    #     raise ValueError('Data must be defined on edges to be ploted as an edge overlay')
+
     x, y, z, elements, is2d, isplanet = processmesh(md, [], options)
+    Edges = md.mesh.edges - 1
 
-    flags = datain > 6.7e-5  # this is appropriate for channel Area (perhaps)
-    hide = np.invert(flags)
+    #First we mask values under a given value define by quantiles
+    if options.exist('edgemin'):
+        minval = options.getfieldvalue('edgemin', 0)
+        minquant = percentileofscore(datain, minval, kind='weak')
+        minquant = minquant / 100
+    else:
+        minquant = 0.7
+        minval = np.quantile(datain, minquant)
+    print("For the overlay we plot values above {:.4g} wich corresponds to the {}% percentile".format(minval, minquant * 100))
 
-    NodeMask = np.zeros(np.shape(md.mesh.x), dtype=bool)
-    HideElt = md.mesh.edges[np.where(hide), 0] - 1
-    NodeMask[HideElt] = True
-    MaskX = np.ma.array(x, mask=NodeMask)
-    MaskY = np.ma.array(y, mask=NodeMask)
+    flags = datain > minval  #6.7e-5  # this is appropriate for channel Area (perhaps)
 
-    EdgeEnd = md.mesh.edges[:, 1] - 1
-    EdgeStart = md.mesh.edges[:, 0] - 1
-    quiverU = MaskX[EdgeEnd] - MaskX[EdgeStart]
-    quiverV = MaskY[EdgeEnd] - MaskY[EdgeStart]
+    edgetype = options.getfieldvalue('edgetype', 'thickness')
 
-    Masked = np.ma.masked_array(datain, mask=hide)
-    if options.exist('cedgelim'):
-        lims = options.getfieldvalue('cedgelim', [Masked.min(), Masked.max()])
-        edgenorm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
+    if edgetype == "color":
+        #create an nodewise dataset from edges
+        NodeMask = np.ones(np.shape(md.mesh.x), dtype=bool)
+        #We grab all the starts from the unmasked edges
+        Starters = Edges[np.where(flags), 0]
+        NodeMask[Starters] = False
+        Xstart = np.ma.array(x, mask=NodeMask)
+        Ystart = np.ma.array(y, mask=NodeMask)
+
+        NodeMask = np.ones(np.shape(md.mesh.x), dtype=bool)
+        Enders = Edges[np.where(flags), 1]
+        NodeMask[Enders] = False
+        Xend = np.ma.array(x, mask=NodeMask)
+        Yend = np.ma.array(y, mask=NodeMask)
+
+        #define vectors from the start and end point of the unmasked edges
+        EdgeEnd = Edges[:, 1]
+        EdgeStart = Edges[:, 0]
+        quiverU = np.ma.array(Xend[EdgeEnd] - Xstart[EdgeStart], mask=np.invert(flags))
+        quiverV = np.ma.array(Yend[EdgeEnd] - Ystart[EdgeStart], mask=np.invert(flags))
+
+        #mask out the values from the data and create colorscale
+        Masked = np.ma.masked_array(datain, mask=np.invert(flags))
+        cbar_extend = 0
+        edgemap = plt.cm.get_cmap('inferno')
+        if options.exist('cedgeaxis'):
+            lims = options.getfieldvalue('cedgeaxis', [Masked.min(), Masked.max()])
+            edgenorm = mpl.colors.Normalize(vmin=lims[0], vmax=lims[1])
+            if lims[0] > Masked.min():
+                edgemap.set_under('r')
+                cbar_extend += 2
+            if lims[1] < Masked.max():
+                edgemap.set_over('k')
+                cbar_extend += 1
+        else:
+            edgenorm = mpl.colors.Normalize(vmin=Masked.min(), vmax=Masked.max())
+
+        if cbar_extend == 0:
+            extend = 'neither'
+        elif cbar_extend == 1:
+            extend = 'max'
+        elif cbar_extend == 2:
+            extend = 'min'
+        elif cbar_extend == 3:
+            extend = 'both'
+
+        ax.quiver(Xstart[EdgeStart], Ystart[EdgeStart], quiverU, quiverV, datain,
+                  units="xy", angles="xy", scale_units="xy", scale=1,
+                  headwidth=0, headlength=0, width=100, headaxislength=0,
+                  norm=edgenorm, cmap=edgemap)
+
+        #plt.colorbar(plt.cm.ScalarMappable(norm=edgenorm, cmap=edgemap), ax=ax, extend=extend, orientation='horizontal', anchor=(1, 0))
+        cbarax = ax.inset_axes([0.02, 0.02, 0.96, 0.05])
+        #inset_axes(ax, width="100%", height="5%", loc='lower center', borderpad=-5)
+        mpl.colorbar.ColorbarBase(cbarax, norm=edgenorm, cmap=edgemap, extend=extend, orientation='horizontal', ticklocation='top')
+
+    elif edgetype == 'thickness':
+        #First we classify a range
+        if options.exist('edgeranges'):
+            edgeranges = options.getfieldvalue('edgeranges', 2)
+        else:
+            edgeranges = 2
+
+        quantrange = np.linspace(minquant, 1, edgeranges + 1)[:-1]
+        flags = []
+        legtext = []
+        for Qindex, quantile in enumerate(quantrange):
+            if quantile < quantrange[-1]:
+                lowquant = np.quantile(datain, quantile)
+                highquant = np.quantile(datain, quantrange[Qindex + 1])
+                flags.append(np.logical_and(datain >= lowquant, datain <= highquant))
+                legtext.append('From  {:.2g} to {:.2g}'. format(lowquant, highquant))
+            else:
+                lowquant = np.quantile(datain, quantile)
+                flags.append((datain > lowquant))
+                legtext.append('More than {:.2g}'.format(lowquant))
+
+        flags = np.asarray(np.squeeze(flags))
+        EdgeEnd = Edges[:, 1]
+        EdgeStart = Edges[:, 0]
+
+        for index in range(edgeranges):
+            linecol = []
+            #We loop on the result and aonly add to a linecollection the flaged edges.
+            for datindex, datval in enumerate(datain):
+                if flags[index, datindex]:
+                    linecol.append([(x[EdgeStart[datindex]], y[EdgeStart[datindex]]),
+                                    (x[EdgeEnd[datindex]], y[EdgeEnd[datindex]])])
+
+            lc = collections.LineCollection(linecol, color='k', linewidth=0.5 + index,
+                                            label=legtext[index])
+            ax.add_collection(lc)
+
+        ax.legend(title='Overlay')
+
     else:
-        edgenorm = mpl.colors.Normalize(vmin=Masked.min(), vmax=Masked.max())
-    edgemap = plt.cm.get_cmap('inferno')
-
-    ax.quiver(MaskX[EdgeStart], MaskY[EdgeStart], quiverU, quiverV, datain,
-              units="xy", angles="xy", scale_units="xy", scale=1,
-              headwidth=0, headlength=0, width=100, headaxislength=0,
-              norm=edgenorm, cmap=edgemap)
-
-    plt.colorbar(plt.cm.ScalarMappable(norm=edgenorm, cmap=edgemap), ax=ax)
+        print("Only 'color', and 'thickness' are accepted as edgeoverlay types")
Index: /issm/trunk-jpl/src/m/plot/plot_manager.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 27197)
+++ /issm/trunk-jpl/src/m/plot/plot_manager.py	(revision 27198)
@@ -27,5 +27,5 @@
     currently generated using matplotlib's AxesGrid toolkit.
 
-    'gridindex' is passed through to specialized plot* functions and 
+    'gridindex' is passed through to specialized plot* functions and
     applyoptions.
 
@@ -89,5 +89,4 @@
     # }}}
 
-    #process data and model
     x, y, z, elements, is2d, isplanet = processmesh(md, data, options)
     data2, datatype = processdata(md, data, options)
Index: /issm/trunk-jpl/src/m/plot/plot_unit.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 27197)
+++ /issm/trunk-jpl/src/m/plot/plot_unit.py	(revision 27198)
@@ -33,13 +33,9 @@
 
     # colormap
-    # {{{ give number of colorlevels and transparency
+    # 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
+    # }}}
+    # define wich colormap to use {{{
     try:
         defaultmap = plt.cm.get_cmap('viridis', colorlevels)
@@ -51,13 +47,7 @@
     else:
         cmap = getcolormap(options)
-    if options.exist('cmap_set_over'):
-        over = options.getfieldvalue('cmap_set_over', 'k')
-        cmap.set_over(over)
-    if options.exist('cmap_set_under'):
-        under = options.getfieldvalue('cmap_set_under', 'k')
-        cmap.set_under(under)
     options.addfield('colormap', cmap)
     # }}}
-    # {{{ if plotting only one of several layers reduce dataset, same for surface
+    # if plotting only one of several layers reduce dataset, same for surface {{{
     if options.getfieldvalue('layer', 0) >= 1:
         plotlayer = options.getfieldvalue('layer', 0)
@@ -68,5 +58,5 @@
         data = data[(plotlayer - 1) * slicesize:plotlayer * slicesize]
     # }}}
-    # {{{ Get the colormap limits
+    # Get the colormap limits {{{
     dataspread = np.nanmax(data) - np.nanmin(data)
     if dataspread != 0.:
@@ -77,4 +67,8 @@
     if options.exist('caxis'):
         lims = options.getfieldvalue('caxis', [np.nanmin(data), np.nanmax(data)])
+        if lims[0] > np.nanmin(data):
+            options.addfielddefault('cmap_set_under', 'r')
+        if lims[1] < np.nanmax(data):
+            options.addfielddefault('cmap_set_over', 'k')
     else:
         if limextent == 0.:
@@ -85,6 +79,28 @@
         else:
             lims = [np.nanmin(data), np.nanmax(data)]
-    # }}}
-    # {{{ Set the spread of the colormap (default is normal
+
+    cbar_extend = 0
+    if options.exist('cmap_set_over'):
+        over = options.getfieldvalue('cmap_set_over', 'k')
+        cmap.set_over(over)
+        cbar_extend += 1
+    if options.exist('cmap_set_under'):
+        under = options.getfieldvalue('cmap_set_under', 'r')
+        cmap.set_under(under)
+        cbar_extend += 2
+    # }}}
+
+    # colorbar extension {{{
+    if cbar_extend == 0:
+        extend = 'neither'
+    elif cbar_extend == 1:
+        extend = 'max'
+    elif cbar_extend == 2:
+        extend = 'min'
+    elif cbar_extend == 3:
+        extend = 'both'
+    options.addfielddefault('cbar_extend', extend)
+    # }}}
+    # Set the spread of the colormap (default is normal) {{{
     if options.exist('log'):
         norm = mpl.colors.LogNorm(vmin=lims[0], vmax=lims[1])
@@ -95,5 +111,5 @@
 
     # Plot depending on the datatype
-    # {{{ data are on elements
+    # data are on elements {{{
     if datatype == 1:
         if is2d:
@@ -102,4 +118,5 @@
             else:
                 triangles = mpl.tri.Triangulation(x, y, elements)
+
             tri = ax.tripcolor(triangles, data, colorlevels, cmap=cmap, norm=norm, alpha=alpha, edgecolors=edgecolor)
         else:
@@ -153,6 +170,5 @@
 
     # }}}
-    # {{{ data are on nodes
-
+    # data are on nodes {{{
     elif datatype == 2:
         if is2d:
@@ -162,15 +178,8 @@
             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)
+            if edgecolor == 'None':
+                tri = ax.tripcolor(triangles, data, cmap=cmap, norm=norm, alpha=alpha, shading='gouraud')
             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)
+                tri = ax.tripcolor(triangles, data, cmap=cmap, norm=norm, alpha=alpha, edgecolors=edgecolor)
         else:
             #first deal with the colormap
@@ -178,5 +187,4 @@
             loccmap.set_array(lims)
             loccmap.set_clim(vmin=lims[0], vmax=lims[1])
-
     #deal with prism sides
             recface = np.vstack((elements[:, 0], elements[:, 1], elements[:, 4], elements[:, 3])).T
@@ -215,5 +223,5 @@
 
     # }}}
-    # {{{ plotting quiver
+    # plotting quiver {{{
     elif datatype == 3:
         if is2d:
@@ -223,14 +231,24 @@
         return
     # }}}
-    # {{{ plotting P1 Patch (TODO)
+    # plotting P1 Patch (TODO) {{{
     elif datatype == 4:
         print('plot_unit message: P1 patch plot not implemented yet')
         return
     # }}}
-    # {{{ plotting P0 Patch (TODO)
+    # plotting P0 Patch (TODO) {{{
     elif datatype == 5:
         print('plot_unit message: P0 patch plot not implemented yet')
         return
     # }}}
+    # plotting edges  {{{
+    elif datatype == 6:
+        if is2d:
+            triangles = mpl.tri.Triangulation(x, y, elements)
+            tri = ax.tripcolor(triangles, data, cmap=cmap, norm=norm, alpha=alpha, edgecolors=edgecolor)
+
+        else:
+            print("edge plotting is not implemented for 3D")
+        return
+    # }}}
     else:
         raise ValueError('datatype = %d not supported' % datatype)
Index: /issm/trunk-jpl/src/m/plot/plotmodel.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 27197)
+++ /issm/trunk-jpl/src/m/plot/plotmodel.py	(revision 27198)
@@ -28,33 +28,27 @@
     numberofplots = options.numberofplots
 
-    #get number of subplots
-    subplotwidth = ceil(sqrt(numberofplots))
+    #get the "optimal" number of subfigures in a row/col
+    if (np.nanmax(md.mesh.x) - np.nanmin(md.mesh.x)) > (np.nanmax(md.mesh.y) - np.nanmin(md.mesh.y)):
+        maxrow = ceil(sqrt(numberofplots))
+        maxcol = ceil(numberofplots / maxrow)
+    else:
+        maxcol = ceil(sqrt(numberofplots))
+        maxrow = ceil(numberofplots / maxcol)
 
-    # TODO: Check that commenting this out is correct; we do not need a hold
-    # under matplotlib, right?
-    #
-    # #get hold
-    # hold = options.list[0].getfieldvalue('hold', False)
-
-    #if nrows and ncols specified, then bypass
+    #If any  of nrows or ncols are given we use that
     if options.list[0].exist('nrows'):
         nrows = options.list[0].getfieldvalue('nrows')
-        nr = True
+        if options.list[0].exist('ncols'):
+            ncols = options.list[0].getfieldvalue('ncols')
+        else:
+            ncols = ceil(numberofplots / nrows)
+    elif options.list[0].exist('ncols'):
+        ncols = options.list[0].getfieldvalue('ncols')
+        nrows = ceil(numberofplots / ncols)
     else:
-        nrows = np.ceil(numberofplots / subplotwidth)
-        nr = False
-
-    if options.list[0].exist('ncols'):
-        ncols = options.list[0].getfieldvalue('ncols')
-        nc = True
-    else:
-        ncols = int(subplotwidth)
-        nc = False
+        nrows = maxrow
+        ncols = maxcol
     ncols = int(ncols)
     nrows = int(nrows)
-
-    #check that nrows and ncols were given at the same time!
-    if nr != nc:
-        raise Exception('plotmodel error message: nrows and ncols need to be specified together, or not at all')
 
     # Go through plots
Index: /issm/trunk-jpl/src/m/plot/processdata.py
===================================================================
--- /issm/trunk-jpl/src/m/plot/processdata.py	(revision 27197)
+++ /issm/trunk-jpl/src/m/plot/processdata.py	(revision 27198)
@@ -8,5 +8,7 @@
     datatype = 2 -> nodes
     datatype = 3 -> node quivers
-    datatype = 4 -> patch
+    datatype = 4 -> P1 patch
+    datatype = 5 -> P0 patch
+    datatype = 6 -> edges
 
     Usage:
@@ -24,4 +26,9 @@
     else:
         numberofvertices2d = np.nan
+
+    try:
+        numberofedges = md.mesh.numberofedges
+    except AttributeError:
+        numberofedges = np.nan
 
     if options.exist('amr'):
@@ -52,5 +59,5 @@
     # }}}
 
-    #  log {{{
+    # log {{{
     if options.exist('log'):
         cutoff = options.getfieldvalue('log', 1)
@@ -58,5 +65,5 @@
     # }}}
 
-    #  quiver plot {{{
+    # quiver plot {{{
     if datasize[1] > 1 and datasize[0] != numberofvertices + 1:
         if datasize[0] == numberofvertices and datasize[1] in [2, 3]:
@@ -80,5 +87,5 @@
     # }}}
 
-    #  element data{{{
+    # element data{{{
     if datasize[0] == numberofelements and datasize[1] == 1:
         #initialize datatype if non patch
@@ -107,6 +114,6 @@
     # }}}
 
-    #  node data {{{
-    if datasize[0] in [numberofvertices, numberofvertices2d] and datasize[1] == 1:
+    # node data {{{
+    elif datasize[0] in [numberofvertices, numberofvertices2d] and datasize[1] == 1:
         datatype = 2
         # AMR {{{
@@ -115,5 +122,5 @@
             procdata = procdata[nonan]
         # }}}
-        #  Mask {{{
+        # Mask {{{
         if options.exist('mask'):
             flags = options.getfieldvalue('mask')
@@ -133,6 +140,19 @@
     # }}}
 
-    #  spc time series {{{
-    if datasize[0] == numberofvertices + 1:
+    # edge data {{{
+    elif datasize[0] in [numberofedges] and datasize[1] == 1:
+        datatype = 6
+        procdata = np.zeros((md.mesh.numberofelements))
+        repeat = np.zeros((md.mesh.numberofelements))
+        for index, edge in enumerate(md.mesh.edges[:, -2:]):
+            procdata[edge - 1] += data[index] * np.asarray(edge - 1 > -1, dtype=int)
+            repeat[edge - 1] += np.asarray(edge - 1 > -1, dtype=int)
+        procdata = procdata / repeat
+
+        # }}}
+    # }}}
+
+    # spc time series {{{
+    elif datasize[0] == numberofvertices + 1:
         datatype = 2
         spccol = options.getfieldvalue('spccol', 0)
@@ -141,12 +161,7 @@
         procdata = procdata[0:-1, spccol]
 
-        #mask?
-
-        #layer projection?
-
-        #control arrow density if quiver plot
     # }}}
 
-    # convert rank - 2 array to rank - 1 {{{
+    # convert rank -2 array to rank -1 {{{
     if np.ndim(procdata) == 2 and np.shape(procdata)[1] == 1:
         procdata = procdata.reshape(-1, )
@@ -169,7 +184,6 @@
         procdata[np.isnan(procdata)] = nanfill
         procdata = np.ma.array(procdata, mask=np.isnan(procdata))
-        print('from nan processing {} and {}'.format(lb, ub))
         if nanfill < lb:
-            options.addfielddefault('cmap_set_under', 'k')
+            options.addfielddefault('cmap_set_under', 'r')
         elif nanfill > ub:
             options.addfielddefault('cmap_set_over', 'k')