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 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 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 test_tri_smooth_gradient():
    # Image comparison based on example trigradient_demo.

    def dipole_potential(x, y):
        """ An electric dipole potential V """
        r_sq = x ** 2 + y ** 2
        theta = np.arctan2(y, x)
        z = np.cos(theta) / r_sq
        return (np.max(z) - z) / (np.max(z) - np.min(z))

    # Creating a Triangulation
    n_angles = 30
    n_radii = 10
    min_radius = 0.2
    radii = np.linspace(min_radius, 0.95, n_radii)
    angles = np.linspace(0, 2 * np.pi, n_angles, endpoint=False)
    angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
    angles[:, 1::2] += np.pi / n_angles
    x = (radii * np.cos(angles)).flatten()
    y = (radii * np.sin(angles)).flatten()
    V = dipole_potential(x, y)
    triang = mtri.Triangulation(x, y)
    xmid = x[triang.triangles].mean(axis=1)
    ymid = y[triang.triangles].mean(axis=1)
    mask = np.where(xmid * xmid + ymid * ymid < min_radius * min_radius, 1, 0)
    triang.set_mask(mask)

    # Refine data - interpolates the electrical potential V
    refiner = mtri.UniformTriRefiner(triang)
    tri_refi, z_test_refi = refiner.refine_field(V, subdiv=3)

    # Computes the electrical field (Ex, Ey) as gradient of -V
    tci = mtri.CubicTriInterpolator(triang, -V)
    (Ex, Ey) = tci.gradient(triang.x, triang.y)
    E_norm = np.sqrt(Ex ** 2 + Ey ** 2)

    # Plot the triangulation, the potential iso-contours and the vector field
    plt.figure()
    plt.gca().set_aspect("equal")
    plt.triplot(triang, color="0.8")

    levels = np.arange(0.0, 1.0, 0.01)
    cmap = cm.get_cmap(name="hot", lut=None)
    plt.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap, linewidths=[2.0, 1.0, 1.0, 1.0])
    # Plots direction of the electrical vector field
    plt.quiver(
        triang.x,
        triang.y,
        Ex / E_norm,
        Ey / E_norm,
        units="xy",
        scale=10.0,
        zorder=3,
        color="blue",
        width=0.007,
        headwidth=3.0,
        headlength=4.0,
    )

