def print_matrix(row, col, row_label, col_label, matrix_1, matrix_2, number, same_row_col):
    f, ax = plt.subplots(2, sharex = True)
    ax[0] = plt.subplot2grid((2,2), (0,0), colspan = 2)
    ax[1] = plt.subplot2grid((2,2), (1,0), colspan = 2)
        
    # Assigns color according to values in fitness_heat arrays
    norm = MidpointNormalize(midpoint=0)
    #im1 = ax[0].imshow(matrix_1, cmap=plt.cm.seismic, interpolation='none')
    im1 = ax[0].imshow(matrix_1, norm=norm, cmap=plt.cm.seismic, interpolation='none')
    divider1 = make_axes_locatable(ax[0])
    cax1 = divider1.append_axes("right", size="2%", pad=0.05)

    im2 = ax[1].imshow(matrix_2, cmap=plt.cm.seismic, interpolation='none')
    #im2 = ax[1].imshow(matrix_2, norm=norm, cmap=plt.cm.seismic, interpolation='none')
    divider2 = make_axes_locatable(ax[1])
    cax2 = divider2.append_axes("right", size="2%", pad=0.05)

    # Sorting out labels for y-axis in correct order and placing in new list 'yax'
    yax1 = sorted(number.items(), key=lambda (k, v): v)
    yax = []
    for item in yax1:
        if item[0] =='STOP':
            yax.append ('*')
        else:
            yax.append (item[0])
    
    if same_row_col == True:
        xax = yax
        x = np.arange(0, row+0.5, 1)
        ax[0].set_xticks(x)
        ax[0].set_xticklabels(xax, rotation='horizontal', fontsize = 10)
        ax[0].set_xlabel(row_label)
        ax[1].set_xticks(x)
        ax[1].set_xticklabels(xax, rotation='horizontal', fontsize = 10)
        ax[1].set_xlabel(row_label)
    else:
        x = np.arange(0.5, col+0.5, 5)
        xlab = np.arange(1, col, 5)
        plt.xticks(x, xlab, rotation='horizontal')
        ax[1].set_xlabel(col_label)
        
    # Setting up y-axis ticks...want AA labels, centered
    y = np.arange(0, row, 1)
    ax[0].set_yticks(y)
    ax[0].set_yticklabels(yax, rotation='horizontal', fontsize = 10)
    ax[0].set_ylabel(row_label)
    ax[1].set_yticks(y)
    ax[1].set_yticklabels(yax, rotation='horizontal', fontsize = 10)
    ax[1].set_ylabel(row_label)

    # Setting up x-axis ticks...want residue number, centered
    

    # Limit axes to appropriate values
    plt.axis([0, col, 0, row])

    plt.colorbar(im1, cax = cax1)
    plt.colorbar(im2, cax = cax2)

    plt.show()

def display_data(data):
    from mpl_toolkits.axes_grid1 import make_axes_locatable

    fig, (trans_ax, f_1_ax, f_2_ax) = plt.subplots(ncols=1, nrows=3, sharex=True, constrained_layout=True, figsize=(6, 4))
    trans_im = trans_ax.imshow(data.trans, extent=data.extent)
    trans_divider = make_axes_locatable(trans_ax)
    trans_ax_cb = trans_divider.new_horizontal(size="3%", pad=0.05)
    trans_fig = trans_ax.get_figure()
    trans_fig.add_axes(trans_ax_cb)
    trans_cb = plt.colorbar(trans_im, cax=trans_ax_cb)
    trans_cb.ax.set_ylabel(r'$P_{trans}$', rotation=0, labelpad=25, size=15)

    f_1_im = f_1_ax.imshow(data.f_1.T, extent=data.extent)
    f_1_divider = make_axes_locatable(f_1_ax)
    f_1_ax_cb = f_1_divider.new_horizontal(size="3%", pad=0.05)
    f_1_fig = f_1_ax.get_figure()
    f_1_fig.add_axes(f_1_ax_cb)
    f_1_cb = plt.colorbar(f_1_im, cax=f_1_ax_cb)
    f_1_cb.ax.set_ylabel(r'$P_{f1}$', rotation=0, labelpad=15, size=15)

    f_2_im = f_2_ax.imshow(data.f_2.T, extent=data.extent)
    f_2_divider = make_axes_locatable(f_2_ax)
    f_2_ax_cb = f_2_divider.new_horizontal(size="3%", pad=0.05)
    f_2_fig = f_2_ax.get_figure()
    f_2_fig.add_axes(f_2_ax_cb)
    f_2_cb = plt.colorbar(f_2_im, cax=f_2_ax_cb)
    f_2_cb.ax.set_ylabel(r'$P_{f2}$', rotation=0, labelpad=13, size=15)

    f_1_ax.set_ylabel(r'Neutron Angle $\phi$ (rad)')
    f_2_ax.set_xlabel(r'Detector Angle $\theta$ (rad)')
    fig.suptitle('Data Sinograms')
    plt.subplots_adjust(top=0.92, bottom=0.125, left=0.100, right=0.9, hspace=0.2, wspace=0.2)

