Foreground modeling/removal and Power Spectra generation
pip install ps-eor==0.19
You will need the following packages:
The main goal of this package is to provide a simple and clean interface for Foreground removal and PS estimation tasks that can be used by other program/scripts. As such it does not provide an end to end script, but provide examples of such scripts (see for example fits_image_to_vvis_cube.py, generate_ps.py and run_gpr.py)
Below we briefly document the main tasks that can be perform with this package:
data_cube_i = datacube.CartDataCube.load_from_fits_image(files, umin, umax, theta_fov)
This will do the following steps:
It is possible (recommended) to save/load to an h5 format with the save()/load() method.
One can also regrid the data and make an image with the regrid() and image() method.
The main FG removal code is GPR, but PCA, GMCA (the python version) and Poly fitting are also implemented.
To run GPR, one do the following:
data_cube_i = datacube.CartDataCube.load_from_fits_image(files_i, umin, umax, theta_fov)
data_cube_v = datacube.CartDataCube.load_from_fits_image(files_v, umin, umax, theta_fov)
eor_bin_list = pspec.EorBinList(data_cube_i.freqs)
# Create an EoR bin 122-134 MHz with a 120-136MHz range for the FG fitting
eor_bin_list.add_freq(1, 122, 134, 120, 136)
eor = eor_bin_list.get(1)
gpr_config = fitutil.GprConfig.load(gpr_config_filename)
gpr_fit = fgfit.GprForegroundFit(gpr_config)
gpr_res = gpr_fit.run(eor.get_slice_fg(data_cube_i), eor.get_slice_fg(data_cube_v))
This return a GprForegroundResult object which have the following attributes:
And the following method:
On can then save/load those CartDataCube as needed for later processing.
The CartDataCube cubes of the GPR model and residual contains the error covariance from the GPR model that need to be taken into account when generating the power spectra (see fitutil.get_ps_err_from_cov_err for an exemple of how this can be done.).
Look at gpr_config.parset and gpr_config_v.parset for examples of GPR configuration.
The PS code into account automatically the error covariance of the GPR model.
It is possible to generate spatial only PS, Cylindrically averaged PS (2D) or spherically averaged PS (3D).
# Create a PS configuration
el = 2 * np.pi * (np.arange(data_cube.ru.min(), data_cube.ru.max(), du))
ps_conf = pspec.PowerSpectraConfig(el)
pb = datacube.LofarHBAPrimaryBeam()
# Create a PS generation object
ps_gen = pspec.PowerSpectraCart(eor, ps_conf, pb)
# Create a Spatial PS, plot it and save it to a file
ps = ps_gen.get_ps(data_cube)
ps.plot(title='Spatial power spectra')
plt.savefig('ps.pdf')
ps.save_to_txt('ps.txt')
# Create a Cylindrically averaged PS
ps2d = ps_gen.get_ps2d(data_cube)
ps2d.plot(title='Cylindrically averaged power spectra')
plt.savefig('ps2d.pdf')
ps2d.save_to_txt('ps2d.txt')
# Create a Spherically averaged PS
kbins = np.logspace(np.log10(ps_gen.kmin), np.log10(0.5), 10)
ps3d = ps_gen.get_ps3d(kbins, data_cube)
ps3d.plot(title='Spherically averaged power spectra')
plt.savefig('ps3d.pdf')
ps3d.save_to_txt('ps3d.txt')