def plot3D(data_3d, name='3D Plot'):

    X = [point[0] for point in data_3d]
    Y = [point[1] for point in data_3d]
    Z = [point[2] for point in data_3d]

    fig = plt.figure(name)
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(X, Y, Z, marker='o', c='b', zdir='y')

    # Create cubic bounding box to simulate equal aspect ratio
    max_range = np.array(
        [max(X) - min(X), max(Y) - min(Y), max(Z) - min(Z)]).max()
    Xb = 0.5 * max_range * \
        np.mgrid[-1:2:2, -1:2:2, -1:2:2][0].flatten() + 0.5 * (max(X) + min(X))
    Yb = 0.5 * max_range * \
        np.mgrid[-1:2:2, -1:2:2, -1:2:2][1].flatten() + 0.5 * (max(Y) + min(Y))
    Zb = 0.5 * max_range * \
        np.mgrid[-1:2:2, -1:2:2, -1:2:2][2].flatten() + 0.5 * (max(Z) + min(Z))

    # Comment or uncomment following both lines to test the fake bounding box:
    for xb, yb, zb in zip(Xb, Yb, Zb):
        ax.plot([xb], [yb], [zb], 'w')

    ax.set_xlabel('X (m)')
    ax.set_ylabel('Z (m)')
    ax.set_zlabel('Y (m)')

    plt.locator_params(nbins=4)
    plt.show()

    def plot_diff(self, x=None, ax=None, kwargs={}, offset = 0.):
        """
        Plots the difference between the Feature fit and the original Spectrum
        on an axis.
        If x is undefined, x is taken as inside the window region.
        If ax is undefined, a new axis is generated
        Kwargs must be a dictionary of legal matplotlib keywords
        """
        if x == None:
            x = self.x
        if not ax:
            fig = plt.figure()
            ax = fig.add_subplot(111)
            
        kill_label = False
        if 'label' not in kwargs.keys():
            kwargs['label'] = 'Difference'
            kill_label = True
        

        y = self(x) + offset

        window = self.window
        self.set_window([x.min(),x.max()])
        
        y -= self.y
        ax.plot(x,y,**kwargs)
        plt.locator_params(axis = 'y', nbins = 4)
        self.set_window(window)
        
        if kill_label:
            del kwargs['label']

def show_trand_chart(stats):
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        sys.exit(0)

    # prepare data
    net_lines = 0
    xlabels, axisx, axisy = [], [], []
    for stat in stats.items():
        weeks_ago, added, deleted = stat[0], stat[1]['add'], stat[1]['del']
        net_lines += added - deleted
        if net_lines == 0:
            continue

        xlabels.append(date_of_weeks_ago(weeks_ago))
        axisx.append(weeks_ago)
        axisy.append(net_lines)

    # start to plot
    plt.locator_params(axis = 'x', nbins = 4)
    plt.grid(True)
    plt.xlabel('weeks ago')
    plt.ylabel('LoC')
    plt.title('LineOfCode trend')
    plt.xticks(axisx, xlabels)
    plt.plot(axisx, axisy)
    plt.show()

