Flag out NaN weights in coadd & allow for weights to have different WCS than data #474

keflavich · 2024-10-01T03:15:57Z

I encountered a severe error case in which the weights array became NaN after reprojection. Flagging out locations where weights is NaN solved the issue.

I cannot come up with an MWE; this occurs deep in the guts of a complicated mosaicking attempt, and it happened for only 1 field out of 270. I think it might be the consequence of a numerical error that occurs only under unique conditions, but it also looks logically like this fix should work.

codecov · 2024-10-01T03:20:48Z

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 91.46%. Comparing base (9f4ad1c) to head (c83e4f1).
Report is 3 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #474      +/-   ##
==========================================
+ Coverage   91.00%   91.46%   +0.45%     
==========================================
  Files          25       25              
  Lines        1078     1089      +11     
==========================================
+ Hits          981      996      +15     
+ Misses         97       93       -4

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

keflavich · 2025-01-25T16:03:13Z

Currently broken:

FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[interp-False-arrays] - TypeError: input_data should either be an HDU object or a tuple of (array, WCS) or (array, Header)
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[interp-False-filenames] - ValueError: Output WCS has celestial components but input WCS does not
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[interp-False-hdus] - ValueError: Output WCS has celestial components but input WCS does not
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[interp-True-arrays] - TypeError: input_data should either be an HDU object or a tuple of (array, WCS) or (array, Header)
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[interp-True-filenames] - ValueError: Output WCS has celestial components but input WCS does not
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[interp-True-hdus] - ValueError: Output WCS has celestial components but input WCS does not
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[exact-False-arrays] - TypeError: input_data should either be an HDU object or a tuple of (array, WCS) or (array, Header)
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[exact-False-filenames] - NotImplementedError: Currently only data with a 2-d celestial WCS can be reprojected using flux-conserving algorithm
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[exact-False-hdus] - NotImplementedError: Currently only data with a 2-d celestial WCS can be reprojected using flux-conserving algorithm
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[exact-True-arrays] - TypeError: input_data should either be an HDU object or a tuple of (array, WCS) or (array, Header)
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[exact-True-filenames] - NotImplementedError: Currently only data with a 2-d celestial WCS can be reprojected using flux-conserving algorithm
FAILED reproject/mosaicking/tests/test_coadd.py::TestReprojectAndCoAdd::test_coadd_with_weights[exact-True-hdus] - NotImplementedError: Currently only data with a 2-d celestial WCS can be reprojected using flux-conserving algorithm
FAILED reproject/mosaicking/tests/test_coadd.py::test_coadd_solar_map - TypeError: input_data should either be an HDU object or a tuple of (array, WCS) or (array, Header)

keflavich · 2025-01-25T16:32:01Z

Tests now pass locally (w/o sunpy installed...)

keflavich · 2025-01-26T15:40:39Z

@astrofrog could you review this? There are several critical fixes for spectral-cube's cube mosaicking.

Major changes are:

Weights are no longer assumed to be on the same grid as the data, and they can bring their own WCS
input and output array dtypes are not assumed to be the same; input is swapped to match output (under the assumption that that's cheaper - but perhaps this needs to be made optional / documented)
bad weights are removed, which was the original point of this PR, but might have been a side-effect of the first

keflavich · 2025-01-26T18:30:24Z

For coverage checks: Do we have a dask fits reader, or could we add that in? The byte swapping does get tested in spectral-cube tests, since those sometimes load different endian arrays, but I don't see a simple way to do that here.

astrofrog

This looks mostly good to me, but one thing I'm not sure about is whether the byteswapping is happening too early here. Could you explain what issues you have been seeing here related to this? The byteswapping should probably ideally happen inside each individual reprojection chunk in case the input isn't a dask array and to avoid doubling memory usage?

keflavich · 2025-01-30T12:39:15Z

in case the input isn't a dask array

Technically I think this might be needed for numpy arrays too, but the FITS reader for numpy arrays already handles that?

In any case, I'll report the issues and see if we can move the byteswapping. I haven't worked up a minimal example for the dask byteswapped data; iirc, it was any dask array loaded with a fits reader

astrofrog · 2025-01-30T12:48:52Z

@keflavich - just to check though, why is the byteswapping needed at all? Where does it crash if we don't byteswap?

keflavich · 2025-01-30T13:00:17Z

it doesn't crash. It interprets the data as having the wrong endianness and all values end up as 2e-205 and similar.

astrofrog · 2025-01-30T14:16:59Z

Ok interesting, is that making use of reproject_interp?

keflavich · 2025-01-30T14:26:58Z

yes

keflavich · 2025-01-31T20:40:47Z

so, I didn't figure out how to build a MWE within reproject, but if you comment out the _byteordermatch call and run the spectral cube tests, you'll get errors like this:

spectral_cube/tests/test_regrid.py ...................................................F.FFF............F.F.F.F.......                                                                                                                   [100%]

================================================================================================================== FAILURES ===================================================================================================================
_____________________________________________________________________________________________________ test_mosaic_cubes[True-100-False0] ______________________________________________________________________________________________________

use_memmap = False, data_adv = PosixPath('/private/var/folders/k_/7qh4l0nn72b7qgq15pkd4hw40000gt/T/pytest-of-adam/pytest-171/test_mosaic_cubes_True_100_Fal0/adv.fits'), use_dask = True, spectral_block_size = 100

    @pytest.mark.parametrize('spectral_block_size,use_memmap', ((None, False),
                                                                (100, False),
                                                                (None, True),
                                                                (100, False),
                                                                (1, True),
                                                                (1, False),
                                                                ))
    def test_mosaic_cubes(use_memmap, data_adv, use_dask, spectral_block_size):

        pytest.importorskip('reproject')

        # Read in data to use
        cube, data = cube_and_raw(data_adv, use_dask=use_dask)

        # cube is doppler-optical by default, which uses the rest wavelength,
        # which isn't auto-computed, resulting in nan pixels in the WCS transform
        cube._wcs.wcs.restwav = constants.c.to(u.m/u.s).value / cube.wcs.wcs.restfrq

        expected_header = combine_headers(cube.header, cube.header)
        expected_wcs = WCS(expected_header).celestial

        # Make two overlapping cubes of the data
        part1 = cube[:, :round(cube.shape[1]*2./3.), :]
        part2 = cube[:, round(cube.shape[1]/3.):, :]

        assert part1.wcs.wcs.restwav != 0
        assert part2.wcs.wcs.restwav != 0

        # Mosaic give the expected header.
        result = mosaic_cubes([part1, part2],
                              # order is not a mosaic_cubes argument order='nearest-neighbor',
                              target_header=cube.header,
                              # not used roundtrip_coords=False,
                              spectral_block_size=spectral_block_size,
                              save_to_tmp_dir=False,
                              verbose=True,
                              use_memmap=use_memmap,
                              method='cube')

        # Check that the shapes are the same
        assert result.shape == cube.shape

        assert repr(cube.wcs.celestial) == repr(result.wcs.celestial)

        # Check WCS in reprojected matches wcs_out
        # (comparing WCS failed for no reason we could discern)
        assert repr(expected_wcs) == repr(result.wcs.celestial)

        # Check that values of original and result are comparable
>       np.testing.assert_almost_equal(result.filled_data[:].value,
                                       cube.filled_data[:].value, decimal=9)

spectral_cube/tests/test_regrid.py:651:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
../../mambaforge/envs/py312/lib/python3.12/contextlib.py:81: in inner
    return func(*args, **kwds)
../../mambaforge/envs/py312/lib/python3.12/contextlib.py:81: in inner
    return func(*args, **kwds)
../../mambaforge/envs/py312/lib/python3.12/site-packages/numpy/_utils/__init__.py:85: in wrapper
    return fun(*args, **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _

args = (<function assert_array_almost_equal.<locals>.compare at 0x284f77920>, array([[[ 2.68308185e+185, -2.51563734e+252],
 ...0.39923534]],

       [[0.74684203, 0.85496775],
        [0.65468046, 0.96629004],
        [0.64365533, 0.74426717]]]))
kwds = {'err_msg': '', 'header': 'Arrays are not almost equal to 9 decimals', 'precision': 9, 'verbose': True}

    @wraps(func)
    def inner(*args, **kwds):
        with self._recreate_cm():
>           return func(*args, **kwds)
E           AssertionError:
E           Arrays are not almost equal to 9 decimals
E
E           Mismatched elements: 24 / 24 (100%)
E           Max absolute difference among violations: 3.89622506e+287
E           Max relative difference among violations: 4.55716028e+287
E            ACTUAL: array([[[ 2.683081851e+185, -2.515637336e+252],
E                   [-9.039787586e+039,  3.043705986e+222],
E                   [-1.223304808e+051,  2.728961237e+167]],...
E            DESIRED: array([[[0.214686889, 0.909502371],
E                   [0.817235373, 0.903880475],
E                   [0.70687927 , 0.118076201]],...

../../mambaforge/envs/py312/lib/python3.12/contextlib.py:81: AssertionError
------------------------------------------------------------------------------------------------------------ Captured stdout call -------------------------------------------------------------------------------------------------------------
INFO: Using memory [spectral_cube.cube_utils]
INFO: Using Cube method [spectral_cube.cube_utils]
-------------------------------------------------------------------------------------------------------------- Captured log call --------------------------------------------------------------------------------------------------------------
INFO     astropy:cube_utils.py:1004 Using memory
INFO     astropy:cube_utils.py:1004 Using Cube method

astrofrog · 2025-02-01T09:06:27Z

Ok perfect, thanks!

astrofrog · 2025-02-01T10:07:18Z

Ok I have a MWE with just reproject_interp, so I'll see how we can fix it:

import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.utils.data import get_pkg_data_filename
from reproject import reproject_interp
from dask import array as da

hdu1 = fits.open(get_pkg_data_filename("galactic_center/gc_2mass_k.fits"))[0]
hdu2 = fits.open(get_pkg_data_filename("galactic_center/gc_msx_e.fits"))[0]

output_array = np.zeros(
    shape=(hdu1.header["NAXIS2"], hdu1.header["NAXIS1"]), dtype=">f8"
)

data = da.from_array(hdu2.data)
wcs = WCS(hdu2.header)

reproject_interp(
    (data, wcs), hdu1.header, output_array=output_array, block_size=(50, 50)
)

print(np.nanmin(output_array), np.nanmax(output_array))

and it happens for reproject_exact too, so this is something that needs fixing in the generic dispatcher.

astrofrog · 2025-02-01T13:42:45Z

Ok I've spotted the bug, I will open a separate PR shortly

astrofrog · 2025-02-01T13:53:25Z

@keflavich - see #487

astrofrog · 2025-02-01T20:02:42Z

You should be able to rebase against main now as I have merged the other PR

astrofrog · 2025-02-02T00:55:24Z

(Once you rebase this I will review and merge ASAP)

…e very rare, but it did happen.

…ped data

keflavich · 2025-02-02T03:08:34Z

rebase done

astrofrog

This looks mostly good but a couple of comments - in particular we shouldn't break APE-14 support. I don't think you need to check has_celestial, if the weight's WCS isn't compatible with the mosaic WCS then an error will happen during reprojection anyway. We shouldn't ignore the WCS just because it's not celestial?

astrofrog · 2025-02-03T14:07:18Z

reproject/mosaicking/coadd.py

+                weights_in, weights_wcs = parse_input_weights(
+                    input_weights[idata], hdu_weights=hdu_weights, return_wcs=True
+                )
+                if weights_wcs is None or not weights_wcs.has_celestial:


weights_wcs.has_celestial breaks support for APE 14 WCS

astrofrog · 2025-02-03T14:07:44Z

reproject/mosaicking/coadd.py

@@ -211,7 +211,12 @@ def reproject_and_coadd(
            if input_weights is None:
                weights_in = None
            else:
-                weights_in = parse_input_weights(input_weights[idata], hdu_weights=hdu_weights)
+                weights_in, weights_wcs = parse_input_weights(
+                    input_weights[idata], hdu_weights=hdu_weights, return_wcs=True


I think this is the only place parse_input_weights is used so I don't think we need the kwarg, we can just always return the WCS

fine by me, i can remove the kwarg

astrofrog · 2025-02-03T14:10:18Z

reproject/utils.py

-        return input_weights.data
+        if return_wcs:
+            ww = WCS(input_weights.header)
+            ww = ww if ww.has_celestial else None


Here we again break compatibility with APE 14 WCS - why is it needed?

This check is serving a critical role and cannot be removed, so we have to think of a replacement.

This check is asking whether the FITS header contains a WCS that should be interpreted as the WCS for the weights. Since WCS() (i.e., WCS with an empty dict as input) is still a valid WCS object but contains no transforms, we need some check here whether the WCS exists. Annoyingly, you can't just check that naxis matches ndim, since naxis defaults to 2.

Maybe...

ww = ww if ww.array_shape == input_weights.data.shape

...I'll check this...

no, I don't think that works because the shape comes from the data when reading.

what is the APE14 way to check if there is a transformation associated with an axis?

I don't think there's any general way to get at this. All WCSes contain valid transforms, so the only solution is to use the full WCS validation. I've added that now. It's less clear for specific cases, but at least it solves the general case.

keflavich · 2025-02-03T14:14:02Z

has_celestial isn't necessary, but we need something to check that there is a valid WCS. The default is just WCS(), and it is possible (likely) for a user to pass a FITS HDU with a weight array and no WCS information in the header. What's the appropriate way to do this check under APE14?

keflavich · 2025-02-03T14:15:16Z

Your suggestion seems to be to just let it break downstream if there is an invalid WCS? That's fine, I guess. I'm not a fan of that strategy in general; I like to catch obvious errors where I know they'll occur.

astrofrog · 2025-02-03T14:45:27Z

@keflavich - I think reproject_and_coadd can be used for any arbitrary WCS that might not even be celestial? Or is that incorrect? In that sense, even ignoring APE 14 it does not make sense to restrict this?

keflavich · 2025-02-03T14:54:14Z

right.... the cube modifications made that possible. OK, I'll get rid of has_celestial, but I think the error messages may be opaque

keflavich · 2025-02-03T14:59:50Z

This fails now with:

E           NotImplementedError: Currently only data with a 2-d celestial WCS can be reprojected using flux-conserving algorithm

should we just skip these tests?

keflavich · 2025-02-03T15:01:52Z

Also, ValueError: Output WCS has celestial components but input WCS does not. That's the expected error; I haven't looked more closely yet but at least the error is clear.

astrofrog · 2025-02-05T13:25:02Z

pre-commit.ci autofix

for more information, see https://pre-commit.ci

keflavich mentioned this pull request Oct 1, 2024

Giant Cube: HCNO (12m + 7m + TP) ACES-CMZ/reduction_ACES#386

Open

keflavich force-pushed the flag_out_nan_weights branch from 2f67e7d to 280472f Compare January 24, 2025 20:11

This was referenced Jan 24, 2025

For weighted reproject_and_coadd, don't assume weight WCS is same as input #486

Open

Cube mosaicing: Expand capabilities substantially radio-astro-tools/spectral-cube#868

Draft

keflavich mentioned this pull request Jan 26, 2025

Continuum mosaic: Edge/noise problems between fields ACES-CMZ/reduction_ACES#411

Open

keflavich mentioned this pull request Jan 26, 2025

WIP: Cubewcs mosaic and dask reproject radio-astro-tools/spectral-cube#894

Open

astrofrog reviewed Jan 30, 2025

View reviewed changes

keflavich mentioned this pull request Feb 1, 2025

Fix a bug that caused big-endian dask arrays to not be reprojected correctly #487

Merged

keflavich force-pushed the flag_out_nan_weights branch from a28250a to 0201830 Compare February 1, 2025 18:16

keflavich changed the title ~~Flag out NaN weights in coadd~~ Flag out NaN weights in coadd & allow for weights to have different WCS than data Feb 1, 2025

keflavich added 3 commits February 1, 2025 21:19

Remove NaN weights. I do not know how this happens, and it seems to b…

4e57318

…e very rare, but it did happen.

move "reset" earlier

5249326

issue 486: get the WCS from the weights

d4008bf

keflavich and others added 8 commits February 1, 2025 21:19

halfway to fixing the errors

3d73a61

fixed tests

e8fc9ac

jump through some insane hoops to make dask arrays work with byteswap…

33e119a

…ped data

add new weight test that isn't on the same grid

0cbe895

add another test for the weights parser

a72f957

Fix codestyle

7e9dd28

Fix codestyle

54e60d3

remove byteordermatch

8e42fba

keflavich force-pushed the flag_out_nan_weights branch from 0201830 to 8e42fba Compare February 2, 2025 02:19

astrofrog requested changes Feb 3, 2025

View reviewed changes

remove has_celestial checks and always return_wcs

21ae34a

check that WCS is valid before using it

3d4dd90

[pre-commit.ci] auto fixes from pre-commit.com hooks

c83e4f1

for more information, see https://pre-commit.ci

Flag out NaN weights in coadd & allow for weights to have different WCS than data #474

Are you sure you want to change the base?

Flag out NaN weights in coadd & allow for weights to have different WCS than data #474

Conversation

keflavich commented Oct 1, 2024

codecov bot commented Oct 1, 2024 • edited Loading

Codecov Report

keflavich commented Jan 25, 2025

keflavich commented Jan 25, 2025

keflavich commented Jan 26, 2025

keflavich commented Jan 26, 2025

astrofrog left a comment

Choose a reason for hiding this comment

keflavich commented Jan 30, 2025

astrofrog commented Jan 30, 2025

keflavich commented Jan 30, 2025

astrofrog commented Jan 30, 2025

keflavich commented Jan 30, 2025

keflavich commented Jan 31, 2025

astrofrog commented Feb 1, 2025

astrofrog commented Feb 1, 2025 • edited Loading

astrofrog commented Feb 1, 2025

astrofrog commented Feb 1, 2025

astrofrog commented Feb 1, 2025

astrofrog commented Feb 2, 2025

keflavich commented Feb 2, 2025

astrofrog left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

keflavich commented Feb 3, 2025

keflavich commented Feb 3, 2025

astrofrog commented Feb 3, 2025 • edited Loading

keflavich commented Feb 3, 2025

keflavich commented Feb 3, 2025

keflavich commented Feb 3, 2025

astrofrog commented Feb 5, 2025

codecov bot commented Oct 1, 2024 •

edited

Loading

astrofrog commented Feb 1, 2025 •

edited

Loading

astrofrog commented Feb 3, 2025 •

edited

Loading