docs: def query(self, coords, mode='random_sample'):
          """
          Returns A0 at the given coordinates. There are several different query
          modes, which handle the probabilistic nature of the map differently.
  
          Args:
              coords (:obj:`astropy.coordinates.SkyCoord`): The coordinates to query.
              mode (Optional[:obj:`str`]): Five different query modes are available:
                  ``'random_sample'``, ``'random_sample_per_pix'``, ``'samples'``,
                  ``'median'`` and ``'mean'``. The ``mode`` determines how the output
                  will reflect the probabilistic nature of the IPHAS dust map.
  
          Returns:
              Monochromatic extinction, A0, at the specified coordinates, in mags.
              The shape of the output depends on the ``mode``, and on whether
              ``coords`` contains distances.
  
              If ``coords`` does not specify distance(s), then the shape of the
              output begins with `coords.shape`. If `coords` does specify
              distance(s), then the shape of the output begins with
              ``coords.shape + ([number of distance bins],)``.
  
              If ``mode`` is ``'random_sample'``, then at each coordinate/distance, a
              random sample of reddening is given.
  
              If ``mode`` is ``'random_sample_per_pix'``, then the sample chosen for
              each angular pixel of the map will be consistent. For example, if
              two query coordinates lie in the same map pixel, then the same
              random sample will be chosen from the map for both query
              coordinates.
  
              If ``mode`` is ``'median'``, then at each coordinate/distance, the
              median reddening is returned.
  
              If ``mode`` is ``'mean'``, then at each coordinate/distance, the mean
              reddening is returned.
  
              Finally, if ``mode`` is ``'samples'``, then all at each
              coordinate/distance, all samples are returned.
          """
  
          # Check that the query mode is supported
          valid_modes = [
              'random_sample',
              'random_sample_per_pix',
              'samples',
              'median',
              'mean']
  
          if mode not in valid_modes:
              raise ValueError(
                  '"{}" is not a valid `mode`. Valid modes are:\n'
                  '  {}'.format(mode, valid_modes))
  
          n_coords_ret = coords.shape[0]
  
          # Determine if distance has been requested
          has_dist = hasattr(coords.distance, 'kpc')
          d = coords.distance.kpc if has_dist else None
  
          # Convert coordinates to pixel indices
          pix_idx = self._coords2idx(coords)
  
          # Determine which coordinates are out of bounds
          mask_idx = (pix_idx == self._n_pix)
          if np.any(mask_idx):
              pix_idx[mask_idx] = 0
  
          # Which samples to extract
          if mode == 'random_sample':
              samp_idx = np.random.randint(0, self._n_samples, pix_idx.size)
              n_samp_ret = 1
          elif mode == 'random_sample_per_pix':
              samp_idx = np.random.randint(0, self._n_samples, self._n_pix)[pix_idx]
              n_sample_ret = 1
          else:
              samp_idx = slice(None)
              n_samp_ret = self._n_samples
  
          # Which distances to extract
          if has_dist:
              d = coords.distance.pc
              dist_idx_ceil = np.searchsorted(self._dists, d)
  
              if isinstance(samp_idx, slice):
                  ret = np.empty((n_coords_ret, n_samp_ret), dtype='f4')
              else:
                  ret = np.empty((n_coords_ret,), dtype='f4')
  
              # d < d(nearest distance slice)
              idx_near = (dist_idx_ceil == 0)
              if np.any(idx_near):
                  a = d[idx_near] / self._dists[0]
                  if isinstance(samp_idx, slice):
                      ret[idx_near] = a[:,None] * self._data['A0'][pix_idx[idx_near], 0, samp_idx]
                  else:
                      ret[idx_near] = a[:] * self._data['A0'][pix_idx[idx_near], 0, samp_idx[idx_near]]
  
              # d > d(farthest distance slice)
              idx_far = (dist_idx_ceil == self._n_dists)
              if np.any(idx_far):
                  if isinstance(samp_idx, slice):
                      ret[idx_far] = self._data['A0'][pix_idx[idx_far], -1, samp_idx]
                  else:
                      ret[idx_far] = self._data['A0'][pix_idx[idx_far], -1, samp_idx[idx_far]]
  
              # d(nearest distance slice) < d < d(farthest distance slice)
              idx_btw = ~idx_near & ~idx_far
  
              if np.any(idx_btw):
                  d_ceil = self._dists[dist_idx_ceil[idx_btw]]
                  d_floor = self._dists[dist_idx_ceil[idx_btw]-1]
                  a = (d_ceil - d[idx_btw]) / (d_ceil - d_floor)
                  if isinstance(samp_idx, slice):
                      ret[idx_btw] = (
                          (1.-a[:,None]) * self._data['A0'][pix_idx[idx_btw], dist_idx_ceil[idx_btw], samp_idx]
                          +    a[:,None] * self._data['A0'][pix_idx[idx_btw], dist_idx_ceil[idx_btw]-1, samp_idx])
                  else:
                      ret[idx_btw] = (
                          (1.-a[:]) * self._data['A0'][pix_idx[idx_btw], dist_idx_ceil[idx_btw], samp_idx[idx_btw]]
                          +    a[:] * self._data['A0'][pix_idx[idx_btw], dist_idx_ceil[idx_btw]-1, samp_idx[idx_btw]])
          else:
              # TODO: Harmonize order of distances & samples with Bayestar.
              ret = self._data['A0'][pix_idx, :, samp_idx]
  
          # Reduce the samples in the requested manner
          samp_axis = 1 if has_dist else 2
  
          if mode == 'median':
              ret = np.median(ret, axis=samp_axis)
          elif mode == 'mean':
              ret = np.mean(ret, axis=samp_axis)
  
          if np.any(mask_idx):
              ret[mask_idx] = np.nan
  
          return ret