reproject_adaptive collapses for very large shape_out: returns all-NaN + zero footprint

[Joshuaalbert]

I’m seeing a deterministic size-dependent failure in reproject_adaptive when shape_out is very large (~1.8e9 pixels). Even with a trivial source plane (all ones) and clearly overlapping WCS, the function returns either:

1. plane all-NaN and footprint all-false (sum=0), or
2. plane all-zero and footprint has exactly 1 true pixel (sum=1)

No exception is raised.

This makes it hard to detect upstream because the output looks “valid” from an API perspective but is physically wrong.

Expected behavior

Either, A) a correct reprojection with nonzero footprint and finite plane values, or B) a raised exception indicating internal limitation (memory / indexing / shape limit)

Ideally A

Minimal verifiable complete example

This reproducer constructs a target WCS that clearly overlaps with source WCS, and uses an all-ones source plane.

import sys
import numpy as np
import astropy
from astropy import wcs
from reproject import reproject_adaptive
import reproject


def make_target_wcs(num_l: int, num_m: int, dl_deg: float, dm_deg: float, ra0_deg: float, dec0_deg: float):
    # 4D WCS like our pipeline, then use .celestial
    target_w = wcs.WCS(naxis=4)
    target_w.wcs.ctype = ["RA---SIN", "DEC--SIN", "FREQ", "STOKES"]
    target_w.wcs.cunit = ["deg", "deg", "Hz", ""]
    # NOTE: reference pixel placement is intentional and matches our pipeline:
    target_w.wcs.crpix = [num_l / 2.0, num_m / 2.0 + 1.0, 1.0, 1.0]
    target_w.wcs.cdelt = [-abs(dl_deg), +abs(dm_deg), 1.0, 1.0]
    target_w.wcs.crval = [ra0_deg, dec0_deg, 1.0, 1.0]  # STOKES=I code is 1
    return target_w.celestial


def make_source_wcs(num_l_src: int, num_m_src: int, dl_deg: float, dm_deg: float, ra0_deg: float, dec0_deg: float):
    src_w = wcs.WCS(naxis=2)
    src_w.wcs.ctype = ["RA---SIN", "DEC--SIN"]
    src_w.wcs.cunit = ["deg", "deg"]
    src_w.wcs.crpix = [(num_l_src + 1) / 2.0, (num_m_src + 1) / 2.0]
    src_w.wcs.cdelt = [-abs(dl_deg), +abs(dm_deg)]
    src_w.wcs.crval = [ra0_deg, dec0_deg]
    return src_w


def run_case(num_l_out: int, num_m_out: int, label: str):
    ra0_deg = 60.0
    dec0_deg = 0.0

    # ~1 arcsec/pixel
    dl_deg = (1.009918463060501 / 3600.0)
    dm_deg = (1.009918463060501 / 3600.0)

    # Big-enough source plane so target overlap is guaranteed
    num_l_src = 81920
    num_m_src = 81920

    src_wcs = make_source_wcs(num_l_src, num_m_src, dl_deg, dm_deg, ra0_deg, dec0_deg)
    tgt_wcs = make_target_wcs(num_l_out, num_m_out, dl_deg, dm_deg, ra0_deg, dec0_deg)

    # All ones input (dec, ra)
    src_plane = np.ones((num_m_src, num_l_src), dtype=np.float32)

    out_plane, footprint = reproject_adaptive(
        (src_plane, src_wcs),
        tgt_wcs,
        kernel="gaussian",
        conserve_flux=True,
        shape_out=(num_m_out, num_l_out),
        return_footprint=True,
        parallel=False,
    )

    fp_sum = float(np.sum(footprint))
    fp_size = int(footprint.size)
    out_nan = int(np.isnan(out_plane).sum())
    out_max = np.nanmax(out_plane)
    out_min = np.nanmin(out_plane)

    print(f"\n[{label}] shape_out=(num_m,num_l)=({num_m_out},{num_l_out})")
    print(f"  footprint sum/size: {fp_sum}/{fp_size}")
    print(f"  out NaN count: {out_nan} / {out_plane.size}")
    print(f"  out min/max: {out_min} / {out_max}")


def main():
    print("Python:", sys.version)
    print("astropy:", astropy.__version__)
    print("reproject:", reproject.__version__)

    # Control case: smaller (works)
    run_case(17856, 17856, "small-square")

    # Large cases: failures observed (either fp_sum==0 and NaNs, or fp_sum==1 and zeros)
    run_case(43274, 42739, "large-nonsquare-odd")
    run_case(43274, 42738, "large-nonsquare-even")


