Python matplotlib.tri.Triangulation() Examples

def test_triplot_return():
    # Check that triplot returns the artists it adds
    from matplotlib.figure import Figure
    ax = Figure().add_axes([0.1, 0.1, 0.7, 0.7])
    triang = mtri.Triangulation(
        [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0],
        triangles=[[0, 1, 3], [3, 2, 0]])
    assert ax.triplot(triang, "b-") is not None, \
        'triplot should return the artist it adds' 

def _get_tri(self):
        import matplotlib.tri as tri
        
        _x,_y = [],[]
        # ~ df = 1
        for n in self.get_nodes():
            _x.append(n.x)
            # ~ _x.append(n.x + n.ux*df)
            _y.append(n.y)
            # ~ _y.append(n.y + n.uy*df)
            
        tg = []
        for e in self.get_elements():
            ni,nj,nm = e.get_nodes()
            tg.append([ni.label, nj.label, nm.label])
            
        tr = tri.Triangulation(_x,_y, triangles=tg)
        return tr 

def test_tripcolor():
    x = np.asarray([0, 0.5, 1, 0,   0.5, 1,   0, 0.5, 1, 0.75])
    y = np.asarray([0, 0,   0, 0.5, 0.5, 0.5, 1, 1,   1, 0.75])
    triangles = np.asarray([
        [0, 1, 3], [1, 4, 3],
        [1, 2, 4], [2, 5, 4],
        [3, 4, 6], [4, 7, 6],
        [4, 5, 9], [7, 4, 9], [8, 7, 9], [5, 8, 9]])

    # Triangulation with same number of points and triangles.
    triang = mtri.Triangulation(x, y, triangles)

    Cpoints = x + 0.5*y

    xmid = x[triang.triangles].mean(axis=1)
    ymid = y[triang.triangles].mean(axis=1)
    Cfaces = 0.5*xmid + ymid

    plt.subplot(121)
    plt.tripcolor(triang, Cpoints, edgecolors='k')
    plt.title('point colors')

    plt.subplot(122)
    plt.tripcolor(triang, facecolors=Cfaces, edgecolors='k')
    plt.title('facecolors') 

def triangulate_bary(bary):
    """
    Triangulate a single barycentric triangle using matplotlib.

    Argument
    --------
    bary: barycentric matrix obtained using get_barymat.

    Return
    ------
    dely.edges: array (nedges, 2) that contains the indices of the two vertices
                 that form each edge after the triangulation.
    dely.triangles:array (ntriangles, 3) that contains the indices of the three
                   vertices that form each triangle after the triangulation. 
    """

    x = numpy.cos(-numpy.pi/4.)*bary[:, 0] + numpy.sin(-numpy.pi/4.)*bary[:, 1]
    y = bary[:, 2]

    dely = Triang.Triangulation(x, y)
    return dely.edges, dely.triangles 

def test_triinterpcubic_geom_weights():
    # Tests to check computation of weights for _DOF_estimator_geom:
    # The weight sum per triangle can be 1. (in case all angles < 90 degrees)
    # or (2*w_i) where w_i = 1-alpha_i/np.pi is the weight of apex i ; alpha_i
    # is the apex angle > 90 degrees.
    (ax, ay) = (0., 1.687)
    x = np.array([ax, 0.5*ax, 0., 1.])
    y = np.array([ay, -ay, 0., 0.])
    z = np.zeros(4, dtype=np.float64)
    triangles = [[0, 2, 3], [1, 3, 2]]
    sum_w = np.zeros([4, 2])  # 4 possibilities ; 2 triangles
    for theta in np.linspace(0., 2*np.pi, 14):  # rotating the figure...
        x_rot = np.cos(theta)*x + np.sin(theta)*y
        y_rot = -np.sin(theta)*x + np.cos(theta)*y
        triang = mtri.Triangulation(x_rot, y_rot, triangles)
        cubic_geom = mtri.CubicTriInterpolator(triang, z, kind='geom')
        dof_estimator = mtri.triinterpolate._DOF_estimator_geom(cubic_geom)
        weights = dof_estimator.compute_geom_weights()
        # Testing for the 4 possibilities...
        sum_w[0, :] = np.sum(weights, 1) - 1
        for itri in range(3):
            sum_w[itri+1, :] = np.sum(weights, 1) - 2*weights[:, itri]
        assert_array_almost_equal(np.min(np.abs(sum_w), axis=0),
                                  np.array([0., 0.], dtype=np.float64)) 

