Returns E(B-V) at the specified location(s) on the sky.
        
        Args:
            coords (`astropy.coordinates.SkyCoord`): The coordinates to query.

        Returns:
            A float array of reddening, in units of E(B-V), at the given
            coordinates. The shape of the output is the same as the shape of the
            coordinates stored by `coords`.

def query(self, coords):
        """
        Returns E(B-V) at the specified location(s) on the sky.
        
        Args:
            coords (`astropy.coordinates.SkyCoord`): The coordinates to query.

        Returns:
            A float array of reddening, in units of E(B-V), at the given
            coordinates. The shape of the output is the same as the shape of the
            coordinates stored by `coords`.
        """
        # gal = coords.transform_to('galactic')
        gal = coords
        l = gal.l.deg
        b = gal.b.deg

        # Detect scalar input
        scalar_input = not hasattr(l, '__len__')
        if scalar_input:
            l = np.array([l])
            b = np.array([b])

        # Fill return array with NaNs
        ebv = np.empty(l.shape, dtype='f8')
        ebv[:] = np.nan

        # Fill northern cap
        idx = (b >= 65.) & (b <= 90.)
        ebv[idx] = self._lb2ebv_northcap(l[idx], b[idx])

        # Fill southern cap
        idx = (b <= -65.) & (b >= -90.)
        ebv[idx] = self._lb2ebv_southcap(l[idx], b[idx])

        # Fill northern midplane
        idx = (b < 65.) & (b >= 10.)
        ebv[idx] = self._lb2ebv_midnorth(l[idx], b[idx])

        # Fill southern midplane
        idx = (b > -65.) & (b <= -10.)
        ebv[idx] = self._lb2ebv_midsouth(l[idx], b[idx])

        if scalar_input:
            ebv = ebv[0]

        return ebv