Consider marray as a uniform masked-array representation across array backends

[mwcraig]

Follow-up to a question raised on #979 (https://github.com/astropy/ccdproc/pull/979/changes#r3837117088): should ccdproc use marray, the array-API-compatible masked-array package, as a uniform "data + mask" representation across backends?

Current state

Masked data is not part of the array API standard, so ccdproc carries masks around by hand, in several different shapes:

• numpy.ma.MaskedArray inputs, detected with hasattr(x, "mask") / isinstance checks and split into .data / .mask (cosmicray_median, core.py; block_reduce/block_average comment at core.py:~2283).
• CCDData.mask, which astropy keeps as a numpy array even when .data lives in another namespace, so every consumer coerces it with xp.asarray(ccd.mask, dtype=xp.bool, device=...) (Combiner.__init__, combine() in three places, flat_correct, transform_image, wcs_project, rebin, cosmicray_lacosmic, cosmicray_median).
• Internal (data, mask) pairs in Combiner (_data_arr, _data_arr_mask) with the mask applied via xp.where / NaN substitution before reductions.

This works but the attribute juggling is repeated, easy to get subtly wrong (#932 was exactly a dropped mask; #965 was the mask stacking), and mask device/namespace placement has to be re-established at every site.

What marray would and would not give

• Would: one object carrying data + mask in any array-API namespace, with mask-aware reductions (mean, sum, std, ...) and elementwise ops, so Combiner could drop much of its NaN/xp.where bookkeeping, and functions could accept/return one masked type instead of branching on MaskedArray vs CCDData vs bare array.
• Would not: help with anything that leaves the array API. scipy.ndimage.median_filter (cosmicray_median), astroscrappy, reproject and scipy.ndimage in transform_image all see only the underlying data, so the fill-then-filter pattern around them stays regardless (see the Honor the input mask in cosmicray_median #979 discussion).

Open questions

1. Is a new required dependency acceptable, or would it be optional with a fallback (which would mean keeping the manual paths anyway)?
2. Maturity / maintenance status of marray, and whether its masked reductions match the semantics ccdproc currently documents (e.g. Combiner weighted averages exclude masked and clipped pixels, Exclude masked weights from average combinations #952).
3. Boundary with CCDData: astropy's NDData.mask setter forces numpy, so a marray would have to be unpacked at the CCDData boundary in both directions unless astropy grows array-API-aware masks.
4. Performance on dask / jax / CuPy compared to the existing xp.where approach.

Suggested first step if this is pursued: a prototype that replaces Combiner._data_arr + _data_arr_mask with a single marray and measures the change in combine() on dask and numpy, before touching the public API.

This is a design discussion, not a bug; no change is proposed for #979 or #932.