def plot_diagnostics(image,data_filtered,pl_image,smooth_pl,inter_image,inter_pl,final_image,vgridsize,lg,ilng,flng,IV,output_file):
    fig3 = plt.figure(figsize=(8,vgridsize*2))
    gs3 = gridspec.GridSpec(vgridsize,2)
    for i in range(0,vgridsize):
        pl = np.power(lg,abs(IV[i,1]))    
        ax1 = fig3.add_subplot(gs3[i,0])
        [a,b,c,d] = ax1.plot(lg,np.multiply(pl,image[i,:]),lg,np.multiply(pl,data_filtered[i,:]),lg,np.multiply(pl,pl_image[i,:]),lg,np.multiply(pl,smooth_pl[i,:]))
        plt.setp(ax1.get_xticklabels(), visible=False)
        ax1.set_ylim([min(np.multiply(pl,image[i,:])),max(np.multiply(pl,image[i,:]))]) 
        plt.locator_params(axis='y',nbins=5)    
        if i < 1:
            ax1.set_title('Raw Images')
            plt.legend([a,b,c,d], ['Data', 'Data boxcar', 'Power Law','Power Law boxcar'],bbox_to_anchor=(0, 2.0, 1., .102),loc=9,
               borderaxespad=0.)
        ax2 = fig3.add_subplot(gs3[i,1],sharex=ax1)
        [e,f,g] = ax2.plot(lg,inter_image[i,:],lg,inter_pl[i,:],lg,final_image[i,:])
        plt.setp(ax2.get_xticklabels(), visible=False) 
        ax2.set_xlim([ilng+3,flng-3])
        ax2.set_ylim([min(inter_image[i,:]),max(final_image[i,:])])     
        ax2.hlines(0,ilng,flng,linestyles='dashed',color='0.75')
        #plt.locator_params(axis='y',nbins=5)     
        if i <1 :
            ax2.set_title('Subtracted Images')
            plt.legend([e,f,g], ['Subtracted Image', 'Subtracted Power Law', 'Final Image'],bbox_to_anchor=(0, 1.8, 1., .102),loc=9,
               borderaxespad=0.)
    plt.setp(ax1.get_xticklabels(),visible=True)
    plt.setp(ax2.get_xticklabels(),visible=True)
    plt.savefig(output_file+'_diagnostics.png')
    plt.close()

def main():
     #The url is made up of the prefix, day, and suffix:
     prefix = "http://www.wunderground.com/history/airport/KLGA/2016/1/"
     suffix = "/DailyHistory"
     days = []
     depths = []
     mins = [36,34,36,15,13,26,32,34,41,42,28,26,24,24,34,42,31,20,18,30,27,22,25,24,30,35,34,29,32,30,37] #min temps for January 2016
 #    days = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]          #Sets up a list to store days
 #    depths = []           #Sets up a list so store snow depths
     for day in range(1,32): #For each day
          days.append(day)       #Add the day to the list
          url = prefix+str(day)+suffix      #Make the url
          M = getDepthFromWeb(url)     #Call the function to extract Snow Depth
          
          depths.append(M) #Add the temp to the list
          print day, M 
 #    depths = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,7,27,20,15,9,6,4,3,2]
     depthScale = [(8)*((n*(n*.09))+1) for n in depths]
     plt.scatter(days, mins, s = depthScale, color = 'r', alpha=0.5)
     plt.title("Minimum Temps in January with Snow Depths as Point Size", y=1.02)       #Title for plot
     plt.xlabel('Days', fontsize=16).set_color("white")                     #Label for x-axis
     plt.ylabel('Minimum Temperature of the Day', fontsize=16).set_color("white")                   #Label for the y-axis
     plt.axis([0,32,10,50])
     plt.locator_params(axis = 'x', nbins = 16)
     plt.locator_params(axis = 'y', nbins = 20)
     l1 = plt.scatter([],[], s=8, edgecolors='none', color = 'r')
     labels = [" = 0 in. or Trace Amounts of Snow"]
     leg = plt.legend([l1], labels, ncol=4, frameon=True, fontsize=12,
     handlelength=1, loc = 'upper left', borderpad = .5,
     handletextpad=1, title='Size of Points Corresponds to Snow Depth', scatterpoints = 1)
     plt.show()

    def plot_kmeans_components(self, fig1, gs, kmeans_clusters, clrs, plot_title='Hb', num_subplots=1,
                               flag_separate=1, gridspecs=[0, 0]):
        with sns.axes_style('darkgrid'):
            if flag_separate:
                ax1 = fig1.add_subplot(2, 1, num_subplots)
            else:
                ax1 = eval('fig1.add_subplot(gs' + gridspecs + ')')

            plt.gca().set_color_cycle(clrs)
            for ii in xrange(0, size(kmeans_clusters, 1)):
                plt.plot(kmeans_clusters[:, ii], lw=4, label=str(ii))

            plt.locator_params(axis='y', nbins=4)
            sns.axlabel("Time (seconds)", "a.u")
            ax1.legend(prop={'size': 14}, loc='center left', bbox_to_anchor=(1, 0.5), ncol=1, fancybox=True,
                       shadow=True)

            plt.title(plot_title, fontsize=14)

            plt.ylim((min(kmeans_clusters) - 0.0001,
                      max(kmeans_clusters) + 0.0001))
            plt.xlim((0, size(kmeans_clusters, 0)))
            plt.axhline(y=0, linestyle='-', color='k', linewidth=1)
            self.plot_vertical_lines_onset()
            self.plot_vertical_lines_offset()

