#!/usr/bin/env python
#
# bh.py
# Reads the Burstein & Heiles (1982; BH) dust reddening map.
#
# Copyright (C) 2016 Gregory M. Green
#
# 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 numpy as np
import astropy.coordinates as coordinates
import h5py
import os
from .std_paths import *
from .map_base import DustMap, ensure_flat_galactic
[docs]def ascii2h5(bh_dir=None):
"""
Convert the Burstein & Heiles (1982) dust map from ASCII to HDF5.
"""
if bh_dir is None:
bh_dir = os.path.join(data_dir_default, 'bh')
fname = os.path.join(bh_dir, '{}.ascii')
f = h5py.File('bh.h5', 'w')
for region in ('hinorth', 'hisouth'):
data = np.loadtxt(fname.format(region), dtype='f4')
# Reshape and clip
data.shape = (210, 201) # (R, N)
data = data[:201] # Last 9 records are empty
# Use NaNs where no data
data[data < -9000] = np.nan
dset = f.create_dataset(
region,
data=data,
chunks=True,
compression='gzip',
compression_opts=3
)
dset.attrs['axes'] = ('R', 'N')
dset.attrs['description'] = (
'HI 21cm column densities, in units of 10*NHYD. '
'R = 100 + [(90^o-|b|) sin(l)]/[0.3 degrees]. '
'N = 100 + [(90^o-|b|) cos (l)]/[0.3 degrees].'
)
for region in ('rednorth', 'redsouth'):
data = np.loadtxt(fname.format(region), dtype='f4')
# Reshape and clip
data.shape = (94, 1200) # (R, N)
data = data[:93] # Last record is empty
# Use NaNs where no data
data[data < -9000] = np.nan
dset = f.create_dataset(
region,
data=data,
chunks=True,
compression='gzip',
compression_opts=3
)
dset.attrs['axes'] = ('R', 'N')
dset.attrs['description'] = (
'E(B-V), in units of 0.001 mag. '
'R = (|b| - 10) / (0.6 degrees). '
'N = (l + 0.15) / 0.3 - 1.'
)
f.attrs['description'] = (
'The Burstein & Heiles (1982) dust map.'
)
f.close()
[docs]class BHQuery(DustMap):
"""
Queries the Burstein & Heiles (1982) reddening map.
"""
[docs] def __init__(self, bh_dir=None):
"""
Args:
bh_dir (Optional[str]): The directory containing the Burstein &
Heiles dust map. Defaults to `None`, meaning that the default
directory is used.
"""
if bh_dir is None:
bh_dir = os.path.join(data_dir_default, 'bh')
f = h5py.File(os.path.join(bh_dir, 'bh.h5'), 'r')
self._hinorth = f['hinorth'][:]
self._hisouth = f['hisouth'][:]
self._rednorth = f['rednorth'][:]
self._redsouth = f['redsouth'][:]
f.close()
def _lb2RN_northcap(self, l, b):
R = 100. + (90. - b) * np.sin(np.radians(l)) / 0.3
N = 100. + (90. - b) * np.cos(np.radians(l)) / 0.3
return np.round(R).astype('i4'), np.round(N).astype('i4')
def _lb2RN_southcap(self, l, b):
R = 100. + (90. + b) * np.sin(np.radians(l)) / 0.3
N = 100. + (90. + b) * np.cos(np.radians(l)) / 0.3
return np.round(R).astype('i4'), np.round(N).astype('i4')
def _lb2RN_mid(self, l, b):
R = (np.abs(b) - 10.) / 0.6
N = (np.mod(l, 360.) + 0.15) / 0.3 - 1
return np.round(R).astype('i4'), np.round(N).astype('i4')
def _lb2ebv_northcap(self, l, b):
R, N = self._lb2RN_northcap(l, b)
return -0.0372 + self._hinorth[R,N] * 0.0000357
def _lb2ebv_southcap(self, l, b):
R, N = self._lb2RN_southcap(l, b)
return -0.0372 + self._hisouth[R,N] * 0.0000357
def _lb2ebv_midnorth(self, l, b):
R, N = self._lb2RN_mid(l, b)
return self._rednorth[R,N] * 0.001
def _lb2ebv_midsouth(self, l, b):
R, N = self._lb2RN_mid(l, b)
return self._redsouth[R,N] * 0.001
[docs] @ensure_flat_galactic
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
def main():
#ascii2h5()
bh = BHQuery()
# Calculate E(B-V) on a grid
l = np.arange(-180, 180., 0.1)
b = np.arange(-90., 90.01, 0.1)
l, b = np.meshgrid(l, b)
c = coordinates.SkyCoord(l, b, frame='galactic', unit='deg')
ebv = bh.query(c)
# Apply gamma stretch
gamma = 0.8
img = np.power(np.abs(ebv), gamma) * np.sign(ebv)
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(12,6), dpi=300)
ax = fig.add_subplot(1,1,1, axisbg='blue')
ax.imshow(
img,
origin='lower',
interpolation='none',
cmap='Greys',
aspect='equal',
extent=(-180, 180., -90., 90.),
vmin=-np.power(0.1, gamma),
vmax=np.power(0.5, gamma),
rasterized=True
)
ax.set_xlim(ax.get_xlim()[::-1])
ax.set_xlabel(r'$\ell$', fontsize=18)
ax.set_ylabel(r'$b$', fontsize=18)
ax.set_title(r'$\mathrm{Burstein - Heiles \ \left( 1982 \right)}$', fontsize=22)
fig.savefig(os.path.join(output_dir, 'bh.svg'), dpi=300, bbox_inches='tight')
#plt.show()
return 0
if __name__ == '__main__':
main()