def plot_thsections(dbz,vvel,ht,dayt,**kwargs):

	fig,ax = plt.subplots(2,1,sharex=True)

	'colormap for vvel'
	orig_cmap = cm.bwr
	shifted_cmap = shiftedColorMap(orig_cmap, midpoint=0.7, name='shifted')

	' add images and colorbar'
	im0=ax[0].imshow(dbz,interpolation='none',cmap='nipy_spectral',vmin=-20,vmax=60,aspect='auto',origin='lower')
	im1=ax[1].imshow(vvel,interpolation='none',cmap=shifted_cmap,vmin=-10,vmax=4,aspect='auto',origin='lower')
	divider0 = make_axes_locatable(ax[0])
	cax0 = divider0.append_axes("right", size="2%", pad=0.05)
	cbar0 = plt.colorbar(im0, cax=cax0)
	divider1 = make_axes_locatable(ax[1])
	cax1 = divider1.append_axes("right", size="2%", pad=0.05)
	cbar1 = plt.colorbar(im1, cax=cax1)

	if 'echotop' in kwargs:
		echot=kwargs['echotop']
		ax[0].plot(echot,color='k')

	format_yaxis(ax[0],ht)
	format_yaxis(ax[1],ht)

	format_xaxis(ax[1], dayt, minutes_tick=30, labels=True)
	ax[1].invert_xaxis()

	ax[0].set_ylabel('Hgt MSL [km]')
	ax[1].set_ylabel('Hgt MSL [km]')
	ax[1].set_xlabel(r'$\Leftarrow$'+' Time [UTC]')

	plt.subplots_adjust(hspace=0.05)
	plt.suptitle('SPROF observations. Date: '+dayt[0].strftime('%Y-%b'))
	plt.draw()

def display_images(mu_im, mu_f_im, p_im, title='Images'):
    from mpl_toolkits.axes_grid1 import make_axes_locatable

    fig, (mu_ax, mu_f_ax, p_ax) = plt.subplots(ncols=1, nrows=3, sharex=True, constrained_layout=True, figsize=(5, 6))

    mu_im = mu_ax.imshow(mu_im.data, extent=mu_im.extent)
    mu_divider = make_axes_locatable(mu_ax)
    mu_ax_cb = mu_divider.new_horizontal(size="5%", pad=0.05)
    mu_fig = mu_ax.get_figure()
    mu_fig.add_axes(mu_ax_cb)
    mu_cb = plt.colorbar(mu_im, cax=mu_ax_cb)
    mu_cb.ax.set_ylabel(r'$\mu$', rotation=0, labelpad=10, size=15)

    mu_f_im = mu_f_ax.imshow(mu_f_im.data, extent=mu_f_im.extent)
    mu_f_divider = make_axes_locatable(mu_f_ax)
    mu_f_ax_cb = mu_f_divider.new_horizontal(size="5%", pad=0.05)
    mu_f_fig = mu_f_ax.get_figure()
    mu_f_fig.add_axes(mu_f_ax_cb)
    mu_f_cb = plt.colorbar(mu_f_im, cax=mu_f_ax_cb)
    mu_f_cb.ax.set_ylabel(r'$\frac{\mu_f}{\mu}$', rotation=0, labelpad=10, size=15)

    p_im = p_ax.imshow(p_im.data, extent=p_im.extent)
    p_divider = make_axes_locatable(p_ax)
    p_ax_cb = p_divider.new_horizontal(size="5%", pad=0.05)
    p_fig = p_ax.get_figure()
    p_fig.add_axes(p_ax_cb)
    p_cb = plt.colorbar(p_im, cax=p_ax_cb)
    p_cb.ax.set_ylabel(r'$p$', rotation=0, labelpad=10, size=15)

    mu_f_ax.set_ylabel('Y (cm)')
    p_ax.set_xlabel('X (cm)')
    fig.suptitle(title)
    plt.show()