def main():
  global algorithms
  global datadir
  global history
  global legend_location
  args = argparser.parse_args()
  algorithms = set(algorithms) - set(args.exclude_algorithm)
  datadir = args.datadir
  history = args.history
  legend_location = args.legend_location
  plt.rc('font', size=5)
  plt.rc('axes', color_cycle=['r','g','b','y'])
  plt.figure(figsize=(2,2))
  plt.locator_params(axis='x', nbins=5)
  if args.plot_all or args.steps_frequency:
    plot_steps_frequency()
  if args.plot_all or args.stepssum_vs_opcount:
    plot_stepssum_vs_opcount()
  if args.plot_all or args.nodesmax_vs_opcount:
    plot_nodesmax_vs_opcount()
  if args.plot_all or args.stepssum_vs_history:
    plot_stepssum_vs_history()
  if args.plot_all or args.nodesmax_vs_history:
    plot_nodesmax_vs_history()
  if args.plot_all or args.stepsavg_vs_history:
    plot_stepsavg_vs_history()
  if args.plot_all or args.stepsmed_vs_history:
    plot_stepsmed_vs_history()
  if args.plot_all or args.stepsavgdev_vs_history:
    plot_stepsavgdev_vs_history()
  if args.plot_all or args.stepsmedmax_vs_history:
    plot_stepsmedmax_vs_history()

def init_plotting():

    plt.rcParams['figure.figsize']       = (16, 9)
    plt.rcParams['figure.dpi']           = 75
    plt.rcParams['font.size']            = 16
    plt.rcParams['font.family']          = 'Times New Roman'
    plt.rcParams['axes.labelsize']       = 18
    plt.rcParams['axes.titlesize']       = 20
    plt.rcParams['legend.fontsize']      = plt.rcParams['font.size']
    plt.rcParams['xtick.labelsize']      = plt.rcParams['font.size']
    plt.rcParams['ytick.labelsize']      = plt.rcParams['font.size']
    plt.rcParams['xtick.major.size']     = 3
    plt.rcParams['xtick.minor.size']     = 3
    plt.rcParams['xtick.major.width']    = 1
    plt.rcParams['xtick.minor.width']    = 1
    plt.rcParams['ytick.major.size']     = 3
    plt.rcParams['ytick.minor.size']     = 3
    plt.rcParams['ytick.major.width']    = 1
    plt.rcParams['ytick.minor.width']    = 1
    plt.rcParams['legend.frameon']       = True
    plt.rcParams['legend.shadow']        = True
    plt.rcParams['legend.loc']           = 'lower left'
    plt.rcParams['legend.numpoints']     = 1
    plt.rcParams['legend.scatterpoints'] = 1
    plt.rcParams['axes.linewidth']       = 1
    plt.rcParams['savefig.dpi']          = 300
    plt.rcParams['xtick.minor.visible']  = 'False'
    plt.rcParams['ytick.minor.visible']  = 'False'
    plt.gca().xaxis.set_ticks_position('bottom')
    plt.gca().yaxis.set_ticks_position('left')
    plt.locator_params(nticks=4)

