docs: def _krige(X, y, coords, variogram_function,
             variogram_model_parameters, coordinates_type):
      """Sets up and solves the ordinary kriging system for the given
      coordinate pair. This function is only used for the statistics calculations.
  
      Parameters
      ----------
      X: ndarray
          float array [n_samples, n_dim], the input array of coordinates
      y: ndarray
          float array [n_samples], the input array of measurement values
      coords: ndarray
          float array [1, n_dim], point at which to evaluate the kriging system
      variogram_function: callable
          function that will be called to evaluate variogram model
      variogram_model_parameters: list
          user-specified parameters for variogram model
      coordinates_type: str
          type of coordinates in X array, can be 'euclidean' for standard
          rectangular coordinates or 'geographic' if the coordinates are lat/lon
  
      Returns
      -------
      zinterp: float
          kriging estimate at the specified point
      sigmasq: float
          mean square error of the kriging estimate
      """
  
      zero_index = None
      zero_value = False
  
      # calculate distance between points... need a square distance matrix
      # of inter-measurement-point distances and a vector of distances between
      # measurement points (X) and the kriging point (coords)
      if coordinates_type == 'euclidean':
          d = squareform(pdist(X, metric='euclidean'))
          bd = np.squeeze(cdist(X, coords[None, :], metric='euclidean'))
  
      # geographic coordinate distances still calculated in the old way...
      # assume X[:, 0] ('x') => lon, X[:, 1] ('y') => lat
      # also assume problem is 2D; check done earlier in initializing variogram
      elif coordinates_type == 'geographic':
          x1, x2 = np.meshgrid(X[:, 0], X[:, 0], sparse=True)
          y1, y2 = np.meshgrid(X[:, 1], X[:, 1], sparse=True)
          d = great_circle_distance(x1, y1, x2, y2)
          bd = great_circle_distance(X[:, 0], X[:, 1],
                                     coords[0] * np.ones(X.shape[0]),
                                     coords[1] * np.ones(X.shape[0]))
  
      # this check is done when initializing variogram, but kept here anyways...
      else:
          raise ValueError("Specified coordinate type '%s' "
                           "is not supported." % coordinates_type)
  
      # check if kriging point overlaps with measurement point
      if np.any(np.absolute(bd) <= 1e-10):
          zero_value = True
          zero_index = np.where(bd <= 1e-10)[0][0]
  
      # set up kriging matrix
      n = X.shape[0]
      a = np.zeros((n+1, n+1))
      a[:n, :n] = - variogram_function(variogram_model_parameters, d)
      np.fill_diagonal(a, 0.0)
      a[n, :] = 1.0
      a[:, n] = 1.0
      a[n, n] = 0.0
  
      # set up RHS
      b = np.zeros((n+1, 1))
      b[:n, 0] = - variogram_function(variogram_model_parameters, bd)
      if zero_value:
          b[zero_index, 0] = 0.0
      b[n, 0] = 1.0
  
      # solve
      res = np.linalg.solve(a, b)
      zinterp = np.sum(res[:n, 0] * y)
      sigmasq = np.sum(res[:, 0] * -b[:, 0])
  
      return zinterp, sigmasq