docs: def execute(self, style, xpoints, ypoints, mask=None, backend='vectorized',
                  n_closest_points=None):
          """Calculates a kriged grid and the associated variance.
  
          This is now the method that performs the main kriging calculation.
          Note that currently measurements (i.e., z values) are considered
          'exact'. This means that, when a specified coordinate for interpolation
          is exactly the same as one of the data points, the variogram evaluated
          at the point is forced to be zero. Also, the diagonal of the kriging
          matrix is also always forced to be zero. In forcing the variogram
          evaluated at data points to be zero, we are effectively saying that
          there is no variance at that point (no uncertainty,
          so the value is 'exact').
  
          In the future, the code may include an extra 'exact_values' boolean
          flag that can be adjusted to specify whether to treat the measurements
          as 'exact'. Setting the flag to false would indicate that the
          variogram should not be forced to be zero at zero distance
          (i.e., when evaluated at data points). Instead, the uncertainty in
          the point will be equal to the nugget. This would mean that the
          diagonal of the kriging matrix would be set to
          the nugget instead of to zero.
  
          Parameters
          ----------
          style : str
              Specifies how to treat input kriging points. Specifying 'grid'
              treats xpoints and ypoints as two arrays of x and y coordinates
              that define a rectangular grid. Specifying 'points' treats
              xpoints and ypoints as two arrays that provide coordinate pairs
              at which to solve the kriging system. Specifying 'masked'
              treats xpoints and ypoints as two arrays of x and y coordinates
              that define a rectangular grid and uses mask to only evaluate
              specific points in the grid.
          xpoints : array_like, shape (N,) or (N, 1)
              If style is specific as 'grid' or 'masked',
              x-coordinates of MxN grid. If style is specified as 'points',
              x-coordinates of specific points at which to solve
              kriging system.
          ypoints : array_like, shape (M,) or (M, 1)
              If style is specified as 'grid' or 'masked',
              y-coordinates of MxN grid. If style is specified as 'points',
              y-coordinates of specific points at which to solve kriging
              system. Note that in this case, xpoints and ypoints must have
              the same dimensions (i.e., M = N).
          mask : bool, array_like, shape (M, N), optional
              Specifies the points in the rectangular grid defined
              by xpoints and ypoints that are to be excluded in the
              kriging calculations. Must be provided if style is specified
              as 'masked'. False indicates that the point should not be
              masked, so the kriging system will be solved at the point.
              True indicates that the point should be masked, so the kriging
              system should will not be solved at the point.
          backend : str, optional
              Specifies which approach to use in kriging.
              Specifying 'vectorized' will solve the entire kriging problem
              at once in a vectorized operation. This approach is faster but
              also can consume a significant amount of memory for large grids
              and/or large datasets. Specifying 'loop' will loop through each
              point at which the kriging system is to be solved.
              This approach is slower but also less memory-intensive.
              Specifying 'C' will utilize a loop in Cython.
              Default is 'vectorized'.
          n_closest_points : int, optional
              For kriging with a moving window, specifies the number of
              nearby points to use in the calculation. This can speed up the
              calculation for large datasets, but should be used
              with caution. As Kitanidis notes, kriging with a moving window
              can produce unexpected oddities if the variogram model
              is not carefully chosen.
  
          Returns
          -------
          zvalues : ndarray, shape (M, N) or (N, 1)
              Z-values of specified grid or at the specified set of points.
              If style was specified as 'masked', zvalues will
              be a numpy masked array.
          sigmasq : ndarray, shape (M, N) or (N, 1)
              Variance at specified grid points or at the specified
              set of points. If style was specified as 'masked', sigmasq
              will be a numpy masked array.
          """
  
          if self.verbose:
              print("Executing Ordinary Kriging...\n")
  
          if style != 'grid' and style != 'masked' and style != 'points':
              raise ValueError("style argument must be 'grid', "
                               "'points', or 'masked'")
  
          if n_closest_points is not None and n_closest_points <= 1:
              # If this is not checked, nondescriptive errors emerge
              # later in the code.
              raise ValueError("n_closest_points has to be at least two!")
  
          xpts = np.atleast_1d(np.squeeze(np.array(xpoints, copy=True)))
          ypts = np.atleast_1d(np.squeeze(np.array(ypoints, copy=True)))
          n = self.X_ADJUSTED.shape[0]
          nx = xpts.size
          ny = ypts.size
          a = self._get_kriging_matrix(n)
  
          if style in ['grid', 'masked']:
              if style == 'masked':
                  if mask is None:
                      raise IOError("Must specify boolean masking array "
                                    "when style is 'masked'.")
                  if mask.shape[0] != ny or mask.shape[1] != nx:
                      if mask.shape[0] == nx and mask.shape[1] == ny:
                          mask = mask.T
                      else:
                          raise ValueError("Mask dimensions do not match "
                                           "specified grid dimensions.")
                  mask = mask.flatten()
              npt = ny*nx
              grid_x, grid_y = np.meshgrid(xpts, ypts)
              xpts = grid_x.flatten()
              ypts = grid_y.flatten()
  
          elif style == 'points':
              if xpts.size != ypts.size:
                  raise ValueError("xpoints and ypoints must have "
                                   "same dimensions when treated as "
                                   "listing discrete points.")
              npt = nx
          else:
              raise ValueError("style argument must be 'grid', "
                               "'points', or 'masked'")
  
          if self.coordinates_type == 'euclidean':
              xpts, ypts = _adjust_for_anisotropy(np.vstack((xpts, ypts)).T,
                                                  [self.XCENTER, self.YCENTER],
                                                  [self.anisotropy_scaling],
                                                  [self.anisotropy_angle]).T
              # Prepare for cdist:
              xy_data = np.concatenate((self.X_ADJUSTED[:, np.newaxis],
                                        self.Y_ADJUSTED[:, np.newaxis]), axis=1)
              xy_points = np.concatenate((xpts[:, np.newaxis],
                                          ypts[:, np.newaxis]), axis=1)
          elif self.coordinates_type == 'geographic':
              # In spherical coordinates, we do not correct for anisotropy.
              # Also, we don't use scipy.spatial.cdist, so we do not have to
              # format the input data accordingly.
              pass
  
          if style != 'masked':
              mask = np.zeros(npt, dtype='bool')
  
          c_pars = None
          if backend == 'C':
              try:
                  from .lib.cok import _c_exec_loop, _c_exec_loop_moving_window
              except ImportError:
                  print('Warning: failed to load Cython extensions.\n'
                        '   See https://github.com/bsmurphy/PyKrige/issues/8 \n'
                        '   Falling back to a pure python backend...')
                  backend = 'loop'
              except:
                  raise RuntimeError("Unknown error in trying to "
                                     "load Cython extension.")
  
              c_pars = {key: getattr(self, key) for key in ['Z', 'eps', 'variogram_model_parameters', 'variogram_function']}
  
          if n_closest_points is not None:
              if self.coordinates_type == 'geographic':
                  # To make use of the KDTree, we have to convert the
                  # spherical coordinates into three dimensional Euclidean
                  # coordinates, since the standard KDTree cannot handle
                  # the periodicity.
                  # Do the conversion just for the step involving the KDTree:
                  lon_d = self.X_ADJUSTED[:, np.newaxis] * np.pi / 180.0
                  lat_d = self.Y_ADJUSTED[:, np.newaxis] * np.pi / 180.0
                  xy_data = np.concatenate((np.cos(lon_d) * np.cos(lat_d),
                                            np.sin(lon_d) * np.cos(lat_d),
                                            np.sin(lat_d)), axis=1)
                  lon_p = xpts[:, np.newaxis] * np.pi / 180.0
                  lat_p = ypts[:, np.newaxis] * np.pi / 180.0
                  xy_points = np.concatenate((np.cos(lon_p) * np.cos(lat_p),
                                              np.sin(lon_p) * np.cos(lat_p),
                                              np.sin(lat_p)), axis=1)
  
              from scipy.spatial import cKDTree
              tree = cKDTree(xy_data)
              bd, bd_idx = tree.query(xy_points, k=n_closest_points, eps=0.0)
  
              if self.coordinates_type == 'geographic':
                  # Between the nearest neighbours from Euclidean search,
                  # calculate the great circle distance using the standard method:
                  x_points = np.tile(xpts[:,np.newaxis],(1,n_closest_points))
                  y_points = np.tile(ypts[:,np.newaxis],(1,n_closest_points))
                  bd = core.great_circle_distance(x_points, y_points,
                                                  self.X_ADJUSTED[bd_idx],
                                                  self.Y_ADJUSTED[bd_idx])
  
              if backend == 'loop':
                  zvalues, sigmasq = \
                      self._exec_loop_moving_window(a, bd, mask, bd_idx)
              elif backend == 'C':
                  zvalues, sigmasq = \
                      _c_exec_loop_moving_window(a, bd, mask.astype('int8'),
                                                 bd_idx, self.X_ADJUSTED.shape[0],
                                                 c_pars)
              else:
                  raise ValueError('Specified backend {} for a moving window '
                                   'is not supported.'.format(backend))
          else:
              if self.coordinates_type == 'euclidean':
                  bd = cdist(xy_points,  xy_data, 'euclidean')
              elif self.coordinates_type == 'geographic':
                  bd = core.great_circle_distance(xpts[:,np.newaxis],
                                                  ypts[:,np.newaxis],
                                                  self.X_ADJUSTED, self.Y_ADJUSTED)
  
              if backend == 'vectorized':
                  zvalues, sigmasq = self._exec_vector(a, bd, mask)
              elif backend == 'loop':
                  zvalues, sigmasq = self._exec_loop(a, bd, mask)
              elif backend == 'C':
                  zvalues, sigmasq = _c_exec_loop(a, bd, mask.astype('int8'),
                                                  self.X_ADJUSTED.shape[0],
                                                  c_pars)
              else:
                  raise ValueError('Specified backend {} is not supported for '
                                   '2D ordinary kriging.'.format(backend))
  
          if style == 'masked':
              zvalues = np.ma.array(zvalues, mask=mask)
              sigmasq = np.ma.array(sigmasq, mask=mask)
  
          if style in ['masked', 'grid']:
              zvalues = zvalues.reshape((ny, nx))
              sigmasq = sigmasq.reshape((ny, nx))
  
          return zvalues, sigmasq