Query the cKDtree for the nearest neighbours and Euclidean
        distance from x,y points.

        Returns 0, 0 if a cKDtree has not been constructed
        (switch tree=True if you need this routine)

        Parameters
        ----------
         lon : 1D array of longitudinal coordinates in radians
         lat : 1D array of latitudinal coordinates in radians
         k   : number of nearest neighbours to return
             (default: 1)
         max_distance
             : maximum Euclidean distance to search
               for neighbours (default: 2.0)

        Returns
        -------
         d    : Euclidean distance between each point and their
                nearest neighbour(s)
         vert : vertices of the nearest neighbour(s)

def nearest_vertices(self, lon, lat, k=1, max_distance=2.0 ):
        """
        Query the cKDtree for the nearest neighbours and Euclidean
        distance from x,y points.

        Returns 0, 0 if a cKDtree has not been constructed
        (switch tree=True if you need this routine)

        Parameters
        ----------
         lon : 1D array of longitudinal coordinates in radians
         lat : 1D array of latitudinal coordinates in radians
         k   : number of nearest neighbours to return
             (default: 1)
         max_distance
             : maximum Euclidean distance to search
               for neighbours (default: 2.0)

        Returns
        -------
         d    : Euclidean distance between each point and their
                nearest neighbour(s)
         vert : vertices of the nearest neighbour(s)
        """

        if self.tree == False or self.tree == None:
            return 0, 0

        lons = np.array(lon).reshape(-1,1)
        lats = np.array(lat).reshape(-1,1)

        xyz = np.empty((lons.shape[0],3))
        x,y,z = lonlat2xyz(lons, lats)

        xyz[:,0] = x[:].reshape(-1)
        xyz[:,1] = y[:].reshape(-1)
        xyz[:,2] = z[:].reshape(-1)

        dxyz, vertices = self._cKDtree.query(xyz, k=k, distance_upper_bound=max_distance)


        if k == 1:   # force this to be a 2D array
            vertices = np.reshape(vertices, (-1, 1))

        ## Now find the angular separation / great circle distance: dlatlon


        vertxyz = self.points[vertices].transpose(0,2,1)
        extxyz  = np.repeat(xyz, k, axis=1).reshape(vertxyz.shape)

        angles = np.arccos((extxyz * vertxyz).sum(axis=1))

        return angles, vertices