Base class to handle nearest neighbour, linear, and cubic interpolation.
        Given a triangulation of a set of nodes and values at the nodes,
        this method interpolates the value at the given xi,yi coordinates.

        Parameters
        ----------
         xi    : float / array of floats, shape (l,)
                 x Cartesian coordinate(s)
         yi    : float / array of floats, shape (l,)
                 y Cartesian coordinate(s)
         zdata : array of floats, shape (n,)
                 value at each point in the triangulation
                 must be the same size of the mesh
         order : int (default=1)
                 order of the interpolatory function used
                  0 = nearest-neighbour
                  1 = linear
                  3 = cubic

        Returns
        -------
         zi    : float / array of floats, shape (l,)
                interpolates value(s) at (xi, yi)
         err   : int / array of ints, shape (l,)
                whether interpolation (0), extrapolation (1) or error (other)

def interpolate(self, xi, yi, zdata, order=1):
        """
        Base class to handle nearest neighbour, linear, and cubic interpolation.
        Given a triangulation of a set of nodes and values at the nodes,
        this method interpolates the value at the given xi,yi coordinates.

        Parameters
        ----------
         xi    : float / array of floats, shape (l,)
                 x Cartesian coordinate(s)
         yi    : float / array of floats, shape (l,)
                 y Cartesian coordinate(s)
         zdata : array of floats, shape (n,)
                 value at each point in the triangulation
                 must be the same size of the mesh
         order : int (default=1)
                 order of the interpolatory function used
                  0 = nearest-neighbour
                  1 = linear
                  3 = cubic

        Returns
        -------
         zi    : float / array of floats, shape (l,)
                interpolates value(s) at (xi, yi)
         err   : int / array of ints, shape (l,)
                whether interpolation (0), extrapolation (1) or error (other)
        """

        if order == 0:
            zierr = np.zeros_like(xi, dtype=np.int)
            return self.interpolate_nearest(xi, yi, zdata), zierr
        elif order == 1:
            return self.interpolate_linear(xi, yi, zdata)
        elif order == 3:
            return self.interpolate_cubic(xi, yi, zdata)
        else:
            raise ValueError("order must be 0, 1, or 3")