def test_extendednorm():

    a = np.zeros((4, 5))
    a[0,0] = -9999
    a[1,1] = 1.1
    a[2,2] = 2.2
    a[2,4] = 1.9
    a[3,3] = 9999999

    cm = mpl.cm.get_cmap('jet')
    bounds = [0,1,2,3]
    norm = cleo.colors.ExtendedNorm(bounds, cm.N, extend='both')

    #fig, (ax1, ax2) = plt.subplots(1, 2)
    fig = plt.figure()
    ax1 = fig.add_subplot(1, 2, 1)
    ax2 = fig.add_subplot(1, 2, 2)
    imax = ax1.imshow(a, interpolation='None', norm=norm, cmap=cm,
                     origin='lower');
    divider = make_axes_locatable(ax1)
    cax = divider.append_axes("right", size="5%", pad=0.2)
    plt.colorbar(imax, cax=cax, extend='both')

    ti = cm(norm(a))
    ax2.imshow(ti, interpolation='None', origin='lower')
    divider = make_axes_locatable(ax2)
    cax = divider.append_axes("right", size="5%", pad=0.2)
    cbar = mpl.colorbar.ColorbarBase(cax, extend='both', cmap=cm,
                                     norm=norm)
    fig.tight_layout()

    def show(self):
        
        nclasses = len(self._landclass)
        gs=plt.GridSpec(nclasses, 2)
        
        row = 0
        fig=plt.figure()
        for c in self._landclass:
            ax=fig.add_subplot(gs[row,0], title=c)
            im=ax.imshow(self._landclass[c].get_raster())
            
            from mpl_toolkits.axes_grid1 import make_axes_locatable
            divider = make_axes_locatable(plt.gca())
            cax = divider.append_axes("right", "5%", pad="3%")            
            cb=fig.colorbar(im,cax=cax)
            
            ax=fig.add_subplot(gs[row,1],title='Classified '+c)
            im=ax.imshow(self._landclass[c].get_classraster())
            
            
            divider = make_axes_locatable(plt.gca())
            cax = divider.append_axes("right", "5%", pad="3%")            
            cb=fig.colorbar(im,cax=cax)
            
            cb.set_ticks( list(range(1,self._landclass[c].get_nclasses()+1)))
            cb.set_ticklabels(self._landclass[c].get_classes_str())            

            row +=1
        
        fig.tight_layout()

def plot_residuals(directory):
    residuals = []
    for directory in _get_cc_dirs(directory):
        residual_file = directory + os.sep + 'results/residuals.dat'
        residuals.append(np.loadtxt(residual_file))
    all_residuals = np.array(residuals)
    nr_f = all_residuals.shape[1] / 2
    fig, axes = plt.subplots(1, 2, figsize=(6, 2.5))

    ax = axes[0]
    im = ax.imshow(all_residuals[:, 0:nr_f], interpolation='none')
    divider = make_axes_locatable(ax)
    ax_cb = divider.new_horizontal(size="5%", pad=0.05)
    fig.add_axes(ax_cb)
    fig.colorbar(im, cax=ax_cb)
    ax.set_xlabel('frequency')
    ax.set_title(r'Magnitude $\Delta(log(R))$')

    ax = axes[1]
    im = ax.imshow(all_residuals[:, nr_f:], interpolation='none')
    divider = make_axes_locatable(ax)
    ax_cb = divider.new_horizontal(size="5%", pad=0.05)
    fig.add_axes(ax_cb)
    fig.colorbar(im, cax=ax_cb)
    ax.set_xlabel('frequency')
    ax.set_title(r'Phase $\Delta \phi$')
    fig.tight_layout()
    fig.savefig('residuals.png', dpi=300)

def createColorbar(patches, cMin=None, cMax=None, nLevs=5,
                   label=None, orientation='horizontal', *args, **kwargs):
    cbarTarget = plt
    cax = None
    divider = None
    #if hasattr(patches, 'figure'):
        #cbarTarget = patches.figure

    #print( patches)

    if hasattr(patches, 'ax'):
        divider = make_axes_locatable(patches.ax)
    elif hasattr(patches, 'get_axes'):
        divider = make_axes_locatable(patches.get_axes())
    
    if divider:
        if orientation == 'horizontal':
            cax = divider.append_axes("bottom", size=0.25, pad=0.65)
        else:
            cax = divider.append_axes("right", size=0.25, pad=0.05)
            #cax = divider.append_axes("right", size="5%", pad=0.05)
            #cbar3 = plt.colorbar(im3, cax=cax3)
        
    #print(patches,  cax)
    cbar = cbarTarget.colorbar(patches, cax=cax,
                               orientation=orientation)

    setCbarLevels(cbar, cMin, cMax, nLevs)

    if label is not None:
        cbar.set_label(label)

    return cbar

