Build ioda and ioda-converters with Python 3.13

repos/spack_stack/spack_repo/spack_stack/packages/ioda/package.py

@@ -71,7 +71,8 @@ class Ioda(CMakePackage):
     depends_on("oops@1.10.0.20260331", when="@2.9.0.20260326")
     depends_on("oops@1.10.0.20250827", when="@2.9.0.20250826")
     depends_on("python")
-    depends_on("python@3.9:3.11", when="@2.9:")
+    # https://github.com/JCSDA/spack-stack/issues/2116
+    depends_on("python@3.9:3.13", when="@2.9:")
     depends_on("py-pybind11")
     depends_on("py-pycodestyle", type=("build", "test"))
     depends_on("py-netcdf4", type=("build", "test"))

repos/spack_stack/spack_repo/spack_stack/packages/ioda_converters/ioda-converters-numpy2-api-and-precision.patch

@@ -0,0 +1,107 @@
+diff --git a/src/compo/airnow2ioda_nc.py b/src/compo/airnow2ioda_nc.py
+index acf0fbf..b784b6d 100755
+--- a/src/compo/airnow2ioda_nc.py
++++ b/src/compo/airnow2ioda_nc.py
+@@ -97,7 +97,8 @@ def read_monitor_file(sitefile, is_epa):
+
+
+ def filter_bad_values(df):
+-    df.loc[(df.obs > 3000) | (df.obs < 0), 'obs'] = np.NaN
++    # np.NaN was removed in numpy 2.0; use np.nan
++    df.loc[(df.obs > 3000) | (df.obs < 0), 'obs'] = np.nan
+     return df
+
+
+diff --git a/src/gsi_ncdiag/gsi_ncdiag.py b/src/gsi_ncdiag/gsi_ncdiag.py
+index 18d1791..f1b0b34 100755
+--- a/src/gsi_ncdiag/gsi_ncdiag.py
++++ b/src/gsi_ncdiag/gsi_ncdiag.py
+@@ -1190,7 +1190,9 @@ def grabobsidx(obsdata, platform, var):
+         codes = uv_bufrtypes[platform]
+     else:
+         codes = conv_bufrtypes[platform]
+-    idx = np.logical_and(np.in1d(code, codes), idx2)
++    # np.in1d was removed in numpy 2.0; np.isin is the drop-in replacement
++    # (identical for the 1-D inputs used here).
++    idx = np.logical_and(np.isin(code, codes), idx2)
+
+     return idx
+
+diff --git a/src/marine/ndbc_hfradar2ioda.py b/src/marine/ndbc_hfradar2ioda.py
+index 3345b5d..e24e6f2 100755
+--- a/src/marine/ndbc_hfradar2ioda.py
++++ b/src/marine/ndbc_hfradar2ioda.py
+@@ -65,7 +65,10 @@ class Observation(object):
+                 qcKey = vName[j], iconv.OqcName()
+                 if vals_u[i] != '--':
+                     count += 1
+-                    obs_date = int(time[i])
++                    # numpy 2 no longer converts a size-1, non-0-d array with
++                    # int(); .item() is the explicit form and still raises if
++                    # the slice holds more than one element.
++                    obs_date = int(time[i].item())
+                     locKey = lats[i], lons[i], obs_date
+                     if j == 0:
+                         self.data[locKey][valKey] = vals_u[i]
+diff --git a/src/pyiodaconv/meteo_utils.py b/src/pyiodaconv/meteo_utils.py
+index 750a07b..2ee02ae 100644
+--- a/src/pyiodaconv/meteo_utils.py
++++ b/src/pyiodaconv/meteo_utils.py
+@@ -61,6 +61,14 @@ class meteo_utils(object):
+         temp_K  - temperature (k)
+         '''
+
++        # Evaluate in double precision. Under numpy < 2, combining a numpy
++        # float32 scalar with a Python float promoted the result to float64;
++        # NEP 50 (numpy >= 2) keeps it float32, which loses several digits in
++        # the polynomial fit below. Coerce explicitly so the result no longer
++        # depends on the dtype the caller happens to pass in.
++        pres_Pa = float(pres_Pa)
++        temp_K = float(temp_K)
++
+         es = self.e_sub_s(temp_K)
+
+         # Even at P=1050hPa and T=55C, sat. vap. pres only contributes to ~15% of total pressure.
+@@ -81,6 +89,13 @@ class meteo_utils(object):
+         polynomial fit of Goff-Gratch (1946) formulation. (Walko, 1991)
+         '''
+
++        # Evaluate in double precision. Under numpy < 2, combining a numpy
++        # float32 scalar with a Python float promoted the result to float64;
++        # NEP 50 (numpy >= 2) keeps it float32, which loses several digits in
++        # the polynomial fit below. Coerce explicitly so the result no longer
++        # depends on the dtype the caller happens to pass in.
++        temp_K = float(temp_K)
++
+         c = [610.5851, 44.40316, 1.430341, 0.2641412e-1, 0.2995057e-3, 0.2031998e-5, 0.6936113e-8, 0.2564861e-11, -0.3704404e-13]
+         x = max(-80., temp_K-self.C_2_K)
+         es = c[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*(c[4]+x*(c[5]+x*(c[6]+x*(c[7]+x*c[8])))))))
+@@ -113,6 +128,14 @@ class meteo_utils(object):
+         temp_K  - temperature (k)
+         '''
+
++        # Evaluate in double precision. Under numpy < 2, combining a numpy
++        # float32 scalar with a Python float promoted the result to float64;
++        # NEP 50 (numpy >= 2) keeps it float32, which loses several digits in
++        # the polynomial fit below. Coerce explicitly so the result no longer
++        # depends on the dtype the caller happens to pass in.
++        pres_Pa = float(pres_Pa)
++        temp_K = float(temp_K)
++
+         esi = self.e_sub_i(temp_K)
+
+         # Even at P=1050hPa and T=55C, sat. vap. pres only contributes to ~15% of total pressure.
+@@ -133,6 +156,13 @@ class meteo_utils(object):
+         polynomial fit of Goff-Gratch (1946) formulation. (Walko, 1991)
+         '''
+
++        # Evaluate in double precision. Under numpy < 2, combining a numpy
++        # float32 scalar with a Python float promoted the result to float64;
++        # NEP 50 (numpy >= 2) keeps it float32, which loses several digits in
++        # the polynomial fit below. Coerce explicitly so the result no longer
++        # depends on the dtype the caller happens to pass in.
++        temp_K = float(temp_K)
++
+         c = [.609868993E03, .499320233E02, .184672631E01, .402737184E-1, .565392987E-3, .521693933E-5, .307839583E-7, .105785160E-9, .161444444E-12]
+         x = max(-80., temp_K-self.C_2_K)
+         esi = c[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*(c[4]+x*(c[5]+x*(c[6]+x*(c[7]+x*c[8])))))))