def test_qhull_triangle_orientation():
    # github issue 4437.
    xi = np.linspace(-2, 2, 100)
    x, y = map(np.ravel, np.meshgrid(xi, xi))
    w = np.logical_and(x > y - 1, np.logical_and(x < -1.95, y > -1.2))
    x, y = x[w], y[w]
    theta = np.radians(25)
    x1 = x*np.cos(theta) - y*np.sin(theta)
    y1 = x*np.sin(theta) + y*np.cos(theta)

    # Calculate Delaunay triangulation using Qhull.
    triang = mtri.Triangulation(x1, y1)

    # Neighbors returned by Qhull.
    qhull_neighbors = triang.neighbors

    # Obtain neighbors using own C++ calculation.
    triang._neighbors = None
    own_neighbors = triang.neighbors

    assert_array_equal(qhull_neighbors, own_neighbors) 

def test_delaunay_duplicate_points():
    # x[duplicate] == x[duplicate_of]
    # y[duplicate] == y[duplicate_of]
    npoints = 10
    duplicate = 7
    duplicate_of = 3

    np.random.seed(23)
    x = np.random.random((npoints))
    y = np.random.random((npoints))
    x[duplicate] = x[duplicate_of]
    y[duplicate] = y[duplicate_of]

    # Create delaunay triangulation.
    triang = mtri.Triangulation(x, y)

    # Duplicate points should be ignored, so the index of the duplicate points
    # should not appear in any triangle.
    assert_array_equal(np.unique(triang.triangles),
                       np.delete(np.arange(npoints), duplicate)) 

def test_trirefiner_fortran_contiguous_triangles():
    # github issue 4180.  Test requires two arrays of triangles that are
    # identical except that one is C-contiguous and one is fortran-contiguous.
    triangles1 = np.array([[2, 0, 3], [2, 1, 0]])
    assert not np.isfortran(triangles1)

    triangles2 = np.array(triangles1, copy=True, order='F')
    assert np.isfortran(triangles2)

    x = np.array([0.39, 0.59, 0.43, 0.32])
    y = np.array([33.99, 34.01, 34.19, 34.18])
    triang1 = mtri.Triangulation(x, y, triangles1)
    triang2 = mtri.Triangulation(x, y, triangles2)

    refiner1 = mtri.UniformTriRefiner(triang1)
    refiner2 = mtri.UniformTriRefiner(triang2)

    fine_triang1 = refiner1.refine_triangulation(subdiv=1)
    fine_triang2 = refiner2.refine_triangulation(subdiv=1)

    assert_array_equal(fine_triang1.triangles, fine_triang2.triangles) 

def test_tripcolor():
    x = np.asarray([0, 0.5, 1, 0,   0.5, 1,   0, 0.5, 1, 0.75])
    y = np.asarray([0, 0,   0, 0.5, 0.5, 0.5, 1, 1,   1, 0.75])
    triangles = np.asarray([
        [0, 1, 3], [1, 4, 3],
        [1, 2, 4], [2, 5, 4],
        [3, 4, 6], [4, 7, 6],
        [4, 5, 9], [7, 4, 9], [8, 7, 9], [5, 8, 9]])

    # Triangulation with same number of points and triangles.
    triang = mtri.Triangulation(x, y, triangles)

    Cpoints = x + 0.5*y

    xmid = x[triang.triangles].mean(axis=1)
    ymid = y[triang.triangles].mean(axis=1)
    Cfaces = 0.5*xmid + ymid

    plt.subplot(121)
    plt.tripcolor(triang, Cpoints, edgecolors='k')
    plt.title('point colors')

    plt.subplot(122)
    plt.tripcolor(triang, facecolors=Cfaces, edgecolors='k')
    plt.title('facecolors') 