def compare(data, model, res, chi, filename):

    head = np.percentile(data, 99.9)
    toe = np.percentile(data, 0.1)
    chirange = 5. # head*0.05 # np.percentile(data, 9.3)
    resrange = (head-toe)*0.5 # np.percentile(data, 50.0)

    width, height = np.shape(data)
    figw, figh = (width/dippy)*2., (height/dippy)*2.

    fig = plt.figure(figsize=(figw, figh))
#    fig = plt.figure(dpi=dippy, figsize=(figw, figh))

    gs = gridspec.GridSpec(2, 2)

#    axr = fig.add_subplot(221)
    axr = plt.subplot(gs[0])
    axr.set_axis_off()
    axr.set_title(r"Data", fontsize=40)

    caxr = axr.imshow(data, cmap=plt.cm.gray, vmin=toe, vmax=head)
    divr = make_axes_locatable(plt.gca())
    cbarr = divr.append_axes("right", "5%", pad="3%")
    fig.colorbar(caxr, cax=cbarr)

#    axf = fig.add_subplot(222)
    axf = plt.subplot(gs[1])
    axf.set_axis_off()
    axf.set_title(r"Model Image", fontsize=40)

    caxf = axf.imshow(model, cmap=plt.cm.gray, vmin=toe, vmax=head)
    divf = make_axes_locatable(plt.gca())
    cbarf = divf.append_axes("right", "5%", pad="3%")
    fig.colorbar(caxf, cax=cbarf)

#    axres = fig.add_subplot(223)
    axres = plt.subplot(gs[2])
    axres.set_axis_off()
    axres.set_title(r"Residual Image",fontsize=40)

    caxres = axres.imshow(res, cmap=plt.cm.gray, vmin=-resrange, vmax=resrange)
    divres = make_axes_locatable(plt.gca())
    cbarres = divres.append_axes("right", "5%", pad="3%")
    fig.colorbar(caxres, cax=cbarres)

#    axchi = fig.add_subplot(224)
    axchi = plt.subplot(gs[3])
    axchi.set_axis_off()
    axchi.set_title(r"Chi Image", fontsize=40)

    caxchi = axchi.imshow(chi, cmap=plt.cm.gray, vmin=-chirange, vmax=chirange)
    divchi = make_axes_locatable(plt.gca())
    cbarchi = divchi.append_axes("right", "5%", pad="3%")
    fig.colorbar(caxchi, cax=cbarchi)

#    fig.subplots_adjust(hspace=0.001)
    fig.tight_layout()

#    fig.savefig(filename, dpi=60., bbox_inches='tight')
    fig.savefig(filename)

def compare_to_model(srcs, img, fig=None, axarr=None): 
    """ plots true image, model image, and difference (much like above) 
        Input: 
            srcs: python list of PointSrcParams
            img : FitsImage object

        Output:
            fig, axarr
    """
    if fig is None or axarr is None:
        fig, axarr = plt.subplots(1, 3)

    # generate model image
    model_img = gen_model_image(srcs, img)
    vmin = min(img.nelec.min(), model_img.min())
    vmax = max(img.nelec.max(), model_img.max())

    im1 = axarr[0].imshow(np.log(img.nelec), interpolation='none', origin='lower', vmin=np.log(vmin), vmax=np.log(vmax))
    axarr[0].set_title('log data ($\log(x_{n,m})$')
    im2 = axarr[1].imshow(np.log(model_img), interpolation='none', origin='lower', vmin=np.log(vmin), vmax=np.log(vmax))
    axarr[1].set_title('log model ($\log(F_{n,m})$)')
    divider2 = make_axes_locatable(axarr[1])
    cax2 = divider2.append_axes('right', size='10%', pad=.05)
    cbar2 = fig.colorbar(im2, cax=cax2)

    im3 = axarr[2].imshow(img.nelec-model_img, interpolation='none', origin='lower')
    axarr[2].set_title('Difference: Data - Model')
    divider3 = make_axes_locatable(axarr[2])
    cax3 = divider3.append_axes('right', size='10%', pad=.05)
    cbar3 = fig.colorbar(im3, cax=cax3)
    return fig, axarr

