How To Deconvolve a Map¶
LUKE CURRENTLY UPDATING THIS DOC PAGE.
The deconvolution wrapper has lots of input options with more functionality likely to be added in the Future. Not every possible combination of inputs has been tested and so some errors may get spat out. If you experience such an error please contact Luke Jew or Richard Grumitt who can help you to commit a fix.
For a description of the routine see the overleaf document.
deconvolveMap function¶
The cbasspy.mapManipulations.deconvolveMap.deconvolveMap() function is a wrapper that calls the necessary functions in the correct order. It has lots of input options that are discussed in various categories below. This document concludes with an example of how to call the function.
Necessary parameters¶
| Parameter | Description |
|---|---|
| inMap | Input healpix map to be smoothed/deconvolved. Should be a HEALPix map with shape either (Npix,) or (3,Npix). |
| theta | Angles at which the beam profile has been sampled, note this is the beam profile that the map is being deconvolved by. |
| beam | The beam profile to deconvolve the map by. |
| fwhm | The fwhm (in radians) of the 2D Gaussian that the map will be smoothed to. |
| NSIDE_out | The Nside to downgrade the output map to |
Miscellanious parameters¶
| Parameter | Description |
|---|---|
| niter = 3 | Number of iterations to take when estimating alms from maps. You will not likely need to change this. |
| returnLots = False | Whether to return lots of outputs. You likely will want to set this to True. If False, only returns outMap, outVarMap, outMap_dg, outVarMap_dg, outVarMapScaleFactor. If True, returns outMap, outVarMap, outMap_dg, outVarMap_dg, outVarMapScaleFactor, ell, inBeamBl_I, inBeamBl_P, outBeamBl_I, outBeamBl_P, ratio_I, ratio_P, omega_beam, omega_out |
| inCov = None | The covariance map to use. If this is set to None then no noise information is used. If given then it is used during the degrading of Nside. A smoothed covariance map is also returned, along with a scaling factor that must be applied to it. If not None then shape should be either the same as inMap or (6,Npix). |
| inMapCoord = ‘c’. | If rotating the map then the funciton needs the coordinate system of the input map. |
| numOfNoiseSims = 3. | How many noise simulations to use when estimating the variance map correction factor. |
Window function parameters¶
It takes ~1-2 hours to estimate the window function for high Nside maps. For speed there are options to save/load in pre-calculated window functions. The window fucntions are saved with names [windowFileName]I.txt and [windowFileName]P.txt and the columns are ell, inBeamBl, outBeamBl, ratio.
| Parameter | Description |
|---|---|
| saveWindowFns = False | Whether to save the window functions as a txt file. |
| loadWindowFns = False | Whether to load in a window function. |
| windowFileName =’windowFunctions’ | The suffix of the window function file names. |
Declination mask parameters¶
The deconvolution/smoothing/downgrading can have strange effects at the edge of the maps. Due to the C-BASS scan strategy the edge of the map is a line of constant declination. Therefore mask pixels in the output map between some minimum and maximum declination. This is particularly important when estimating the variance map scale factors. The default is to cut out the southern hemisphere.
| Parameter | Description |
|---|---|
| useDecMask = False | Whether to apply masking. |
| outMapCoord = ‘g’ | The output coordinate system of the maps. They will be rotated to this frame if not already in it. |
| cutMin = -90. | The minimum declination (in degrees) to cut out. |
| cutMax = 0. | The maximum declination (in degrees) to cut out. |
Tapering parameters¶
There is option to apply a cosine taper to the ratio of the window functions.
| Parameter | Description |
|---|---|
| doTapering = False | Whether to apply tapering to the ratio of the window functions |
| ell_min = 200 | The multipole moment to begin the tapering. |
| ell_max = 2000 | The mulitpole moment to end the tapering. |
Output description¶
If returnLots is False
| Output | Description |
|---|---|
| outMap | The deconvolved and smoothed map. |
| outVarMap | The deconvolved and smoothed variance map(s) (or None is inCov=None). |
| outMap_dg | The deconvolved, smoothed and degraded map. |
| outVarMap_dg | The deconvolved, smoothed and degraded variance map(s) (or None is inCov=None). |
| outVarMapScaleFactor | Scaling factor(s) to apply to outVarMap_dg. |
If returnLots is True then the following is also returned
| Output | Description |
|---|---|
| ell | Multipole moments. |
| inBeamBl_I | In beam total intensity legendre transform. |
| inBeamBl_P | In beam polarization legendre transform. |
| outBeamBl_I | Gaussian beam total intensity legendre transform (approximate). |
| outBeamBl_P | Gaussian beam polarization legendre transform (approximate). |
| ratio_I | Ratio of the total intensity transfer functions. |
| ratio_P | Ratio of the polarization transfer functions. |
| omega_beam | Solid angle of the input beam. |
| omega_out | Solid angle of ouput beam. |
Example function call¶
First you need to import some python modules:
import numpy as np
import healpy as h
from scipy import constants
from cbasspy.mapManipulations.deconvolveMap import deconvolveMap
import matplotlib.pyplot as plt
Here is an example loading in the relevant data:
inmap_name = '/users/luke/CBASS/cbass_analysis/IDL/NIGHT_v20_allFiles_masked5pc_C_1024_ol500_c_map.fits'
inCov_name = '/users/luke/CBASS/cbass_analysis/IDL/NIGHT_v20_allFiles_masked5pc_C_1024_ol500_c_cov.fits'
inMap = h.read_map(inmap_name,(0,1,2))
inCov = h.read_map(inCov_name,(0,1,2,3,4,5))
# Beam
deconvBeam_name = '/users/luke/CBASS/cbass_analysis/IDL/cbass-beam.txt'
beamFileData = np.loadtxt(deconvBeam_name)
theta = beamFileData[:,0]
beam = beamFileData[:,1]
# output options
fwhm = 1.0*np.pi/180.
NSIDE_out = 256
The function would then be called:
outMap, outVarMap, outMap_dg, outVarMap_dg, outVarMapScaleFactor,\
ell, inBeamBl_I, inBeamBl_P, outBeamBl_I, outBeamBl_P, ratio_I, ratio_P,\
omega_beam, omega_out = deconvolveMap(inMap,theta,beam,fwhm,NSIDE_out,\
niter=3,inCov=inCov,returnLots=True,inMapCoord='c',numOfNoiseSims=3,\
saveWindowFns=False,loadWindowFns=False,windowFileName='windowFunctions',\
useDecMask=True,outMapCoord='g',cutMin=-90,cutMax=0,\
doTapering=True,ell_min=200,ell_max=2000)
The output maps could be saved for future use, after applying the correction factor:
# Save the total intensity map
h.write_map('NIGHT_v20_allFiles_masked5pc_C_1024_ol500_c_map_g_1deg_0256.fits',
outMap_dg,coord='g',column_units='K_antenna',overwrite=True)
# Make a copy of the downgraded/etc... variance map
correctedVarOut = np.copy(outVarMap_dg)
# Apply correction factor to copy
correctedVarOut[0] = correctedVarOut[0]*outVarMapScaleFactor[0]
correctedVarOut[1] = correctedVarOut[1]*outVarMapScaleFactor[1]
correctedVarOut[2] = correctedVarOut[2]*outVarMapScaleFactor[2]
# Set bad pixels back to healpix unseen value
correctedVarOut[0][outVarMap_dg[0]==h.UNSEEN] = h.UNSEEN
correctedVarOut[1][outVarMap_dg[1]==h.UNSEEN] = h.UNSEEN
correctedVarOut[2][outVarMap_dg[2]==h.UNSEEN] = h.UNSEEN
# Save the covariance map
h.write_map('NIGHT_v20_allFiles_masked5pc_C_1024_ol500_c_cov_g_1deg_0256.fits',
correctedVarOut,coord='g',column_units='K_antenna^2',overwrite=True)
You could plot the window functions:
# Plot window functions
plt.figure()
plt.loglog(ell,ratio_I**2,label='I Ratio',c='DarkOrange',ls='-')
plt.loglog(ell,inBeamBl_I**2,label='I C-BASS',c='purple',ls='--')
plt.loglog(ell,outBeamBl_I**2,label='I Gaussian',c='DarkGreen',ls=':')
plt.xlabel(r'$\ell$',fontsize=15)
plt.ylabel(r'$w=b_\ell^2$',fontsize=15)
plt.tight_layout()
plt.grid()
plt.ylim([1e-27,5])
plt.legend()
plt.savefig('{0}_I.png'.format(windowFileName))
plt.show()
plt.figure()
plt.loglog(ell,ratio_P**2,label='P Ratio',c='DarkOrange',ls='-')
plt.loglog(ell,inBeamBl_P**2,label='P C-BASS',c='purple',ls='--')
plt.loglog(ell,outBeamBl_P**2,label='P Gaussian',c='DarkGreen',ls=':')
plt.xlabel(r'$\ell$',fontsize=15)
plt.ylabel(r'$w=b_\ell^2$',fontsize=15)
plt.tight_layout()
plt.grid()
plt.ylim([1e-27,5])
plt.legend()
plt.savefig('{0}_P.png'.format(windowFileName))
Assuming that your input maps were in correctly calibrated antenna temperatures then you can rescale the outputs to Jy in beam:
nu = 4.76e9
inMap_S = 2.*constants.k*omega_beam/(constants.c/nu)**2*1e26*np.array(inMap)
outMap_S = 2.*constants.k*omega_out/(constants.c/nu)**2*1e26*np.array(outMap)