def test_triinterpcubic_geom_weights():
    # Tests to check computation of weights for _DOF_estimator_geom:
    # The weight sum per triangle can be 1. (in case all angles < 90 degrees)
    # or (2*w_i) where w_i = 1-alpha_i/np.pi is the weight of apex i; alpha_i
    # is the apex angle > 90 degrees.
    (ax, ay) = (0., 1.687)
    x = np.array([ax, 0.5*ax, 0., 1.])
    y = np.array([ay, -ay, 0., 0.])
    z = np.zeros(4, dtype=np.float64)
    triangles = [[0, 2, 3], [1, 3, 2]]
    sum_w = np.zeros([4, 2])  # 4 possibilities; 2 triangles
    for theta in np.linspace(0., 2*np.pi, 14):  # rotating the figure...
        x_rot = np.cos(theta)*x + np.sin(theta)*y
        y_rot = -np.sin(theta)*x + np.cos(theta)*y
        triang = mtri.Triangulation(x_rot, y_rot, triangles)
        cubic_geom = mtri.CubicTriInterpolator(triang, z, kind='geom')
        dof_estimator = mtri.triinterpolate._DOF_estimator_geom(cubic_geom)
        weights = dof_estimator.compute_geom_weights()
        # Testing for the 4 possibilities...
        sum_w[0, :] = np.sum(weights, 1) - 1
        for itri in range(3):
            sum_w[itri+1, :] = np.sum(weights, 1) - 2*weights[:, itri]
        assert_array_almost_equal(np.min(np.abs(sum_w), axis=0),
                                  np.array([0., 0.], dtype=np.float64)) 

def test_triplot_return():
    # Check that triplot returns the artists it adds
    from matplotlib.figure import Figure
    ax = Figure().add_axes([0.1, 0.1, 0.7, 0.7])
    triang = mtri.Triangulation(
        [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0],
        triangles=[[0, 1, 3], [3, 2, 0]])
    if ax.triplot(triang, "b-") is None:
        raise AssertionError("triplot should return the artist it adds") 

def test_triplot_return():
    # Check that triplot returns the artists it adds
    from matplotlib.figure import Figure
    ax = Figure().add_axes([0.1, 0.1, 0.7, 0.7])
    triang = mtri.Triangulation(
        [0.0, 1.0, 0.0, 1.0], [0.0, 0.0, 1.0, 1.0],
        triangles=[[0, 1, 3], [3, 2, 0]])
    assert ax.triplot(triang, "b-") is not None, \
        'triplot should return the artist it adds' 

def test_qhull_triangle_orientation():
    # github issue 4437.
    xi = np.linspace(-2, 2, 100)
    x, y = map(np.ravel, np.meshgrid(xi, xi))
    w = np.logical_and(x > y - 1, np.logical_and(x < -1.95, y > -1.2))
    x, y = x[w], y[w]
    theta = np.radians(25)
    x1 = x*np.cos(theta) - y*np.sin(theta)
    y1 = x*np.sin(theta) + y*np.cos(theta)

    # Calculate Delaunay triangulation using Qhull.
    triang = mtri.Triangulation(x1, y1)

    # Neighbors returned by Qhull.
    qhull_neighbors = triang.neighbors

    # Obtain neighbors using own C++ calculation.
    triang._neighbors = None
    own_neighbors = triang.neighbors

    assert_array_equal(qhull_neighbors, own_neighbors) 