def plotarray(zonearray,Z,dirname,runid):
    
    """
    Plots HY properties
    """
    
    assert(len(Z.shape)==2)
    assert(len(zonearray.shape)==2)
    
    fig =plt.figure(figsize=(11,8.5))
    fig.subplots_adjust(hspace=0.45, wspace=0.3)
    
    """Compute the log of HY """
    nrows,ncols=Z.shape
    logZ =np.log10(Z)
    logZ1d=np.reshape(logZ, (nrows*ncols,))
   
  
    """Histogram of log HY """  
    ax3=fig.add_subplot(2,1,1)
    mybins=np.arange(-3.0,4.0,0.333333)
    ax3.hist(logZ1d, bins=mybins, normed=0, facecolor='green', alpha=0.75)
    ax3.set_xlabel('Log10 HY')
    ax3.set_ylabel('Frequency')
    ax3.set_ylim(0,30000)
    ax3.grid(True)    
    
    """Plot of HY Zone Array """
    ax1=fig.add_subplot(2,2,3)
    im1 =ax1.imshow(zonearray,interpolation='bilinear')
    """
    create an axes on the right side of ax1. The width of cax will be 5%
    of ax and the padding between cax1 and ax1 will be fixed at 0.05 inch.
    """
    divider= make_axes_locatable(ax1)
    cax1=divider.append_axes("right",size="5%",pad=0.05)
    cbar1=plt.colorbar(im1, cax=cax1,ticks=[0,1])
    cbar1.ax.set_yticklabels(['Low','High']) #Vertically Oriented Colorbar
    ax1.set_title('HYZones'+"{0:04d}".format(runid))
       
    """Plot of HY Array"""
    ax =fig.add_subplot(2,2,4)
    im =ax.imshow(logZ,interpolation='bilinear',vmin=-2,vmax=2)
    """
    create an axes on the right side of ax. The width of cax will be 5%
    of ax and the padding between cax and ax will be fixed at 0.05 inch.
    """
    divider= make_axes_locatable(ax)
    cax=divider.append_axes("right",size="5%",pad=0.05)    
    #Add colorbar, make sure to specify tick locations to match desired ticklabels
    cbar=plt.colorbar(im, cax=cax,ticks=[-2,-1,0,1,2])
    cbar.ax.set_yticklabels(['0.01','0.1','1','10','>100']) #Vertically Oriented Colorbar
    ax.set_title('HYField'+"{0:04d}".format(runid))
    
    """Save the figure """
    fout =dirname +"/HY" + "{0:04d}".format(runid)+".png"
    plt.tight_layout()
    plt.savefig(fout)
    plt.clf()

