|
| 1 | +#!/usr/bin/env python |
| 2 | +# |
| 3 | +# csfd.py |
| 4 | +# Reads the "Corrected SFD" dust reddening map of Chiang (2023). |
| 5 | +# |
| 6 | +# Copyright (C) 2023 Gregory M. Green |
| 7 | +# |
| 8 | +# This program is free software; you can redistribute it and/or modify |
| 9 | +# it under the terms of the GNU General Public License as published by |
| 10 | +# the Free Software Foundation; either version 2 of the License, or |
| 11 | +# (at your option) any later version. |
| 12 | +# |
| 13 | +# This program is distributed in the hope that it will be useful, |
| 14 | +# but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 15 | +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 16 | +# GNU General Public License for more details. |
| 17 | +# |
| 18 | +# You should have received a copy of the GNU General Public License along |
| 19 | +# with this program; if not, write to the Free Software Foundation, Inc., |
| 20 | +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
| 21 | +# |
| 22 | + |
| 23 | +from __future__ import print_function, division |
| 24 | + |
| 25 | +import os |
| 26 | +import numpy as np |
| 27 | +import healpy as hp |
| 28 | +import astropy.io.fits as fits |
| 29 | +import astropy.units as units |
| 30 | + |
| 31 | +from .std_paths import * |
| 32 | +from .healpix_map import HEALPixQuery |
| 33 | +from . import fetch_utils |
| 34 | +from . import dustexceptions |
| 35 | + |
| 36 | + |
| 37 | +class CSFDQuery(HEALPixQuery): |
| 38 | + """ |
| 39 | + Queries the Corrected SFD dust map of Chiang (2023). This map is based |
| 40 | + on SFD, but contains a correction to remove contamination from |
| 41 | + large-scale structure (i.e., external galaxies). |
| 42 | + """ |
| 43 | + |
| 44 | + def __init__(self, map_fname=None, mask_fname=None): |
| 45 | + """ |
| 46 | + Args: |
| 47 | + map_fname (Optional[:obj:`str`]): Filename of the CSFD EBV map. |
| 48 | + Defaults to ```None``, meaning that the default location is |
| 49 | + used. |
| 50 | + mask_fname (Optional[:obj:`str`]): Filename of the CSFD mask map. |
| 51 | + Defaults to ```None``, meaning that the default location is |
| 52 | + used. |
| 53 | + """ |
| 54 | + |
| 55 | + if map_fname is None: |
| 56 | + map_fname = os.path.join(data_dir(), 'csfd', 'csfd_ebv.fits') |
| 57 | + if mask_fname is None: |
| 58 | + mask_fname = os.path.join(data_dir(), 'csfd', 'mask.fits') |
| 59 | + |
| 60 | + try: |
| 61 | + with fits.open(map_fname) as hdulist: |
| 62 | + ebv_data = hdulist['xtension'].data[:]['T'].flatten() |
| 63 | + with fits.open(mask_fname) as hdulist: |
| 64 | + mask_data = hdulist['xtension'].data[:]['T'].flatten() |
| 65 | + except IOError as error: |
| 66 | + print(dustexceptions.data_missing_message('csfd', |
| 67 | + 'CSFD (Chiang 2023)')) |
| 68 | + raise error |
| 69 | + |
| 70 | + super(CSFDQuery, self).__init__(ebv_data, False, 'galactic', |
| 71 | + flags=mask_data) |
| 72 | + |
| 73 | + def query(self, coords, **kwargs): |
| 74 | + """ |
| 75 | + Returns CSFD reddening on the same scale as SFD (similar to E(B-V)) at |
| 76 | + the specified location(s) on the sky. Also optionally returns a |
| 77 | + bit mask, where the bits (ordered from least to most significant) have |
| 78 | + the following meanings:: |
| 79 | + |
| 80 | + Bit 0: 'LSS_corr' - This bit is set in the footprint within which |
| 81 | + the LSS is reconstructed, and CSFD = SFD - LSS (otherwise |
| 82 | + CSFD = SFD). |
| 83 | + Bit 1: 'no_IRAS' - Set in the area with no IRAS data (DIRBE data |
| 84 | + filled in SFD); LSS removal in CSFD is done using a 1 deg |
| 85 | + smoothed LSS. |
| 86 | + Bit 2: 'cosmology' - Set in the area where both the LSS and CSFD |
| 87 | + are most reliable for precision cosmology analyses. |
| 88 | +
|
| 89 | + Args: |
| 90 | + coords (:obj:`astropy.coordinates.SkyCoord`): The coordinates to |
| 91 | + query. |
| 92 | + return_flags (Optional[:obj:`bool`]): If ``True``, then a |
| 93 | + bit mask is returned as well, indicating where CSFD |
| 94 | + has been corrected for large-scale structure, where IRAS data |
| 95 | + was used, and where the map is suitable for cosmology. See |
| 96 | + above description of bits. Defaults to ``False``. |
| 97 | +
|
| 98 | + Returns: |
| 99 | + A float array of the reddening, at the given coordinates. The |
| 100 | + shape of the output is the same as the shape of the input |
| 101 | + coordinate array, ``coords``. If ``return_flags`` is ``True``, |
| 102 | + a second array (a bit mask) of the same shape is returned. See |
| 103 | + above description of the meaning of each bit. |
| 104 | + """ |
| 105 | + return super(CSFDQuery, self).query(coords, **kwargs) |
| 106 | + |
| 107 | + |
| 108 | +def fetch(clobber=False): |
| 109 | + """ |
| 110 | + Downloads the Corrected SFD dust map of Chiang (2023). |
| 111 | +
|
| 112 | + Args: |
| 113 | + clobber (Optional[bool]): If ``True``, any existing file will be |
| 114 | + overwritten, even if it appears to match. If ``False`` (the |
| 115 | + default), ``fetch()`` will attempt to determine if the dataset |
| 116 | + already exists. This determination is not 100\% robust against data |
| 117 | + corruption. |
| 118 | + """ |
| 119 | + file_spec = [ |
| 120 | + ('csfd_ebv.fits', '31cd2eec51bcb5f106af84a610ced53c'), |
| 121 | + ('mask.fits', '9142f5a5d184125836a68b6f48d1113f') |
| 122 | + ] |
| 123 | + for fn,md5sum in file_spec: |
| 124 | + fname = os.path.join(data_dir(), 'csfd', fn) |
| 125 | + # Download from Zenodo |
| 126 | + url = 'https://zenodo.org/record/8207175/files/{}'.format(fn) |
| 127 | + fetch_utils.download_and_verify(url, md5sum, fname, clobber=clobber) |
| 128 | + |
0 commit comments