docs: def _bilinear_interp(xyref, zref, xi, yi):
          """
          Perform bilinear interpolation of four 2D arrays located at
          points on a regular grid.
  
          Parameters
          ----------
          xyref : list of 4 (x, y) pairs
              A list of 4 ``(x, y)`` pairs that form a rectangle.
  
          refdata : 3D `~numpy.ndarray`
              A 3D `~numpy.ndarray` of shape ``(4, nx, ny)``.  The first
              axis corresponds to ``xyref``, i.e. ``refdata[0, :, :]`` is
              the 2D array located at ``xyref[0]``.
  
          xi, yi : float
              The ``(xi, yi)`` point at which to perform the
              interpolation.  The ``(xi, yi)`` point must lie within the
              rectangle defined by ``xyref``.
  
          Returns
          -------
          result : 2D `~numpy.ndarray`
              The 2D interpolated array.
          """
  
          if len(xyref) != 4:
              raise ValueError('xyref must contain only 4 (x, y) pairs')
  
          if zref.shape[0] != 4:
              raise ValueError('zref must have a length of 4 on the first '
                               'axis.')
  
          xyref = [tuple(i) for i in xyref]
          idx = sorted(range(len(xyref)), key=xyref.__getitem__)
          xyref = sorted(xyref)   # sort by x, then y
          (x0, y0), (_x0, y1), (x1, _y0), (_x1, _y1) = xyref
  
          if x0 != _x0 or x1 != _x1 or y0 != _y0 or y1 != _y1:
              raise ValueError('The refxy points do not form a rectangle.')
  
          if not np.isscalar(xi):
              xi = xi[0]
          if not np.isscalar(yi):
              yi = yi[0]
  
          if not x0 <= xi <= x1 or not y0 <= yi <= y1:
              raise ValueError('The (x, y) input is not within the rectangle '
                               'defined by xyref.')
  
          data = np.asarray(zref)[idx]
          weights = np.array([(x1 - xi) * (y1 - yi), (x1 - xi) * (yi - y0),
                              (xi - x0) * (y1 - yi), (xi - x0) * (yi - y0)])
          norm = (x1 - x0) * (y1 - y0)
  
          return np.sum(data * weights[:, None, None], axis=0) / norm