def _add_LR_contour(x,y,color,label,columns,data,theta_star,threshold):
    ax = plt.gca()
    xvar, yvar, loc = _get_variables(ax,columns)
    
    X, Y, Z = _get_data_slice(xvar,yvar,columns,data,theta_star)
    
    triang = tri.Triangulation(X, Y)
    
    plt.tricontour(triang,Z,[threshold], colors='r')

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

    contour(X, Y, C, N, **kwargs)
    contour(X, Y, C, V, **kwargs)
    
    Create a contour 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 = 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.tricontour(triang, data, **kwargs)
        elif arglen == 4:
            h = plt.tricontour(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 plotResult(self, **k):
     x = []
     y = []
     z = []
     for n in self.nodes:
         x.append(n.x)
         y.append(n.y)
         z.append(k['r'][n.i-1])
     t = tri.Triangulation(x, y)
     plt.tricontour(t, z, 15, linewidths=0.5)
     #plt.tricontourf(t, z, 15, cmap=plt.cm.rainbow)
     plt.colorbar()
     #plt.plot(x, y, 'ko', ms=3)
     

 def stab_plot(self):
     x = self.x
     y = self.y
     z = self.z
     # initialize the delauney triangulation:
     triang = tri.Triangulation(x, y)
     # do the plots:
     plt.tricontourf(triang, z, colors=self.color, alpha=0.3, levels=[1.0, 2.0])
     plt.tricontour(triang, z, colors=self.color, levels=[1.0, 2.0])
     # title of the plot
     plt.title("Stability of " + self.model_name, fontsize=16)
     # labels
     plt.xlabel(self.param_names[0], fontsize=18)
     plt.ylabel(self.param_names[1], fontsize=18)

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 contour_paths(self, epsilon):
     '''Extract the polygon patches for the provided epsilon'''
     figure = pyplot.figure()
     contours = pyplot.tricontour(self.triang, self.vals, levels=[epsilon])
     paths = Paths()
     if len(contours.collections) == 0:
         return paths
     for path in contours.collections[0].get_paths():
         paths.append(Path(path.vertices[:, 0] + 1j*path.vertices[:, 1]))
     pyplot.close(figure)
     return paths

    def plot_contour(self, **kwargs):
        """Contour plot of the xy plane

        Parameters
        ----------
        **kwargs
            Forwarded to `plt.tricontour()`.
        """
        x, y, _ = self.pos
        contour = plt.tricontour(x, y, self.data, **kwargs)
        self._decorate_plot()
        return contour

def test_tri_smooth_contouring():
    # Image comparison based on example tricontour_smooth_user.
    n_angles = 20
    n_radii = 10
    min_radius = 0.15

    def z(x, y):
        r1 = np.sqrt((0.5 - x) ** 2 + (0.5 - y) ** 2)
        theta1 = np.arctan2(0.5 - x, 0.5 - y)
        r2 = np.sqrt((-x - 0.2) ** 2 + (-y - 0.2) ** 2)
        theta2 = np.arctan2(-x - 0.2, -y - 0.2)
        z = -(
            2 * (np.exp((r1 / 10) ** 2) - 1) * 30.0 * np.cos(7.0 * theta1)
            + (np.exp((r2 / 10) ** 2) - 1) * 30.0 * np.cos(11.0 * theta2)
            + 0.7 * (x ** 2 + y ** 2)
        )
        return (np.max(z) - z) / (np.max(z) - np.min(z))

    # First create the x and y coordinates of the points.
    radii = np.linspace(min_radius, 0.95, n_radii)
    angles = np.linspace(0 + n_angles, 2 * np.pi + n_angles, n_angles, endpoint=False)
    angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
    angles[:, 1::2] += np.pi / n_angles
    x0 = (radii * np.cos(angles)).flatten()
    y0 = (radii * np.sin(angles)).flatten()
    triang0 = mtri.Triangulation(x0, y0)  # Delaunay triangulation
    z0 = z(x0, y0)
    xmid = x0[triang0.triangles].mean(axis=1)
    ymid = y0[triang0.triangles].mean(axis=1)
    mask = np.where(xmid * xmid + ymid * ymid < min_radius * min_radius, 1, 0)
    triang0.set_mask(mask)

    # Then the plot
    refiner = mtri.UniformTriRefiner(triang0)
    tri_refi, z_test_refi = refiner.refine_field(z0, subdiv=4)
    levels = np.arange(0.0, 1.0, 0.025)
    plt.triplot(triang0, lw=0.5, color="0.5")
    plt.tricontour(tri_refi, z_test_refi, levels=levels, colors="black")

def tricontourf_deformed(a, mesh, d,  name, cmap, dpi, magf):
    """Plot contour with the tricoutour function and the boundary line with
    the boundary node.

    """
    fig = plt.figure(name, dpi=dpi)
    ax1 = fig.add_subplot(1, 1, 1, aspect='equal')
    c = mesh.nodes_coord
    bn = mesh.boundary_nodes

    xx, yy, zz = c[:, 0]+d[::2]*magf, c[:, 1]+d[1::2]*magf, a

    dx = d[::2]*magf
    dy = d[1::2]*magf

    ccx = np.append(c[bn[:, 1], 0]+dx[bn[:, 1]], c[bn[0, 1], 0]+dx[bn[0, 1]])
    ccy = np.append(c[bn[:, 1], 1]+dy[bn[:, 1]], c[bn[0, 1], 1]+dy[bn[0, 1]])
    plt.plot(ccx, ccy, '-k')

    triangles = []
    for n1, n2, n3, n4 in mesh.ele_conn:
        triangles.append([n1, n2, n3])
        triangles.append([n1, n3, n4])

    triangles = np.asarray(triangles)

    CS2 = plt.tricontourf(xx, yy, triangles, zz, N=10, origin='lower',
                          cmap=cmap)

    CS3 = plt.tricontour(xx, yy, triangles, zz, N=10, origin='lower',colors='k')


    plt.xlabel(r'$x$', fontsize=14)
    plt.ylabel(r'$y$', fontsize=14)

    divider = make_axes_locatable(ax1)
    cax = divider.append_axes("right", size=0.3, pad=0.1)

    cbar = plt.colorbar(CS2, cax=cax)
    cbar.ax.set_label(name, fontsize=12)

    plt.clabel(CS3, fontsize=8, colors='k', fmt='%1.1f')

    limits=plt.axis('off')
    # plt.savefig('1.png', transparent=True, dpi=300)
    #plt.axes().set_aspect('equal')
    #plt.axes().autoscale_view(True, True, True)

    plt.draw()

 def plot(self, epsilons, **kwargs):
     epsilons = list(numpy.sort(epsilons))
     padepsilons = [epsilons[0]*0.9] + epsilons + [epsilons[-1]*1.1]
     X = []
     Y = []
     Z = []
     for epsilon in padepsilons:
         paths = self.contour_paths(epsilon)
         for path in paths:
             X += list(numpy.real(path.vertices[:-1]))
             Y += list(numpy.imag(path.vertices[:-1]))
             Z += [epsilon] * (len(path.vertices) - 1)
     contours = pyplot.tricontour(X, Y, Z, levels=epsilons,
                                  colors=pyplot.rcParams['axes.color_cycle']
                                  )
     plot_finish(contours, **kwargs)
     return contours

  def contour(self, prop=None):
    if(prop=='Tdust'):
        data=self.T[0:len(self.x)]['1']
        label="Dust temperature [K]"
    elif(prop=='rho'):
        data=np.log(self.rho[0:len(self.x)]['1'])
        label="Log some density [?]"
    else:
        try:
            data=self.prop
            label='Unknown property'
        except:
            print "Please specify property to plot. Valid properties are [Tdust, rho]"
            return

    x_au = self.x/1.49e13
    y_au = self.y/1.49e13

    tri=Triangulation(x_au,y_au)

    plt.close('all')
    plt.figure()
    ax=plt.subplot(111)
    ax.set_xscale("log",nonposx='clip')
    ax.set_yscale("log",nonposy='clip')
    ax.set_xlim([0.1,200])
    ax.set_ylim([0.1,200])
    plt.xlabel('r [AU]')
    plt.ylabel('z [AU]')

    nmax=data.max()
    nmin=data.min()
    levels=np.arange(12) * (nmax-nmin)/12. + nmin
    plt.tricontourf(tri, data, levels)
    cbar=plt.colorbar()
    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)

 def computeContours(self):
     extend = str(self.uExtend.itemText(self.uExtend.currentIndex()))
     x, y, z = self._data
     levels = self.getLevels()
     if self._isMPLOk()==True: # If so, we can use the new tricontour fonction
         cs = plt.tricontour(x, y, z, levels, extend=extend)
     else:
         gx = x.reshape(self._nr,self._nc)
         gy = y.reshape(self._nr,self._nc)
         gz = z.reshape(self._nr,self._nc)
         cs = plt.contour(gx, gy, gz, levels, extend=extend)
     lines = list()
     levels = [float(l) for l in cs.levels]
     self.progressBar.setRange(0, len(cs.collections))
     for i, line in enumerate(cs.collections):
         linestrings = []
         for path in line.get_paths():
             if len(path.vertices) > 1:
                 linestrings.append(path.vertices)
         lines.append([ i, levels[i], linestrings])
         self.progressBar.setValue(i+1)
     self.progressBar.setValue(0)
     return lines

    def lineContourFeatures(self):
        x,y,z=self.data()
        levels = self.levels()
        usegrid=self.isGridded() and self._useGrid
        try:
            if usegrid:
                gx,gy,gz=self.gridContourData()
                cs = contour(gx, gy, gz, levels )
            else:
                trig,z=self.trigContourData()
                cs = tricontour(trig, z, levels )
        except:
            raise ContourGenerationError.fromException(sys.exc_info())

        fields = self.fields()
        zfield=self.zFieldName()
        dx,dy=self._origin
        for i, line in enumerate(cs.collections):
            level=float(cs.levels[i])
            glines = []
            try:
                for path in line.get_paths():
                    if len(path.vertices) > 1:
                        points=[QgsPointXY(x,y) for x,y in path.vertices]
                        glines.append(points)
                geom=QgsGeometry.fromMultiPolylineXY(glines)
                geom.translate(dx,dy)
                feat = QgsFeature(fields)
                feat.setGeometry(geom)
                feat['index']=i
                feat[zfield]=level
                feat['label']=self._levelLabel(level)
                yield feat
            except:
                message=sys.exc_info()[1]
                self._feedback.reportError(message)

