Add Source Catalog Step

[bmorris3]

(work in progress)

Resolves RCAL-764
Closes #1083

This PR ports the Source Catalog Step from jwst. We're currently discussing what level of support we'll implement in this PR for AB/Vega mags.

Related PRs:

• Adding schema for Source Catalog rad#374
• Add Source Catalog Model roman_datamodels#317

Checklist

• added entry in CHANGES.rst under the corresponding subsection
• updated relevant tests
• updated relevant documentation
• updated relevant milestone(s)
• added relevant label(s)
• ran regression tests, post a link to the Jenkins job below. How to run regression tests on a PR

[codecov]

Codecov Report

Attention: Patch coverage is 96.61836% with 21 lines in your changes are missing coverage. Please review.

Project coverage is 76.78%. Comparing base (2df177e) to head (1898ca5).
Report is 16 commits behind head on main.

Files	Patch %	Lines
romancal/source_catalog/source_catalog.py	96.28%	17 Missing ⚠️
romancal/source_catalog/source_catalog_step.py	94.11%	3 Missing ⚠️
romancal/source_catalog/detection.py	96.66%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1102      +/-   ##
==========================================
- Coverage   77.15%   76.78%   -0.38%     
==========================================
  Files         107      113       +6     
  Lines        7101     7933     +832     
==========================================
+ Hits         5479     6091     +612     
- Misses       1622     1842     +220     

Flag	Coverage Δ
nightly	?

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[schlafly]

The code looks reasonable to me here. I would propose that we record things in native units (i.e., I guess MJy for the L3 images, DN/s / e/s for the L2 images), and just provide the photometry metadata block to allow conversions.

I was not able to completely understand what the aperture correction block is intended to do. We should enlist someone on Webb to provide information there. I would guess that we can compute the values from the PSF and should just use those for the moment.

Do you have a sense for what the units coming out of the aperture and PSF fluxes are? I am getting confused conceptually about whether we need to be multiplying by a factor like the area of a pixel or if that's already being done. i.e., for an L2 image, the images have been flattened so that one gets uniform response for a flat surface brightness. For an L2 file, then, the total flux is just sum(flux * pixel_area / nominal_pixel_area). For an L3 file, the fluxes are sum(flux * pixel_area_in_sr). Are we doing those conversions?

[schlafly]

Re units, yes, this is annoying.

I want the answer to be "let's report the fluxes in the units on the images." DN/s is fine for the L2 catalogs. We should include the conversion to Jy in the metadata; that will be the meta.photometry.mjysr conversion factor, times the pixel area in sr, divided by 1e6. We should also include the zero point, so that AB magnitudes are -2.5log10(DN/s) + ZPT; that is -2.5log10(mjysr * pixel area in sr / 1e6) + 2.5*log10(3631) (3631 Jy defines the AB zero point). Most of that code is already in your PR; please check to see if I've flopped a sign somewhere.

For the L3 images, reporting the answer in MJy/sr is too horrible. I'd propose that if the image is in MJy/sr you convert to Jy are report the fluxes in those units; that's the existing code in the PR. And then in the metadata include the conversion to Jy (1) and the AB zero point (2.5 * log10(3631)).

Does that cover it?

[larrybradley]

This PR is awaiting an update the source catalog schema to include ref_files. @nden

This PR includes a commit that bypasses that check. DO NOT MERGE. I'll remove that last commit when the schemas are updated.

[nden]

@larrybradley We discussed today what needs to go in that section and I am going to update the PR tomorrow.

[larrybradley]

After talking with @nden, I've updated this PR to assign CRDS version only if reference files are used (and an associated test is marked xfail).

