From 702c804120cadfd0c77f3dc6fb1a67a5ad7aecf5 Mon Sep 17 00:00:00 2001 From: Thomas Robitaille Date: Tue, 24 Mar 2026 11:08:19 +0000 Subject: [PATCH] Add a new Montage accuracy test in ICRS since inaccuracies in existing test were due to conversion to/from Galactic --- .../tests/test_reproject.py | 71 ++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/reproject/spherical_intersect/tests/test_reproject.py b/reproject/spherical_intersect/tests/test_reproject.py index e3a05c883..b26d24d5e 100644 --- a/reproject/spherical_intersect/tests/test_reproject.py +++ b/reproject/spherical_intersect/tests/test_reproject.py @@ -66,6 +66,54 @@ def test_reproject_celestial_slices_2d(): ] ) +# Same geometry as INPUT_HDR/OUTPUT_HDR but entirely in ICRS (no Galactic conversion). +# This isolates the spherical polygon intersection accuracy from any differences +# in Galactic coordinate conversions between astropy and Montage. + +INPUT_HDR_ICRS = """ +WCSAXES = 2 / Number of coordinate axes +CRPIX1 = 299.628 / Pixel coordinate of reference point +CRPIX2 = 299.394 / Pixel coordinate of reference point +CDELT1 = -0.001666666 / [deg] Coordinate increment at reference point +CDELT2 = 0.001666666 / [deg] Coordinate increment at reference point +CUNIT1 = 'deg' / Units of coordinate increment and value +CUNIT2 = 'deg' / Units of coordinate increment and value +CTYPE1 = 'RA---CAR' / Right ascension, plate caree projection +CTYPE2 = 'DEC--CAR' / Declination, plate caree projection +CRVAL1 = 0.0 / [deg] Coordinate value at reference point +CRVAL2 = 0.0 / [deg] Coordinate value at reference point +LONPOLE = 180.0 / [deg] Native longitude of celestial pole +LATPOLE = 0.0 / [deg] Native latitude of celestial pole +EQUINOX = 2000.0 / [yr] Equinox of equatorial coordinates +""" + +OUTPUT_HDR_ICRS = """ +WCSAXES = 2 / Number of coordinate axes +CRPIX1 = 2.5 / Pixel coordinate of reference point +CRPIX2 = 2.5 / Pixel coordinate of reference point +CDELT1 = -0.001500000 / [deg] Coordinate increment at reference point +CDELT2 = 0.001500000 / [deg] Coordinate increment at reference point +CUNIT1 = 'deg' / Units of coordinate increment and value +CUNIT2 = 'deg' / Units of coordinate increment and value +CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection +CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection +CRVAL1 = 0.4960464682000 / [deg] Coordinate value at reference point +CRVAL2 = -0.4956564684000 / [deg] Coordinate value at reference point +LONPOLE = 180.0 / [deg] Native longitude of celestial pole +LATPOLE = -0.4956564684000 / [deg] Native latitude of celestial pole +EQUINOX = 2000.0 / [yr] Equinox of equatorial coordinates +""" + +# Reference values from Montage mProject (run with -f -x 1.0 flags) +MONTAGE_ICRS_REF = np.array( + [ + [1.0, 1.5555337416634756, 2.0, np.nan], + [2.111112839105374, 2.6666450687397223, 3.1111069641586675, np.nan], + [3.0, 3.5555337849185467, 4.0, np.nan], + [np.nan, np.nan, np.nan, np.nan], + ] +) + @pytest.mark.parametrize("wcsapi", (False, True)) def test_reproject_celestial_montage(wcsapi): @@ -79,10 +127,31 @@ def test_reproject_celestial_montage(wcsapi): array, footprint = _reproject_celestial(DATA, wcs_in, wcs_out, (4, 4)) - # TODO: improve agreement with Montage - at the moment agreement is ~10% + # The ~10% disagreement with Montage is due to differences in the + # Galactic-to-ICRS coordinate transformation between astropy and Montage, + # not the polygon intersection algorithm (see test_reproject_celestial_icrs_montage + # which agrees to ~1e-6 when no Galactic conversion is involved). np.testing.assert_allclose(array, MONTAGE_REF, rtol=0.09) +@pytest.mark.parametrize("wcsapi", (False, True)) +def test_reproject_celestial_icrs_montage(wcsapi): + # Accuracy compared to Montage, using only ICRS coordinates (no Galactic + # conversion). This isolates the polygon intersection algorithm from + # any differences in Galactic-to-ICRS coordinate transformations between + # astropy and Montage. + + wcs_in = WCS(fits.Header.fromstring(INPUT_HDR_ICRS, sep="\n")) + wcs_out = WCS(fits.Header.fromstring(OUTPUT_HDR_ICRS, sep="\n")) + + if wcsapi: + wcs_in, wcs_out = as_high_level_wcs(wcs_in), as_high_level_wcs(wcs_out) + + array, footprint = _reproject_celestial(DATA, wcs_in, wcs_out, (4, 4)) + + np.testing.assert_allclose(array, MONTAGE_ICRS_REF, rtol=1.0e-5) + + def test_reproject_flipping(): # Regression test for a bug that caused issues when the WCS was oriented # in a way that meant polygon vertices were clockwise.