Xmax, Ymax = X[area_max].flat[amax], Y[area_max].flat[amax]
lsteps_ = int(max(refi_magU[area_max])/lmax*lsteps)
lx = Xmin + np.sign(Xmax-Xmin)*np.linspace(0, abs(Xmax-Xmin)**(1./args.lpow), 1 + lsteps_)**args.lpow
ly = Ymin + np.sign(Ymax-Ymin)*np.linspace(0, abs(Ymax-Ymin)**(1./args.lpow), 1 + lsteps_)**args.lpow
#plt.plot(lx, ly, '*')
lx, ly = (lx[1:] + lx[:-1])/2, (ly[1:] + ly[:-1])/2
#plt.plot(lx, ly, 'D')
labelpos = tuple(map(tuple, np.vstack((lx,ly)).T))
linewidth, fontsize = .5, 6

### Plot detailed contours
refi_magU = np.maximum(refi_magU, lmin)
dl, ds = ((lmax-lmin)/lsteps/2,1) if np.max(magU) > lmax else (0,0)
levels = np.linspace(lmin, lmax + dl, 1 + 2*lsteps + ds)
if grayscale:
    plt.tricontour(refi_triang, refi_magU, levels=levels, colors='lightgray', linewidths=0.05)
else:
    cmap = cm.get_cmap('coolwarm')
    plt.tricontourf(refi_triang, refi_magU, cmap=cmap, levels=levels)

