ioda and ioda-converters with Python 3.13 - ioda over constrained; ioda-converters: numpy 2 failures

[climbfuji]

Summary

Disclaimer: Prepared with the help of Anthropic Claude. This would have taken me a week at least otherwise.

Building ioda@2.9.0.20260326 (JCSDA/ioda 9e0eb39) and ioda-converters@0.0.1.20250830(JCSDA-internal/ioda-convertersa91f432) in a spack-stack environment on Python 3.13 pulls py-numpy@2andpy-pandas@3`.

ioda builds without problems and all ctests pass, including the pyioda tests. The existing constraint in the ioda recipe (`depends_on("python@3.9:3.11") can therefore be relaxed.

Different story for ioda-converters:

pandas 3 breaks first: its dedicated string dtype is an extension dtype, not an np.dtype, so
passing it to numpy raises TypeError: Cannot interpret '<StringDtype(storage='python', na_value=nan)>' as a data type. We therefore pin py-pandas@2.3.3. Everything below is with that pin in
place — 24 of 375 ctests fail, 18 of them from numpy 2 API and type-promotion changes.

Environment: gcc@14.2.1, python@3.13, py-numpy@2.4.6, py-pandas@2.3.3.

spack-stack will carry two patches in the ioda-converters recipe, applied when='^py-numpy@2:' (PR). @eap @srherbener — please check if these are needed in develop. They are applied here specifically to make commit a91f432 pass the ctests in Spack. One alternative is to recreate the reference solutions when updating to py-numpy@2, with the potential drawback that the bug in numpy will get fixed and the reference solution then depends on the exact version of numpy. I am sure there are other alternatives, too.

1. numpy.ma does not implement NEP 50 — 11 tests

numpy >= 2 applies NEP 50 weak-scalar promotion to ndarray but not inside numpy.ma
(numpy/numpy#27029), so an operation between a masked array and a Python scalar widens the result:

                       ndarray              masked array
numpy 1.26     0.0*f32 -> float32     0.0*f32 -> float32,  0*i32 -> int32
numpy 2.4      0.0*f32 -> float32     0.0*f32 -> float64,  0*i32 -> int64

The converters read inputs with netCDF4.Dataset(...)[:], which returns masked arrays, and
IodaWriter.WriteObsVars passes Vvals.dtype straight to create_var — so the widened dtype is
written out and nccmp reports DOUBLE <> FLOAT / INT64 <> INT.

Test	Site (at a91f432)	Expression
178, 181 imsfv3_scf, imsfv3grid_scf	src/land/imsfv3_scf2ioda.py:102-105	0*sncv.astype('int32'), 0.0*sncv
32 glider	src/marine/glider2ioda.py:66-67	ncd.variables['temperature_qc'][:]-1
43 pace_oc_l2	src/marine/pace_oc2ioda.py:134	data_in['chlor_a']*0.0
57, 58 viirs_jpss1_oc_l2, modis_aqua_oc_l2	src/marine/viirs_modis_l2_oc2ioda.py:138,142	data_in['poc']*0.0
84 mopitt_co	src/compo/mopitt_co_nc2ioda.py:127-130,172,191	xa_gd*vmr2col/u_conv, hPa2Pa*pr_gd[flg]
91 tropomi_co_total	src/compo/tropomi_no2_co_nc2ioda.py:148	0.5*(albedo1+albedo2)
174 saber_timed_v2_0	src/conventional/saber2ioda.py:200	obs_data[('longitude', ...)] - 360
108, 126 prepbufr_ncep_api_sfcshp2ioda, bufr_ncep_prepbufr_adpupa_api	src/ncepbufr/bufr_ncep_prepbufr_adpupa.py:77, prepbufr_sfcshp_api.py:63	(r.get('timeOffset')*3600).astype(np.int64)

Two do not present as type diffs:

• 174 raises ValueError: No missing value defined for dtype float64 — map_missing_values()
  matches against numpy.float32, and longitude - 360 is now float64.
• 108, 126 report dateTime off by one second: the multiply happens in float64, so the
  following astype(np.int64) truncates differently. src/ncepbufr/prepbufr_adpsfc_api.py:60 has
  the identical expression; its test passes today but should be fixed too.

Casting the result is not always enough. Legacy value-based casting also narrowed the Python
scalar to the array's dtype before operating, so xa_gd * vmr2col multiplied by
float32(2.12e13), not by the double. Where the constant is not exactly representable in float32
(vmr2col = 2.12e13, u_conv = 6.02214076e19 in mopitt_co_nc2ioda.py), doing the arithmetic in
float64 and rounding at the end still leaves aprioriTerm ~1.2e-7 off, above
IODA_CONV_COMP_TOL. The scalar has to be narrowed:

xa_gd = xa_gd * xa_gd.dtype.type(vmr2col) / xa_gd.dtype.type(u_conv)

Constants that are exact in float32 (0.0, 0.5, -1, 3600, 1e2) are unaffected.

2. numpy 2 API removals — 4 tests

• 66, 67 gsidiag_conv_uv, gsidiag_conv — src/gsi_ncdiag/gsi_ncdiag.py:1193, np.in1d -> np.isin
• 89, 90 airnow, airnow_epalist — src/compo/airnow2ioda_nc.py:100, np.NaN -> np.nan

A sweep of the tree found no other removed aliases.

3. int() on a size-1, ndim>0 array — 1 test

33 hfradar, src/marine/ndbc_hfradar2ioda.py:68 — TypeError: only 0-dimensional arrays can be converted to Python scalars. Deprecated in numpy 1.25, an error in numpy 2.
int(time[i].item()) is the exact equivalent and still raises on a multi-element slice.

4. NEP 50 for numpy scalars: meteo_utils runs in single precision — 2 tests

150 generic_bufr_raob and 153 ship report specificHumidity differences of ~1e-5 relative.
This is NEP 50 proper, opposite in direction to section 1: under numpy 1 np.float32(x) - 273.15
promoted to float64, so the Goff-Gratch polynomial in src/pyiodaconv/meteo_utils.py was evaluated
in double precision; under numpy 2 it stays float32 and loses digits to cancellation.

The reference files were generated in double precision, so coerce explicitly rather than depend on
the caller's dtype: pres_Pa = float(pres_Pa) / temp_K = float(temp_K) at the top of r_sub_s
(:55), e_sub_s (:77), r_sub_i (:107), e_sub_i (:129). specific_humidity is the only entry
point called anywhere in the repo.

Not numpy

• 26, 27, 28, 148, 149 — ModuleNotFoundError: No module named 'scipy' (godae_*2ioda.py,
  apply_BG.py). Missing py-scipy dependency, fixed in the recipe.
• 2 iodaconv_gsi_ncdiag_coding_norms — pre-existing pycodestyle E502, unrelated.
• src/compo/airnow2ioda_nc.py:132 uses pd.to_timedelta(..., unit='H'), deprecated in pandas 2
  and removed in 3. Only a warning under the pin, but it blocks unpinning pandas later.

Verification

Checked against numpy 1.26.4 rather than by eye:

• meteo_utils patched on numpy 2.4.4 returns bit-identical doubles to the unpatched module on
  numpy 1.26.4 over the Td/P range these converters screen for.
• mopitt aprioriTerm, full chain replayed on 1500 synthetic profiles: 0/1500 differences from
  numpy 1.26.4, versus 281/1500 unpatched.
• imsfv3_scf2ioda run end to end against a stub pyioda: all six output variables recover
  float32/int32.

With both patches plus py-scipy, 23 of the 24 failures clear; only the pre-existing test 2 remains (already excluded in the Spack recipe).