From b43c0820d08a78be1e10cdbbaed797be0c10a793 Mon Sep 17 00:00:00 2001 From: mansurjisan Date: Sat, 24 Oct 2026 08:56:33 -0400 Subject: [PATCH 1/2] fix(pacific): restore ATL/SECOFS/AK behaviour regressed by the Pacific ops features ADT datum offset applied per day again (ATL -0.45 via adt_datum_offset), ops ADT ROI kept for -180/180 systems, 23 RTOFS ops slots (obc_min_slots), OPS_OBC_DT NameError, 3-hourly SECOFS/AK nudging (nudging_dt default 10800), rerun archive extras opt-in (obc_archive_extended), surface restore off and refusing to run unless explicitly forced, Pacific forecast_hours 48. Adds guard tests. --- nos_utils/config.py | 30 +++++++++---- nos_utils/forcing/adt.py | 27 +++++++++--- nos_utils/forcing/nudging.py | 6 +-- nos_utils/forcing/rtofs.py | 10 ++--- nos_utils/forcing/surface_restore.py | 15 +++++++ nos_utils/orchestrator.py | 50 ++++++++++++++-------- tests/test_stofs_config.py | 64 ++++++++++++++++++++++++++-- 7 files changed, 159 insertions(+), 43 deletions(-) diff --git a/nos_utils/config.py b/nos_utils/config.py index e9ef202..70b265c 100644 --- a/nos_utils/config.py +++ b/nos_utils/config.py @@ -238,6 +238,8 @@ class ForcingConfig: # ADT satellite SSH blending (STOFS-3D-ATL uses CMEMS ADT to correct boundary SSH) adt_enabled: bool = False adt_datum_offset: float = -0.45 + # Rerun archive/tar also carries temnu/salnu/surface_restore with the Pacific names; off keeps main's set. MJ (10/24/26) + obc_archive_extended: bool = False # Ops ESMF map for the ADT regrid (yaml ocean.adt.weight_file); a bare name resolves against # FIXofs. Unset or missing falls back to the bilinear regrid. MJ (10/02/26) adt_weight_file: Optional[Path] = None @@ -245,7 +247,7 @@ class ForcingConfig: nudging_enabled: bool = False nudging_timescale_seconds: float = 86400.0 # Temporal resolution of the nudging files (seconds). - nudging_dt: float = 21600.0 + nudging_dt: float = 10800.0 # Apply 2-step temporal moving average to nudging fields. nudging_time_smoothing: bool = False # Surface restore enabled and smoothing @@ -436,6 +438,8 @@ class ForcingConfig: # OFSes have shorter OBC windows (typically ~10 records for 54h runs) # and would be wrongly flagged if this defaulted to 21. obc_min_timesteps: int = 0 + # Minimum RTOFS slots in the ops OBC timeline (ops: 23 6-hourly records, n012 to f120). MJ (10/24/26) + obc_min_slots: int = 23 # --- Tidal settings --- # Tidal constituents to use @@ -707,7 +711,7 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig": lon_min=93.1098, lon_max=289.6750, lat_min=-33.0388, lat_max=66.3254, pdy=pdy, cyc=cyc, - nowcast_hours=24, forecast_hours=108, + nowcast_hours=24, forecast_hours=48, gfs_resolution="0p25", met_num=2, n_levels=84, model_dt=90.0, @@ -724,6 +728,8 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig": nudging_timescale_seconds=86400.0, nwm_product="medium_range_mem1", obc_min_timesteps=4, + obc_min_slots=4, + obc_archive_extended=True, obc_dt=86400.0, nudging_dt=86400.0, obc_time_smoothing=True, @@ -732,8 +738,8 @@ def for_stofs_3d_pac(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig": obc_roi_2d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2}, obc_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2}, nudge_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2}, - surface_restore_enabled=True, - surface_restore_time_smoothing=True, + # Surface restore is untested (reads T and S from one file, restore nodes unverified): off. MJ (10/24/26) + surface_restore_enabled=False, ) defaults.update(overrides) return cls(**defaults) @@ -750,7 +756,7 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig lon_min=93.1098, lon_max=289.6750, lat_min=-33.0388, lat_max=66.3254, pdy=pdy, cyc=cyc, - nowcast_hours=24, forecast_hours=108, + nowcast_hours=24, forecast_hours=48, gfs_resolution="0p25", met_num=2, nws=4, n_levels=84, model_dt=90.0, @@ -764,6 +770,8 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig nudging_timescale_seconds=86400.0, nwm_product="medium_range_mem1", obc_min_timesteps=4, + obc_min_slots=4, + obc_archive_extended=True, obc_dt=86400.0, nudging_dt=86400.0, obc_time_smoothing=True, @@ -772,8 +780,8 @@ def for_stofs_3d_pac_ufs(cls, pdy: str, cyc: int, **overrides) -> "ForcingConfig obc_roi_2d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2}, obc_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2}, nudge_roi_3d={"x1": 222, "x2": 2703, "y1": 1051, "y2": 2641, "x_stride": 2, "y_stride": 2}, - surface_restore_enabled=True, - surface_restore_time_smoothing=True, + # Surface restore is untested (reads T and S from one file, restore nodes unverified): off. MJ (10/24/26) + surface_restore_enabled=False, # DATM grid — STOFS3D_PAC preset, verified from weight_gfs.nc # nx = round((290.0-92.5)/0.025)+1 = 7901; ny = round((67.0-(-33.5))/0.025)+1 = 4021 datm_lon_min=92.5, datm_lon_max=290.0, @@ -995,7 +1003,7 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig": if _adt_wt is None and str(_sys.get("name", "")).startswith("stofs_3d_atl"): _adt_wt = "stofs_3d_atl_ufs.adt_weight.nc" _adt_wt = Path(_adt_wt) if _adt_wt else None - _adt_off = float(adt.get("datum_offset", 0.0)) if isinstance(adt, dict) else 0.0 + _adt_off = float(adt.get("datum_offset", -0.45)) if isinstance(adt, dict) else -0.45 # Ops river inputs (static FIX source_sink.in / msource.th / vsink.th, ops source # order): STOFS-3D-ATL by name (the coupled executable now has PREC_EVAP), yaml @@ -1077,7 +1085,7 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig": kwargs["rtofs_daily_3z"] = _strict_bool(ocean.get("rtofs_daily_3z"), "ocean.rtofs_daily_3z", False) kwargs["obc_dt"] = float(obc.get("dt", 21600.0)) kwargs["obc_time_smoothing"] = _strict_bool(obc.get("time_smoothing"), "obc.time_smoothing", False) - kwargs["nudging_dt"] = float(nudge.get("dt", kwargs["obc_dt"])) + kwargs["nudging_dt"] = float(nudge.get("dt", 10800.0)) kwargs["nudging_time_smoothing"] = _strict_bool(nudge.get("time_smoothing"), "nudging.time_smoothing", False) surf_res = ocean.get("surface_restore", {}) if isinstance(ocean, dict) else {} @@ -1182,6 +1190,10 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig": kwargs["obc_min_timesteps"] = int(yaml_min_ts) elif is_stofs_like: kwargs["obc_min_timesteps"] = 21 + if obc.get("min_slots") is not None: + kwargs["obc_min_slots"] = int(obc["min_slots"]) + kwargs["obc_archive_extended"] = _strict_bool( + obc.get("archive_extended"), "obc.archive_extended", False) # Declarative prep extras override. When a root-level `prep.extras` # block is present it is authoritative and supersedes the diff --git a/nos_utils/forcing/adt.py b/nos_utils/forcing/adt.py index d422365..4016dc2 100644 --- a/nos_utils/forcing/adt.py +++ b/nos_utils/forcing/adt.py @@ -50,6 +50,10 @@ # Ops ADT subset ROI, used only when the SSH_1 destination grid cannot be read. # MJ (10/01/26) +ADT_LON_MIN = -62.5 +ADT_LON_MAX = -51.5 +ADT_LAT_MIN = 7.0 +ADT_LAT_MAX = 54.0 ADT_PAD = 1.0 # Ops holds ADT-missing cells as the CMEMS int fill carried through float32. MJ (10/02/26) @@ -155,9 +159,21 @@ def _dest_bounds(self, ssh_path: Path): float(np.nanmin(lat)) - ADT_PAD, float(np.nanmax(lat)) + ADT_PAD) except Exception as e: log.warning(f"Cannot read SSH_1 grid extent ({e}) — using config domain bounds") + if self._pm180(): + return ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX lon_min, lon_max, lat_min, lat_max = self.config.domain return lon_min - ADT_PAD, lon_max + ADT_PAD, lat_min - ADT_PAD, lat_max + ADT_PAD + def _pm180(self) -> bool: + from ..coords import lon_convention + return lon_convention(self.config) != "0360" + + def _roi(self, ssh_path): + """Ops fixed ROI for -180/180 systems (STOFS-3D-ATL); SSH_1 extent for the 0-360 Pacific. MJ (10/24/26)""" + if self._pm180(): + return ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX + return self._dest_bounds(ssh_path) + def _find_adt_data(self) -> Optional[Path]: """Newest available CMEMS ADT file (valid day, else previous day).""" files = self._find_adt_files() @@ -309,7 +325,7 @@ def dlon_diff(a, b): if row.min() < 0 or row.max() >= n_b or col.min() < 0 or col.max() >= n_a: raise ValueError("map row/col out of range") - bounds = self._dest_bounds(ssh_path) + bounds = self._roi(ssh_path) days = [] for f in adt_files: a, lons, lats = self._subset_ops_roi(f, bounds) @@ -323,7 +339,8 @@ def dlon_diff(a, b): or np.abs(yca - LA.ravel()).max() > ADT_COORD_TOL): raise ValueError("ADT ROI coordinates do not match map xc_a/yc_a") out = self._apply_esmf_map(S, row, col, n_b, a.ravel()) - days.append(out.astype(np.float32) if self.ops_numerics else out) + off = self.config.adt_datum_offset + days.append((out + off).astype(np.float32) if self.ops_numerics else out + off) if self.ops_numerics: stack = np.stack(days) ok = np.isfinite(stack) @@ -429,7 +446,7 @@ def _read_adt(self, adt_path: Path, bounds=None): bounds is (lon_min, lon_max, lat_min, lat_max) in -180/180; default is the ops ROI. Returns (adt[lat, lon], lons, lats) or None. Fill values become NaN. """ - lon_min, lon_max, lat_min, lat_max = bounds or self._dest_bounds(None) + lon_min, lon_max, lat_min, lat_max = bounds or self._roi(None) try: ds = Dataset(str(adt_path)) try: @@ -439,7 +456,7 @@ def _read_adt(self, adt_path: Path, bounds=None): lons = np.array(ds.variables[lon_name][:], dtype=np.float64) lats = np.array(ds.variables[lat_name][:], dtype=np.float64) adt_var = "adt" if "adt" in ds.variables else "surf_el" - adt_data = ds.variables[adt_var][...] + adt_data = ds.variables[adt_var][...] + self.config.adt_datum_offset from ..coords import lon_convention, shift_grid_0360 if lon_convention(self.config) == "0360": @@ -577,7 +594,7 @@ def _apply_adt_blend(self, ssh_path: Path, adt_field, work_dir: Path, shutil.copy2(ssh_path, output) ds = Dataset(str(output), "r+") - ssh = ds.variables["ssh"][...] + (self.config.adt_datum_offset or 0) + ssh = ds.variables["ssh"] nt, ny, nx = ssh.shape dst_lon = np.array(ds.variables["xlon"][:], dtype=np.float64) diff --git a/nos_utils/forcing/nudging.py b/nos_utils/forcing/nudging.py index 8e0f1bc..d1586ea 100644 --- a/nos_utils/forcing/nudging.py +++ b/nos_utils/forcing/nudging.py @@ -742,7 +742,7 @@ def _process_python(self) -> ForcingResult: # time records (rare, but the reader supports it), they are # treated as consecutive 6h steps starting at the file's stamped # forecast hour. - rtofs_dt_input = self.config.nudging_dt + rtofs_dt_input = self.config.nudging_dt if getattr(self.config, "rtofs_daily_3z", False) else 21600.0 is_split = getattr(self.config, "rtofs_daily_3z", False) if is_split: @@ -899,8 +899,8 @@ def get_time_key(fp): return fp.name.split('_daily')[0] # the model t=0 anchor (= cycle - nowcast_hours). target_dt = getattr(self.config, "nudging_dt", 10800.0) if self.config.obc_ops_timeline: - # Ops phase-anchored records, time 0 = phase start. - target_dt = self.config.nudging_dt + # step_nu_tr = 21600 s: phase-anchored 6-hourly records, time 0 = phase start. MJ (10/02/26) + target_dt = 21600.0 # The output time axis covers [0, sim_duration] at ``target_dt`` # cadence — i.e. 0, 3h, 6h, ..., nowcast+forecast hours from the diff --git a/nos_utils/forcing/rtofs.py b/nos_utils/forcing/rtofs.py index 623dd48..7778d3e 100644 --- a/nos_utils/forcing/rtofs.py +++ b/nos_utils/forcing/rtofs.py @@ -1296,11 +1296,11 @@ def _slice_phase_records(self, path: Path) -> None: model_t0, _, dur = self._get_output_window(self.phase) cycle = datetime.strptime(self.config.pdy, "%Y%m%d") + timedelta(hours=self.config.cyc) start = (model_t0 - (cycle - timedelta(hours=self.config.nowcast_hours))).total_seconds() - n = int(np.ceil(dur / OPS_OBC_DT)) + 1 + n = int(np.ceil(dur / self.config.obc_dt)) + 1 tmp = path.with_name(path.name + ".phase") with Dataset(str(path)) as src: t = np.asarray(src.variables["time"][:], dtype=np.float64) - dt = float(np.diff(t).min()) if len(t) > 1 else OPS_OBC_DT + dt = float(np.diff(t).min()) if len(t) > 1 else self.config.obc_dt idx = np.clip(np.rint(start / dt).astype(int) + np.arange(n), 0, len(t) - 1) if start / dt + n - 1 > len(t) - 1: self._ops_note(f"{path.name}: the exe series ends before the {n}-record " @@ -1586,7 +1586,7 @@ def _ops_slot_times(self) -> List[datetime]: t0 = cycle_dt - timedelta(hours=int(self.config.nowcast_hours)) span = (cycle_dt - t0).total_seconds() + ( int(self.config.forecast_hours) + self.buffer_hours) * 3600.0 - n = max(self.config.obc_min_timesteps, int(np.ceil(span / self.config.obc_dt)) + 1) + n = max(self.config.obc_min_slots, int(np.ceil(span / self.config.obc_dt)) + 1) return [t0 + timedelta(seconds=self.config.obc_dt * k) for k in range(n)] def _ops_candidates(self, date: datetime, kind: str): @@ -1655,10 +1655,10 @@ def _find_ops_files(self) -> Tuple[List[Path], List[Path]]: if have and n_prev[kind]: self._select_notes.append( f"RTOFS {kind}: {n_prev[kind]} of {len(have)} slots filled from the previous cycle") - miss = [vt for vt in slots[:self.config.obc_min_timesteps] if vt not in picked[kind]] if have else [] + miss = [vt for vt in slots[:self.config.obc_min_slots] if vt not in picked[kind]] if have else [] if miss: self._select_notes.append( - f"RTOFS {kind}: {len(miss)} of the first {self.config.obc_min_timesteps} slots missing " + f"RTOFS {kind}: {len(miss)} of the first {self.config.obc_min_slots} slots missing " f"in both cycles (first {miss[0]:%Y-%m-%d %Hz})") for m in self._select_notes: log.warning(m) diff --git a/nos_utils/forcing/surface_restore.py b/nos_utils/forcing/surface_restore.py index 48869b9..a9c3a95 100755 --- a/nos_utils/forcing/surface_restore.py +++ b/nos_utils/forcing/surface_restore.py @@ -43,6 +43,12 @@ def __init__( self.weight_file = weight_file self.phase = phase + def find_input_files(self) -> List[Path]: + """Abstract in the base class and was never implemented, so the class could not be instantiated. MJ (10/24/26)""" + from .rtofs import RTOFSProcessor + return RTOFSProcessor(self.config, self.input_path, self.output_path, + phase=self.phase).find_input_files_by_type()[1] + def process(self) -> ForcingResult: if not HAS_NETCDF4: return ForcingResult( @@ -56,6 +62,15 @@ def process(self) -> ForcingResult: warnings=["Surface restore disabled in config"], ) + import os + if os.environ.get("NOS_SURFACE_RESTORE_EXPERIMENTAL") != "1": + msg = ("Surface restore is enabled but this processor is untested and known incomplete " + "(T and S read from one file, boundary instead of restore nodes, placeholder time axis); " + "set NOS_SURFACE_RESTORE_EXPERIMENTAL=1 to run it anyway. MJ (10/24/26)") + log.error(msg) + return ForcingResult(success=False, source=self.SOURCE_NAME, errors=[msg]) + log.warning("Surface restore running in experimental mode: output is unvalidated. MJ (10/24/26)") + self.create_output_dir() # Reuse RTOFS 3D discovery to get the raw daily/6-hourly files diff --git a/nos_utils/orchestrator.py b/nos_utils/orchestrator.py index 4004bdd..43456be 100644 --- a/nos_utils/orchestrator.py +++ b/nos_utils/orchestrator.py @@ -679,26 +679,40 @@ def _prev_cycle_dirs(self, max_days: int = 3) -> List[Path]: dirs += [r / f"{self.run_name}.{day}", r / day] return list(dict.fromkeys(dirs)) + # Mapping from the active OBC filename (what SCHISM reads) to the + # standard archive basename under $COMOUT_PREV/rerun/ (see + # stofs_3d_atl_create_obc_3d_th_non_adjust.sh:641-642). The archive + # file is stored as `...nc`. + _OBC_ARCHIVE_BASE: Dict[str, str] = { + "elev2D.th.nc": "elev2dth_non_adj", + "TEM_3D.th.nc": "tem3dth", + "SAL_3D.th.nc": "sal3dth", + "uv3D.th.nc": "uv3dth", + } + _OBC_TAR_BASE = ("elev2D.th.nc", "TEM_3D.th.nc", "SAL_3D.th.nc", + "uv3D.th.nc", "TEM_nu.nc", "SAL_nu.nc") + @property def _OBC_ARCHIVE_NAMES(self) -> Dict[str, str]: - """Mapping of active OBC files to their archive basenames. - - Dynamically includes nudging and surface restore if enabled, and handles - the difference between ATL (elev2dth_non_adj) and PAC (elev2dth). - """ - names = { - "elev2D.th.nc": "elev2dth_non_adj" if getattr(self.config, "dynamic_adjust_enabled", False) else "elev2dth", - "TEM_3D.th.nc": "tem3dth", - "SAL_3D.th.nc": "sal3dth", - "uv3D.th.nc": "uv3dth", - } - if getattr(self.config, "nudging_enabled", False): - names["TEM_nu.nc"] = "temnu" - names["SAL_nu.nc"] = "salnu" - if getattr(self.config, "surface_restore_enabled", False): - names["surface_restore.nc"] = "surface_restore" + """Rerun-archive map: as nos-utils main unless obc_archive_extended is set (Pacific). MJ (10/24/26)""" + names = dict(self._OBC_ARCHIVE_BASE) + if getattr(self.config, "obc_archive_extended", False): + names["elev2D.th.nc"] = "elev2dth_non_adj" if getattr(self.config, "dynamic_adjust_enabled", False) else "elev2dth" + if getattr(self.config, "nudging_enabled", False): + names["TEM_nu.nc"] = "temnu" + names["SAL_nu.nc"] = "salnu" + if getattr(self.config, "surface_restore_enabled", False): + names["surface_restore.nc"] = "surface_restore" return names + @property + def _obc_tar_files(self) -> List[str]: + files = list(self._OBC_TAR_BASE) + if getattr(self.config, "obc_archive_extended", False) and getattr( + self.config, "surface_restore_enabled", False): + files.append("surface_restore.nc") + return files + def _obc_archive_name(self, archive_base: str, phase: str) -> str: """``...nc``; the forecast-phase file gets ``.fcst`` before the suffix. @@ -1256,7 +1270,7 @@ def _build_archive_manifest( # --- COMMON: OBC (phase-specific, backward-compat) --- { "tar_name": f"{prefix}.{cycle}.{pdy}.obc.{phase}.tar", - "payload_files": list(self._OBC_ARCHIVE_NAMES.keys()), + "payload_files": self._obc_tar_files, "condition": True, "label": "OBC", }, @@ -1589,7 +1603,7 @@ def archive_to_comout(self, result: PrepResult, comout: Path) -> List[Path]: # Tar OBC files (phase-specific): obc.nowcast.tar / obc.forecast.tar, # each with all active files (boundary + nudging + surface restore). - obc_files = list(self._OBC_ARCHIVE_NAMES.keys()) + obc_files = self._obc_tar_files existing_obc = [work_dir / f for f in obc_files if (work_dir / f).exists()] if existing_obc: # Phase-specific tar (backward compat): obc.nowcast.tar / obc.forecast.tar diff --git a/tests/test_stofs_config.py b/tests/test_stofs_config.py index 5cbb328..5cd2401 100644 --- a/tests/test_stofs_config.py +++ b/tests/test_stofs_config.py @@ -316,7 +316,7 @@ def test_domain_bounds_0_360(self): def test_run_window(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) assert cfg.nowcast_hours == 24 - assert cfg.forecast_hours == 108 + assert cfg.forecast_hours == 48 def test_model_dt_is_90(self): """PAC model runs at dt=90s (verified from param.nml_6globaloutput).""" @@ -338,7 +338,8 @@ def test_no_dynamic_adjust(self): def test_no_adt(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) - assert cfg.adt_enabled is False + assert cfg.adt_enabled is True + assert cfg.adt_datum_offset == pytest.approx(-0.4715517) def test_nudging_enabled(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) @@ -356,7 +357,8 @@ def test_nwm_product(self): def test_obc_min_timesteps(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) - assert cfg.obc_min_timesteps == 21 + assert cfg.obc_min_timesteps == 4 + assert cfg.obc_min_slots == 4 def test_gfs_resolution(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) @@ -508,3 +510,59 @@ def test_ufs_tasks_parsed(self, tmp_path): assert cfg.ufs_datm_tasks == 120 assert cfg.ufs_total_tasks == 4072 + + +class TestPacificBranchKeepsExistingSystems: + """Pacific (#74) options must not change ATL / SECOFS / AK. MJ (10/24/26)""" + + @staticmethod + def _from_yaml(tmp_path, ocean, run=None): + import yaml + doc = { + "grid": {"domain": {"lon_min": -98.5, "lon_max": -52.5, + "lat_min": 7.3, "lat_max": 52.6}, "n_levels": 51}, + "model": {"physics": {"nws": 2}, "run": run or {"nowcast_hours": 24, "forecast_hours": 96}}, + "forcing": {"ocean": ocean}, + } + f = tmp_path / "x.yaml" + f.write_text(yaml.dump(doc)) + return ForcingConfig.from_yaml(f, pdy="20260401", cyc=12) + + def test_atl_adt_offset_with_and_without_yaml_key(self, tmp_path): + stofs = {"obc": {"roi_2ds": {"x1": 1, "x2": 2, "y1": 1, "y2": 2}, + "roi_3dz": {"x1": 1, "x2": 2, "y1": 1, "y2": 2}}, + "adt": {"enabled": True}} + assert self._from_yaml(tmp_path, stofs).adt_datum_offset == -0.45 + assert ForcingConfig.for_stofs_3d_atl(pdy="20260401", cyc=12).adt_datum_offset == -0.45 + stofs["adt"]["datum_offset"] = -0.4715517 + assert self._from_yaml(tmp_path, stofs).adt_datum_offset == pytest.approx(-0.4715517) + + def test_atl_keeps_ops_rtofs_slots_and_qc(self, tmp_path): + cfg = ForcingConfig.for_stofs_3d_atl(pdy="20260401", cyc=12) + assert cfg.obc_min_slots == 23 + assert cfg.obc_min_timesteps == 21 + assert cfg.obc_dt == 21600.0 + assert cfg.obc_archive_extended is False + assert cfg.surface_restore_enabled is False + + def test_secofs_and_unspecified_yaml_keep_3h_nudging(self, tmp_path): + assert ForcingConfig.for_secofs(pdy="20260401", cyc=12).nudging_dt == 10800.0 + assert ForcingConfig.for_secofs_ufs(pdy="20260401", cyc=12).nudging_dt == 10800.0 + cfg = self._from_yaml(tmp_path, {"nudging": {"enabled": True}}) + assert cfg.nudging_dt == 10800.0 + assert cfg.obc_archive_extended is False + + def test_pacific_sets_its_own_values(self): + cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) + assert cfg.nudging_dt == 86400.0 + assert cfg.obc_dt == 86400.0 + assert cfg.obc_archive_extended is True + assert cfg.surface_restore_enabled is False + + def test_surface_restore_refuses_without_experimental_env(self, tmp_path, monkeypatch): + from nos_utils.forcing.surface_restore import SurfaceRestoreProcessor + monkeypatch.delenv("NOS_SURFACE_RESTORE_EXPERIMENTAL", raising=False) + cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12, surface_restore_enabled=True) + res = SurfaceRestoreProcessor(cfg, tmp_path, tmp_path / "out").process() + assert not res.success + assert "untested" in res.errors[0] From 177f8eed4129217ffe19868d410f2ce7e5bb502a Mon Sep 17 00:00:00 2001 From: mansurjisan Date: Sat, 24 Oct 2026 09:10:49 -0400 Subject: [PATCH 2/2] fix(rtofs): keep depth in TSUV_1.nc and close datasets; tolerate null yaml ocean blocks; more guard tests --- nos_utils/config.py | 6 ++++- nos_utils/forcing/adt.py | 6 ++--- nos_utils/forcing/rtofs.py | 11 +++++---- tests/test_stofs_config.py | 48 +++++++++++++++++++++++++++++++++++++- 4 files changed, 61 insertions(+), 10 deletions(-) diff --git a/nos_utils/config.py b/nos_utils/config.py index 70b265c..d537c6d 100644 --- a/nos_utils/config.py +++ b/nos_utils/config.py @@ -1082,13 +1082,17 @@ def from_yaml(cls, yaml_path, **overrides) -> "ForcingConfig": if rtofs_3d_region: kwargs["rtofs_3d_region"] = str(rtofs_3d_region) + ocean = ocean if isinstance(ocean, dict) else {} + obc = obc if isinstance(obc, dict) else {} + nudge = nudge if isinstance(nudge, dict) else {} kwargs["rtofs_daily_3z"] = _strict_bool(ocean.get("rtofs_daily_3z"), "ocean.rtofs_daily_3z", False) kwargs["obc_dt"] = float(obc.get("dt", 21600.0)) kwargs["obc_time_smoothing"] = _strict_bool(obc.get("time_smoothing"), "obc.time_smoothing", False) kwargs["nudging_dt"] = float(nudge.get("dt", 10800.0)) kwargs["nudging_time_smoothing"] = _strict_bool(nudge.get("time_smoothing"), "nudging.time_smoothing", False) - surf_res = ocean.get("surface_restore", {}) if isinstance(ocean, dict) else {} + surf_res = ocean.get("surface_restore") + surf_res = surf_res if isinstance(surf_res, dict) else {} kwargs["surface_restore_enabled"] = _strict_bool(surf_res.get("enabled"), "ocean.surface_restore.enabled", False) kwargs["surface_restore_time_smoothing"] = _strict_bool(surf_res.get("time_smoothing"), "ocean.surface_restore.time_smoothing", False) if surf_res.get("weight_file"): diff --git a/nos_utils/forcing/adt.py b/nos_utils/forcing/adt.py index 4016dc2..fe508f4 100644 --- a/nos_utils/forcing/adt.py +++ b/nos_utils/forcing/adt.py @@ -149,7 +149,7 @@ def blend_ssh(self, ssh_path: Path, work_dir: Path) -> Optional[Path]: def _dest_bounds(self, ssh_path: Path): """ADT subset (lon_min, lon_max, lat_min, lat_max): SSH_1 grid extent padded - by ADT_PAD, else the config domain bounds. MJ (10/01/26) + by ADT_PAD, else the ops ROI (-180/180 systems) or the padded config domain (0-360). MJ (10/01/26) """ try: with Dataset(str(ssh_path)) as ds: @@ -158,7 +158,7 @@ def _dest_bounds(self, ssh_path: Path): return (float(np.nanmin(lon)) - ADT_PAD, float(np.nanmax(lon)) + ADT_PAD, float(np.nanmin(lat)) - ADT_PAD, float(np.nanmax(lat)) + ADT_PAD) except Exception as e: - log.warning(f"Cannot read SSH_1 grid extent ({e}) — using config domain bounds") + log.warning(f"Cannot read SSH_1 grid extent ({e}) — using the fallback ROI") if self._pm180(): return ADT_LON_MIN, ADT_LON_MAX, ADT_LAT_MIN, ADT_LAT_MAX lon_min, lon_max, lat_min, lat_max = self.config.domain @@ -585,7 +585,7 @@ def _apply_adt_blend(self, ssh_path: Path, adt_field, work_dir: Path, """Apply the ops ADT blend to SSH_1.nc. Formula: SSH_final(t) = SSH_rtofs(t) - SSH_rtofs(t=0) + ADT + datum_offset - (ADT datum offset, e.g. -0.45 is applied here), in float32 as ncap2. Where either + (ADT already carries adt_datum_offset, applied per day on read/regrid), in float32 as ncap2. Where either side is invalid the point is left as the -30000 fill (dry in gen_3Dth). """ try: diff --git a/nos_utils/forcing/rtofs.py b/nos_utils/forcing/rtofs.py index 7778d3e..5d01f95 100644 --- a/nos_utils/forcing/rtofs.py +++ b/nos_utils/forcing/rtofs.py @@ -1019,6 +1019,7 @@ def get_time_key(fp): return fp.name.split('_daily')[0] file_groups = [[f] for f in files_3d] for group in file_groups: + datasets = [] try: datasets = [Dataset(str(f)) for f in group] # Merge variables from multiple files into one dictionary @@ -1030,7 +1031,6 @@ def get_time_key(fp): return fp.name.split('_daily')[0] subset.update(self._stofs_subset_roi(ds, roi, ["temperature", "salinity", "u", "v"])) if "temperature" not in subset: - for ds in datasets: ds.close() continue for t in range(subset["temperature"].shape[0]): @@ -1043,12 +1043,13 @@ def get_time_key(fp): return fp.name.split('_daily')[0] if all_lon is None: all_lon = np.array(subset["Longitude"]) all_lat = np.array(subset["Latitude"]) - if all_depth is None and "Depth" in ds.variables: - all_depth = np.array(ds.variables["Depth"][:]) - - for ds in datasets: ds.close() + if all_depth is None and "Depth" in datasets[0].variables: + all_depth = np.array(datasets[0].variables["Depth"][:]) except Exception as e: log.warning(f"Failed to read 3D group {group[0].name}: {e}") + finally: + for ds in datasets: + ds.close() if not all_temp: return None diff --git a/tests/test_stofs_config.py b/tests/test_stofs_config.py index 5cd2401..126966c 100644 --- a/tests/test_stofs_config.py +++ b/tests/test_stofs_config.py @@ -336,7 +336,7 @@ def test_no_dynamic_adjust(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) assert cfg.dynamic_adjust_enabled is False - def test_no_adt(self): + def test_adt_enabled_with_own_offset(self): cfg = ForcingConfig.for_stofs_3d_pac(pdy="20260401", cyc=12) assert cfg.adt_enabled is True assert cfg.adt_datum_offset == pytest.approx(-0.4715517) @@ -566,3 +566,49 @@ def test_surface_restore_refuses_without_experimental_env(self, tmp_path, monkey res = SurfaceRestoreProcessor(cfg, tmp_path, tmp_path / "out").process() assert not res.success assert "untested" in res.errors[0] + + def test_from_yaml_parses_min_slots_and_archive_extended(self, tmp_path): + base = {"nudging": {"enabled": True}} + cfg = self._from_yaml(tmp_path, base) + assert cfg.obc_min_slots == 23 and cfg.obc_archive_extended is False + cfg = self._from_yaml(tmp_path, {"obc": {"min_slots": 4, "archive_extended": True}}) + assert cfg.obc_min_slots == 4 and cfg.obc_archive_extended is True + + def test_from_yaml_tolerates_null_ocean_blocks(self, tmp_path): + cfg = self._from_yaml(tmp_path, {"nudging": None, "surface_restore": None}) + assert cfg.nudging_dt == 10800.0 + assert self._from_yaml(tmp_path, None).nudging_dt == 10800.0 + + def test_atl_rerun_archive_map_is_the_four_base_entries(self): + from nos_utils.orchestrator import PrepOrchestrator + orch = PrepOrchestrator(ForcingConfig.for_stofs_3d_atl(pdy="20260401", cyc=12), {}) + assert orch._OBC_ARCHIVE_NAMES == { + "elev2D.th.nc": "elev2dth_non_adj", "TEM_3D.th.nc": "tem3dth", + "SAL_3D.th.nc": "sal3dth", "uv3D.th.nc": "uv3dth"} + + def test_tsuv_keeps_depth_and_logs_no_failure(self, tmp_path, caplog): + import logging + import numpy as np + from netCDF4 import Dataset + from nos_utils.forcing.rtofs import RTOFSProcessor + rng = np.random.default_rng(0) + files = [] + for tag in ("n000", "f006"): + p = tmp_path / f"rtofs_glo_3dz_{tag}_6hrly_hvr_US_east.nc" + with Dataset(p, "w") as d: + d.createDimension("MT", 1); d.createDimension("Depth", 3) + d.createDimension("Y", 5); d.createDimension("X", 6) + d.createVariable("Depth", "f4", ("Depth",))[:] = [0, 10, 20] + d.createVariable("Longitude", "f4", ("Y", "X"))[:] = np.tile(np.arange(6.), (5, 1)) - 70 + d.createVariable("Latitude", "f4", ("Y", "X"))[:] = np.tile(np.arange(5.)[:, None], (1, 6)) + 30 + for v in ("temperature", "salinity", "u", "v"): + d.createVariable(v, "f4", ("MT", "Depth", "Y", "X"))[:] = rng.random((1, 3, 5, 6)) + files.append(p) + cfg = ForcingConfig.for_stofs_3d_atl(pdy="20260401", cyc=12) + cfg.obc_roi_3d = {"x1": 0, "x2": 4, "y1": 0, "y2": 3} + proc = RTOFSProcessor(cfg, tmp_path, tmp_path / "out") + with caplog.at_level(logging.WARNING): + out = proc._stofs_prepare_tsuv(files, tmp_path) + assert "Failed to read 3D group" not in caplog.text + with Dataset(out) as d: + assert "depth" in d.variables