In [2]:
%pylab inline
# combine_dustmaps.py: make a composite extinction map (Marshall,Green,Drimmel)
###############################################################################
import sys
import os, os.path
import pickle
import numpy
import h5py
import healpy
from galpy.util import save_pickles
import mwdust
_DEGTORAD= numpy.pi/180.
_ERASESTR= " "
_DRYISHRUN= False
def dist2distmod(dist):
"""dist in kpc to distance modulus"""
return 5.*numpy.log10(dist)+10.
def distmod2dist(distmod):
"""distance modulus to distance in kpc"""
return 10.**(distmod/5.-2.)
_greendir= os.path.join(os.getenv('DUST_DIR'),'green19')
_GREEN15DISTMODS= numpy.linspace(4.,19.,31)
_GREEN15DISTS= distmod2dist(_GREEN15DISTMODS)
In [2]:
greendata = h5py.File(os.path.join(_greendir,'bayestar2019.h5'),'r')
In [3]:
greendata['/pixel_info'].attrs['DM_bin_edges']
Out[3]:
In [12]:
numpy.linspace(4,18.875,120)
Out[12]:
In [191]:
numpy.linspace(4,19,31)
Out[191]:
In [22]:
drimmel= mwdust.Drimmel03()
In [23]:
dists = _GREEN15DISTS
In [188]:
dists
Out[188]:
In [7]:
dists = _GREEN15DISTS
plt.plot(dists,drimmel(270,0,dists))
plt.xlabel('dist in kpc')
plt.ylabel('E(B-V) in mag')
plt.xscale('log')
plt.yscale('linear')
In [8]:
healpy.pix2ang(256,0,nest = True, lonlat = True)
Out[8]:
In [145]:
def get_nside(healpix_number = 1):
"""
returns the nside for specific healpix level
"""
return(np.power(2,healpix_number))
def number_of_healpixels(healpix_number = 1):
"""
returns the number of pixels for a specific level
"""
return(np.power(4,healpix_number)*12)
In [147]:
healpy.pix2ang(256,0,nest = True, lonlat = True)
Out[147]:
In [10]:
hpx_level = 8
nside = get_nside(hpx_level)
print(nside)
nHpx = number_of_healpixels(hpx_level)
ext = np.zeros(shape = (nHpx,len(_GREEN15DISTS)))
for i in range(nHpx):
if i%1000 == 0:
print(i, nHpx)
lb = healpy.pix2ang(nside,i,nest = True, lonlat = True)
ext[i] = drimmel(lb[0],lb[1],_GREEN15DISTS)
In [11]:
from astropy.io import fits
fits.writeto('drimmel_map',ext)
In [15]:
ext.shape
Out[15]:
In [17]:
pixinfo = numpy.recarray(len(ext),dtype = [("nside","uint32"),("healpix_index","uint64")])
pixinfo["nside"] = np.ones(len(ext))*256
pixinfo["healpix_index"] = np.arange(len(ext))
In [18]:
outfile = h5py.File("drimmel.h5","w")
outfile.create_dataset("pixel_info", data = pixinfo)
outfile.create_dataset("best_fit", data = ext)
outfile.close()
In [19]:
drimmel = h5py.File("drimmel.h5",'r')
In [4]:
import dustmaps
In [13]:
from dustmaps.drimmel import DrimmelQuery
from astropy.coordinates import SkyCoord
import astropy.units as u
In [ ]:
drimmel = DrimmelQuery(map_fname='/home/rybizki/programme/dust_green/drimmel/drimmel.h5')
dr = mwdust.Drimmel03()
In [177]:
from healpy import ang2pix
In [179]:
ang2pix(256, 0, 0, nest=True, lonlat=True)
Out[179]:
In [180]:
healpy.pix2ang(256,311296,nest = True, lonlat = True)
Out[180]:
In [182]:
hpx_level = 8
nside = get_nside(hpx_level)
print(nside)
nHpx = number_of_healpixels(hpx_level)
i = 311297
lb = healpy.pix2ang(nside,i,nest = True, lonlat = True)
l = lb[0]
b = lb[1]
d = dists
coords = SkyCoord(l*np.ones_like(d)*u.deg, b*np.ones_like(d)*u.deg, d*u.kpc, frame='galactic')
ext_1 = drimmel(coords)
ext_2 = dr(l,b,d)
for i,item in enumerate(d):
print(i,"dist = %.2f" %item, "green=%.2f" %ext_1[i], "bovy=%.2f" %ext_2[i]) #ext_1[i]-ext_2[i])
In [183]:
from dustmaps import bayestar
In [184]:
bs = bayestar.BayestarQuery()
In [186]:
x = h5py.File("/home/rybizki/programme/dust_green/bayestar/bayestar2019.h5","r")
In [187]:
In [ ]: