""" Planck inputs management.
Component separation inputs are downloaded from magique3 under /redtruck/ashdown/repository.
A subset of this repository is needed to run the smica comp sep pipelines.
"""
import os.path as op
import os
import glob
import healpy
import numpy as np
import subprocess
import spherelib.astro as astro
import spherelib
import pyfits
REPO = os.getenv("REDTRUCK") # local mirror of m3 /redtruck repo
if REPO is None:
print "WARNING: empty enviroment variable REDTRUCK"
RELEASE = os.getenv("RELEASE") # data set name in [dx11dr2, ffp8]
if RELEASE is None:
#os.environ["RELEASE"] = "dx11dr2"
RELEASE = "dx11dr2"
print "WARNING: using default data set dx11dr2"
else:
print "INFO: using %s data set"%RELEASE
[docs]def ensure_arbo(filename):
""" Ensure that filename's subdirectory exists in the local REDTRUCK repository
"""
if not op.exists(op.dirname(filename)):
os.makedirs(op.dirname(filename))
return
[docs]def check_file_is_present(filename, pla=False, product="MAP.MAP_ID"):
""" Check that filename is already present, and download it if not.
"""
if op.exists(filename):
return
else:
basename = filename.split(REPO)[-1]
dirname = op.dirname(filename)
ensure_arbo(filename)
if pla:
source = "http://pla.esac.esa.int/pla/aio/product-action?%s=%s"%(product, op.basename(filename))
subprocess.check_call(["wget", "-c" , "-O", op.join(dirname,source.split("=")[-1]), source])
else:
if 'extra' in op.dirname(filename):
source = op.join("/data/lejeune/planck", basename[1:])
subprocess.check_call(["rsync", "-avz" , "--copy-links", "lejeune@apccluster1.in2p3.fr:"+source, filename])
elif 'delouis' in op.dirname(filename):
source = op.join("/redtruck", basename[1:])
subprocess.check_call(["rsync", "-avz" , "--copy-links", "lejeune@magique3.iap.fr:"+source, filename])
else:
source = op.join("/redtruck/ashdown/repository", basename[1:])
subprocess.check_call(["rsync", "-avz" , "--copy-links", "lejeune@magique3.iap.fr:"+source, filename])
return
[docs]def parse_unit(inmap):
"""
"""
unit_fr = spherelib.read_key(inmap, "TUNIT1").strip()
if unit_fr=='Kcmb':
unit_fr = 'K_CMB'
return unit_fr
[docs]def is_lfi(name):
return get_freq(name)<100
[docs]def is_hfi(name):
return get_freq(name)>=100
[docs]def is_polarized(name, release=None):
if release is None:
release = RELEASE
if release=="R2.00":
return get_freq(name)<100
else:
return get_freq(name)<545
[docs]def get_fwhm_beam(name):
""" Return approximate FWHM.
"""
name = "%03d"%get_freq(name)
dfwhm = dict({"030":32.29, "044":27.66, "070":13.21, "100":9.69, "143":7.30, "217":5.02, "353":4.94, "545":4.83, "857":4.64})
return dfwhm[name]
[docs]def RIMO2Bl(fn, name):
""" Read beam window function from RIMO and write them to cl file.
"""
hdulist = pyfits.open(fn)
index = hdulist.index_of("BEAMWF_%sx%s"%(name,name))
if is_lfi(name):
bl = hdulist[index].data["BL"]
bl = np.hstack((bl, np.zeros((4001-len(bl)))))
else:
bl = hdulist[index].data["NOMINAL"]
fn = op.join(op.dirname(fn), "BEAMS_%sx%s.fits"%(name,name))
healpy.write_cl(fn, np.squeeze(bl))
return fn
[docs]def get_bl(name, release=None):
""" Return transfer function from file.
"""
if release is None:
release = RELEASE
if '_' in name:
name = name.replace("_", '-')
if 'ds' in name:
name = name.replace("-", "")
if release=="R2.00":
dirmap = op.join(REPO, "pla", release, "beams")
if is_lfi(name):
RIMOfile = op.join(dirmap, "LFI_RIMO_R2.50.fits")
else:
RIMOfile = op.join(dirmap, "HFI_RIMO_Beams-100pc_R2.00.fits")
check_file_is_present(RIMOfile, pla=True, product="DOCUMENT_MAP.DOCUMENT_ID")
blfile = RIMO2Bl(RIMOfile, name)
if release in ['dx11dr2', 'rd11c']:
dirbl = op.join(REPO, "exchanges", "dx11", "beams")
hfidir = op.join(dirbl, "hfi_extended")
lfidir = op.join(dirbl, "lfi_extended")
if is_lfi(name):
blfile = op.join(lfidir, "bls_GB_%s_T_extended.fits"%name)
elif name=="857":
blfile = op.join(hfidir, "qb_bl_ds_%s-2_x_%s-2_extended.fits"%(name,name))
else:
blfile = op.join(hfidir, "qb_bl_ds_%s_x_%s_extended.fits"%(name,name))
elif 'ffp8' in release:
dirbl = op.join(REPO, "simulations", "ffp8", "beams")
hfidir = op.join(dirbl, "hfi_extended")
lfidir = op.join(dirbl, "lfi_extended")
if is_lfi(name):
blfile = op.join(lfidir, "bls_GB_%s_T_extended.fits"%name)
else:
blfile = op.join(hfidir, "qb_bl_ds_%s_x_%s_extended.fits"%(name,name))
check_file_is_present(blfile)
return healpy.read_cl(blfile)
[docs]def get_freq(name):
""" Return central frequency as int
"""
return int(name.split("_")[0].split('-')[0].split('ds')[0])
[docs]def get_unit_factor(name, fr='K_CMB', to='K_CMB', release=None):
""" Return extra conversion factor MJy/sr -> KCMB
"""
if release is None:
release = RELEASE
freq = get_freq(name)
cv_base = astro.convfact(freq=freq*1e9, fr=fr, to=to)
cv1 = astro.convfact(freq=freq*1e9, fr=fr, to="MJy/sr")
cv2 = astro.convfact(freq=freq*1e9, fr="K_CMB", to=to)
if release in ["R2.00", 'dx11dr2', 'rd11c']:
dcv = dict({"545":58.062295, "545_1":57.083062, "545_2":58.882429, "545_4":58.059482,
"857":2.2703657, "857_1":2.1890447, "857_2":2.3455933, "857_3":2.2132231,"857_4":2.4021213})
if freq in [545,857]:
cv = cv1/dcv[name]*cv2
else:
cv = cv_base
elif 'ffp8' in release:
if freq==545:
cv = cv1/57.895428*cv2
elif freq==857:
cv = cv1/2.2291155*cv2
else:
cv = cv_base
return cv
[docs]def get_mapfile(name, set, release=None):
""" Return input map for a given name in [030,044,...,857] and dataset in [full, hr1, hr2, ...]
Default release is read from RELEASE environment variable
"""
product = "MAP.MAP_ID"
pla = False
if release is None:
release = RELEASE
freq = get_freq(name)
if release=="R2.00":
pla = True
dirmap = op.join(REPO, "pla", release, "maps")
hfidset = dict({"full":"full", "hr1":"full-ringhalf-1", "hr2":"full-ringhalf-2", "yr1":"year-1", "yr2":"year-2", "hm1":"halfmission-1", "hm2":"halfmission-2"})
lfidset = dict({"full":"full", "hr1":"full-ringhalf-1", "hr2":"full-ringhalf-2", "yr1":"year-1-3", "yr2":"year-2-4", "hm1":"halfmission-1", "hm2":"halfmission-2"})
if is_lfi(name):
release = "R2.01"
inmap = op.join(dirmap, "LFI_SkyMap_%03d_1024_%s_%s.fits"%(freq,release,lfidset[set]))
else:
inmap = op.join(dirmap, "HFI_SkyMap_%03d_2048_%s_%s.fits"%(freq,release,hfidset[set]))
if release=='rd11c':
if is_lfi(name) or freq>353:
release = 'dx11dr2'
else:
dirmap = op.join(REPO, "delouis", "m3data", "PROD_RD11c")
hfidset = dict({"full":"ful", "hr1":"hr1", "hr2":"hr2", "yr1":"yr1", "yr2":"yr2"})
inmap = op.join(dirmap, "%03dGHz_%s.fits"%(freq,hfidset[set]))
check_file_is_present(inmap)
return inmap
if release=='dx11dr2':
dirmap = op.join(REPO, "exchanges", "dx11", "maps")
hfidir = op.join(dirmap, "hfi_zodi_removed_bpleak_corrected")
lfidir = op.join(dirmap, "lfi_bpleak_corrected")
lfidset = dict({"full":"full", "hr1":"full_hr1", "hr2":"full_hr2", "yr1":"yr1+yr3", "yr2":"yr2+yr4", "hm1":"yr1+yr3", "hm2":"yr2+yr4"})
hfidset = dict({"full":"full", "hr1":"full_hr1", "hr2":"full_hr2", "yr1":"yr1", "yr2":"yr2", "hm1":"hm1", "hm2":"hm2"})
if set[0:2]=='ds':
inmap = op.join(hfidir, "hfi_dx11d_%03d_%s_full_nozodi.fits"%(freq,set[-1]))
elif is_lfi(name):
inmap = op.join(lfidir, "lfi_dx11d_%s_%s_bplcorrected.fits"%(name,lfidset[set]))
elif name=="857":
inmap = op.join(hfidir, "hfi_dx11d_%03d_2_%s_nozodi_bplcorrected_dust_ground.fits"%(freq,hfidset[set]))
else:
inmap = op.join(hfidir, "hfi_dx11d_%s_%s_nozodi_bplcorrected_dust_ground.fits"%(name,hfidset[set]))
elif 'ffp8' in release:
dirmap = op.join(REPO, "simulations", "ffp8", "maps_preprocessed")
if len(release)>4:
dirmap = op.join(dirmap, release)
lfidset = dict({"full":"full_map", "hr1":"full_map_hr1", "hr2":"full_map_hr2", "yr1":"y13_map", "yr2":"y24_map", "hm1":"y13_map", "hm2":"y24_map"})
hfidset = dict({"full":"full_map", "hr1":"full_map_hr1", "hr2":"full_map_hr2", "yr1":"y1_map", "yr2":"y2_map", "hm1":"hm1_map", "hm2":"hm2_map"})
if is_lfi(name) and set in lfidset.keys():
inmap = op.join(dirmap, "ffp8_total_nobpm_%s_%s.fits"%(name,lfidset[set]))
elif set in hfidset.keys():
inmap = op.join(dirmap, "ffp8_total_nobpm_%s_%s.fits"%(name,hfidset[set]))
else:
dirmap = op.join(dirmap, "components")
inmap = op.join(dirmap, "ffp8_%s_%s_full_map.fits"%(set,name))
check_file_is_present(inmap, pla=pla, product=product)
return inmap
[docs]def get_psmaskfile(name, release=None):
""" Return point source mask provided by the MHW team.
"""
if release is None:
release = RELEASE
freq = get_freq(name)
name = "%03d"%freq
fwhm = np.int(get_fwhm_beam(name))
product = "MAP.MAP_ID"
pla = False
if release=="R2.00":
pla = True
dirmap = op.join(REPO, "pla", release, "masks")
if is_lfi(name):
offmap = op.join(dirmap, "LFI_Mask_PointSrc_2048_R2.00.fits")
else:
offmap = op.join(dirmap, "HFI_Mask_PointSrc_2048_R2.00.fits")
if release in ['dx11dr2', 'rd11c']:
dirmask = op.join(REPO, "exchanges", "dx11", "masks")
hfidir = op.join(dirmask, "sources_hfi")
lfidir = op.join(dirmask, "sources_lfi")
if is_lfi(name):
offmap = op.join(lfidir, 'mask_ps_%dGHz_beam%damin_nside2048.00_DX11D_full_v1.0_eff_beams.fits.gz'%(freq,fwhm))
else:
offmap = op.join(hfidir, 'dx11d_%s_mask_sources_int_snr_5_radius_3sigma.fits.gz'%(name))
elif 'ffp8' in release:
dirmask = op.join(REPO, "simulations", "ffp8", "masks")
hfidir = op.join(dirmask, "sources_hfi")
lfidir = op.join(dirmask, "sources_lfi")
if is_lfi(name):
offmap = op.join(lfidir, 'mask_ps_%dGHz_beam%damin_nside2048.00_FFP8_nonblind_variable_holesize.fits.gz'%(freq,fwhm))
else:
offmap = op.join(hfidir, 'ffp8_%s_mask_sources_int_snr_5_radius_3sigma.fits.gz'%(name))
check_file_is_present(offmap, product=product, pla=pla)
return offmap
[docs]def get_catalogfile(name, release=None):
""" Return point source catalog provided by the MHW team.
"""
if release is None:
release = RELEASE
freq = get_freq(name)
name = "%03d"%freq
product = "SOURCE_LIST.NAME"
pla = False
if release=="R2.00":
pla = True
dirmap = op.join(REPO, "pla", release, "catalogs")
if is_lfi(name):
fn = op.join(dirmap, "COM_PCCS_%s_R1.30.fits"%name)
else:
fn = op.join(dirmap, "COM_PCCS_%s_R1.20.fits"%name)
if release in ['dx11dr2', 'rd11c']:
dircat = op.join(REPO, "exchanges", "dx11", "catalogues")
hfidir = op.join(dircat, "hfi")
lfidir = op.join(dircat, "lfi")
if is_lfi(name):
fn = op.join(lfidir, "cat_%dNB_4col.dat"%freq)
else:
fn = op.join(hfidir, "COM_PCCS_%s_DX11d.fits"%name)
elif 'ffp8' in release:
dircat = op.join(REPO, "simulations", "ffp8", "catalogues")
hfidir = op.join(dircat, "hfi")
lfidir = op.join(dircat, "lfi")
if is_lfi(name):
fn = op.join(lfidir, "cat_%dNB_4col.dat"%freq)
else:
fn = op.join(hfidir, "ffp8_%s_catalogue.fits"%name)
check_file_is_present(fn, product=product, pla=pla)
return fn
[docs]def load_hfi_catalog(fn, name, snr_t, nside=2048):
""" Load HFI MHW catalog from COM PCCP usual fits format
(glon, glat in columns 2 and 3 ; flux and error in columns 6 and 7)
"""
beam = get_fwhm_beam(name)
dtr = np.pi/180.0
ratio = 2.0*np.sqrt(2.0*np.log(2.0))
sigma = (beam/ratio)/60.0*dtr
beam_area = 2*np.pi*(sigma**2)
freq = get_freq(name)
if 'ffp8' in fn:
flux = spherelib.read_col(fn, col=3)
err = spherelib.read_col(fn, col=4)
glon = spherelib.read_col(fn, col=1)
glat = spherelib.read_col(fn, col=2)
else:
flux = spherelib.read_col(fn, col=6)
err = spherelib.read_col(fn, col=7)
glon = spherelib.read_col(fn, col=2)
glat = spherelib.read_col(fn, col=3)
ns = len(flux)
cv = astro.convfact (freq=freq*1e9, fr='mJy/sr', to='K_CMB')
theta= astro.glat2theta(glat)
phi = astro.glon2phi(glon)
cat = np.zeros((ns, 5))
cat[:,0] = healpy.ang2pix(nside, theta, phi)
cat[:,2] = glat
cat[:,1] = glon
cat[:,3] = flux*cv/beam_area
cat[:,4] = err*cv/beam_area
snr = cat[:,3]/cat[:,4]
return cat,theta,phi,snr
[docs]def load_lfi_catalog(fn, name, snr_t, nside=1024):
""" Load LFI MHW catalog from 4 columns txt format
(glon, gcolat, flux and snr)
"""
freq = get_freq(name)
tmpcat = np.loadtxt(fn)
ns = np.shape(tmpcat)[0]
glat = 90-tmpcat[:,0]
glon = tmpcat[:,1]
theta= astro.glat2theta(glat)
phi = astro.glon2phi(glon)
cat = np.zeros((ns, 5))
idxsort = np.argsort(glat)
cat[:,0] = healpy.ang2pix(nside, theta[idxsort], phi[idxsort])
cat[:,2] = glat[idxsort]
cat[:,1] = glon[idxsort]
cat[:,3] = tmpcat[[idxsort],2]
snr = np.squeeze(tmpcat[[idxsort],3])
cat[:,4] = cat[:,3]/snr
return cat,theta,phi,snr
[docs]def load_catalog(name, snr_t):
""" Load MHW catalog and applied SNR cut
The output catalog includes [pixnum, glon, glat, flux, flux_err]
"""
fn = get_catalogfile(name)
if get_freq(name) > 70 or 'COM' in fn:
cat,theta,phi,snr = load_hfi_catalog(fn, name, snr_t)
else:
cat,theta,phi,snr = load_lfi_catalog(fn, name, snr_t)
# apply snr cut
cat = cat[snr>snr_t,:]
ns = np.shape(cat)[0]
phi = phi[snr>snr_t]
theta = theta[snr>snr_t]
snr = snr[snr>snr_t]
return cat,theta,phi,snr
[docs]def get_ps2galmask(radius=0, release=None):
""" Return minimum galactic mask.
TODO : to be replaced by a function which computes it !
"""
if release is None:
release = RELEASE
if release=='rd11c' or "R2.00":
release = 'dx11dr2'
if radius>0:
strapo = "_apo%d"%radius
else:
strapo = ""
fn = op.join(REPO, "extra", "masks", "ps2galmask_%s%s.fits.gz"%(release,strapo))
check_file_is_present(fn)
return fn
[docs]def get_galmask(fsky=60, radius=0):
""" Return a galactic mask file.
"""
fn = op.join(REPO, "exchanges", "dx9", "masks", "galactic", "combined_mask_0.%d_sky_fraction.fits.gz"%fsky)
check_file_is_present(fn)
if radius>0:
strapo = "apo%d"%radius
fnapo = op.join(REPO, "extra", "masks", "combined_mask_0.%d_sky_fraction_%s.fits.gz"%(fsky,strapo))
if not op.exists(fnapo):
mask = healpy.read_map(fn)
spherelib.apodize_mask(mask, radius)
healpy.write_map(fnapo[:-3], mask)
subprocess.check_call(["gzip", fnapo[:-3]])
return fnapo
return fn
[docs]def get_binfile():
""" Return xpol bin file
"""
fn = op.join(REPO, "extra", "planck_xpol_binned_tt_v47_2048_0240_16148_RingCal.txt")
check_file_is_present(fn)
return fn
[docs]def get_lambdaCDMfile(year='2015'):
""" Return likelihood best fit file
"""
if year=='2013':
fn = op.join(REPO, "extra", "base_planck_lowl_lowLike.bestfit_cl")
else:
fn = op.join(REPO, "extra", "base_plikHM_TT_lowTEB.minimum.theory_cl")
check_file_is_present(fn)
return fn
[docs]def get_PlanckClfile(year='2015'):
""" Return planck empirical cmb spectra file
"""
fn = op.join(REPO, "extra", "COM_PowerSpect_CMB_R2.00.fits")
check_file_is_present(fn)
return fn
[docs]def download_all():
""" Download data a priori
"""
names = ["030", "044", "070", "100","143","217","353", "545", "857"]
for name in names:
get_bl(name)
#get_psmaskfile(name)
#get_catalogfile(name)
#for dset in ["full", "hr1", "hr2", "yr1", "yr2", "hm1", "hm2"]:
# get_mapfile(name, dset)
#get_ps2galmask()
#get_ps2galmask(radius=60)
#get_binfile()
#get_lambdaCDMfile()
[docs]def get_calibfile(release=None):
"""
"""
if release is None:
release = RELEASE
if release=='dx11dr2':
return "/workdir/lejeune/planck/outputs/dx11d/calib-mean.txt"
elif release=='rd11c':
return "/workdir/lejeune/planck/outputs/rd11c_calib.txt"
else:
return None
[docs]def get_commonmask(field='int', apo=30, release=None, reso="005a_2048"):
"""
"""
strapo = "_apo_%03da"%apo if apo is not None else ""
if release is None:
release = RELEASE
if 'ffp8' in release:
fn = op.join(REPO, "comparison", "ffp8", "ffp8_common_%s_mask_%s%s.fits"%(field, reso,strapo))
elif release=='dx11dr2':
if field=='int':
fn = op.join(REPO, "comparison", "dx11_v2", "dx11_v2_common_%s_mask_%s%s.fits"%(field, reso,strapo))
else:
fn = op.join(REPO, "comparison", "dx11_v2", "dx11_v2_common_%s_mask_extended2_%s%s.fits"%(field, reso,strapo))
check_file_is_present(fn)
return fn
[docs]def get_misspixmask(field='int', apo=30, release=None, reso="005a_2048"):
strapo = "_apo_%03da"%apo if apo is not None else ""
if release is None:
release = RELEASE
if release=='dx11dr2':
fn = op.join(REPO, "comparison", "dx11_v2", "dx11_v2_common_%s_mask_misspix_hm_%s%s.fits"%(field, reso,strapo))
check_file_is_present(fn)
return fn
[docs]def get_commonmask_mll(field='int', apo=30, release=None):
"""
"""
fn = get_commonmask(field=field, apo=apo, release=release)
fn_mll = fn.replace(".fits.gz", "_mll")
if op.exists(fn_mll+".npy"):
mll = np.load(fn_mll+".npy")
else:
mask = healpy.read_map(fn)
mll = spherelib.mask2mll(mask,3000)
np.save(fn_mll, mll)
return mll
[docs]def get_outputcmb(team='smica', field='int', release=None, reso="005a_2048", dset="full"):
"""
"""
if release is None:
release = RELEASE
if field=='pol' and release=='dx11dr2':
field = 'pol_case1'
ddset = dict({"full":"cmb", "hrhd":"cmb_hrhd", "hrhs":"cmb_hrhs",
"hmhd":"cmb_hmhd", "hmhs":"cmb_hmhs","yrhd":"cmb_yrhd", "yrhs":"cmb_yrhs", })
if 'ffp8' in release:
fn = op.join(REPO, "comparison", "ffp8", "ffp8_%s_%s_%s_%s.fits"%(team, field, ddset[dset],reso))
elif release=='dx11dr2':
fn = op.join(REPO, "comparison", "dx11_v2", "dx11_v2_%s_%s_%s_%s.fits"%(team, field, ddset[dset],reso))
check_file_is_present(fn)
return fn