Major improvements and fixes for companion fitting

[kammerje]

This is a big PR comprising of multiple improvements and fixes:

1. Get pixel area in steradian directly from FITS file header instead of computing it from pySIAF pixel scale. The value is slightly different and I'm assuming that the JWST pipeline used what is written to the FITS file header so this will make things more consistent.
2. Enable contrast curve calculation for MIRI 4QPM (simply copying Aarynn's implementation), MIRI LYOT is still not implemented!
3. Bugfixes when injecting and recovering companions into multiple datasets at once, especially if using high-pass filtering.
4. Increased size of automatic blurring to 2 times the theoretical value, based on testing with beta Pic data this is necessary to not obtain Fourier ripples.
5. Added background subtraction step in companion fitting routine. Before, it could happen that the KLIP-subtracted data is negative in the background if there is a bright source such as an extended disk in the images. This lead to very poor companion fitting results and wrong photometry. There is an option to activate a dedicated background removal before fitting the PSFs. Improvement for beta Pic data are from contrast within ~15% to contrast within ~3%.
6. It is now possible to not save the stage 2 files when injecting and recovering companion, to limit the number of disk space required.
7. Centering has been changed to (shape-1)/2 after a discussion with @AarynnCarter and consistently works with the injection and recovery module.

[kammerje]

As an IMPORTANT note, the following code needs to be added to pyklip.rdi.py in order to make RDI work with high-pass filtering (after line 56):

if isinstance(highpass, bool):
if highpass:
data = high_pass_filter_imgs(data, numthreads=None)
else:
# should be a number
if isinstance(highpass, (float, int)):
highpass = float(highpass)
fourier_sigma_size = (data.shape[1]/(highpass)) / (2np.sqrt(2np.log(2)))
data = high_pass_filter_imgs(data, numthreads=None, filtersize=fourier_sigma_size)

[kammerje]

This plot shows the new background removal step and how it improves the measured contrast (injected contrast = 1e-4):

[mperrin]

For point 1 about pixel scale, we ought to be able to do more than assume what the pipeline is doing. Let’s ask the developers.

The subtle part is that there is not a single unique correct pixel area; it varies slightly over the field of view due to optical distortions. So we should also double check which field point is being referred to. These are going to be tiny differences I expect but we should quantify.

[kammerje]

For point 1 about pixel scale, we ought to be able to do more than assume what the pipeline is doing. Let’s ask the developers.

The subtle part is that there is not a single unique correct pixel area; it varies slightly over the field of view due to optical distortions. So we should also double check which field point is being referred to. These are going to be tiny differences I expect but we should quantify.

Yes, please do so! The main interest for me is to reverse exactly what the JWST pipeline did when it converted Jy to Jy/sr.

[mperrin]

For point 1 about pixel scale, we ought to be able to do more than assume what the pipeline is doing. Let’s ask the developers.
The subtle part is that there is not a single unique correct pixel area; it varies slightly over the field of view due to optical distortions. So we should also double check which field point is being referred to. These are going to be tiny differences I expect but we should quantify.

Yes, please do so! The main interest for me is to reverse exactly what the JWST pipeline did when it converted Jy to Jy/sr.

It turns out the pipeline is getting the average pixel area for each mode from the PHOTOM calibration files in CRDS. See docs here for that calibration step: https://jwst-pipeline.readthedocs.io/en/stable/jwst/photom/main.html
Specifically the subsection entitled “Pixel area data”

[maxwellmb]

As an IMPORTANT note, the following code needs to be added to pyklip.rdi.py in order to make RDI work with high-pass filtering (after line 56):

if isinstance(highpass, bool): if highpass: data = high_pass_filter_imgs(data, numthreads=None) else: # should be a number if isinstance(highpass, (float, int)): highpass = float(highpass) fourier_sigma_size = (data.shape[1]/(highpass)) / (2_np.sqrt(2_np.log(2))) data = high_pass_filter_imgs(data, numthreads=None, filtersize=fourier_sigma_size)

Yesterday I put in a pull request with a similar fix within pyklip!

[kammerje]

@AarynnCarter @semaphoreP @mperrin (tagging all relevant people)

After some time of silence, I finally found the time to finish preparing this big PR from dev/jk into develop. In preparation for this PR, I have already merged develop into dev/jk a few days ago, resolved all the conflicts, tested things with the ERS data, and included a couple of additional fixes and improvements. Below, you can find a list of all the updates that this PR includes. Items in bold are new compared to the list posted at the top of this PR (several months ago). I recommend merging this PR as soon as possible before develop diverts too far from dev/jk again and new merging conflicts arise!

As an additional note for @semaphoreP, this PR is fully compatible with the pyKLIP jwst branch which I also updated this week. I reverted back to the previous version which does not need pysiaf as a dependency and also implemented a fallback option for when the SVO FPS is unavailable (see here).

1. Get pixel area in steradian directly from FITS file header instead of computing it from pysiaf pixel scale and propagate it through the spaceKLIP database.
2. Enable contrast curve computation for MIRI 4QPM (simply copying @AarynnCarter's implementation), MIRI LYOT is still not implemented!
3. Bugfixes when injecting and recovering companions into multiple datasets at once, especially if using high-pass filtering.
4. Bugfixes when computing automatic blurring factor. As noted by @mperrin, the previous implementation was wrong. The automatic blurring factor is now set to 1.5 times the theoretically required value based on experience with the beta Pic data.
5. Add background subtraction step in companion fitting routine. Before, it could happen that the KLIP-subtracted data is negative in the background if there is a bright source such as an extended disk in the images. This lead to very poor companion fitting results and wrong photometry. There is now an option to activate a dedicated background removal step before fitting the PSFs.
6. Enabling possibility to not save stage 2 files when injecting and recovering companions to limit the required disk space.
7. Change centering to (shape - 1) / 2 after a discussion with @AarynnCarter. Consistently works with the injection and recovery module.
8. Ensuring compatibility with the pyKLIP JWST class and enabling processing of data without PIXAR_A2 SCI header keyword while not invoking additional dependencies on pyKLIP. This is required to process data for which the JWST pipeline photometry step has been skipped for instance.
9. Update installation instructions on GitHub (not on Read The Docs!).
10. Add fallback options for all SVO FPS calls for the case that the service is temporarily unavailable.
11. Pull PSF masks directly from CRDS instead of computing them manually in make_psfmasks.py.