def test_triinterpcubic_geom_weights():
    # Tests to check computation of weights for _DOF_estimator_geom:
    # The weight sum per triangle can be 1. (in case all angles < 90 degrees)
    # or (2*w_i) where w_i = 1-alpha_i/np.pi is the weight of apex i; alpha_i
    # is the apex angle > 90 degrees.
    (ax, ay) = (0., 1.687)
    x = np.array([ax, 0.5*ax, 0., 1.])
    y = np.array([ay, -ay, 0., 0.])
    z = np.zeros(4, dtype=np.float64)
    triangles = [[0, 2, 3], [1, 3, 2]]
    sum_w = np.zeros([4, 2])  # 4 possibilities; 2 triangles
    for theta in np.linspace(0., 2*np.pi, 14):  # rotating the figure...
        x_rot = np.cos(theta)*x + np.sin(theta)*y
        y_rot = -np.sin(theta)*x + np.cos(theta)*y
        triang = mtri.Triangulation(x_rot, y_rot, triangles)
        cubic_geom = mtri.CubicTriInterpolator(triang, z, kind='geom')
        dof_estimator = mtri.triinterpolate._DOF_estimator_geom(cubic_geom)
        weights = dof_estimator.compute_geom_weights()
        # Testing for the 4 possibilities...
        sum_w[0, :] = np.sum(weights, 1) - 1
        for itri in range(3):
            sum_w[itri+1, :] = np.sum(weights, 1) - 2*weights[:, itri]
        assert_array_almost_equal(np.min(np.abs(sum_w), axis=0),
                                  np.array([0., 0.], dtype=np.float64)) 

def test_tripcolor():
    x = np.asarray([0, 0.5, 1, 0,   0.5, 1,   0, 0.5, 1, 0.75])
    y = np.asarray([0, 0,   0, 0.5, 0.5, 0.5, 1, 1,   1, 0.75])
    triangles = np.asarray([
        [0, 1, 3], [1, 4, 3],
        [1, 2, 4], [2, 5, 4],
        [3, 4, 6], [4, 7, 6],
        [4, 5, 9], [7, 4, 9], [8, 7, 9], [5, 8, 9]])

    # Triangulation with same number of points and triangles.
    triang = mtri.Triangulation(x, y, triangles)

    Cpoints = x + 0.5*y

    xmid = x[triang.triangles].mean(axis=1)
    ymid = y[triang.triangles].mean(axis=1)
    Cfaces = 0.5*xmid + ymid

    plt.subplot(121)
    plt.tripcolor(triang, Cpoints, edgecolors='k')
    plt.title('point colors')

    plt.subplot(122)
    plt.tripcolor(triang, facecolors=Cfaces, edgecolors='k')
    plt.title('facecolors') 

def test_delaunay_duplicate_points():
    # x[duplicate] == x[duplicate_of]
    # y[duplicate] == y[duplicate_of]
    npoints = 10
    duplicate = 7
    duplicate_of = 3

    np.random.seed(23)
    x = np.random.random((npoints))
    y = np.random.random((npoints))
    x[duplicate] = x[duplicate_of]
    y[duplicate] = y[duplicate_of]

    # Create delaunay triangulation.
    triang = mtri.Triangulation(x, y)

    # Duplicate points should be ignored, so the index of the duplicate points
    # should not appear in any triangle.
    assert_array_equal(np.unique(triang.triangles),
                       np.delete(np.arange(npoints), duplicate)) 

def test_delaunay_insufficient_points(x, y):
    with pytest.raises(ValueError):
        mtri.Triangulation(x, y) 

def test_qhull_large_offset():
    # github issue 8682.
    x = np.asarray([0, 1, 0, 1, 0.5])
    y = np.asarray([0, 0, 1, 1, 0.5])

    offset = 1e10
    triang = mtri.Triangulation(x, y)
    triang_offset = mtri.Triangulation(x + offset, y + offset)
    assert len(triang.triangles) == len(triang_offset.triangles) 

def test_no_modify():
    # Test that Triangulation does not modify triangles array passed to it.
    triangles = np.array([[3, 2, 0], [3, 1, 0]], dtype=np.int32)
    points = np.array([(0, 0), (0, 1.1), (1, 0), (1, 1)])

    old_triangles = triangles.copy()
    tri = mtri.Triangulation(points[:, 0], points[:, 1], triangles)
    edges = tri.edges
    assert_array_equal(old_triangles, triangles) 

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. * np.cos(7.*theta1) +
              (np.exp((r2/10)**2)-1)*30. * np.cos(11.*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)
    triang0.set_mask(np.hypot(x0[triang0.triangles].mean(axis=1),
                              y0[triang0.triangles].mean(axis=1))
                     < min_radius)

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

