Base class to handle nearest neighbour, linear, and cubic interpolation.
        Given a triangulation of a set of nodes on the unit sphere, along with data
        values at the nodes, this method interpolates (or extrapolates) the value
        at a given longitude and latitude.

        Parameters
        ----------
         lons : float / array of floats, shape (l,)
                longitudinal coordinate(s) on the sphere
         lats : float / array of floats, shape (l,)
                latitudinal coordinate(s) on the sphere
         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,)
            interpolated value(s) at (lons, lats)
         err : int / array of ints, shape (l,)
            whether interpolation (0), extrapolation (1) or error (other)

def interpolate(self, lons, lats, zdata, order=1):
        """
        Base class to handle nearest neighbour, linear, and cubic interpolation.
        Given a triangulation of a set of nodes on the unit sphere, along with data
        values at the nodes, this method interpolates (or extrapolates) the value
        at a given longitude and latitude.

        Parameters
        ----------
         lons : float / array of floats, shape (l,)
                longitudinal coordinate(s) on the sphere
         lats : float / array of floats, shape (l,)
                latitudinal coordinate(s) on the sphere
         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,)
            interpolated value(s) at (lons, lats)
         err : int / array of ints, shape (l,)
            whether interpolation (0), extrapolation (1) or error (other)

        """


        shape = lons.shape

        lons, lats = self._check_integrity(lons, lats)

        if order not in [0,1,3]:
            raise ValueError("order must be 0, 1, or 3")
        if zdata.size != self.npoints:
            raise ValueError("data must be of size {}".format(self.npoints))

        zdata = self._shuffle_field(zdata)

        zi, zierr, ierr = _stripack.interp_n(order, lats, lons,\
                                      self._x, self._y, self._z, zdata,\
                                      self.lst, self.lptr, self.lend)

        if ierr != 0:
            print('Warning some points may have errors - check error array\n'.format(ierr))
            zi[zierr < 0] = np.nan

        return zi.reshape(shape), zierr.reshape(shape)