This PR is ready for review and can be merged (actually, this PR should go in after #1201).

[larrybradley]

Updates:

• Added support for L2 ImageModel input
• Added a flags columns
  • L2: returns 1 if source centroid pixel has DQ of DO_NOT_USE (dq=1)
  • L3: returns 1 if source centroid pixel has weight=0
  • otherwise, returns 0
  • Added DMS387 regtest
• Computes pixel scale from WCS if model.meta.photometry.pixelarea_steradians < 0
• Flux columns have units of μJy
• Added SourceCatalogStep to the High-level pipeline -- removed

@schlafly

[schlafly]

Looks good. I left a few more comments in line.

One more bigger question: I can derive decent aperture corrections. That would look something like a dictionary with the different filters, each one containing an aperture correction corresponding to a different number of pixels.

I went ahead and computed those numbers. For something like F062, a 30% encircled energy aperture would be significantly smaller than 1 pixel, which doesn't sound like a good idea to me. I'm not sure if it's worth doing anything with these values since obviously apertures of 1, 2, and 3 pixels are just placeholders and we shouldn't treat them as too real. But it's straightforward to get the numbers we need.

apcorr script:

bandpasses = ['F062', 'F087', 'F106', 'F129', 'F158', 'F184', 'F213', 'F146']
apcorr = dict()
aps = [CircularAperture((22, 22), i) for i in range(1, 21)]
wfi = webbpsf.WFI()
wfi.detector = 'SCA07'
wfi.detector_position = (2000, 2000)
for b in bandpasses:
    wfi.filter = b
    res = wfi.calc_psf(oversample=8)
    apfluxes = aperture.aperture_photometry(res[1].data, aps)
    apcorr[b] = dict({k+1: 1/apfluxes[f'aperture_sum_{k}'][0] for k in range(len(aps))})

which gives the following results

{'F062': {1: 1.6700499921933796,
  2: 1.2406223550413196,
  3: 1.171500438628268,
  4: 1.1321238506514975,
  5: 1.1066169334397933,
  6: 1.0883394570637264,
  7: 1.0743343843212518,
  8: 1.0634687494542465,
  9: 1.0551564574694712,
  10: 1.048747672133765,
  11: 1.0437767353435445,
  12: 1.039853878649696,
  13: 1.036751531375191,
  14: 1.034221885295224,
  15: 1.0321223150468193,
  16: 1.0302963477379874,
  17: 1.0286271949235084,
  18: 1.0270375912043688,
  19: 1.025473215962012,
  20: 1.0239242249852876},
 'F087': {1: 1.7548200976882269,
  2: 1.2667229981824595,
  3: 1.1963852768292245,
  4: 1.1551348920967295,
  5: 1.1287836950817338,
  6: 1.1081258286825466,
  7: 1.0926845875008024,
  8: 1.0801798473590478,
  9: 1.069732235873413,
  10: 1.0610535206343705,
  11: 1.0538453947418234,
  12: 1.0480115125787612,
  13: 1.0432119490471066,
  14: 1.0393805001873908,
  15: 1.0362337578398295,
  16: 1.033705178714667,
  17: 1.0316134614036887,
  18: 1.0299071369672412,
  19: 1.0284973359015745,
  20: 1.0273222469201642},
 'F106': {1: 1.8909117430203468,
  2: 1.3072985206344572,
  3: 1.216454734735456,
  4: 1.1748730832612557,
  5: 1.1447836188112714,
  6: 1.124533137673792,
  7: 1.1073241782151584,
  8: 1.094017891875251,
  9: 1.0830610638570812,
  10: 1.0735439599204273,
  11: 1.0653329916840983,
  12: 1.0583372352653677,
  13: 1.0523158561842867,
  14: 1.0472571372093753,
  15: 1.0429353442392868,
  16: 1.039328195105576,
  17: 1.0363559105274074,
  18: 1.0338205548515904,
  19: 1.0317387247450083,
  20: 1.0299969334686752},
 'F129': {1: 2.142876975388637,
  2: 1.4046201136368834,
  3: 1.237525895143953,
  4: 1.1962644228659125,
  5: 1.1651149571628274,
  6: 1.1419804985425677,
  7: 1.125601014350914,
  8: 1.110728671199529,
  9: 1.098618783076443,
  10: 1.0886780741048823,
  11: 1.0800047149070524,
  12: 1.0721828630119183,
  13: 1.065274711117077,
  14: 1.059236921077065,
  15: 1.0539552793931504,
  16: 1.0492331573188327,
  17: 1.0452044353140784,
  18: 1.0416321128068269,
  19: 1.0385241941493994,
  20: 1.035926244759529},
 'F158': {1: 2.4666216623476296,
  2: 1.523356356062162,
  3: 1.264653117360973,
  4: 1.2197884392502485,
  5: 1.1885541668945767,
  6: 1.163136552349395,
  7: 1.143475912233356,
  8: 1.1299539532258107,
  9: 1.1171599697145103,
  10: 1.1056676475523504,
  11: 1.0961919618892066,
  12: 1.088090151346644,
  13: 1.080910052799954,
  14: 1.0741432657206869,
  15: 1.0681455617597457,
  16: 1.062634974536718,
  17: 1.0578072442036803,
  18: 1.0534387206017906,
  19: 1.049443831159602,
  20: 1.0459112058352351},
 'F184': {1: 3.067122486201166,
  2: 1.7675213544222386,
  3: 1.3682676945587233,
  4: 1.262244308521063,
  5: 1.2262167797490469,
  6: 1.204927786865037,
  7: 1.179300954103238,
  8: 1.1551190213900593,
  9: 1.1390836170954974,
  10: 1.126922949845646,
  11: 1.1165722839432453,
  12: 1.1055274356456966,
  13: 1.0952123575918813,
  14: 1.0872476419828878,
  15: 1.0810661667902073,
  16: 1.0749468699061622,
  17: 1.0684176178889564,
  18: 1.062803804168304,
  19: 1.0583392473604245,
  20: 1.0545922469331852},
 'F213': {1: 3.426840472053474,
  2: 1.8731314306436042,
  3: 1.4882021387933424,
  4: 1.2846913937870525,
  5: 1.2466026113506612,
  6: 1.2200548671695983,
  7: 1.2026352311485449,
  8: 1.179817727011062,
  9: 1.157580688061812,
  10: 1.1429582994845708,
  11: 1.1309823353702886,
  12: 1.1221053124348357,
  13: 1.1129320634254833,
  14: 1.1032499688043347,
  15: 1.0944127839055557,
  16: 1.0873949002140282,
  17: 1.0819818612048009,
  18: 1.0768932573673755,
  19: 1.071238620642699,
  20: 1.0657072976695716},
 'F146': {1: 2.1837919773088257,
  2: 1.42172678986078,
  3: 1.24842595877169,
  4: 1.2000070805482999,
  5: 1.1698412393628825,
  6: 1.1466918475334875,
  7: 1.1286469788639037,
  8: 1.1141888841891825,
  9: 1.10246416309668,
  10: 1.0921833209092404,
  11: 1.0831989222002902,
  12: 1.075421802239501,
  13: 1.068620828182733,
  14: 1.0627115055554754,
  15: 1.057457822960673,
  16: 1.0528242617316381,
  17: 1.0487527941306776,
  18: 1.045172162510871,
  19: 1.0420520185280877,
  20: 1.0393071064618822}}

[schlafly]

There's this schema for the L2 catalogs:
https://github.com/spacetelescope/rad/blob/main/src/rad/resources/schemas/source_catalog-1.0.0.yaml

I think we should be able to fill out the metadata by something like:

for key in source_catalog_model.meta.keys():
    source_catalog_model.meta[key] = input_model.meta[key]

which may work transparently for both L2 and L3 files? I think right now the metadata is being filled with filler values?

[larrybradley]

I've added this. optical_element is not in the L2 SourceCatalog schema, so those values are not copied.

[larrybradley]

Ready for final review

[schlafly]

Looks good to me. @ddavis-stsci , want to review?

[ddavis-stsci]

This looks fine to me for the basic catalog product.
At some point we'll need a discussion of how to trigger this for the production catalog products. It appears to me the trigger for the single and multiband products will be quite different.
Will the input parameters for all the catalogs be the same or will they be custom based on the product?

[schlafly]

I agree that we should discuss when to trigger this. These single band products are computationally cheap (much less than the rest of the ELP / HLP) and my proposal is that we make them whenever the ELP / HLP runs on an image.
As you point out, the multi-band products are another story. Those need to wait on all of the individual bands. We're only doing them as part of DRs, so the triggering, as you say, will look a lot different.