def draw_pdf_contours(dist, border=False, nlevels=200, subdiv=8, **kwargs):
    '''Draws pdf contours over an equilateral triangle (2-simplex).

    Arguments:

        `dist`: A distribution instance with a `pdf` method.

        `border` (bool): If True, the simplex border is drawn.

        `nlevels` (int): Number of contours to draw.

        `subdiv` (int): Number of recursive mesh subdivisions to create.

        kwargs: Keyword args passed on to `plt.triplot`.
    '''
    from matplotlib import ticker, cm
    import math
 
    refiner = tri.UniformTriRefiner(_triangle)
    trimesh = refiner.refine_triangulation(subdiv=subdiv)
    pvals = [dist.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]
 
    plt.tricontourf(trimesh, pvals, nlevels, **kwargs)
    plt.axis('equal')
    plt.xlim(0, 1)
    plt.ylim(0, 0.75**0.5)
    plt.axis('off')
    if border is True:
        plt.hold(1)
        plt.triplot(_triangle, linewidth=1)

def plot(filename):
    import os
    from matplotlib.pyplot import clf, tricontour, tricontourf, \
        gca, savefig, rc, minorticks_on

    if not os.path.exists(filename):
        return -1

    rc('text', usetex=True)
    clf()
    x, y, tri, ux, uy = load_velocity(filename)
    tricontourf(x, y, tri, ux, 16)
    tricontour(x, y, tri, ux, 16, linestyles='-',
               colors='black', linewidths=0.5)
    minorticks_on()
    gca().set_aspect('equal')
    gca().tick_params(direction='out', which='both')
    gca().set_xticklabels([])
    gca().set_yticklabels([])

    name, _ = os.path.splitext(filename)
    name = os.path.basename(name)

    savefig('{0}.png'.format(name), dpi=300, bbox_inches='tight')
    savefig('{0}.pdf'.format(name), bbox_inches='tight')

    def graph_grid(self, debug=False):
        ''' 2D xy plot of bathymetry and mesh.
            No inputs required so far'''
        if debug or self._debug:
            print 'Plotting grid...'
        nv = self._var.nv.T -1
        h = self._var.h
        tri = Tri.Triangulation(self._var.lon, self._var.lat, triangles=nv)

        levels=np.arange(-100,-4,1)   # depth contours to plot

        fig = plt.figure(figsize=(18,10))
        plt.rc('font',size='22')
        ax = fig.add_subplot(111,aspect=(1.0/np.cos(np.mean(self._var.lat)*np.pi/180.0)))
        plt.tricontourf(tri, -h,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
        plt.triplot(tri)
        plt.ylabel('Latitude')
        plt.xlabel('Longitude')
        plt.gca().patch.set_facecolor('0.5')
        cbar=plt.colorbar()
        cbar.set_label('Water Depth (m)', rotation=-90,labelpad=30)

        scale = 1
        ticks = ticker.FuncFormatter(lambda lon, pos: '{0:g}'.format(lon/scale))
        ax.xaxis.set_major_formatter(ticks)
        ax.yaxis.set_major_formatter(ticks)
        plt.grid()
        plt.show()

        if debug or self._debug:
            print '...Passed'

 def plot(self, N=6, cm=plt.cm.jet):
     plt.figure()
     plt.gca().set_aspect('equal')
     plt.tricontourf(self.triang, self.density, N, cm=cm)
     plt.colorbar()
     plt.tricontour(self.triang, self.density, N, colors='k')
     plt.show()

def fcontour_plot_dataset(file_path,hdf_file_name,x_variable,y_variable,z_variable,
                          clims=None,num_levels=41):
    '''Script to make a contour plot of a dataset from an HDF5 files of several
    scans combined.
    '''
    #Create a new figure
    plt.figure()
    plt.suptitle(hdf_file_name)
    #Open file
    hdf_file = h5py.File(file_path + hdf_file_name,'r')
    #Mask off any NAN entries is x; indicates scan wasn't as wide as widest scan
    mask = np.isfinite(hdf_file[x_variable])
    #Make triangulation for Delauney triangulation plot
    triang = tri.Triangulation(hdf_file[x_variable][mask],
                               hdf_file[y_variable][mask])
    #Create contour plot
    if clims:
        contour_levels = np.linspace(clims[0],clims[1],num_levels)
        plt.tricontourf(triang,hdf_file[z_variable][mask],contour_levels,extend='both')
    else:
        plt.tricontourf(triang,hdf_file[z_variable][mask],num_levels,extend='both')
    plt.xlabel(x_variable)
    plt.ylabel(y_variable)
    cb = plt.colorbar()
    cb.ax.set_ylabel(z_variable)
    plt.show()

def isosurf(G,Z,titre):
    """ trace isosurface de Z sur le maillage G"""
    triang=tri.Triangulation(G.X[:,0],G.X[:,1],triangles=G.Tbc)
    plt.tricontourf(triang, Z)
    plt.colorbar()
    plt.title(titre)
    return

def my_plot(u,v,t,daystr,levels):
    #boston light swim
    ax= [-71.10, -70.10, 41.70, 42.70] # region to plot
    vel_arrow = 0.2 # velocity arrow scale
    subsample = 8  # subsampling of velocity vectors

    # find velocity points in bounding box
    ind = np.argwhere((lonc >= ax[0]) & (lonc <= ax[1]) & (latc >= ax[2]) & (latc <= ax[3]))

    np.random.shuffle(ind)
    Nvec = int(len(ind) / subsample)
    idv = ind[:Nvec]
    # tricontourf plot of water depth with vectors on top
    plt.figure(figsize=(20,10))
    plt.subplot(111,aspect=(1.0/np.cos(lat[:].mean()*np.pi/180.0)))
    #tricontourf(tri, t,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
    plt.tricontourf(tri, t,levels=levels,shading='faceted')
    plt.axis(ax)
    plt.gca().patch.set_facecolor('0.5')
    cbar=plt.colorbar()
    cbar.set_label('Forecast Surface Temperature (C)', rotation=-90)
    plt.tricontour(tri, t,levels=[0])
    Q = plt.quiver(lonc[idv],latc[idv],u[idv],v[idv],scale=10)
    maxstr='%3.1f m/s' % vel_arrow
    qk = plt.quiverkey(Q,0.92,0.08,vel_arrow,maxstr,labelpos='W')
    plt.title('NECOFS Surface Velocity, Layer %d, %s UTC' % (ilayer, daystr))
    plt.plot(lon_track,lat_track,'m-o')
    plt.plot(lon_buoy,lat_buoy,'y-o')

def contour(data, x, y, label=None, log=False):

    tri=Triangulation(x,y)

    plt.close('all')
    plt.figure()
    ax=plt.subplot(111)
    ax.minorticks_on()
    if(log):
        ax.set_xscale("log",nonposx='clip')
        ax.set_yscale("log",nonposy='clip')

    ax.set_xlim([min(x.min(),y.min()),max(x.max(),y.max())])
    ax.set_ylim([min(x.min(),y.min()),max(x.max(),y.max())])
    plt.xlabel('r [AU]')
    plt.ylabel('z [AU]')

    nmax=data.max()
    nmin=data.min()
    levels=np.logspace(np.log10(nmin),np.log10(nmax),num=12)

    plt.tricontourf(tri, data, levels, norm=colors.LogNorm(vmin=nmin, vmax=nmax))
    cbar=plt.colorbar(format='%.2e')
    cbar.set_label(label)

    CS=plt.tricontour(tri, data, levels, colors='black', linewidths=1.5)
    plt.clabel(CS, fontsize=8, inline=1)
    cbar.add_lines(CS)

    plt.triplot(tri, color='black', alpha=0.2)

    plt.show(block=False)

def test_tricontourf_decreasing_levels():
    # github issue 5477.
    x = [0.0, 1.0, 1.0]
    y = [0.0, 0.0, 1.0]
    z = [0.2, 0.4, 0.6]
    plt.figure()
    with pytest.raises(ValueError):
        plt.tricontourf(x, y, z, [1.0, 0.0])

def testSquare():

    ptlist = {
        "vertices": np.array(
            ((0.0, 0.0), (0.5, 0.0), (1.0, 0.0), (0.0, 0.5), (0.5, 0.5), (1.0, 0.5), (0.0, 1.0), (0.5, 1.0), (1.0, 1.0))
        )
    }
    t = triangle.triangulate(ptlist)
    t1 = triangle.triangulate(ptlist, "qa0.001")

    triangle.plot.compare(plt, t, t1)
    #    plt.show()

    L, M = FE.assembleMatrices(t)
    #    print L
    #    print '\n\n'
    #    print M

    np.savetxt("textL", L)
    np.savetxt("textM", M)

    eig = FE.eigenvalues(L, M)
    elist = eig[0]
    efunc = eig[1]
    print elist[0]
    print elist[1]
    print elist[2]

    #    vertices = np.asarray(t['vertices'])
    #    faces = np.asarray(t['triangles'])
    #    x = vertices[:,0]
    #    y = vertices[:,1]

    #    z = efunc[1]

    #    plt.figure()
    #    plt.tricontourf(x,y,faces,z,cmap='afmhot')
    #    plt.show()
    print "****************************"

    L, M = FE.assembleMatrices(t1)
    eig = FE.eigenvalues(L, M)
    elist = eig[0]
    efunc = eig[1]
    for j in range(10):
        print elist[j]

    vertices = np.asarray(t1["vertices"])
    faces = np.asarray(t1["triangles"])
    x = vertices[:, 0]
    y = vertices[:, 1]
    z = efunc[:, 5]

    plt.figure()
    plt.tricontourf(x, y, z, 100, cmap="afmhot")
    plt.show()

    print "***************************\n\n\n\n\n"

def show_data_domain_2D(samples, data, Q_ref, ref_markers=None,
        ref_colors=None, xlabel=r'$q_1$', ylabel=r'$q_2$',
        triangles=None, save=True, interactive=False, filenames=None):
    r"""
    Plot the data domain D using a triangulation based on the generating
    samples with a marker for various :math:`Q_{ref}`. Assumes that the first
    dimension of data is :math:`q_1`.

    :param samples: Samples to plot
    :type samples: :class:`~numpy.ndarray` of shape (num_samples, ndim)
    :param data: Data associated with ``samples``
    :type data: :class:`numpy.ndarray`
    :param Q_ref: reference data value
    :type Q_ref: :class:`numpy.ndarray` of shape (M, 2)
    :param list ref_markers: list of marker types for :math:`Q_{ref}`
    :param list ref_colors: list of colors for :math:`Q_{ref}`
    :param string xlabel: x-axis label
    :param string ylabel: y-axis label
    :param triangles: triangulation defined by ``samples``
    :type triangles: :class:`tri.Triuangulation.triangles`
    :param bool save: flag whether or not to save the figure
    :param bool interactive: flag whether or not to show the figure
    :param list filenames: file names for the unmarked and marked domain plots

    """
    if ref_markers == None:
        ref_markers = markers
    if ref_colors == None:
        ref_colors = colors
    if type(triangles) == type(None):
        triangulation = tri.Triangulation(samples[:, 0], samples[:, 1])
        triangles = triangulation.triangles
    if filenames == None:
        filenames = ['domain_q1_q2_cs.eps', 'q1_q2_domain_Q_cs.eps']

    Q_ref = util.fix_dimensions_data(Q_ref, 2)

    # Create figure
    plt.tricontourf(data[:, 0], data[:, 1], np.zeros((data.shape[0],)),
            triangles=triangles, colors='grey') 
    plt.autoscale(tight=True)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.savefig(filenames[0], bbox_inches='tight', transparent=True,
            pad_inches=0)
    # Add truth markers
    for i in xrange(Q_ref.shape[0]):
        plt.scatter(Q_ref[i, 0], Q_ref[i, 1], s=60, c=ref_colors[i],
                marker=ref_markers[i])
    if save:
        plt.savefig(filenames[1], bbox_inches='tight', transparent=True,
            pad_inches=0)
    if interactive:
        plt.show()
    else:
        plt.close()

    def graphGrid(self,narrowGrid=False, plot=False):
        #nx.draw(self.graph, self.pointIDXY)
        #plt.show()

        #lat = self.data.variables['lat'][:]
        #lon = self.data.variables['lon'][:]
        #nv = self.data.variables['nv'][:].T -1
        #h = self.data.variables['h'][:]
        lat = self.lat
        lon = self.lon
        nv = self.nv.T - 1
        h = self.h

        tri = Tri.Triangulation(lon, lat, triangles=nv)  # xy or latlon based on how you are #Grand Passage

        levels=np.arange(-38,6,1)   # depth contours to plot

        fig = plt.figure(figsize=(18,10))
        plt.rc('font',size='22')
        ax = fig.add_subplot(111,aspect=(1.0/np.cos(np.mean(lat)*np.pi/180.0)))
        plt.tricontourf(tri, -h,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
        plt.triplot(tri)
        plt.ylabel('Latitude')
        plt.xlabel('Longitude')
        plt.gca().patch.set_facecolor('0.5')
        cbar=plt.colorbar()
        cbar.set_label('Water Depth (m)', rotation=-90,labelpad=30)

        scale = 1
        ticks = ticker.FuncFormatter(lambda lon, pos: '{0:g}'.format(lon/scale))
        ax.xaxis.set_major_formatter(ticks)
        ax.yaxis.set_major_formatter(ticks)
        plt.grid()

        maxlat, maxlon = np.max(self.maxcoordinates,axis=0)
        minlat, minlon = np.min(self.mincoordinates,axis=0)
        if narrowGrid:
            ax.set_xlim(minlon,maxlon)
            ax.set_ylim(minlat,maxlat)


        zz = len(self.elements)
        for i,v in enumerate(self.elements):
            source = self.pointIDXY[v[0]]
            target = self.pointIDXY[v[-1]]
            lab = '({:.6},{:.6})-({:.6},{:.6})'.format(source[0], source[1],
                                                       target[0], target[1])

            plt.scatter(self.lonc[v], self.latc[v],
                        s=80, label=lab, c=plt.cm.Set1(i/zz))

        #plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=2, ncol=3,fontsize='14', borderaxespad=0.)
        plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=2, ncol=3)
        #plt.legend()
        if plot:
            plt.show()

    def get_SMC_plot(age):
        sfr = np.array([])
        for i in np.arange(len(smc_coor)):
            sfr = np.append(sfr, get_SFH(smc_coor["ra"][i], \
                            smc_coor["dec"][i], age, smc_coor, smc_sfh))

        plt.tricontourf(smc_coor["ra"], smc_coor["dec"], sfr)
        plt.title(str(int(age)) + ' Myr')

        return plt

def contourf(*arguments, **kwargs):
    """Call signatures::

    contourf(X, Y, C, N, **kwargs)
    contourf(X, Y, C, V, **kwargs)
    
    Create a contourf plot of a 2-D llc array (with tricontour).
    
    *C* is the array of color values.

    *N* is the number of levels

    *V* is a list of levels
    
    *X* and *Y*, specify the (*x*, *y*) coordinates of
    the grid points

    **kwargs are passed to tricontour.
    
    """

    arglen = len(arguments)
    h = []
    if arglen >= 3:
        data = np.copy(arguments[2].flatten())
        x = arguments[0].flatten()
        y = arguments[1].flatten()

        # Create the Triangulation; 
        # no triangles so Delaunay triangulation created. 
        triang = tri.Triangulation(x, y)
        ntri = triang.triangles.shape[0]

        # Mask off unwanted triangles.
        mask = np.where(data[triang.triangles].prod(axis=1)==0., 1, 0)
        triang.set_mask(mask)
            
        if arglen == 3:
            h = plt.tricontourf(triang, data, **kwargs)
        elif arglen == 4:
            h = plt.tricontourf(triang, data, arguments[3], **kwargs)
        else:
            print("wrong number of arguments")
            print("need at least 3 or 4 arguments")
            sys.exit(__doc__)

        # show the triangles for debugging
        #plt.triplot(triang, color='0.7')

    else:
        print("wrong number of arguments")
        print("need at least x,y,fld")
        sys.exit(__doc__)

    return h

 def isosurf(self,Z,titre,front=True):
     """ trace isosurface de Z sur le maillage G"""
     triang=tri.Triangulation(self.X[:,0],self.X[:,1],triangles=self.Tbc)
     plt.tricontourf(triang, Z)
     if front:
         ARF=self.arfront()
         for L in ARF:
             plt.plot(self.X[L,0],self.X[L,1],lw=2,color='k')
     plt.colorbar()
     plt.title(titre)
     return

    def get_LMC_plot(age):
        sfr = np.array([])
        for i in np.arange(len(lmc_coor)):
            sfr = np.append(sfr, get_SFH(lmc_coor["ra"][i], \
                            lmc_coor["dec"][i], age, lmc_coor, lmc_sfh))

        plt.tricontourf(lmc_coor["ra"], lmc_coor["dec"], sfr)
        plt.title(str(int(age)) + ' Myr')
        plt.ylim(-73, -64)

        return plt

def draw_pdf_contours(dist, nlevels=200, subdiv=8, **kwargs):
    import math

    refiner = tri.UniformTriRefiner(triangle)
    trimesh = refiner.refine_triangulation(subdiv=subdiv)
    pvals = [dist.pdf(xy2bc(xy)) for xy in zip(trimesh.x, trimesh.y)]

    plt.tricontourf(trimesh, pvals, nlevels, **kwargs)
    plt.axis('equal')
    plt.xlim(0, 1)
    plt.ylim(0, 0.75**0.5)
    plt.axis('off')

def _add_obj_contour(x,y,color,label,columns,data,theta_star):
    ax = plt.gca()
    xvar, yvar, loc = _get_variables(ax,columns)
    try:
        X, Y, Z = _get_data_slice(xvar,yvar,columns,data,theta_star)
        
        triang = tri.Triangulation(X, Y)
        cmap = plt.cm.get_cmap('Greys')
        
        plt.tricontourf(triang,Z,cmap=cmap)
    except:
        print('Objective contour plot for', xvar, yvar,'slice failed')

def PlotAutocorrelations(system,shifts=[],update_datafiles=True,all_windows=True):
  outdir=_OutputDir(system)
  if update_datafiles:system.UpdateDataFiles(new_only=False)
  for w in system.windows:
    if (not all_windows) and hasattr(w,"auto_correlation_times"):continue
    t,cvs=w.ReadDataFile()
    ndata=len(t)
    if not shifts:
      dts=[]
      for i in range(ndata):
        dt=2**i
        if dt<ndata:dts.append(dt)
        else:break
    else:
      dts=[el for el in shifts if el<ndata]
    auto_corr_times=[]
    for i,cv in enumerate(cvs):
      cv=npy.array(cv)
      cl=[]
      for dt in dts:cl.append(npy.corrcoef(cv[:-dt],cv[dt:])[0,1])
      bools=npy.array(cl)<0.1
      auto_corr_time=dts[-1]
      for j in range(len(dts)-1,-1,-1):
        if not bools[j]:break
        auto_corr_time=dts[j]
      plt.figure()
      plt.plot(dts,cl,'-o')
      plt.xlim([1,dts[-1]+1])
      plt.xscale("log")
      plt.xlabel("Time shift")
      plt.ylabel("Autocorrelation")   
      plt.title("Autocorr for {0} for window {1}:{2}".format(system.cv_list[i].name,w.name,auto_corr_time))
      plt.savefig(os.path.join(outdir,"{0}_{1}_autocorr.png".format(w.name,system.cv_list[i].name)))
      plt.close()
      auto_corr_times.append(auto_corr_time)
    w.auto_correlation_times=auto_corr_times
  cv1=[]
  cv2=[]
  auto_corr=[]
  for w in system.windows:
    cv1.append(w.cv_values[0])
    cv2.append(w.cv_values[1])
    auto_corr.append(npy.log10(npy.max(w.auto_correlation_times)))
  plt.figure()
  plt.tricontourf(cv1,cv2,auto_corr, 20)
  plt.xlabel(system.cv_list[0].name+" "+system.cv_list[0].units)
  plt.ylabel(system.cv_list[1].name+" "+system.cv_list[1].units)
  cbar=plt.colorbar()
  cbar.set_label("log10 decorrelation time", rotation=270)
  plt.title("Decorrelation times")
  plt.savefig(os.path.join(outdir,"decorrelation_times.png"))
  plt.close()
  return

    def graphGrid(self,narrowGrid=False, plot=False):
        ''' A method to graph the grid with the shortest path plotted on the
        grid. narrowGrid will limit the field of view down to only show the
        paths drawn and not the entire grid. If only one path is drawn, this
        can skew the way the grid looks, and is sometime better to view the
        whole grid and zoom in. The plot option is in cause you want to choose
        when to plot the graph in a different part of the code.'''
        lat = self.lat
        lon = self.lon
        nv = self.nv.T - 1
        h = self.h

        tri = Tri.Triangulation(lon, lat, triangles=nv)  # xy or latlon based on how you are #Grand Passage

        levels=np.arange(-38,6,1)   # depth contours to plot

        fig = plt.figure(figsize=(18,10))
        plt.rc('font',size='22')
        ax = fig.add_subplot(111,aspect=(1.0/np.cos(np.mean(lat)*np.pi/180.0)))
        plt.tricontourf(tri, -h,levels=levels,shading='faceted',cmap=plt.cm.gist_earth)
        plt.triplot(tri)
        plt.ylabel('Latitude')
        plt.xlabel('Longitude')
        plt.gca().patch.set_facecolor('0.5')
        cbar=plt.colorbar()
        cbar.set_label('Water Depth (m)', rotation=-90,labelpad=30)

        scale = 1
        ticks = ticker.FuncFormatter(lambda lon, pos: '{0:g}'.format(lon/scale))
        ax.xaxis.set_major_formatter(ticks)
        ax.yaxis.set_major_formatter(ticks)
        plt.grid()

        maxlat, maxlon = np.max(self.maxcoordinates,axis=0)
        minlat, minlon = np.min(self.mincoordinates,axis=0)
        if narrowGrid:
            ax.set_xlim(minlon,maxlon)
            ax.set_ylim(minlat,maxlat)


        zz = len(self.elements)
        for i,v in enumerate(self.elements):
            source = self.pointIDXY[v[0]]
            target = self.pointIDXY[v[-1]]
            lab = '({:.6},{:.6})-({:.6},{:.6})'.format(source[0], source[1],
                                                       target[0], target[1])

            plt.scatter(self.lonc[v], self.latc[v],
                        s=80, label=lab, c=plt.cm.Set1(i/zz))

        plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=2, ncol=3)
        if plot:
            plt.show()