PSF fitting methods

[bmorris3]

Description

This PR introduces methods for fitting model PSFs computed with WebbPSF to ImageModels in romancal.

Resolves RCAL-573

Closes #781

New dependencies

The model PSF is computed by webbpsf, and the model PSF is fit to Roman ImageModels using the new photutils.psf module refactor in astropy/photutils#1558 (and futher PRs that followed it). I've pinned photutils to the dev version until the next release (expected within about a week), and pinned to the latest version of webbpsf.

Like CRDS, WebbPSF depends on several reference files which must be configured on a given machine before tests can be run. I've updated the CI setup so that GitHub Actions has access to those files. One limitation is that these WebbPSF reference files are versioned along with the software, so there are no dev files easily available to test against. For this reason, I've left WebbPSF out of the devdeps tests.

The last new dependency is scikit-learn. For crowded groups of sources, photutils has a SourceGrouper class, which uses a clustering algorithm implemented in scikit-learn to determine which sources are close enough to fit simultaneously with a joint model.

Implementation

There are two public methods introduced in this PR:

• create_gridded_psf_model uses WebbPSF to construct an ePSF model that is compatible with the photutils psf module. The function primarily serves to provide the user with reasonable defaults based on previous testing.

• fit_psf_to_image_model uses fits the PSF model to a Roman ImageModel. At this time, defaults have been optimized to give PSF fits with a good tradeoff between runtime and astrometric accuracy. Further tests should be done in the future to test photometric accuracy.

Tests are included to ensure that the inferred centroids of bright, high S/N, synthetic point sources are within 3-sigma of their injected values.

Checklist

• added entry in CHANGES.rst under the corresponding subsection
• updated relevant tests
• updated relevant documentation
• updated relevant milestone(s)
• added relevant label(s)

[codecov]

Codecov Report

Patch coverage is 85.10% of modified lines.

Files Changed	Coverage
romancal/lib/psf.py	85.10%

📢 Thoughts on this report? Let us know!.

[bmorris3]

@zacharyburnett @WilliamJamieson – it appears the failure on py39-oldestdeps-cov comes from the photutils minpin for numpy >= 1.22, while the others usually require >= 1.20. Any tips on what to do there?

[WilliamJamieson]

@zacharyburnett @WilliamJamieson – it appears the failure on py39-oldestdeps-cov comes from the photutils minpin for numpy >= 1.22, while the others usually require >= 1.20. Any tips on what to do there?

I don't know why romancal has not updated its numpy lower pin. NumPy <= 1.21 are no longer supported according to NEP 29. @zacharyburnett should probably look into this.

[zacharyburnett]

pinned in #800

[bmorris3]

I'll note here that this file on Box that is shared between several folks at INS. If changes are necessary in the future, contact Marcio Melendez and/or Brad Sappington.

[schlafly]

No changes needed here, but for romanisim it's annoying that we have to update this each time a new webbpsf release happens. I haven't understood why this file can't have a name like webbpsf-latest-release-test.gz or something, and when a new version comes out the new test data file gets updated. That would mean that each new version doesn't require new CI updates.

[bmorris3]

The next release of photutils may be out within a week or so, and I will update this pin when it's available.

[schlafly]

This looks great to me.

Some things to think about:

• If we didn't use IterativePSFPhotometry, could we avoid the scikit-learn dependency? It is pretty big and people might push back. On the other hand, it would be nice to be able to use the IterativePSFPhotometry eventually. I don't actually see it in the code, but I had thought that source grouping was a crowded field feature and we could get by with single star fits here.
• Re romancal vs. stcal vs. photutils: do you think fit_psf_to_image_model(...) would be sensible as part of photutils?
• dq_to_boolean_mask might fit better in romancal.lib.dqflags, but I don't feel strongly.
• There's a world where we precompute the gridded PSFs for each filter / SCA combination and save them somewhere (e.g., in CRDS). How long is that computation?
• Does your wiki page use this code? Can you link to it as part of this PR?
• Can you provide a little bit of run time information?
• I do think we need to understand the flipping of the stamps, though not as part of this PR. But let's at least make an issue.
• This probably more naturally belongs under source_detection/psf.py rather than lib/psf.py?

[schlafly]

I wonder if there is a more natural place for these constants, but I don't have one in mind. I think the filters are fairly fixed now, but I don't actually know that either. If not, we need to update these for consistency's sake as the filters change.

Some day someone will ask to to real integrals over these filters in create_gridded_psf_model (i.e., not monochromatic), but I don't think we should go there for this PR.

[nden]

Is this info a good candidate for a reference file in CRDS? Then it can be changed by end users or INS without changing the code.

[schlafly]

It's not crazy. I think if we wanted to go this route, we should instead get the whole filter curves into CRDS. i.e. what Goddard has here:
https://roman.gsfc.nasa.gov/science/RRI/Roman_effarea_20201130.xlsx
If we had that kind of thing we could also put other summary information about the bandpasses in those files; e.g., central wavelengths. The existing 'photom' file is some part of the way toward that end; one option would be to add more attributes there. Obviously, let's not make that part of this PR, though!