def __init__(self, triangulation):
        if not isinstance(triangulation, Triangulation):
            raise ValueError("Expected a Triangulation object")
        self._triangulation = triangulation 

def __init__(self, triangulation):
        if not isinstance(triangulation, Triangulation):
            raise ValueError("Expected a Triangulation object")
        self._triangulation = triangulation 

def test_delaunay():
    # No duplicate points, regular grid.
    nx = 5
    ny = 4
    x, y = np.meshgrid(np.linspace(0.0, 1.0, nx), np.linspace(0.0, 1.0, ny))
    x = x.ravel()
    y = y.ravel()
    npoints = nx*ny
    ntriangles = 2 * (nx-1) * (ny-1)
    nedges = 3*nx*ny - 2*nx - 2*ny + 1

    # Create delaunay triangulation.
    triang = mtri.Triangulation(x, y)

    # The tests in the remainder of this function should be passed by any
    # triangulation that does not contain duplicate points.

    # Points - floating point.
    assert_array_almost_equal(triang.x, x)
    assert_array_almost_equal(triang.y, y)

    # Triangles - integers.
    assert len(triang.triangles) == ntriangles
    assert np.min(triang.triangles) == 0
    assert np.max(triang.triangles) == npoints-1

    # Edges - integers.
    assert len(triang.edges) == nedges
    assert np.min(triang.edges) == 0
    assert np.max(triang.edges) == npoints-1

    # Neighbors - integers.
    # Check that neighbors calculated by C++ triangulation class are the same
    # as those returned from delaunay routine.
    neighbors = triang.neighbors
    triang._neighbors = None
    assert_array_equal(triang.neighbors, neighbors)

    # Is each point used in at least one triangle?
    assert_array_equal(np.unique(triang.triangles), np.arange(npoints)) 

def test_trianalyzer_mismatched_indices():
    # github issue 4999.
    x = np.array([0., 1., 0.5, 0., 2.])
    y = np.array([0., 0., 0.5*np.sqrt(3.), -1., 1.])
    triangles = np.array([[0, 1, 2], [0, 1, 3], [1, 2, 4]], dtype=np.int32)
    mask = np.array([False, False, True], dtype=bool)
    triang = mtri.Triangulation(x, y, triangles, mask=mask)
    analyser = mtri.TriAnalyzer(triang)
    # numpy >= 1.10 raises a VisibleDeprecationWarning in the following line
    # prior to the fix.
    triang2 = analyser._get_compressed_triangulation() 

def threshold_map(self, metric):
        lin = np.linspace(0, 1, 64)

        triang = tri.Triangulation(self.results.threshold1.values, self.results.threshold2.values)
        interpolator = tri.LinearTriInterpolator(triang, self.results[metric])
        Xi, Yi = np.meshgrid(lin, lin)
        zi = interpolator(Xi, Yi)
        plt.figure(figsize=(6, 6))
        img = plt.imshow(zi[::-1], extent=[0, 1, 0, 1])
        plt.colorbar(img)
        plt.xlabel("threshold1")
        plt.ylabel("threshold2") 