def plot_posterior_hist(ax, variable, samples, validation_data=None):
    pyplot.locator_params(axis = 'x', nbins = 4)
    histogram = ax.hist(samples, bins=12)
    if validation_data:
        axvline(x=validation_data[variable])
    ax.set_xlabel(variable)
    ax.set_ylim(0, max(histogram[0])*1.1)

    def plot(self, fig=None):
        import matplotlib.pyplot as plt
        if fig == None and not isinstance(self.fig, plt.Figure):
            self.fig = plt.figure()
        else:
            self.fig = fig
        # plt.gcf(self.fig)
        plt.figure(self.fig.number)
        plt.hold(True)

        from numpy import max as npmax

        mode = "logpow"
        if mode == "logpow":
            gdata = self.get_logpow()
            plt.ylabel("power [db]")

        plt.plot(self.get_xdata() / 1000., gdata)

        ymax = int(npmax(gdata[1:]) / 10.) * 10 + 10
        ymin = ymax - 80

        plt.ylim(ymin, ymax)

        plt.xlabel('Frequency [kHz]')
        plt.locator_params(nbins=20, axis='x', tight=True)
        # plt.locator_params(nbins=15, axis='y', tight=True, fontsize=1)
        plt.grid()
        plt.title(self.name)
        return self

def graph_drugs_line(items,nic):

	months = Config.filenames
	months = [x.strip('.csv') for x in months]
	months.reverse()

	for drug in items.keys():
		items_list = []
		nics=[]
		for month in months:
			try:
				items_list.append(items[drug][month])
				nics.append(nic[drug][month])
			except KeyError:
				print drug + ' not all information available'				
				break
		else:						 
			
			index = numpy.arange(len(months))
			graph = plt.plot(index, items_list, 'r.-', index, nics, 'g.-')
			
			lessMonths = [months[int(i)] for i in numpy.linspace(0,len(months)-1,num=6)]
			ax = plt.gca()
			plt.locator_params(nbins=len(lessMonths))
			ax.set_xticklabels(lessMonths)
			
			ax.set_ylabel('Branded/Generic')
			ax.set_title('Percent branded for chemical: '+drug)

			ax.legend( ('items', 'nic') )

			plt.savefig('Time_ChemPercents_figures/' + drug)
			plt.clf()

def main():
     #The url is made up of the prefix, day, and suffix:
     prefix = "http://www.wunderground.com/history/airport/KLGA/2016/1/"
     suffix = "/DailyHistory"
     days = []          #Sets up a list to store days
     mins = []           #Sets up a list so store min values
     for day in range(1,32): #For each day
          days.append(day)       #Add the day to the list
          url = prefix+str(day)+suffix      #Make the url
          M = getTempFromWeb("Min",url)     #Call the function to extract temp
          mins.append(M) #Add the temp to the list
          print day, M 
 ##    days = [i for i in range(1,32)]   Used for debugging programin
 ##    mins = [36,34,36,15,13,26,32,34,41,42,28,26,24,24,34,42,31,20,18,30,27,22,25,24,30,35,34,29,32,30,37]
     ave = float(sum(mins))/ len(mins)
     print ave
     print len(mins)
     

     
     scaled= [i*100/ave-100 for i in mins]
     
     
   
     plt.plot(days, scaled, color='r', label="Variation of Temp")       #Plot max as red
     plt.axhline(y=0,color='b')
     plt.legend(['% Difference From Avg'], loc='upper left',prop={'size':10})
     plt.title("Variation of January Min Temps from Average Min", y=1.02)       #Title for plot
     plt.xlabel('Days', fontsize=16).set_color("white")                     #Label for x-axis
     plt.ylabel('% Fluctuation from Average', fontsize=16).set_color("white")                   #Label for the y-axis
     plt.axis([0,32,-60,60])
     plt.locator_params(axis = 'x', nbins = 16)
     plt.locator_params(axis = 'y', nbins = 20)

     plt.show()