def plot_class(dist_m, shape_port, dat_port, dat_star, dat_class, ft, humfile, sonpath, base, p):

   if len(shape_port)>2:
      Zdist = dist_m[shape_port[-1]*p:shape_port[-1]*(p+1)]
      extent = shape_port[1]
   else:
      Zdist = dist_m
      extent = shape_port[0]

   print("Plotting ... ")
   # create fig 1
   fig = plt.figure()
   fig.subplots_adjust(wspace = 0, hspace=0.075)
   plt.subplot(2,1,1)
   ax = plt.gca()
   im = ax.imshow(np.vstack((np.flipud(dat_port),dat_star)),cmap='gray',
                  extent=[min(Zdist), max(Zdist), -extent*(1/ft), extent*(1/ft)],origin='upper')
   plt.ylabel('Horizontal distance (m)'); 
   plt.axis('tight')

   plt.tick_params(\
   axis='x',          # changes apply to the x-axis
   which='both',      # both major and minor ticks are affected
   bottom='off',      # ticks along the bottom edge are off
   labelbottom='off') # labels along the bottom edge are off

   try:
      divider = make_axes_locatable(ax)
      cax = divider.append_axes("right", size="5%", pad=0.05)
      plt.colorbar(im, cax=cax)
   except:
      plt.colorbar()

   plt.subplot(2,1,2)
   ax = plt.gca()
   plt.imshow(np.vstack((np.flipud(dat_port), dat_star)),cmap='gray',
              extent=[min(Zdist), max(Zdist), -extent*(1/ft), extent*(1/ft)],origin='upper')
   im = ax.imshow(dat_class, alpha=0.5,extent=[min(Zdist), max(Zdist), -extent*(1/ft), extent*(1/ft)],
                  origin='upper', cmap='YlOrRd', vmin=0, vmax=3)
   plt.ylabel('Horizontal distance (m)'); 
   plt.xlabel('Distance along track (m)')
   plt.axis('tight')

   try:
      divider = make_axes_locatable(ax)
      cax = divider.append_axes("right", size="5%", pad=0.05)
      plt.colorbar(im, cax=cax, extend='max')
   except:
      plt.colorbar(extend='max')

   if len(shape_port)>2:
      custom_save(sonpath,base+'class'+str(p))
   else:
      custom_save(sonpath,base+'class')
   del fig

    def plot_surface(self, save_name=None, show_bicoh=True,
                     cmap='viridis', contour_color='k'):
        '''
        Plot the bispectrum amplitude and (optionally) the bicoherence.

        Parameters
        ----------
        save_name : str, optional
            Save name for the figure. Enables saving the plot.
        show_bicoh : bool, optional
            Plot the bicoherence surface. Enabled by default.
        cmap : {str, matplotlib color map}, optional
            Colormap to use in the plots. Default is viridis.
        contour_color : {str, RGB tuple}, optional
            Color of the amplitude contours.
        '''

        import matplotlib.pyplot as plt
        from mpl_toolkits.axes_grid1 import make_axes_locatable

        if show_bicoh:
            ax1 = plt.subplot(1, 2, 1)
        else:
            ax1 = plt.subplot(1, 1, 1)
        im1 = ax1.imshow(self.bispectrum_logamp, origin="lower",
                         interpolation="nearest", cmap=cmap)
        divider = make_axes_locatable(ax1)
        cax = divider.append_axes("right", size="5%", pad=0.05)
        cbar1 = plt.colorbar(im1, cax=cax)
        cbar1.set_label(r"log$_{10}$ Bispectrum Amplitude")
        ax1.contour(self.bispectrum_logamp,
                    colors=contour_color)
        ax1.set_xlabel(r"$k_1$")
        ax1.set_ylabel(r"$k_2$")

        if show_bicoh:
            ax2 = plt.subplot(1, 2, 2)
            im2 = ax2.imshow(self.bicoherence, origin="lower",
                             interpolation="nearest",
                             cmap=cmap)
            divider = make_axes_locatable(ax2)
            cax2 = divider.append_axes("right", size="5%", pad=0.05)
            cbar2 = plt.colorbar(im2, cax=cax2)
            cbar2.set_label("Bicoherence")
            ax2.set_xlabel(r"$k_1$")
            ax2.set_ylabel(r"$k_2$")

        plt.tight_layout()

        if save_name is not None:
            plt.savefig(save_name)
            plt.close()
        else:
            plt.show()

def make_axr(ax, overlaycolor='darkred', size='23%'):
    axr = ax.twinx()
    div = make_axes_locatable(axr)
    div.append_axes('right', size=size, add_to_figure=False)
    ax.spines['right'].set_color(overlaycolor)
    axr.tick_params(axis='y', colors=overlaycolor)
    axr.yaxis.label.set_color(overlaycolor)
    axr.yaxis.get_offset_text().set_color(overlaycolor)

    div = make_axes_locatable(ax)
    div.append_axes('right', size=size, add_to_figure=False)

    return axr

  def __init__(self, current_surf_dem, glacier_mask, bed_dem, date, \
    output_trace_files, temp_files_path, glacier_thickness_threshold):
    self.output_trace_files = output_trace_files
    self.temp_files_path = temp_files_path
    self.glacier_thickness_threshold = glacier_thickness_threshold

    self.fig = plt.figure(figsize=(16,12))
    self.surf_dem_plot_title = 'Surface DEM ' + date
    self.sub1 = self.fig.add_subplot(131)
    self.sub1.set_title(self.surf_dem_plot_title)
    self.sub1.set_xticks([])
    self.sub1.set_yticks([])
    self.sub1.get_axes().set_frame_on(True)
    self.img1 = self.sub1.imshow(current_surf_dem)

    self.sub1.invert_yaxis() # makes North up
    self.divider1 = make_axes_locatable(plt.gca())
    self.cax1 = self.divider1.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(self.img1, cax=self.cax1)

    self.glacier_mask_plot_title = 'Glacier Mask ' + date
    self.sub2 = self.fig.add_subplot(132)
    self.sub2.set_title(self.glacier_mask_plot_title)
    self.sub2.set_xticks([])
    self.sub2.set_yticks([])
    self.sub2.get_axes().set_frame_on(True)
    self.img2 = self.sub2.imshow(glacier_mask)
    self.sub2.invert_yaxis()
    self.divider2 = make_axes_locatable(plt.gca())
    self.cax2 = self.divider2.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(self.img2, cax=self.cax2)

    self.glacier_thickness_plot_title = 'Glacier Thickness ' + date
    self.sub3 = self.fig.add_subplot(133)
    self.sub3.set_title(self.glacier_thickness_plot_title)
    self.sub3.set_xticks([])
    self.sub3.set_yticks([])
    self.sub3.get_axes().set_frame_on(True)
    self.img3 = self.sub3.imshow( \
      current_surf_dem - glacier_thickness_threshold - bed_dem)
    self.sub3.invert_yaxis()
    self.divider3 = make_axes_locatable(plt.gca())
    self.cax3 = self.divider3.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(self.img3, cax=self.cax3)

    self.fig.tight_layout()
    plt.show(block=False)

    if self.output_trace_files:
      self.fig.savefig(temp_files_path + 'dem_+_glacier_mask_+_thickness_' + \
        date, dpi=self.fig.dpi)

