docs: def _calculate_data_point_zscalars(self, x, y, type_='array'):
          """Determines the Z-scalar values at the specified coordinates
          for use when setting up the kriging matrix. Uses bilinear
          interpolation.
          Currently, the Z scalar values are extracted from the input Z grid
          exactly at the specified coordinates. This means that if the Z grid
          resolution is finer than the resolution of the desired kriged grid,
          there is no averaging of the scalar values to return an average
          Z value for that cell in the kriged grid. Rather, the exact Z value
          right at the coordinate is used."""
  
          if type_ == 'scalar':
              nx = 1
              ny = 1
              z_scalars = None
          else:
              if x.ndim == 1:
                  nx = x.shape[0]
                  ny = 1
              else:
                  ny = x.shape[0]
                  nx = x.shape[1]
              z_scalars = np.zeros(x.shape)
  
          for m in range(ny):
              for n in range(nx):
  
                  if type_ == 'scalar':
                      xn = x
                      yn = y
                  else:
                      if x.ndim == 1:
                          xn = x[n]
                          yn = y[n]
                      else:
                          xn = x[m, n]
                          yn = y[m, n]
  
                  if xn > np.amax(self.external_Z_array_x) or \
                     xn < np.amin(self.external_Z_array_x) or \
                     yn > np.amax(self.external_Z_array_y) or \
                     yn < np.amin(self.external_Z_array_y):
                      raise ValueError("External drift array does not cover "
                                       "specified kriging domain.")
  
                  # bilinear interpolation
                  external_x2_index = \
                      np.amin(np.where(self.external_Z_array_x >= xn)[0])
                  external_x1_index = \
                      np.amax(np.where(self.external_Z_array_x <= xn)[0])
                  external_y2_index = \
                      np.amin(np.where(self.external_Z_array_y >= yn)[0])
                  external_y1_index = \
                      np.amax(np.where(self.external_Z_array_y <= yn)[0])
                  if external_y1_index == external_y2_index:
                      if external_x1_index == external_x2_index:
                          z = self.external_Z_array[external_y1_index, external_x1_index]
                      else:
                          z = (self.external_Z_array[external_y1_index, external_x1_index] *
                               (self.external_Z_array_x[external_x2_index] - xn) +
                               self.external_Z_array[external_y2_index, external_x2_index] *
                               (xn - self.external_Z_array_x[external_x1_index])) / \
                              (self.external_Z_array_x[external_x2_index] -
                               self.external_Z_array_x[external_x1_index])
                  elif external_x1_index == external_x2_index:
                      if external_y1_index == external_y2_index:
                          z = self.external_Z_array[external_y1_index, external_x1_index]
                      else:
                          z = (self.external_Z_array[external_y1_index, external_x1_index] *
                               (self.external_Z_array_y[external_y2_index] - yn) +
                               self.external_Z_array[external_y2_index, external_x2_index] *
                               (yn - self.external_Z_array_y[external_y1_index])) / \
                              (self.external_Z_array_y[external_y2_index] -
                               self.external_Z_array_y[external_y1_index])
                  else:
                      z = (self.external_Z_array[external_y1_index, external_x1_index] *
                           (self.external_Z_array_x[external_x2_index] - xn) *
                           (self.external_Z_array_y[external_y2_index] - yn) +
                           self.external_Z_array[external_y1_index, external_x2_index] *
                           (xn - self.external_Z_array_x[external_x1_index]) *
                           (self.external_Z_array_y[external_y2_index] - yn) +
                           self.external_Z_array[external_y2_index, external_x1_index] *
                           (self.external_Z_array_x[external_x2_index] - xn) *
                           (yn - self.external_Z_array_y[external_y1_index]) +
                           self.external_Z_array[external_y2_index, external_x2_index] *
                           (xn - self.external_Z_array_x[external_x1_index]) *
                           (yn - self.external_Z_array_y[external_y1_index])) / \
                          ((self.external_Z_array_x[external_x2_index] -
                            self.external_Z_array_x[external_x1_index]) *
                           (self.external_Z_array_y[external_y2_index] -
                            self.external_Z_array_y[external_y1_index]))
  
                  if type_ == 'scalar':
                      z_scalars = z
                  else:
                      if z_scalars.ndim == 1:
                          z_scalars[n] = z
                      else:
                          z_scalars[m, n] = z
  
          return z_scalars