def plot_dshift(data, sample, fout, dshift_ind=0, close=True):
    """ 
    PLOT THE POSTERIOR FOR THE DOPPLER SHIFT
    """
    
    fig = plt.figure(figsize=(12,9))
    plt.locator_params(axis = 'x', nbins = 6)

    # Plot a historgram and kernel density estimate
    ax = fig.add_subplot(111)
    sns.distplot(sample[:,11+dshift_ind], hist_kws={"histtype":"stepfilled"}, ax=ax)
    plt.xticks(rotation=45)


    _, ymax = ax.get_ylim()
    ax.vlines(data["dshift"], 0, ymax, lw=3, color="black")

    ax.set_xlabel(r"Doppler shift $d=v/c$")
    ax.set_ylabel("$N_{\mathrm{samples}}$")
    plt.tight_layout()

    plt.savefig(fout+"_dshift%i.png"%dshift_ind, format="png")
    if close:
        plt.close()
    return

def main(plot_type):
    per_size = 5
    nrow, ncol = len(graphs), len(params)
    fig = plt.figure(figsize=(ncol * per_size, nrow * per_size))

    if plot_type.startswith('dist'):
        angle = (10, 45)
    else:
        angle = (15, 210)

    for i, gname in enumerate(graphs):
        for j, param in enumerate(params):
            idx = i * ncol + j + 1
            ax = fig.add_subplot(nrow, ncol, idx, projection='3d')
            plot_surface(gname, param, plot_type,
                         fig, ax=ax,
                         dirname=dirname,
                         angle=angle,
                         use_colorbar=False)
            ax.set_title('{}({})'.format(gname, param))
            plt.locator_params(axis='y', nbins=5)
            plt.locator_params(axis='x', nbins=5)

    fig_dir = 'figs/{}'.format(fig_dirname)
    if not os.path.exists(fig_dir):
        os.makedirs(fig_dir)
    figpath = '{}/{}.pdf'.format(fig_dir, plot_type)
    print(figpath)
    fig.savefig(figpath)

def save_records_plot(file_path, ls_monitors, name, n_train_batches, legend_loc="upper right"):
    """
    Save a plot of a list of monitors' history.
    Args:
        file_path (string): the folder path where to save the plot
        ls_monitors: the list of statistics to plot
        name: name of file to be saved
        n_train_batches: the total number of training batches
    """

    lines = ["--", "-", "-.",":"]
    linecycler = cycle(lines)

    plt.figure()
    for m in ls_monitors:
        X = [i/float(n_train_batches) for i in m.history_minibatch]
        Y = m.history_value
        a, b = zip(*sorted(zip(X, Y)))
        plt.plot(a, b, next(linecycler), label=m.name)

    plt.xlabel('Training epoch')
    plt.ylabel(ls_monitors[0].type)
    plt.legend(loc=legend_loc)
    plt.locator_params(axis='y', nbins=7)
    plt.locator_params(axis='x', nbins=10)
    plt.savefig(file_path + name + ".png")
    tikz_save(file_path + name + ".tikz", figureheight = '\\figureheighttik', figurewidth = '\\figurewidthtik')

