docs: def query(self, coords, return_sigma=False):
          """
          Returns r-band extinction, A_r, at the given coordinates. Can also
          return uncertainties.
  
          Args:
              coords (:obj:`astropy.coordinates.SkyCoord`): The coordinates to query.
              return_sigma (Optional[:obj:`bool`]): If ``True``, returns the uncertainty in
                  extinction as well. Defaults to ``False``.
  
          Returns:
              Extinction in the r-band at the specified coordinates, in mags.
              The shape of the output depends on whether :obj:`coords` contains
              distances.
  
              If :obj:`coords` does not specify distance(s), then the shape of the
              output begins with :obj:`coords.shape`. If :obj:`coords` does specify
              distance(s), then the shape of the output begins with
              ``coords.shape + ([number of distance bins],)``.
          """
          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 distances to extract
          if has_dist:
              d = coords.distance.kpc
              dist_idx_ceil = np.searchsorted(self._dists, d)
  
              ret = np.empty((n_coords_ret,), dtype='f8')
              if return_sigma:
                  sigma_ret = np.empty((n_coords_ret,), dtype='f8')
  
              # d < d(nearest distance slice)
              idx_near = (dist_idx_ceil == 0) & ~mask_idx
              print('d < d(nearest): {:d}'.format(np.sum(idx_near)))
              if np.any(idx_near):
                  a = d[idx_near] / self._dists[0]
                  ret[idx_near] = a[:] * self._A[pix_idx[idx_near], 0]
                  if return_sigma:
                      sigma_ret[idx_near] = a[:] * self._sigma_A[pix_idx[idx_near], 0]
  
              # d > d(farthest distance slice)
              idx_far = (dist_idx_ceil == self._n_dists) & ~mask_idx
              print('d > d(farthest): {:d}'.format(np.sum(idx_far)))
              if np.any(idx_far):
                  ret[idx_far] = self._A[pix_idx[idx_far], -1]
                  if return_sigma:
                      sigma_ret[idx_far] = self._sigma_A[pix_idx[idx_far], -1]
  
              # d(nearest distance slice) < d < d(farthest distance slice)
              idx_btw = ~idx_near & ~idx_far & ~mask_idx
              print('d(nearest) < d < d(farthest): {:d}'.format(np.sum(idx_btw)))
              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)
                  ret[idx_btw] = (
                      (1.-a[:]) * self._A[pix_idx[idx_btw], dist_idx_ceil[idx_btw]]
                      +    a[:] * self._A[pix_idx[idx_btw], dist_idx_ceil[idx_btw]-1])
                  if return_sigma:
                      w0 = (1.-a)**2
                      w1 = a**2
                      norm = 1. / (w0 + w1)
                      w0 *= norm
                      w1 *= norm
                      sigma_ret[idx_btw] = np.sqrt(
                          w0 * self._sigma_A[pix_idx[idx_btw], dist_idx_ceil[idx_btw]]**2
                          + w1 * self._sigma_A[pix_idx[idx_btw], dist_idx_ceil[idx_btw]-1]**2
                      )
          else:
              # TODO: Harmonize order of distances & samples with Bayestar.
              ret = self._A[pix_idx, :]
              if return_sigma:
                  sigma_ret = self._sigma_A[pix_idx, :]
  
          if np.any(mask_idx):
              ret[mask_idx] = np.nan
              if return_sigma:
                  sigma_ret[mask_idx] = np.nan
  
          if return_sigma:
              return ret, sigma_ret
  
          return ret