Source detection step option to fit model PSFs

[bmorris3]

Description

This PR improves astrometry in the Source Detection step by (optionally) making use of PSF fitting methods added in #794.

The centroid accuracy of the PSF fitting methods depends on many factors, like the precision of the PSF model used in the fit (which may be traded off for longer runtimes), the flux of the source, the WFI filter, and to a lesser extent, the detector position, etc. The revised tests for source detection in this PR check the following criteria when using PSF fitting on bright sources in the F087 filter:

• the median centroid residual is smaller than 3 milliarcseconds
• the maximum centroid residual is smaller than 11 milliarcseconds
• the recovered centroids and their errors are all <2-sigma consistent with the injected values

(see here in the diff for these assertions).

These upper-limits in the tests can be revised when we:

1. Precompute PSF models and make them available (e.g. on CRDS). With more accurate (and more expensive to compute) PSF models, we can get more precise centroids, and revise the upper limits in the tests downwards.
2. Have more specific requirements for a given filter and/or a given source flux

This PR supersedes #936.

Resolves RCAL-609

Closes #830
Closes #936

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

Attention: 3 lines in your changes are missing coverage. Please review.

Files	Coverage Δ
romancal/lib/psf.py	76.34% <100.00%> (-8.77%)	⬇️
romancal/source_detection/source_detection_step.py	82.65% <86.95%> (+0.37%)	⬆️

📢 Thoughts on this report? Let us know!

[schlafly]

This looks good. Have you run any romanisim or regtest simulated images? I'm still worried a bit about the DAOStarFinder options, and it would be good to see that we detect most of the real, bright sources.

[schlafly]

Staring a little at the threshold logic, I'm confused. daofind doesn't do it natively, I think, but morally one wants to find all point sources that are more than N sigma above background, for N ~ 5. That looks potentially like what 'calc_threshold' mode is doing, but there it looks like the background is both part of the threshold and subtracted from the data, which seems like doubly subtracting. Meanwhile the mode without calc_threshold doesn't look like it has a notion of RMS and instead only has a notion of background, and the threshold is set exactly at the background?

I'm not sure there are any cases where we don't want to measure and remove a background, so I think I'd recommend just removing the calc_threshold = False branch, and then fixing the calc_threshold = True branch so that doesn't ~doubly subtract the background?

Sorry, I realize that this isn't new code, etc., but we should still fix it.

[bmorris3]

I honestly did not understand what was happening in this block either. I added the else block that you've commented on here to make this PR work, but I agree that the logic needs revision.

[schlafly]

I don't know how hard it is to check, and maybe we defer this to another PR, but we really want to use at least a structured array here eventually, as well as saving all the other columns that may have value.

[bmorris3]

I will ping @PaulHuwe on Slack to ask about this.

[bmorris3]

I confirmed that structured arrays are supported in ASDF. You can test for yourself with:

import numpy as np
import roman_datamodels.datamodels as rdd
from roman_datamodels.maker_utils import mk_level2_image

im = rdd.ImageModel(mk_level2_image())

idx = np.random.randint(0, 100, 10)
x = np.random.uniform(0, 100, 10)
y = np.random.uniform(0, 100, 10)

recarr = np.core.records.fromarrays(
    [idx, x, y], 
    names="idx, x, y", 
    formats=[int, float, float]
)

im.meta['test_recarr'] = recarr

im.to_asdf('test_recarr.asdf')

im_reloaded = rdd.open('test_recarr.asdf')

recarr2 = im_reloaded.meta['test_recarr']

assert np.all(recarr2 == recarr)

I propose that we make a follow-up PR for the move to structured arrays, since that will touch tweakreg a bunch. Thoughts?

[schlafly]

That sounds good, thanks.

[PaulHuwe]

As long as no one cares about accessing this information as text in the metadata, then this is great.

[bmorris3]

By default, this catalog is later deleted by tweakreg, so it shouldn't be useful to anyone unless they stop the pipeline part-way through.

[schlafly]

We definitely don't want a potentially huge catalog converted to text anywhere!

[schlafly]

I think ultimately we want one L2 source catalog for each image with a ~fixed schema. That would argue for merging the PSF & initial finder results together---they're row-by-row matched, right? Could also be part of a separate PR.

[bmorris3]

I agree that we should have all results (source detection and optionally PSF) reported in one table. Since tweakreg currently expects an unstructured array, I won't implement the update in this PR (as discussed in #841 (comment)).

In the follow-up PR, we can have a structured array with columns that tweakreg uses, which will be the PSF fitting results if they're available and the DAO results if not, and also extra columns that tweakreg ignores.

[bmorris3]

@schlafly: All related regression tests are passing (failures are unrelated).

[schlafly]

Great, looks good to me, thanks!

[braingram]

@bmorris3 Do you recall the reasoning for the webbpsf exact pin?

romancal/pyproject.toml

Line 38 in 0ad965c

"webbpsf == 1.2.1",

A new version (1.3.0) is available which appears to be compatible with numpy 2.0 (the version pinned here is not).