[bmorris3]

Great, some options for moving forward: shall we make a ticket to add these quantities to CRDS, reference that ticket in a TODO comment in the file in this PR, and make an RCAL ticket to do a follow-up PR when those quantities are available?

[nden]

sounds good to me

[schlafly]

I think we need to at least file an issue to follow up on this; it worries me.

[schlafly]

I don't see where this noise is added to the data?

[bmorris3]

Good point. So far I'm not yet adding the shot noise to the data array, but I am increasing the uncertainties as the counts change.

[schlafly]

Can we also assert something about the size of x_err / y_err? Maybe x_err < 1 / (amp / (sqrt(noise^2 + amp)) * fudge factor? I just want to be able to catch failure modes like "all of the errors are very large, so we have small delta_x / delta_y but did a poor fit? I'm trying to use the rule-of-thump x_err ~ FWHM / (S/N), which doesn't apply for this undersampled data but isn't terrible, and using 1 pix for the FWHM, and doing a bad approximation to the S/N.

[bmorris3]

RE @schlafly:

If we didn't use IterativePSFPhotometry, could we avoid the scikit-learn dependency? It is pretty big and people might push back. On the other hand, it would be nice to be able to use the IterativePSFPhotometry eventually. I don't actually see it in the code, but I had thought that source grouping was a crowded field feature and we could get by with single star fits here.

Yes! The additional scikit-learn is avoidable if we don't use IterativePSFPhotometry as it's implemented in photutils right now. But it's also possible(/likely?) that other methods like photutils.psf.SourceGrouper, or clustering algorithms within that method, will be implemented that don't depend on the scikit-learn (cc @larrybradley). There's a growing number of options available in scipy.cluster that we can experiment with in photutils to dodge this big dependency without loosing the iterative fitting functionality.

Re romancal vs. stcal vs. photutils: do you think fit_psf_to_image_model(...) would be sensible as part of photutils?

The two methods in this PR are pretty simple – they allow a user fit PSF models to Roman ImageModels. Almost all of the work in this PR was choosing and vetting sensible defaults for the input parameters for PSF fitting specifically for Roman ImageModels. So personally, I think the Roman-specificity of this PR is the primary feature of this PR, making the other repos ill-fitting homes for this PR.

dq_to_boolean_mask might fit better in romancal.lib.dqflags, but I don't feel strongly.

I agree, that's fair. The only reason I've put it in lib.psf here is because the Roman dq flag map defines "GOOD": 0, which needs to be removed before extending the astropy bitmask registry, in order to work correctly for photutils. I suppose that's true for anyone who expects a boolean mask though.

There's a world where we precompute the gridded PSFs for each filter / SCA combination and save them somewhere (e.g., in CRDS). How long is that computation?

This depends on a few parameters, primarily: (1) the number of PSFs to synthesize across each SCA, (2) the width of the FOV of the rendered PSF. The default gridded PSF model from this PR has a FOV of 50 pixels, rendered in 16 locations distributed across the SCA, and takes 10 seconds to run. For all SCAs, runtime will be a few minutes.

Does your wiki page use this code? Can you link to it as part of this PR?

I wrote up a performance evaluation for the code in this repo in the confluence post. The methods in this PR are essentially identical with two exceptions. This PR provides sensible defaults given the results of the analysis, and this PR has comments/docstrings/tests.

Can you provide a little bit of run time information?

With the default settings in this PR, I've run both methods for an exposure the F087 filter on SCA 7. Constructing a gridded PSF model takes 10 seconds. I ran the PSF model fit on a synthetic ImageModel, without supplying initial centroid or flux guesses. Of the 500 sources injected with romanisim, the default finder (DAOStarFind) detected 250 sources and fit their PSFs in 2.7 seconds.

I do think we need to understand the flipping of the stamps, though not as part of this PR. But let's at least make an issue.

Agreed.

This probably more naturally belongs under source_detection/psf.py rather than lib/psf.py?

I'm happy either way. The forward-looking intention for putting it in romancal.lib is to separate it from Source Detection in anticipation of a future photometry step that also uses PSF photometry.

[schlafly]

Did you mean to close spacetelescope/romanisim#63 instead of this one?

[bmorris3]

Did you mean to close spacetelescope/romanisim#63 instead of this one?

Yikes you're right! Thanks @schlafly!

[bmorris3]

Found and worked on a required upstream fix in astropy/astropy#15215. Update: putting off this upstream fix, and adapting the methods in romancal to cope.

[schlafly]

Thanks. Two more minor comments. I don't think we have any regression tests here, so I can't see how this would have an effect there, but we usually require runs of the regression suite with the changes before merging.

[schlafly]

It looks to me like this is likely out now?

[bmorris3]

The last remaining and necessary addition to photutils is from astropy/photutils#1617. I'll ask @larrybradley if the next release is likely to be soon.

[schlafly]

No changes needed here, but for romanisim it's annoying that we have to update this each time a new webbpsf release happens. I haven't understood why this file can't have a name like webbpsf-latest-release-test.gz or something, and when a new version comes out the new test data file gets updated. That would mean that each new version doesn't require new CI updates.