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.

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