def plot_rejection_sampling(thetas, posteriors, x_accepts, bins):
	""" Plot analytical solution and rejection sampling solution on same graph """
	fig, ax = plt.subplots()
	plt.plot(thetas, posteriors, linewidth=3)
	
	#Rejection sampling plot
	hist, bin_edges = numpy.histogram(x_accepts, bins)	
	bin_width = bin_edges[1] - bin_edges[0]
	hist = hist / max(hist)
	ax.bar(bin_edges[:-1], hist, bin_width, color='green')
	ax.tick_params(axis='both', which='major', labelsize=30)

	# Create strings to show numerical mean and standard deviation on graphs
	mean = numpy.mean(x_accepts)
	stdev = numpy.std(x_accepts)
	display_string = ('$\mu_{{MC}} = {0:.3f} $\n$\sigma_{{MC}} = {1:.3f}$').format(mean, stdev)

	plt.xlabel(r'$\theta$', fontsize=30)
	plt.ylabel(r'$\propto P(\theta|x)$', fontsize=30)
	plt.text(0.6, 0.8, display_string, transform=ax.transAxes, fontsize=30)
	plt.locator_params(axis='x', nbins=5)
	plt.savefig('rejection.png', bbox_inches='tight')

	# Plot log
	fig, ax = plt.subplots()
	plt.plot(thetas, -numpy.log(posteriors), linewidth=3)
	ax.bar(bin_edges[:-1], -numpy.log(hist), bin_width, color='green')
	ax.tick_params(axis='both', which='major', labelsize=30)

	plt.xlabel(r'$\theta$', fontsize=30)
	plt.ylabel(r'$\propto log(P(\theta|x))$', fontsize=30)
	plt.locator_params(axis='x', nbins=5)
	plt.text(0.3, 0.5, display_string, transform=ax.transAxes, fontsize=30)

	plt.savefig('rejlog.png', bbox_inches='tight')

def plot_results(dists):
    for i, d in enumerate(dists):
        ax = plt.subplot(3,3,(4*i)+1)
        N, bins, patches = plt.hist(d.data, color="b",ec="k", bins=30, \
                                    range=tuple(d.lims), normed=True, \
                                    edgecolor="k", histtype='bar',linewidth=1.)
        fracs = N.astype(float)/N.max()
        norm = Normalize(-.2* fracs.max(), 1.5 * fracs.max())
        for thisfrac, thispatch in zip(fracs, patches):
            color = cm.gray_r(norm(thisfrac))
            thispatch.set_facecolor(color)
            thispatch.set_edgecolor("w")
        x = np.linspace(d.data.min(), d.data.max(), 100)
        ylim = ax.get_ylim()
        plt.plot(x, d.best.pdf(x), "-r", lw=1.5, alpha=0.7)
        ax.set_ylim(ylim)
        plt.axvline(d.best.MAPP, c="r", ls="--", lw=1.5)
        plt.tick_params(labelright=True, labelleft=False, labelsize=10)
        plt.xlim(d.lims)
        plt.locator_params(axis='x',nbins=10)
        if i < 2:
            plt.setp(ax.get_xticklabels(), visible=False)
        else:
            plt.xlabel(r"[$\mathregular{\alpha}$ / Fe]")
        plt.minorticks_on()
    def hist2D(dist1, dist2):
        """ Plot distribution and confidence contours. """
        X, Y = np.mgrid[dist1.lims[0] : dist1.lims[1] : 20j,
                        dist2.lims[0] : dist2.lims[1] : 20j]
        extent = [dist1.lims[0], dist1.lims[1], dist2.lims[0], dist2.lims[1]]
        positions = np.vstack([X.ravel(), Y.ravel()])
        values = np.vstack([dist1.data, dist2.data])
        kernel = stats.gaussian_kde(values)
        Z = np.reshape(kernel(positions).T, X.shape)
        ax.imshow(np.rot90(Z), cmap="gray_r", extent=extent, aspect="auto",
                  interpolation="spline16")
        plt.axvline(dist1.best.MAPP, c="r", ls="--", lw=1.5)
        plt.axhline(dist2.best.MAPP, c="r", ls="--", lw=1.5)
        plt.tick_params(labelsize=10)
        ax.minorticks_on()
        plt.locator_params(axis='x',nbins=10)
        return
    ax = plt.subplot(3,3,4)
    hist2D(dists[0], dists[1])
    plt.setp(ax.get_xticklabels(), visible=False)
    plt.ylabel("[Z/H]")
    plt.xlim(dists[0].lims)
    plt.ylim(dists[1].lims)
    ax = plt.subplot(3,3,7)
    hist2D(dists[0], dists[2])
    plt.ylabel(r"[$\mathregular{\alpha}$ / Fe]")
    plt.xlabel("log Age (yr)")
    plt.xlim(dists[0].lims)
    plt.ylim(dists[2].lims)
    ax = plt.subplot(3,3,8)
    plt.xlabel("[Z/H]")
    hist2D(dists[1], dists[2])
    plt.xlim(dists[1].lims)
    plt.ylim(dists[2].lims)
    return

