Source code for dustmaps.sfd

#!/usr/bin/env python
#
# sfd.py
# Reads the Schlegel, Finkbeiner & Davis (1998; SFD) dust reddening map.
#
# Copyright (C) 2016  Gregory M. Green, Edward F. Schlafly
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#

from __future__ import print_function, division

import os
import numpy as np

import astropy.wcs as wcs
import astropy.io.fits as fits
from scipy.ndimage import map_coordinates

from .std_paths import *
from .map_base import DustMap, WebDustMap, ensure_flat_galactic
from . import fetch_utils
from . import dustexceptions


[docs]class SFDQuery(DustMap): """ Queries the Schlegel, Finbeiner & Davis (1998) dust reddening map. """
[docs] def __init__(self, map_dir=None): """ Args: map_dir (Optional[str]): The directory containing the SFD map. Defaults to `None`, which means that `dustmaps` will look in its default data directory. """ if map_dir is None: map_dir = os.path.join(data_dir(), 'sfd') self._data = {} base_fname = os.path.join(map_dir, 'SFD_dust_4096') for pole in ['ngp', 'sgp']: fname = '{}_{}.fits'.format(base_fname, pole) try: with fits.open(fname) as hdulist: self._data[pole] = [hdulist[0].data, wcs.WCS(hdulist[0].header)] except IOError as error: print(dustexceptions.data_missing_message('sfd', "SFD'98")) raise error
[docs] @ensure_flat_galactic def query(self, coords, order=1): """ Returns E(B-V) at the specified location(s) on the sky. See Table 6 of Schlafly & Finkbeiner (2011) for instructions on how to convert this quantity to extinction in various passbands. Args: coords (`astropy.coordinates.SkyCoord`): The coordinates to query. order (Optional[int]): Interpolation order to use. Defaults to `1`, for linear interpolation. Returns: A float array containing the SFD E(B-V) at every input coordinate. The shape of the output will be the same as the shape of the coordinates stored by `coords`. """ out = np.zeros(len(coords.l.deg), dtype='f4') for pole in ['ngp', 'sgp']: m = (coords.b.deg >= 0) if pole == 'ngp' else (coords.b.deg < 0) if np.any(m): data, w = self._data[pole] x, y = w.wcs_world2pix(coords.l.deg[m], coords.b.deg[m], 0) out[m] = map_coordinates(data, [y, x], order=order, mode='nearest') return out
[docs]class SFDWebQuery(WebDustMap): """ Remote query over the web for the Schlegel, Finkbeiner & Davis (1998) dust map. This query object does not require a local version of the data, but rather an internet connection to contact the web API. The query functions have the same inputs and outputs as their counterparts in ``SFDQuery``. """
[docs] def __init__(self, api_url=None): super(SFDWebQuery, self).__init__( api_url=api_url, map_name='sfd')
[docs]def fetch(): """ Downloads the Schlegel, Finkbeiner & Davis (1998) dust map, placing it in the data directory for `dustmap`. """ doi = '10.7910/DVN/EWCNL5' for pole in ['ngp', 'sgp']: requirements = {'filename': 'SFD_dust_4096_{}.fits'.format(pole)} local_fname = os.path.join( data_dir(), 'sfd', 'SFD_dust_4096_{}.fits'.format(pole)) print('Downloading SFD data file to {}'.format(local_fname)) fetch_utils.dataverse_download_doi( doi, local_fname, file_requirements=requirements)