code: def nearest_vertices(self, lon, lat, k=1, max_distance=2.0 ):
          
  
          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
docs: 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)