def plotStressStrainFibSet(fiberSet,title,fileName=None,nContours=100,pointSize=50, fiberShape='o'):
  '''Represents graphically the cross-section current stresses and strains.
  The graphics are generated by a triangulation from the x,y coordinates of
  the fibers.

  :param fiberSet: set of fibers to be represented
  :param title:    general title for the graphic
  :param fileName: name of the graphic file (defaults to None: no file generated)
  :param nContours: number of contours to be generated (defaults to 100). If 
         nContours=0 or nContours=None, then each fiber is represented by a
         makred.
  :param pointSize: size of the circles to represent each of the fibers in the 
         set in the case that nContours=0 or nContours=None (defaults to 50)
  :param fiberShape: marker to represent each fiber, in case nContours = 0 or 
                 None.e.g.: "o"->circle, "s"->square, "p"->pentagon (defaults
                 to circle)
  '''
  lsYcoo=list()
  lsZcoo=list()
  lsStrain=list()
  lsStress=list()
  for f in fiberSet:
    fpos=f.getPos()
    lsYcoo.append(fpos.x)
    lsZcoo.append(fpos.y)
    lsStrain.append(f.getStrain())
    lsStress.append(f.getForce()/f.getArea()/1e6)
    
  fig,(ax1,ax2) = plt.subplots(1,2, sharex=True, sharey=True)
  fig.suptitle(title, fontsize=20)
  ax1.set_title('Strain')
  if nContours in [None,0]:
    im1=ax1.scatter(lsYcoo,lsZcoo,s=pointSize,c=lsStrain,marker=fiberShape)
    im2=ax2.scatter(lsYcoo,lsZcoo,s=pointSize,c=lsStress,marker=fiberShape)
  else:
    im1=ax1.tricontourf(lsYcoo,lsZcoo,lsStrain, nContours)
    im2=ax2.tricontourf(lsYcoo,lsZcoo,lsStress, nContours)
  divider1 = make_axes_locatable(ax1)
  cax1 = divider1.append_axes("right", size="20%", pad=0.05)
  cbar1 = plt.colorbar(im1,cax=cax1)
  ax2.set_title('Stress')
#  im2=ax2.tricontourf(lsYcoo,lsZcoo,lsStress, nContours)
#  im2=ax2.scatter(lsYcoo,lsZcoo,50,lsStress)
  divider2 = make_axes_locatable(ax2)
  cax2 = divider2.append_axes("right", size="20%", pad=0.05)
  cbar2 = plt.colorbar(im2,cax=cax2)
  if fileName<>None:
    plt.savefig(fileName)
  plt.show()
  return

def createColorBar(gci, orientation='horizontal', size=0.2, pad=None,
                   **kwargs):
    """Create a Colorbar.

    Shortcut to create a matplotlib colorbar within the ax for a given
    patchset.

    Parameters
    ----------

    gci : matplotlib graphical instance

    orientation : string

    size : float

    pad : float


    **kwargs :
        Forwarded to updateColorBar
    """
    cbarTarget = plt
    cax = None
    divider = None
    #    if hasattr(patches, 'figure'):
    #       cbarTarget = patches.figure

    if hasattr(gci, 'ax'):
        divider = make_axes_locatable(gci.ax)
    if hasattr(gci, 'axes'):
        divider = make_axes_locatable(gci.axes)
    elif hasattr(gci, 'get_axes'):
        divider = make_axes_locatable(gci.get_axes())

    if divider:
        if orientation == 'horizontal':
            if pad is None:
                pad = 0.5
            cax = divider.append_axes("bottom", size=size, pad=pad)
        else:
            if pad is None:
                pad = 0.1
            cax = divider.append_axes("right", size=size, pad=pad)

    cbar = cbarTarget.colorbar(gci, cax=cax, orientation=orientation)

    updateColorBar(cbar, **kwargs)

    return cbar

