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],)``.

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