def select_pairs_delaunay(date_list, tbase_list, pbase_list, norm=True, date12_format='YYMMDD-YYMMDD'):
    """Select Pairs using Delaunay Triangulation based on temporal/perpendicular baselines
    Inputs:
        date_list  : list of date in YYMMDD/YYYYMMDD format
        pbase_list : list of float, perpendicular spatial baseline
        norm       : normalize temporal baseline to perpendicular baseline
    Key points
        1. Define a ratio between perpendicular and temporal baseline axis units (Pepe and Lanari, 2006, TGRS).
        2. Pairs with too large perpendicular / temporal baseline or Doppler centroid difference should be removed
           after this, using a threshold, to avoid strong decorrelations (Zebker and Villasenor, 1992, TGRS).
    Reference:
        Pepe, A., and R. Lanari (2006), On the extension of the minimum cost flow algorithm for phase unwrapping
        of multitemporal differential SAR interferograms, IEEE TGRS, 44(9), 2374-2383.
        Zebker, H. A., and J. Villasenor (1992), Decorrelation in interferometric radar echoes, IEEE TGRS, 30(5), 950-959.
    """
    # Get temporal baseline in days
    date6_list = yymmdd(date_list)
    date8_list = yyyymmdd(date_list)
    #tbase_list = date_list2tbase(date8_list)[0]

    # Normalization (Pepe and Lanari, 2006, TGRS)
    if norm:
        temp2perp_scale = (max(pbase_list)-min(pbase_list)) / (max(tbase_list)-min(tbase_list))
        tbase_list = [tbase*temp2perp_scale for tbase in tbase_list]

    # Generate Delaunay Triangulation
    date12_idx_list = Triangulation(tbase_list, pbase_list).edges.tolist()
    date12_idx_list = [sorted(idx) for idx in sorted(date12_idx_list)]

    # Convert index into date12
    date12_list = [date6_list[idx[0]]+'-'+date6_list[idx[1]]
                   for idx in date12_idx_list]
    if date12_format == 'YYYYMMDD_YYYYMMDD':
        date12_list = yyyymmdd_date12(date12_list)
    return date12_list 

def test_delaunay_insufficient_points(x, y):
    with pytest.raises(ValueError):
        mtri.Triangulation(x, y) 

def test_trianalyzer_mismatched_indices():
    # github issue 4999.
    x = np.array([0., 1., 0.5, 0., 2.])
    y = np.array([0., 0., 0.5*np.sqrt(3.), -1., 1.])
    triangles = np.array([[0, 1, 2], [0, 1, 3], [1, 2, 4]], dtype=np.int32)
    mask = np.array([False, False, True], dtype=bool)
    triang = mtri.Triangulation(x, y, triangles, mask=mask)
    analyser = mtri.TriAnalyzer(triang)
    # numpy >= 1.10 raises a VisibleDeprecationWarning in the following line
    # prior to the fix.
    triang2 = analyser._get_compressed_triangulation() 

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. * np.cos(7.*theta1) +
              (np.exp((r2/10)**2)-1)*30. * np.cos(11.*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)
    triang0.set_mask(np.hypot(x0[triang0.triangles].mean(axis=1),
                              y0[triang0.triangles].mean(axis=1))
                     < min_radius)

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

def __init__(self, triangulation, z, trifinder=None):
        if not isinstance(triangulation, Triangulation):
            raise ValueError("Expected a Triangulation object")
        self._triangulation = triangulation

        self._z = np.asarray(z)
        if self._z.shape != self._triangulation.x.shape:
            raise ValueError("z array must have same length as triangulation x"
                             " and y arrays")

        if trifinder is not None and not isinstance(trifinder, TriFinder):
            raise ValueError("Expected a TriFinder object")
        self._trifinder = trifinder or self._triangulation.get_trifinder()

        # Default scaling factors : 1.0 (= no scaling)
        # Scaling may be used for interpolations for which the order of
        # magnitude of x, y has an impact on the interpolant definition.
        # Please refer to :meth:`_interpolate_multikeys` for details.
        self._unit_x = 1.0
        self._unit_y = 1.0

        # Default triangle renumbering: None (= no renumbering)
        # Renumbering may be used to avoid unnecessary computations
        # if complex calculations are done inside the Interpolator.
        # Please refer to :meth:`_interpolate_multikeys` for details.
        self._tri_renum = None

    # __call__ and gradient docstrings are shared by all subclasses
    # (except, if needed, relevant additions).
    # However these methods are only implemented in subclasses to avoid
    # confusion in the documentation.