HSM: measure PSF/star shapes directly in sky coords (use_sky_coords) - #812
Merged
Conversation
Attach the local (Jacobian) WCS at each object position to the PSF/star stamps and pass use_sky_coords=True to galsim FindAdaptiveMom, so PSFEx and MCCD HSM ellipticities and sizes come out already in world coordinates. This makes the downstream WCS-Jacobian shape rotation (convert_psf_pix2world) redundant. Threaded at every in-repo HSM PSF/star call site: - psfex_interp: multi-epoch (per-CCD WCS log), single-exposure validation and classic (WCS from the galcat FITS_LDAC image header via _galcat_wcs) - mccd_interpolation_script: multi-epoch (per-CCD WCS log) - shapepipe_auxiliary_mccd: interpolation (WCS from galcat LDAC header) local_wcs_list() builds the per-object local WCS: use_sky_coords evaluates the stamp WCS at true_center, so it must be the full WCS linearised at the object's image position, not at the stamp centre. SIGMA_*_HSM (hence HSM_T_* = 2 sigma^2) is now in sky units (arcsec) rather than pixels; this is the intended, physically correct size. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
HSM ellipticities and sizes are measured directly in world coordinates upstream, so the WCS-Jacobian shape rotation (transform_shape, getDecomposition) is redundant. Remove it from both the PSFEx and MCCD branches; the script keeps its remaining job of collating positions (X/Y/ RA/DEC), assigning the MCCD focal-plane CCD id, and merging per-exposure catalogues into per-patch star catalogues. Drops the galsim dependency. Caveat documented in the module docstring: the MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the external mccd fit-validation code, still measured in the pixel frame; they are now passed through unrotated. Update the CANFAR doc to describe the step as collation, not pix->world shape conversion. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
Equivalence test for use_sky_coords=True: draw an elliptical Gaussian on a stamp with a nontrivial local WCS (rotation + shear + 0.187 arcsec/pix scale, positive determinant), and assert that measuring in the pixel frame then applying the removed convert_psf_pix2world.transform_shape rotation (copied in verbatim as reference) agrees with measuring directly in sky coordinates, to numerical precision, for e1/e2 and sigma. A second test pins the deliberate unit change: the sky-frame sigma is the pixel sigma times the WCS scale. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
…ranch The MCCD PSF_MOM_LIST/STAR_MOM_LIST columns come from the external mccd fit-validation code (mccd.auxiliary_fun.mccd_validation), which still measures HSM adaptive moments in the pixel frame. Stripping the WCS-Jacobian rotation there (71d502f) would have handed MergeStarCatMCCD pixel-frame ellipticities. Restore transform_shape + the per-CCD getDecomposition rotation for the MCCD MOM_LIST branch only (and the galsim import it needs); the in-repo PSFEx and MCCD-interpolation paths stay unrotated, as they now measure in world coordinates upstream. Update the module docstring caveat: this branch keeps the rotation until mccd itself adopts use_sky_coords. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
local_wcs_list wraps galsim.AstropyWCS, which follows the 1-indexed FITS pixel convention, but the multi-epoch paths fed it positions from all_world2pix(..., 0) (origin 0) — evaluating the local Jacobian one pixel off. Compute a separate origin=1 position array (gal_pos_wcs) for local_wcs_list; the origin-0 gal_pos stays for the PSF-model interpolation (interpsfex / interp_MCCD), which is its only other use. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w
Contributor
|
We should revisit |
martinkilbinger
requested changes
Jul 14, 2026
cailmdaley
added a commit
that referenced
this pull request
Jul 16, 2026
Martin (PR #812 review): the patch concept is removed in v2.0. For -V v2.0 the handler reads summary/missing_job_32_all.txt and exp_numbers.txt relative to the run root instead of iterating P1..P9; -p restores an explicit P1..P<n> layout when needed. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr
…add patch-less v2.0 Re-scope the script to its actual remaining purpose: collating per-exposure PSF validation catalogues into star catalogues. The WCS-Jacobian shape rotation is gone (HSM moments are measured in sky coordinates upstream), so "convert" no longer describes it. The external-mccd PSF_MOM_LIST/STAR_MOM_LIST rotation branch stays until mccd adopts use_sky_coords. Add a -V/--version option mirroring get_ccds_with_psf. v1.x is unchanged (loop over -P patches, <base>/P<patch>/output, validation_psf_conv-<patch>-<idx>). v2.0 is patch-less: single <base>/output root, patch token dropped from the output filename (validation_psf_conv-<idx>). The downstream combine_runs.bash glob (validation_psf_conv-*) matches both. Path/filename logic extracted into pure collate_paths/output_filename helpers with unit tests; docs collate section updated with both invocations. The script command name follows the file (Dockerfile symlinks by basename). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01WqNFhq9uGpDpprYeuGeb8k
cailmdaley
force-pushed
the
feat/hsm-sky-coords
branch
from
July 20, 2026 19:51
79fafe0 to
baea03a
Compare
Contributor
Author
|
Thanks Martin — both points addressed.
— Fable on behalf of Cail |
cailmdaley
added a commit
that referenced
this pull request
Jul 20, 2026
Align with the collate_star_cat convention (PR #812): v2.0 (patch-less) is now the default catalogue version. Older, patched layouts are selected explicitly with -V. Help text and the pipeline_canfar.md invocations updated; the v2.0 call drops the flag, v1.x shows -V v1.6. Tests set the version explicitly, so none pinned the old default. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01WqNFhq9uGpDpprYeuGeb8k
Contributor
|
Hello, |
This was referenced Jul 21, 2026
Contributor
Author
|
Good point — agreed, and it turns out this was attempted before: #698 (fourth-order moments, PSFHOME-style) was closed unmerged and is now stale. It edits exactly the two functions this PR rewrites. We've reopened #697/#698 to track it and will implement as a follow-up on top of this PR: the spin-2 fourth-moment combinations measured in the sky frame (which removes the pixel→sky rotation the old code needed), plus galsim's spin-0 — Fable on behalf of Cail |
cailmdaley
added a commit
that referenced
this pull request
Sep 5, 2026
* removed tmp files * Moved get_ccds script to pipeline runner and library script * class to download eposure headers (for coverage) * coverage mask: scripts and plotter; using cs_util * coverage code: added Mike Hudson as co-author * documentation (in pipeline_canfar.md) on coverage mask scripts and classes * CCD PSF Handler: Added version for v1.3 * fixed coverage plots * header downloads (coverage): account for existing downloads * canfar-changes * PR: remove non-functioning coverage pipeline function * PR: coverage_map_builder; moved import into try block * PR: field_corners_extractor.py; fixed exp number extraction * PR: added comment * PR: added comment * refactor: consolidate field_corners_extractor helpers; fix regex - Extract _expnum_from_path, _parse_header_to_wcs, _megacam_field_corners as module-level helpers. process_single_header (static, for multiprocessing) and the instance methods now share one implementation. - Remove dead get_wcs_from_header and get_megacam_field (duplicated logic that had already been inlined into process_single_header). - Fix re.search(r'(\d+)\.txt') missing its subject string. The new helper also replaces the hard-coded p[end-6:end] exposure-number extraction in run(), so exposure numbers of any length work. - Declare verbose in params_default so the contract is local instead of relying on cs_util.args.parse_options to inject it. Addresses review feedback on #699. * refactor: coverage_map_builder cleanup - Replace hsp_map + 0 with hsp_map.copy() (healsparse supports .copy since at least 1.x; more explicit intent). - Declare verbose in params_default; currently populated only by cs_util.args.parse_options auto-injection. - Fix typo: "Install te cs_util" -> "Install the cs_util". Addresses review feedback on #699. * docs: declare verbose and annotate survey-specific defaults - coverage_plotter, header_downloader: declare verbose explicitly in params_default. - coverage_plotter: comment acknowledging that FootprintPlotter._regions is a private cs_util attribute; refactor to public accessor when cs_util exposes one. - header_downloader: expand the vos:cfis/pitcairn comment to flag it as UNIONS/CFIS-specific and document override. - build_and_plot_coverage_maps.sh: explain the BUILD_NSIDE=131072 choice (matches UNIONS bit-mask pixel scale, ~0.1"), and note the CoverageMapBuilder default of 2048 for lighter use. Addresses review feedback on #699. * chore: uv lock — register healsparse + skyproj as shapepipe deps The coverage feature adds healsparse (map building) and skyproj (plotting) to pyproject. healsparse was already resolved transitively via cs_util; this pins skyproj 2.5.0 and records both as shapepipe's own direct dependencies so 'uv lock --check' passes. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HvP7uKkfDU2PvWC7YmbuMW * fix: coverage console scripts handle -h; drop dead entry point The five coverage/get-ccds runners fell back to `args = sys.argv` (the full argv, including argv[0]) instead of `sys.argv[1:]`, so under the entry-point smoke test's `python -c "...; run_*()" -h` the parser read the program slot as the integer `-c` option and exited non-zero (`invalid integer value: '-h'`). Slice argv[1:] to match cs_util.args.parse_options' own convention. Also drop the `coverage_pipeline` console script, whose target `coverage_run:run_pipeline` was removed in an earlier commit (dangling entry point -> AttributeError on invocation), and declare `hpgeom` explicitly (imported directly by coverage_map_builder, previously only transitive via healsparse). Together these fix all six entry-point test failures. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HvP7uKkfDU2PvWC7YmbuMW * test: unit + property tests for the coverage-mask feature Cover the pure and lightly-fixtured logic: exposure-number parsing, the CCD-list -> unique-exposure reduction, the MegaCam field-corner convention, the multi-HDU header -> WCS split, power-of-two nside validation, and a regression guard for the -h/argv handling of the console runners. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01HvP7uKkfDU2PvWC7YmbuMW * fix(coverage): subtract missing CCDs in get_ccds_with_psf The v>=1.4 path built the full set of exposure single-HDU IDs but never removed the CCDs whose PSF model is missing, so every candidate CCD was reported as valid. Subtract the missing set so the result is the CCDs that actually have a valid PSF, and report the three counts. * feat(coverage): extract per-CCD corners with optional CCD-list filter Replace the whole-MegaCam-field corner extraction with per-CCD extraction: for each HDU's WCS, build the 4 corners from the CCD pixel bounds (read from pixel_shape, i.e. NAXIS1/NAXIS2) using pixel edges so the quadrilateral covers the full CCD area. Emit one row per CCD keyed on the '<expnum>-<ccd_idx>' ID. Add an optional --ccd_list/-l option: when given, only CCDs in that file (the output of get_ccds_with_psf) are written, so the coverage mask is built from single CCDs with a valid PSF. Resume still keys the done-set on exposure number, now parsed from the CCD-ID first column. Drop the dead _megacam_field_corners whole-field path; keep _parse_header_to_wcs and the multiprocessing structure. * feat(coverage): build maps from per-CCD rows with RA-wrap and pole guards Parse the new per-CCD row format (string CCD ID in column 0, then 8 float corner columns). Add an RA-wrap guard that unwraps polygons straddling the RA=0 seam onto a common branch, and a pole guard that skips polygons reaching |dec| >= 89 deg with a warning (CCD footprints never reach there). Since the CCDs of one exposure do not overlap, accumulating value 1 per CCD polygon still counts exposures-with-valid-PSF per pixel. * fix(coverage): download headers atomically via temp-file rename Copy each header to a '.part' file in the same directory and os.rename on success. An interrupted transfer can no longer leave a partial file that resume treats as complete; the temp file is removed on failure. * docs,build(coverage): align interfaces to the per-CCD flow Fix the console-script name typo (get_ccd_with_psfs -> get_ccds_with_psf) so it matches the docs and script. Harden the build/plot script with 'set -euo pipefail' and document the per-CCD chain (step 0 get_ccds_with_psf feeding extract_field_corners --ccd_list). Update the canfar docs: restore the coverage helpers as shipped, align filenames (ccds_with_psf_<version>.txt, exp_ra_dec_<version>.txt), and describe the per-CCD workflow including step 0. * test(coverage): cover per-CCD extraction, filtering, guards, atomic download Rewrite the coverage tests for the per-CCD design: handler missing-CCD subtraction; per-CCD corner extraction and pixel-shape error; --ccd_list filtering end to end; builder parsing of the new row format; RA-wrap unwrap (a CCD straddling RA=0 fills a small patch, not the complement); pole-guard skip; and atomic-download rename behaviour with a mocked copy. Drop the tests that pinned the removed whole-field convention. * fix(coverage): read true CCD dims from ZNAXIS and resume per-CCD Two correctness fixes in the per-CCD extractor: - Image dimensions: header_downloader fetches fpack tile-compressed headers where NAXIS1/NAXIS2 describe the compressed binary table (byte width, row count), so WCS.pixel_shape reports e.g. (8, 4612) instead of the true (2080, 4612). Read ZNAXIS1/ZNAXIS2 when ZIMAGE is set, falling back to NAXIS1/NAXIS2. _parse_header_to_wcs now returns (wcs, shape) pairs since the WCS drops the Z-keywords; _ccd_corners takes the shape explicitly. - Resume: key the done-set on individual CCD IDs, not exposure numbers. A write interrupted mid-exposure previously left a partial exposure that a -r rerun skipped forever (silent undercount); an expanded --ccd_list rerun skipped exposures whose new CCDs were never written. get_done_ccds reads the present CCD IDs and run() drops only those, so partial exposures are completed and new CCDs added, never duplicating a row. * fix(coverage): single-line CCD list, runner docstrings, plot doc example - header_downloader.get_exposures: wrap the text-format loadtxt in np.atleast_1d so a single-line CCD list (0-d array) does not crash on iteration. - coverage_run: update the extract/build runner docstrings from whole-field 'field corner coordinates' to per-CCD language. - pipeline_canfar docs: add a runnable plot_coverage_map example (the SGC region command) covering the flags build_and_plot_coverage_maps.sh uses. * test(coverage): compressed-header dims, resume path, nexp value contract Extend the coverage tests for the confirmed findings: - fpack-compressed header fixture; assert ZNAXIS dims are recovered and corners span the true CCD width, not the 8-byte compressed NAXIS1. - resume run() path: a new exposure is appended without duplicating existing rows, and a partial exposure (interrupted mid-write) is completed rather than skipped. - the nexp contract on pixel values: two overlapping CCDs from different exposures give value 2 in the overlap and 1 outside. - single-row builder file (atleast_1d/2d guards); single-line get_exposures; RA-wrap strengthened to a positional-parity check (seam CCD vs the same CCD shifted +10 deg agree in pixel count); get_ccds_with_psf subtraction pinned against the real get_all_shdus. * feat(coverage): accept catalogue version v2.0 (9-patch footprint, as v1.6) Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr * feat(coverage): v2.0 defaults to 10 patches; -p overrides the version map The v2 footprint is still growing (>5000 deg2 per the 2026-07 Paris meeting, new Dec<30 region); the newest patch definitions (auxdir/CFIS/tiles_202510) include P10. Default v2.0 to 10 patches and add an explicit -p/--n_patch override so the handler tracks evolving patch layouts without code changes. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr * fix(coverage): v2.0 is patch-less — read summaries from the run root Martin (PR #812 review): the patch concept is removed in v2.0. For -V v2.0 the handler reads summary/missing_job_32_all.txt and exp_numbers.txt relative to the run root instead of iterating P1..P9; -p restores an explicit P1..P<n> layout when needed. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr * docs: state the patch convention per catalogue version (v2.0 is patch-less) Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MqH9wjbJicVmaQdhZ4kBYr * feat(coverage): default get_ccds_with_psf -V to v2.0 (patch-less) Align with the collate_star_cat convention (PR #812): v2.0 (patch-less) is now the default catalogue version. Older, patched layouts are selected explicitly with -V. Help text and the pipeline_canfar.md invocations updated; the v2.0 call drops the flag, v1.x shows -V v1.6. Tests set the version explicitly, so none pinned the old default. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01WqNFhq9uGpDpprYeuGeb8k * chore: remove accidentally committed pytest log artifact Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * coverage: drop the v2.0 patch-less path — v2.0 coverage moves to the Snakemake stack Martin (review, 2026-09-04): the v2.0 pipeline is Snakemake-orchestrated, has no sky patches, and will not use the summary class or its outputs — so summary/missing_job_32_all.txt and exp_numbers.txt simply won't exist at a v2.0 run root. The patch-less path added in July was built on files that are not going to be there. This PR therefore lands as the v1.x tool. Reverts the v2.0 work (15401d9, a3aa4ec, ae10176, 46214c7, 1ffcd4b) and restores the pre-July contract of CcdPsfHandler: - accepted versions are v1.3, v1.4, v1.5, v1.6 again; -V v2.0 now fails loudly with "Invalid version v2.0" - -V default back to v1.6 (was flipped to v2.0 in 1ffcd4b) - the n_patch = 0 / patches = ["."] "patch-less" semantics are gone; patches are always P1..P<n> from the version map (7/7/8/9) - the -p/--n_patch CLI override is dropped: it was introduced in a3aa4ec purely to let the evolving v2.0 footprint be overridden, and the v1.x patch counts are fixed history Docs lose the "Patch convention by catalogue version" note and the patch-less get_ccds_with_psf invocation; the two v2.0 tests (test_v2_is_patchless, test_n_patch_option_overrides_version_default) are removed. The v1.x version->patch-count parametrisation and the invalid-version test stay. Everything else from the July rework is untouched: per-CCD extraction, RA-wrap and pole guards, atomic downloads, resume, ZNAXIS dims. v2.0 coverage will be a separate PR against the Snakemake stack. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com> * chore: relock cs_util to develop 1b15a55 (restores FootprintPlotter.plot_area) --------- Co-authored-by: martinkilbinger <martinkilbinger@cea.fr> Co-authored-by: Cail Daley <cail.daley@cea.fr> Co-authored-by: Claude Fable 5 <noreply@anthropic.com> Co-authored-by: Cail McLean Daley <cailmdaley@gmail.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #775.
ShapePipe measured PSF/star HSM adaptive moments in the pixel frame, then rotated the ellipticities into sky coordinates in a separate post-processing step (
scripts/python/convert_psf_pix2world.py, a WCS-Jacobian decomposition). galsim'sFindAdaptiveMomnow does this at the point of measurement viause_sky_coords=True(available since galsim 2.5; we pin ≥2.8) — per @aguinot's suggestion on #741. This PR moves the WCS knowledge into the HSM call and drops the redundant rotation.What changed
FindAdaptiveMomcall site — PSFEx single-exposure and multi-epoch (psfex_interp.py), MCCD interpolation (mccd_interpolation_script.py,shapepipe_auxiliary_mccd.py) — now attaches the local WCS (linearised at the object's image position) to the stamp and passesuse_sky_coords=True. Multi-epoch paths use the per-CCD WCS log already in scope; single-exposure paths read the WCS from the galcat'sLDAC_IMHEADheader.convert_psf_pix2world.pyno longer rotates the in-repoE1/E2/SIGMA_*_HSMshapes. Per review, it is renamed tocollate_star_cat.pyto reflect its remaining purpose: collate the per-exposure PSF validation catalogues into star catalogues (positions, MCCD CCD id, merge).v2.0support. The script takes-V(defaultv2.0), following theget_ccds_with_psfconvention from Coverage nexp masks #797. the patch-less default reads from the single run root and writesvalidation_psf_conv-<idx>.fits(no patch token; downstreamcombine_runs.bash/config_Ms_psfex_conv.iniglobs match both schemes).-V v1.6 -P …selects the unchanged v1.x behaviour. Invalid versions raise. The collate section ofpipeline_canfar.mdshows both invocations.PSF_MOM_LIST/STAR_MOM_LISTcolumns come from the externalmccdlibrary, which still measures in the pixel frame — that one branch keeps its rotation (with a docstring caveat) untilmccditself adoptsuse_sky_coords.tests/module/test_hsm_sky_coords.py: a synthetic Gaussian stamp under a nontrivial local WCS (0.187″/pix, 28° rotation, shear), measured pixel-frame + oldtransform_shapemath vsuse_sky_coords=True— g1/g2/sigma agree to rtol 1e-6. Plustest_collate_star_cat.pypinning the v1.x/v2.0 path and filename logic. All 14 tests pass in the develop container.Notes for review
moments_sigmanow comes back in world units (arcsec) rather than pixels, soHSM_T_* = 2σ²is arcsec² — which is what the old pipeline produced after the rotation step, so end-of-pipeline catalogues are unchanged. Downstream consumers (merge_starcat, ngmix T-guess) checked: no pixel-scale assumption.getDecompositionto the HSM ellipticities;merge_starcatpasses them through.f153ea28).An independent adversarial review of the full diff ran before opening this PR; both findings it raised (the external-mccd branch above, and the indexing fix) are addressed.
— Claude on behalf of Cail
🤖 Generated with Claude Code
https://claude.ai/code/session_01Gao7Pi5QDBgnTuagV7HT2w