def main():
    dir='/Users/ph290/Public/mo_data/ostia/' # on ELD140
    filename = dir + '*.nc'

    cube = iris.load_cube(filename,'sea_surface_temperature',callback=my_callback)
    #reads in data using a special callback, because it is a nasty netcdf file
    cube.data=cube.data-273.15


    sst_mean = cube.collapsed('time', iris.analysis.MEAN)
    #average all 12 months together
    sst_stdev=cube.collapsed('time', iris.analysis.STD_DEV)

    caribbean = iris.Constraint(
                                    longitude=lambda v: 260 <= v <= 320,
                                    latitude=lambda v: 0 <= v <= 40,
                                    name='sea_surface_temperature'
                                    )

    caribbean_sst_mean = sst_mean.extract(caribbean)
    caribbean_sst_stdev = sst_stdev.extract(caribbean)
    #extract the Caribbean region

    fig = plt.figure()
    ax = plt.axes(projection=ccrs.PlateCarree())
    data=caribbean_sst_mean.data
    data2=caribbean_sst_stdev.data
    lons = caribbean_sst_mean.coord('longitude').points
    lats = caribbean_sst_mean.coord('latitude').points
    lo = np.floor(data.min())
    hi = np.ceil(data.max())
    levels = np.linspace(lo,hi,100)
    lo2 = np.floor(data2.min())
    hi2 = np.ceil(data2.max())
    levels2 = np.linspace(lo2,5,10)
    cube_label = 'latitude: %s' % caribbean_sst_mean.coord('latitude').points
    contour=plt.contourf(lons, lats, data,transform=ccrs.PlateCarree(),levels=levels,xlabel=cube_label)
    #filled contour the annually averaged temperature
    contour2=plt.contour(lons, lats, data2,transform=ccrs.PlateCarree(),levels=levels2,colors='k')
    #contour the standard deviations
    plt.clabel(contour2, inline=0.5, fontsize=12,fmt='%1.1f' )
    ax.add_feature(cartopy.feature.LAND)
    ax.coastlines()
    ax.add_feature(cartopy.feature.RIVERS)
    ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
    #ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
    cbar = plt.colorbar(contour, ticks=np.linspace(lo,hi,7), orientation='horizontal')
    cbar.set_label('Sea Surface Temperature ($^\circ$C)')
    # enable axis ticks
    ax.xaxis.set_visible(True)
    ax.yaxis.set_visible(True)
    # fix 10-degree increments and don't clip range
    plt.locator_params(steps=(1,10), tight=False)
    # add gridlines
    plt.grid(True)
    # add axis labels
    plt.xlabel('longitude')
    plt.ylabel('latitude')
    #plt.show()
    plt.savefig('/home/h04/hador/public_html/twiki_figures/carib_sst_and_stdev.png')