repos/spack_stack/spack_repo/spack_stack/packages/ioda_converters/ioda-converters-numpy2-ma-scalar-dtype.patch

@@ -0,0 +1,235 @@
+diff --git a/src/compo/mopitt_co_nc2ioda.py b/src/compo/mopitt_co_nc2ioda.py
+index 1f60e97..ca86406 100755
+--- a/src/compo/mopitt_co_nc2ioda.py
++++ b/src/compo/mopitt_co_nc2ioda.py
+@@ -124,10 +124,19 @@ class mopitt(object):
+
+             # convert all concentrations and column to correct units to avoid single precision issues
+             u_conv = avogadro / scm2sm
+-            xa_gd = xa_gd * vmr2col / u_conv
+-            xa_tc = xa_tc / u_conv
+-            xr_tc = xr_tc / u_conv
+-            er_tc = er_tc / u_conv
++            # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++            # so an operation between a masked array and a Python scalar widens the
++            # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++            # See https://github.com/numpy/numpy/issues/27029
++            # vmr2col and u_conv are not exactly representable in float32, so it is
++            # not enough to cast the result: numpy < 2 also narrowed the Python
++            # scalar to the array's dtype before operating (value-based casting),
++            # and doing the arithmetic in float64 changes the last bits. Narrow the
++            # scalars explicitly to keep the computation in the array's precision.
++            xa_gd = xa_gd * xa_gd.dtype.type(vmr2col) / xa_gd.dtype.type(u_conv)
++            xa_tc = xa_tc / xa_tc.dtype.type(u_conv)
++            xr_tc = xr_tc / xr_tc.dtype.type(u_conv)
++            er_tc = er_tc / er_tc.dtype.type(u_conv)
+
+             # mopitt number of levels is dependent on surface pressure, for data points with sp<900hPa
+             # nlevs<10. IODA and UFO cannot handle variable nlayers_kernel for a given instrument
+@@ -169,7 +178,8 @@ class mopitt(object):
+
+                 self.outdata[('aprioriTerm', 'RetrievalAncillaryData')] = ap_tc[flg]
+                 self.outdata[('averagingKernel', 'RetrievalAncillaryData')] = ak_tc_dimless[flg]
+-                self.outdata[('pressureVertice', 'RetrievalAncillaryData')] = hPa2Pa * pr_gd[flg]
++                self.outdata[('pressureVertice', 'RetrievalAncillaryData')] = \
++                    (hPa2Pa * pr_gd[flg]).astype(pr_gd.dtype)
+
+                 self.outdata[self.varDict[iodavar]['valKey']] = xr_tc[flg]
+                 self.outdata[self.varDict[iodavar]['errKey']] = er_tc[flg]
+@@ -188,7 +198,8 @@ class mopitt(object):
+                 self.outdata[('averagingKernel', 'RetrievalAncillaryData')] = np.concatenate((
+                     self.outdata[('averagingKernel', 'RetrievalAncillaryData')], ak_tc_dimless[flg]))
+                 self.outdata[('pressureVertice', 'RetrievalAncillaryData')] = np.concatenate((
+-                    self.outdata[('pressureVertice', 'RetrievalAncillaryData')], hPa2Pa * pr_gd[flg]))
++                    self.outdata[('pressureVertice', 'RetrievalAncillaryData')],
++                    (hPa2Pa * pr_gd[flg]).astype(pr_gd.dtype)))
+
+                 self.outdata[self.varDict[iodavar]['valKey']] = np.concatenate(
+                     (self.outdata[self.varDict[iodavar]['valKey']], xr_tc[flg]))
+diff --git a/src/compo/tropomi_no2_co_nc2ioda.py b/src/compo/tropomi_no2_co_nc2ioda.py
+index 0d3f26a..a2b94ea 100755
+--- a/src/compo/tropomi_no2_co_nc2ioda.py
++++ b/src/compo/tropomi_no2_co_nc2ioda.py
+@@ -145,7 +145,12 @@ class tropomi(object):
+                     groups['DETAILED_RESULTS'].variables['surface_albedo_2325'][:].ravel()
+                 albedo2 = ncd.groups['PRODUCT'].groups['SUPPORT_DATA'].\
+                     groups['DETAILED_RESULTS'].variables['surface_albedo_2335'][:].ravel()
+-                albedo = 0.5 * (albedo1 + albedo2)
++                # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++                # so an operation between a masked array and a Python scalar widens the
++                # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++                # Cast back explicitly so the intended netCDF type is written.
++                # See https://github.com/numpy/numpy/issues/27029
++                albedo = (0.5 * (albedo1 + albedo2)).astype(albedo1.dtype)
+
+             # get angles
+             sza = ncd.groups['PRODUCT'].groups['SUPPORT_DATA'].\
+diff --git a/src/conventional/saber2ioda.py b/src/conventional/saber2ioda.py
+index ccb4632..212f711 100755
+--- a/src/conventional/saber2ioda.py
++++ b/src/conventional/saber2ioda.py
+@@ -195,10 +195,16 @@ def get_data_from_file(obs_file_handle):
+     obs_data[('pressure', META_DATA_NAME)] *= 100.0
+
+     # Handle longitudes to be within [-180, 180)
++    # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++    # so an operation between a masked array and a Python scalar widens the
++    # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++    # Cast back explicitly so the intended netCDF type is written.
++    # See https://github.com/numpy/numpy/issues/27029
++    longitudes = obs_data[('longitude', META_DATA_NAME)]
+     obs_data[('longitude', META_DATA_NAME)] = numpy.ma.where(
+-        obs_data[('longitude', META_DATA_NAME)] > 180,
+-        obs_data[('longitude', META_DATA_NAME)] - 360,
+-        obs_data[('longitude', META_DATA_NAME)]
++        longitudes > 180,
++        (longitudes - 360).astype(longitudes.dtype),
++        longitudes
+     )
+
+     # Handle time conversion
+diff --git a/src/land/imsfv3_scf2ioda.py b/src/land/imsfv3_scf2ioda.py
+index ea977a6..d3da70b 100755
+--- a/src/land/imsfv3_scf2ioda.py
++++ b/src/land/imsfv3_scf2ioda.py
+@@ -99,10 +99,15 @@ class imsFV3(object):
+         sncv = sncv.astype('float32')
+         sndv = sndv.astype('float32')
+
+-        qcflg = 0*sncv.astype('int32')
+-        qdflg = 0*sndv.astype('int32')
+-        errsc = 0.0*sncv
+-        errsd = 0.0*sndv
++        # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++        # so an operation between a masked array and a Python scalar widens the
++        # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++        # Cast back explicitly so the intended netCDF type is written.
++        # See https://github.com/numpy/numpy/issues/27029
++        qcflg = (0*sncv.astype('int32')).astype('int32')
++        qdflg = (0*sndv.astype('int32')).astype('int32')
++        errsc = (0.0*sncv).astype(sncv.dtype)
++        errsd = (0.0*sndv).astype(sndv.dtype)
+         errsd[:] = 80.
+
+         times = get_observation_time(self.filename, sncv, ncd)
+diff --git a/src/marine/glider2ioda.py b/src/marine/glider2ioda.py
+index d7486e0..a10abc3 100755
+--- a/src/marine/glider2ioda.py
++++ b/src/marine/glider2ioda.py
+@@ -63,8 +63,13 @@ class Profile(object):
+         with np.errstate(invalid='ignore'):
+             salinity = np.float32(ncd.variables['salinity'][:])
+         errs = np.float32(np.matlib.repmat(0.2, len(lons), 1))
+-        Tqcs = ncd.variables['temperature_qc'][:]-1
+-        Sqcs = ncd.variables['salinity_qc'][:]-1
++        # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++        # so an operation between a masked array and a Python scalar widens the
++        # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++        # Cast back explicitly so the intended netCDF type is written.
++        # See https://github.com/numpy/numpy/issues/27029
++        Tqcs = (ncd.variables['temperature_qc'][:]-1).astype('int32')
++        Sqcs = (ncd.variables['salinity_qc'][:]-1).astype('int32')
+         errs = np.squeeze(errs)
+         ncd.close()
+
+diff --git a/src/marine/pace_oc2ioda.py b/src/marine/pace_oc2ioda.py
+index e6689b1..4342f2a 100755
+--- a/src/marine/pace_oc2ioda.py
++++ b/src/marine/pace_oc2ioda.py
+@@ -131,7 +131,13 @@ def read_input(input_args):
+
+     obs_data[output_var_names[0], global_config['oval_name']] = data_in['chlor_a']
+     # There is not any obs error in the dataset. we need to come up with a reasonable obs error later
+-    obs_data[output_var_names[0], global_config['oerr_name']] = data_in['chlor_a']*0.0
++    # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++    # so an operation between a masked array and a Python scalar widens the
++    # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++    # Cast back explicitly so the intended netCDF type is written.
++    # See https://github.com/numpy/numpy/issues/27029
++    obs_data[output_var_names[0], global_config['oerr_name']] = \
++        (data_in['chlor_a']*0.0).astype(data_in['chlor_a'].dtype)
+     obs_data[output_var_names[0], global_config['opqc_name']] = data_in['l2_flags']
+
+     return (obs_data, basetime, GlobalAttrs)
+diff --git a/src/marine/viirs_modis_l2_oc2ioda.py b/src/marine/viirs_modis_l2_oc2ioda.py
+index 62dfd61..b58d129 100755
+--- a/src/marine/viirs_modis_l2_oc2ioda.py
++++ b/src/marine/viirs_modis_l2_oc2ioda.py
+@@ -135,11 +135,16 @@ def read_input(input_args):
+
+     if global_config['output_poc']:
+         obs_data[output_var_names[0], obsValName] = data_in['poc']
+-        obs_data[output_var_names[0], obsErrName] = data_in['poc']*0.0
++        # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++        # so an operation between a masked array and a Python scalar widens the
++        # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++        # Cast back explicitly so the intended netCDF type is written.
++        # See https://github.com/numpy/numpy/issues/27029
++        obs_data[output_var_names[0], obsErrName] = (data_in['poc']*0.0).astype(data_in['poc'].dtype)
+         obs_data[output_var_names[0], qcName] = data_in['l2_flags']
+     if global_config['output_chl']:
+         obs_data[output_var_names[1], obsValName] = data_in['chlor_a']
+-        obs_data[output_var_names[1], obsErrName] = data_in['chlor_a']*0.0
++        obs_data[output_var_names[1], obsErrName] = (data_in['chlor_a']*0.0).astype(data_in['chlor_a'].dtype)
+         obs_data[output_var_names[1], qcName] = data_in['l2_flags']
+
+     return (obs_data, GlobalAttrs, time_units)
+diff --git a/src/ncepbufr/bufr_ncep_prepbufr_adpupa.py b/src/ncepbufr/bufr_ncep_prepbufr_adpupa.py
+index 19b656e..121dbba 100755
+--- a/src/ncepbufr/bufr_ncep_prepbufr_adpupa.py
++++ b/src/ncepbufr/bufr_ncep_prepbufr_adpupa.py
+@@ -74,7 +74,13 @@ def test_bufr_to_ioda(DATA_PATH, OUTPUT_PATH, date):
+     # to get updated, and then before writing into the output ioda file, the masked array
+     # function filled() needs to be called which will convert the values marked invalid
+     # to the fill value.
+-    hrdr = (r.get('timeOffset') * 3600).astype(np.int64)
++    # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++    # so an operation between a masked array and a Python scalar widens the
++    # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++    # Cast back explicitly so the intended netCDF type is written.
++    # See https://github.com/numpy/numpy/issues/27029
++    timeOffset = r.get('timeOffset')
++    hrdr = (timeOffset * 3600).astype(timeOffset.dtype).astype(np.int64)
+     np.ma.set_fill_value(hrdr, long_missing_value)
+     print("cycleTimeSinceEpoch")
+     cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(date, '%Y%m%d%H%M')))
+diff --git a/src/ncepbufr/prepbufr_adpsfc_api.py b/src/ncepbufr/prepbufr_adpsfc_api.py
+index 7dd4370..77c6e7a 100755
+--- a/src/ncepbufr/prepbufr_adpsfc_api.py
++++ b/src/ncepbufr/prepbufr_adpsfc_api.py
+@@ -57,7 +57,14 @@ def test_bufr_to_ioda(DATA_PATH, OUTPUT_PATH, date):
+    # function filled() needs to be called which will convert the values marked invalid
+    # to the fill value.
+    print("Get time")
+-   dhr = (r.get('obsTimeMinusCycleTime') * 3600).astype(np.int64)  # Needs to be converted to seconds since Epoch time from [-3,3]
++   # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++   # so an operation between a masked array and a Python scalar widens the
++   # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++   # Cast back explicitly so the intended netCDF type is written.
++   # See https://github.com/numpy/numpy/issues/27029
++   obsTimeMinusCycleTime = r.get('obsTimeMinusCycleTime')
++   # Needs to be converted to seconds since Epoch time from [-3,3]
++   dhr = (obsTimeMinusCycleTime * 3600).astype(obsTimeMinusCycleTime.dtype).astype(np.int64)
+    np.ma.set_fill_value(dhr, long_missing_value)
+    print("cycleTimeSinceEpoch") #For now, file time is put in manually 
+    cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(date, '%Y%m%d%H%M')))
+diff --git a/src/ncepbufr/prepbufr_sfcshp_api.py b/src/ncepbufr/prepbufr_sfcshp_api.py
+index 4dc4d99..2a3f811 100755
+--- a/src/ncepbufr/prepbufr_sfcshp_api.py
++++ b/src/ncepbufr/prepbufr_sfcshp_api.py
+@@ -60,7 +60,14 @@ def test_bufr_to_ioda(DATA_PATH, OUTPUT_PATH, date):
+    # function filled() needs to be called which will convert the values marked invalid
+    # to the fill value.
+    print("Get time")
+-   dhr = (r.get('obsTimeMinusCycleTime') * 3600).astype(np.int64)  # Needs to be converted to seconds since Epoch time from [-3,3]
++   # numpy >= 2 does not apply NEP 50 weak-scalar promotion inside numpy.ma,
++   # so an operation between a masked array and a Python scalar widens the
++   # result to float64/int64 (it kept the array's dtype under numpy 1.x).
++   # Cast back explicitly so the intended netCDF type is written.
++   # See https://github.com/numpy/numpy/issues/27029
++   obsTimeMinusCycleTime = r.get('obsTimeMinusCycleTime')
++   # Needs to be converted to seconds since Epoch time from [-3,3]
++   dhr = (obsTimeMinusCycleTime * 3600).astype(obsTimeMinusCycleTime.dtype).astype(np.int64)
+    np.ma.set_fill_value(dhr, long_missing_value)
+    print("cycleTimeSinceEpoch") #For now, file time is put in manually 
+    cycleTimeSinceEpoch = np.int64(calendar.timegm(time.strptime(date, '%Y%m%d%H%M')))