if __name__ == "__main__":
    main()

Output

Python: 3.11.14 (main, Dec  8 2025, 23:47:06) [GCC 12.2.0]
astropy: 7.2.0
reproject: 0.19.0

[small-square] shape_out=(num_m,num_l)=(17856,17856)
  footprint sum/size: 318836736.0/318836736
  out NaN count: 0 / 318836736
  out min/max: 0.9999999998799467 / 1.0000000001127773
/tmp/ipykernel_1664238/1971168938.py:62: RuntimeWarning: All-NaN slice encountered
  out_max = np.nanmax(out_plane)
/tmp/ipykernel_1664238/1971168938.py:63: RuntimeWarning: All-NaN slice encountered
  out_min = np.nanmin(out_plane)

[large-nonsquare-odd] shape_out=(num_m,num_l)=(42739,43274)
  footprint sum/size: 0.0/1849487486
  out NaN count: 1849487486 / 1849487486
  out min/max: nan / nan

[large-nonsquare-even] shape_out=(num_m,num_l)=(42738,43274)
  footprint sum/size: 1.0/1849444212
  out NaN count: 1849444211 / 1849444212
  out min/max: 0.0 / 0.0

Observe

1. small-square: nonzero footprint and finite values.
2. large-nonsquare-odd: footprint sum == 0, output all NaNs.
3. large-nonsquare-even: footprint sum == 1, output max == 0.

[astrofrog]

Thanks for the bug report! @svank would you be happy to take a look at this one?

[svank]

Yes, I'll take care of this next week

[svank]

It looks like this is a problem in the coordinates work (probably wcslib?), not reproject. I'll pull this together nicely later, but the problem occurs in the coordinates work, where we compute a world coordinate for every pixel in the output image, and then convert to a pixel coordinate in the input image. tgt_wcs.pixel_to_world(x, y), when x and y both have shape (43276, 42741) produces an array that has the correct value for the pixel at [0, 0] but is zeros everywhere else. The zeros get processed to the same value in src_wcs.world_to_pixel, so pixel_to_pixel itself gives an array of almost entirely the same constant value. In reproject's coordinate-roundtrip check, those values don't translate back to the original values and so get replaced with NaNs, leading to the broken outputs. (I haven't tried it, but with roundtrip_coords=False, instead of being NaN the coordinate mapping will probably just be wrong.)

The problem seems to depend on the number of coordinates being computed at once---tgt_wcs.pixel_to_world(x, y) fails, but both tgt_wcs.pixel_to_world(x[::2], y[::2]) and tgt_wcs.pixel_to_world(x[:, ::2], y[:, ::2]) succeed (i.e. skipping every other row or every other column), so it's not the range of values in the array or the size of a particular dimension, but the number of values.

The number of pixels (~1.8e9) is less than but suspiciously close to the max value of a signed int32 (~2.1e9). But if in some layer we have one array containing both x and y, or both ra and dec, for every pixel, then that's more values than an int32 can count. Very suspicious.

So in summary, the failure seems to be because of the number of pixels. @Joshuaalbert A probable workaround for now would be to reproject the image in chunks (the full source image, but only a portion of the output image for each chunk)---maybe just two halves, and if you have the halves overlap in the middle of the image (by maybe a 100 pixels should be more than plenty), you can just crop off overlap region of each half and not worry about edge effects.

As an aside, the example code needs > 250 GB of memory to run---the third machine I tried could do it! I'm really curious what data you're processing!

[Joshuaalbert]

The application is radio interferometric forward modelling for the DSA. Forming a single sky model that covers a mosaic area in a single tangent plane at mosaic center. We then simulate different sections of it for each pointing.

[astrofrog]

It should be possible to do the chunking automatically in reproject by the way, will check later!

[astrofrog]

@Joshuaalbert could you try and set for example block_size=(4096,4096) in the call to reproject_adaptive? Even if you disable parallel mode as you have done, operating in chunks will avoid too large arrays and you should be able to find a chunk size that doesn't impact performance negatively.

[astrofrog]

I mean to be clear we do have to fix the underlying issue here, but what I mean is that it would be optimal to operate in blocks/chunks anyway.

[Joshuaalbert]

I did use chunking, though I wrote it myself (hopefully I didn't make a mistake). I didn't realise that the block_size argument would work when parallel=False.

[astrofrog]

Yes it should work even with parallel set to False, I should clarify the documentation. Would you be able to test out the built-in chunking to make sure it works for your use case?