### Plot black contours with labels
levels = np.linspace(lmin, lmax, 1 + lsteps)
CS = plt.tricontour(refi_triang, refi_magU, levels=levels, colors='k', linewidths=.3)
clabels = plt.clabel(CS, levels,
    inline=True,
    use_clabeltext=True,
    fmt='%g',
    manual=labelpos,
    fontsize=6)
[ txt.set_backgroundcolor('white') for txt in clabels ]
[ txt.set_bbox(dict(facecolor='white', edgecolor='none', pad=0.5)) for txt in clabels ]

# Mask off unwanted triangles.
xmid = x[triang.triangles].mean(axis=1)
ymid = y[triang.triangles].mean(axis=1)
mask = np.where(xmid * xmid + ymid * ymid < min_radius * min_radius, 1, 0)
triang.set_mask(mask)

#-----------------------------------------------------------------------------
# Refine data
#-----------------------------------------------------------------------------
refiner = tri.UniformTriRefiner(triang)
tri_refi, z_test_refi = refiner.refine_field(z, subdiv=3)

#-----------------------------------------------------------------------------
# Plot the triangulation and the high-res iso-contours
#-----------------------------------------------------------------------------
plt.figure()
plt.gca().set_aspect('equal')
plt.triplot(triang, lw=0.5, color='white')

levels = np.arange(0., 1., 0.025)
cmap = cm.get_cmap(name='terrain', lut=None)
plt.tricontourf(tri_refi, z_test_refi, levels=levels, cmap=cmap)
plt.tricontour(tri_refi, z_test_refi, levels=levels,
               colors=['0.25', '0.5', '0.5', '0.5', '0.5'],
               linewidths=[1.0, 0.5, 0.5, 0.5, 0.5])

plt.title("High-resolution tricontouring")

plt.show()

start = time.clock()
plt.subplot(211)
xi = np.linspace(-2.1,2.1,ngridx)
yi = np.linspace(-2.1,2.1,ngridy)
zi = griddata(x,y,z,xi,yi,interp='linear')
plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
plt.contourf(xi,yi,zi,15,cmap=plt.cm.rainbow,
             norm=plt.normalize(vmax=abs(zi).max(), vmin=-abs(zi).max()))
plt.colorbar() # draw colorbar
plt.plot(x, y, 'ko', ms=3)
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.title('griddata and contour (%d points, %d grid points)' % (npts, ngridx*ngridy))
print ('griddata and contour seconds: %f' % (time.clock() - start))

# tricontour.
start = time.clock()
plt.subplot(212)
triang = tri.Triangulation(x, y)
plt.tricontour(x, y, z, 15, linewidths=0.5, colors='k')
plt.tricontourf(x, y, z, 15, cmap=plt.cm.rainbow,
                norm=plt.normalize(vmax=abs(zi).max(), vmin=-abs(zi).max()))
plt.colorbar()
plt.plot(x, y, 'ko', ms=3)
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.title('tricontour (%d points)' % npts)
print ('tricontour seconds: %f' % (time.clock() - start))

plt.show()