def plot_FG_border(x1raw, y1raw, x2raw, y2raw, img1, hws_n1 = 35, hws_n2_range = [20, 125]):

    # plot x2GridFSim, y2GridFSim and border-grid
    hws_n1 = 35; hws_n2_range = [20, 125]
    # fitler outlier
    gpi = lstsq_filter(x1raw,y1raw,x2raw,y2raw, 100)
    print len(gpi), len(gpi[gpi])
    x1 = x1raw[gpi]
    y1 = y1raw[gpi]
    x2 = x2raw[gpi]
    y2 = y2raw[gpi]
    # initial drift inter-/extrapolation
    x2GridFSim, y2GridFSim = get_GridFSim(x1, y1, x2, y2, img1)
    # get border_grid
    border_grid = get_border_grid(x1, y1, img1, hws_n2_range)
    
    plt.figure(num=1, figsize=(14, 6), dpi=80, facecolor='w', edgecolor='k')
    plt.subplot(1,3,1)
    ax = plt.gca()
    #im = ax.imshow(np.arange(100).reshape((10,10)))
    im = plt.imshow(x2GridFSim)#;plt.colorbar(orientation='horizontal')
    plt.axis('off')
    # create an axes on the right side of ax. The width of cax will be 5%
    # of ax and the padding between cax and ax will be fixed at 0.05 inch.
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(im, cax=cax)
    plt.subplot(1,3,2)
    ax = plt.gca()
    #im = ax.imshow(np.arange(100).reshape((10,10)))
    im = plt.imshow(y2GridFSim)#;plt.colorbar(orientation='horizontal')
    plt.axis('off')
    # create an axes on the right side of ax. The width of cax will be 5%
    # of ax and the padding between cax and ax will be fixed at 0.05 inch.
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(im, cax=cax)
    plt.subplot(1,3,3)
    ax = plt.gca()
    #im = ax.imshow(np.arange(100).reshape((10,10)))
    im = plt.imshow(border_grid)#;plt.colorbar(orientation='horizontal')
    plt.axis('off')
    # create an axes on the right side of ax. The width of cax will be 5%
    # of ax and the padding between cax and ax will be fixed at 0.05 inch.
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    plt.colorbar(im, cax=cax)
    
    return plt

def test_axes_locatable_position():
    fig, ax = plt.subplots()
    divider = make_axes_locatable(ax)
    cax = divider.append_axes('right', size='5%', pad='2%')
    fig.canvas.draw()
    assert np.isclose(cax.get_position(original=False).width,
                      0.03621495327102808)

def setup_trace_xy(ax, t, y1, y2):
    ax.plot(y1, y2, 'o', color='lightblue')
    ax.set(aspect=1, adjustable='datalim')

    divider = make_axes_locatable(ax)

    hax = divider.append_axes('bottom', size='40%', pad=0, sharex=ax)
    vax = divider.append_axes('left', size='40%', pad=0, sharey=ax)

    for ax in [hax, vax]:
        ax.set(xticks=[], yticks=[], adjustable='datalim')
    ax.figure.add_axes(hax)
    ax.figure.add_axes(vax)

#    hax.set_xlabel('Trace 1 Amplitude')
#    vax.set_ylabel('Trace 2 Amplitude')

    vax.plot(t, y2, marker='o', color='black', mfc='lightblue')
    tt, y2 = scipy.ndimage.zoom(t, 30), scipy.ndimage.zoom(y2, 30)
    vax.fill_between(tt, y2, where=y2 > 0, facecolor='black')

    hax.plot(y1, t, marker='o', color='black', mfc='lightblue')
    tt, y1 = scipy.ndimage.zoom(t, 30), scipy.ndimage.zoom(y1, 30)
    hax.fill_betweenx(tt, y1, where=y1 > 0, facecolor='black')

    hax.margins(0.05)
    vax.margins(0.05)
    return hax, vax