def plot_KDE(sol, var1, var2, fig=None, ax=None, draw=True, save=False, save_as_png=False, dpi=None):
    """
    Like the hexbin plot but a 2D KDE
    Pass mcmcinv object and 2 variable names as strings
    """
    ext = ['png' if save_as_png else 'pdf'][0]
    if fig == None or ax == None:
        fig, ax = plt.subplots(figsize=(3,3))
    MDL = sol.MDL
    if var1 == "R0":
        stoc1 = "R0"
    else:
        stoc1 =  ''.join([i for i in var1 if not i.isdigit()])
        stoc_num1 = [int(i) for i in var1 if i.isdigit()]
    try:
        x = MDL.trace(stoc1)[:,stoc_num1[0]-1]
    except:
        x = MDL.trace(stoc1)[:]
    if var2 == "R0":
        stoc2 = "R0"
    else:
        stoc2 =  ''.join([i for i in var2 if not i.isdigit()])
        stoc_num2 = [int(i) for i in var2 if i.isdigit()]
    try:
        y = MDL.trace(stoc2)[:,stoc_num2[0]-1]
    except:
        y = MDL.trace(stoc2)[:]
    xmin, xmax = min(x), max(x)
    ymin, ymax = min(y), max(y) 
    # Peform the kernel density estimate
    xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
    positions = np.vstack([xx.ravel(), yy.ravel()])
    values = np.vstack([x, y])
    kernel = gaussian_kde(values)
    kernel.set_bandwidth(bw_method='silverman')
#        kernel.set_bandwidth(bw_method=kernel.factor * 2.)
    f = np.reshape(kernel(positions).T, xx.shape)

    ax.set_xlim(xmin, xmax)
    ax.set_ylim(ymin, ymax)
    plt.sca(ax)
    # Contourf plot
    plt.grid(None)
    plt.ticklabel_format(style='sci', axis='both', scilimits=(0,0))
    plt.xticks(rotation=90)
    plt.locator_params(axis = 'y', nbins = 7)
    plt.locator_params(axis = 'x', nbins = 7)
    ax.contourf(xx, yy, f, cmap=plt.cm.viridis, alpha=0.8)
    ax.scatter(x, y, color='k', s=1, zorder=2)

    plt.ylabel("%s" %var2)
    plt.xlabel("%s" %var1)

    if save: 
        fn = 'KDE-%s-%s.%s'%(sol.model_type_str,sol.filename,ext)
        save_figure(fig, subfolder='2D-KDE', fname=fn, dpi=dpi)
    
    plt.close(fig)        
    if draw:    return fig
    else:       return None

def triangular_wave_plot(ax=None):
    """Plot of a triangular wave function.
    """
    import numpy as np
    import matplotlib.pyplot as plt

    if ax is None:
        fig, ax = plt.subplots(1, 1, figsize=(9, 3))
    t = np.array([-3, -2, -1, 0, 1, 2, 3])*np.pi
    x = np.array([0, 1, 0, 1, 0, 1, 0])
    ax.plot(t, x, linewidth=3, label=r'Square wave')
    ax.spines['bottom'].set_position('zero')
    ax.spines['top'].set_color('none')
    ax.spines['left'].set_position('zero')
    ax.spines['right'].set_color('none')
    ax.xaxis.set_ticks_position('bottom')
    ax.yaxis.set_ticks_position('left')
    ax.tick_params(axis='both', direction='inout', which='both', length=5)
    ax.set_xlim((-3*np.pi, 3*np.pi))
    ax.set_ylim((-0.1, 1.1))
    ax.set_xticks(np.linspace(-3*np.pi-0.1, 3*np.pi+0.1, 7))
    ax.set_xticklabels(['$-3\pi$', '$-2\pi$', '$-\pi$', '$0$', '$\pi$', '$2\pi$', '$3\pi$'],
                       fontsize=16)
    plt.locator_params(axis='y', nbins=3)
    ax.annotate(r'$t$', xy=(3*np.pi, 0.1), xycoords = 'data', xytext=(0, 0),
                textcoords = 'offset points', size=18, color='k')
    ax.annotate(r'$x[t]$', xy=(.1, 1.03), xycoords = 'data', xytext=(0, 0),
                textcoords = 'offset points', size=18, color='k')
    ax.grid()
    fig.tight_layout()

    return ax