From 1818558bda7ed03d0bc2d3354116e8a984ce61c6 Mon Sep 17 00:00:00 2001 From: mansurjisan Date: Sun, 11 Oct 2026 20:02:43 -0400 Subject: [PATCH 1/4] adcirc_met sfcf: one hourly 180 h forecast segment, linear interpolation between the 3-hourly GFS records after f120 --- nos_utils/forcing/adcirc_met.py | 65 +++++++++++++++++++++++-------- tests/test_adcirc_met.py | 69 ++++++++++++++++++++------------- 2 files changed, 91 insertions(+), 43 deletions(-) diff --git a/nos_utils/forcing/adcirc_met.py b/nos_utils/forcing/adcirc_met.py index fae7457..36bdd31 100644 --- a/nos_utils/forcing/adcirc_met.py +++ b/nos_utils/forcing/adcirc_met.py @@ -594,7 +594,7 @@ def _write_manifest(path: Path, acquired: List[AcquiredFile], phase: str, # Replicates exstofs_2d_glo_gfs.sh / stofs_2d_glo_surface_forcing.sh / stofs_2d_glo_getges.sh: # one pass per record, ncks -v + ncecat + ncwa -a time as netCDF4 calls. MJ (10/06/26) -SFC_SEGMENTS = ("ncst", "fcst1", "fcst2") +SFC_SEGMENTS = ("ncst", "fcst") SFC_OUTPUTS = (("221", ("pressfc",)), ("222", ("ugrd10m", "vgrd10m")), ("225", ("icec",))) _GETGES_FHEND = 384 _SFCF_SETTLED_S = 900 # an sfcf file this old without a logf marker is a retrospective one. MJ (10/06/26) @@ -622,25 +622,33 @@ def getges_sfcf(comin_gfs, valid: datetime) -> Path: def sfcf_valid_times(cycle: datetime, segment: str, start: Optional[datetime] = None) -> List[datetime]: - """Record times of the ops calls: surface1 ncst (start..cycle) and fcst1 (cycle..+120 h) hourly, - surface3 fcst2 (+120 h..+180 h) every 3 h. ``start`` is the ncst time_beg (multistart).""" + """Hourly record times: ncst (start..cycle) and the single fcst (cycle..+180 h). ``start`` is the + ncst time_beg (multistart). After f120 GFS is 3-hourly, so fcst fills the hours between by linear + interpolation, which is what ADCIRC does between the fcst2 snapshots. MJ (10/11/26)""" if segment not in SFC_SEGMENTS: raise ValueError("segment must be one of {}, got {!r}".format(SFC_SEGMENTS, segment)) if segment == "ncst": - beg, end, step = start or cycle - timedelta(hours=6), cycle, 1 - elif segment == "fcst1": - beg, end, step = cycle, cycle + timedelta(hours=120), 1 + beg, end = start or cycle - timedelta(hours=6), cycle else: - beg, end, step = cycle + timedelta(hours=120), cycle + timedelta(hours=180), 3 + beg, end = cycle, cycle + timedelta(hours=180) if beg > end: raise MetForcingError("sfcf window starts after it ends: {} > {}".format(beg, end)) out = [] while beg <= end: out.append(beg) - beg += timedelta(hours=step) + beg += timedelta(hours=1) return out +def _interp_bracket(cycle: datetime, valid: datetime) -> Optional[Tuple[int, int, float]]: + """(lo, hi, weight) forecast hours when ``valid`` falls between two 3-hourly GFS records, else None.""" + fh = int((valid - cycle).total_seconds() // 3600) + if fh <= 120 or fh % 3 == 0: + return None + lo = fh // 3 * 3 + return lo, lo + 3, (fh - lo) / 3.0 + + def _sfcf_ready(path: Path) -> bool: """Non-empty and, unless old enough to be complete, its logf marker exists (GFS writes sfcf straight into COM).""" if not (path.is_file() and path.stat().st_size > 0): @@ -663,7 +671,7 @@ def _wait_for(path: Path, deadline: float, poll_s: float, sleep, clock) -> bool: def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix: str = "stofs_2d_glo", start: Optional[datetime] = None, wait_s: float = 0.0, poll_s: float = 10.0, sleep=time.sleep, clock=time.monotonic, now=time.time) -> Dict[str, Path]: - """Write ``_.{221,222,225}.nc`` from the GFS sfcf files, as the ops GFS_NCST/FCST1/FCST2 jobs. + """Write ``_.{221,222,225}.nc`` from the GFS sfcf files (ncst, or one hourly 180 h fcst). Valid times at or before ``cycle`` use getges; later ones the current cycle's own sfcf file, waited for up to ``wait_s`` seconds in total. Output files appear under their final names only when all three are @@ -684,7 +692,14 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix done = False try: for rec, valid in enumerate(valids): - if valid <= cycle: + br = _interp_bracket(cycle, valid) if segment == "fcst" else None + if br: + srcs = [sfcf_path(comin_gfs, cycle, h) for h in br[:2]] + for sp in srcs: + if not _wait_for(sp, deadline, poll_s, sleep, clock): + raise FileNotAvailableError("GFS file did not appear before the wait expired: {}".format(sp)) + src = srcs[0] + elif valid <= cycle: if valid == cycle and now() - calendar.timegm(cycle.timetuple()) < _REALTIME_S: _wait_for(sfcf_path(comin_gfs, cycle, 0), deadline, poll_s, sleep, clock) src = getges_sfcf(comin_gfs, valid) @@ -700,13 +715,29 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix for n, names in SFC_OUTPUTS: outs[n] = _open_sfcf_output(nc, part[n], ds, names) lat, lon = ds.variables["lat"][:], ds.variables["lon"][:] - for n, names in SFC_OUTPUTS: - o = outs[n] - if lat.shape != (len(o.dimensions["grid_yt"]), len(o.dimensions["grid_xt"])): # ncecat needs equal shapes. MJ (10/06/26) - raise MetForcingError("{}: grid differs from the first record".format(src)) - o.variables["lat"][rec], o.variables["lon"][rec] = lat, lon - for v in names: - o.variables[v][rec] = ds.variables[v][0] + hi = nc.Dataset(str(srcs[1])) if br else None + try: + if hi is not None: + hi.set_auto_maskandscale(False) + for n, names in SFC_OUTPUTS: + o = outs[n] + if lat.shape != (len(o.dimensions["grid_yt"]), len(o.dimensions["grid_xt"])): # ncecat needs equal shapes. MJ (10/06/26) + raise MetForcingError("{}: grid differs from the first record".format(src)) + o.variables["lat"][rec], o.variables["lon"][rec] = lat, lon + for v in names: + a0 = ds.variables[v][0] + if hi is None: + o.variables[v][rec] = a0 + continue + a1 = hi.variables[v][0] + r = a0.astype(np.float64) * (1.0 - br[2]) + a1.astype(np.float64) * br[2] + fv = ds.variables[v].getncattr("_FillValue") if "_FillValue" in ds.variables[v].ncattrs() else None + if fv is not None: + r = np.where((a0 == fv) | (a1 == fv), fv, r) + o.variables[v][rec] = r.astype(a0.dtype) + finally: + if hi is not None: + hi.close() for o in outs.values(): o.close() outs = {} diff --git a/tests/test_adcirc_met.py b/tests/test_adcirc_met.py index 00cf521..348f128 100644 --- a/tests/test_adcirc_met.py +++ b/tests/test_adcirc_met.py @@ -430,13 +430,11 @@ def test_sfcf_windows(): ncst = am.sfcf_valid_times(CYC, "ncst") assert (len(ncst), ncst[0], ncst[-1]) == (7, CYC - timedelta(hours=6), CYC) assert len(am.sfcf_valid_times(CYC, "ncst", CYC - timedelta(hours=12))) == 13 - f1 = am.sfcf_valid_times(CYC, "fcst1") - assert (len(f1), f1[0], f1[-1]) == (121, CYC, CYC + timedelta(hours=120)) - f2 = am.sfcf_valid_times(CYC, "fcst2") - assert (len(f2), f2[0], f2[1], f2[-1]) == ( - 21, CYC + timedelta(hours=120), CYC + timedelta(hours=123), CYC + timedelta(hours=180)) - with pytest.raises(ValueError): - am.sfcf_valid_times(CYC, "spinup") + f = am.sfcf_valid_times(CYC, "fcst") + assert (len(f), f[0], f[-1]) == (181, CYC, CYC + timedelta(hours=180)) + for old in ("spinup", "fcst1", "fcst2"): + with pytest.raises(ValueError): + am.sfcf_valid_times(CYC, old) with pytest.raises(am.MetForcingError): am.sfcf_valid_times(CYC, "ncst", CYC + timedelta(hours=1)) @@ -488,23 +486,41 @@ def test_ncst_start_override_gives_fewer_records(tmp_path): assert _read(files["225"], "icec").shape[0] == 6 -def test_fcst1_and_fcst2_record_counts_and_stride(tmp_path): +def _fcst_tank(tank, last=180): + _fill_tank(tank, [(CYC, h) for h in range(0, 121)] + [(CYC, h) for h in range(123, last + 1, 3)]) + + +def test_fcst_is_hourly_with_interpolation_after_f120(tmp_path): tank = tmp_path / "g" - _fill_tank(tank, [(CYC, 0)] + [(CYC, h) for h in range(1, 121)]) - f1 = am.build_sfcf_forcing(tank, CYC, "fcst1", tmp_path / "o") - p = _read(f1["221"], "pressfc") - assert p.shape[0] == 121 - assert p[7][0, 0] == _valid_seed(CYC, 7) * 10 + 2 - _fill_tank(tank, [(CYC, h) for h in range(120, 181, 3)]) - f2 = am.build_sfcf_forcing(tank, CYC, "fcst2", tmp_path / "o") - q = _read(f2["222"], "ugrd10m") - assert q.shape[0] == 21 - assert q[1][0, 0] == _valid_seed(CYC, 123) * 10 + 0 + _fcst_tank(tank) + f = am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o") + p = _read(f["221"], "pressfc") + assert p.shape[0] == 181 + s = lambda h, k: _valid_seed(CYC, h) * 10 + k + for h in (0, 7, 120, 123, 126, 180): # native records are exact copies + assert p[h][0, 0] == s(h, 2) + # 121 and 122 are 1/3 and 2/3 of the way from f120 to f123; the seeds step 10 per hour + assert p[121][0, 0] == pytest.approx(s(120, 2) + 10.0) + assert p[122][0, 0] == pytest.approx(s(120, 2) + 20.0) + q = _read(f["222"], "ugrd10m") + assert q[121][1, 2] == pytest.approx(s(120, 0) + 10.0 + 6) + assert _read(f["225"], "icec").shape[0] == 181 + + +def test_fcst_interpolation_keeps_fill_values(tmp_path): + tank = tmp_path / "g" + _fcst_tank(tank) + import netCDF4 as nc + with nc.Dataset(str(am.sfcf_path(tank, CYC, 123)), "a") as d: + d.variables["icec"][0, 0, 0] = np.float32(9.99e20) + f = am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o") + i = _read(f["225"], "icec") + assert i[121][0, 0] == np.float32(9.99e20) and i[120][0, 0] < 1e6 and i[121][0, 1] < 1e6 def test_forecast_waits_then_fails_clearly(tmp_path): tank = tmp_path / "g" - _fill_tank(tank, [(CYC, h) for h in range(120, 178, 3)]) # f180 is missing + _fcst_tank(tank, last=177) # f180 is missing clock = {"t": 0.0} naps = [] @@ -513,7 +529,7 @@ def sleep(s): clock["t"] += s with pytest.raises(am.FileNotAvailableError, match="sfcf180"): - am.build_sfcf_forcing(tank, CYC, "fcst2", tmp_path / "o", wait_s=25, poll_s=10, + am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o", wait_s=25, poll_s=10, sleep=sleep, clock=lambda: clock["t"]) assert naps == [10, 10, 5] assert list((tmp_path / "o").glob("*.nc*")) == [] # no final or .partial files from a failed run @@ -522,17 +538,18 @@ def sleep_and_land(s): _mk_sfcf(tank, CYC, 180, 1) sleep(s) - files = am.build_sfcf_forcing(tank, CYC, "fcst2", tmp_path / "o", wait_s=25, poll_s=10, + files = am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o", wait_s=25, poll_s=10, sleep=sleep_and_land, clock=lambda: clock["t"]) - assert _read(files["221"], "pressfc").shape[0] == 21 + assert _read(files["221"], "pressfc").shape[0] == 181 def test_current_cycle_f000_is_not_waited_for(tmp_path): - # fcst1 record 0 (valid == cycle) goes through getges, so it never waits - _fill_tank(tmp_path / "g", [(CYC - timedelta(hours=6), 6)] + [(CYC, h) for h in range(1, 121)]) - files = am.build_sfcf_forcing(tmp_path / "g", CYC, "fcst1", tmp_path / "o", wait_s=0, + # fcst record 0 (valid == cycle) goes through getges, so it never waits + _fill_tank(tmp_path / "g", [(CYC - timedelta(hours=6), 6)] + [(CYC, h) for h in range(1, 121)] + + [(CYC, h) for h in range(123, 181, 3)]) + files = am.build_sfcf_forcing(tmp_path / "g", CYC, "fcst", tmp_path / "o", wait_s=0, sleep=lambda s: pytest.fail("waited")) - assert _read(files["221"], "pressfc").shape[0] == 121 + assert _read(files["221"], "pressfc").shape[0] == 181 def _age(p, seconds): From f12a42c949a7c66ff71c1a25a0707b62166dd746 Mon Sep 17 00:00:00 2001 From: mansurjisan Date: Sun, 11 Oct 2026 20:27:40 -0400 Subject: [PATCH 2/4] adcirc_met: comment attribution tags on the single-forecast change --- nos_utils/forcing/adcirc_met.py | 4 ++-- tests/test_adcirc_met.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/nos_utils/forcing/adcirc_met.py b/nos_utils/forcing/adcirc_met.py index 36bdd31..6a4103b 100644 --- a/nos_utils/forcing/adcirc_met.py +++ b/nos_utils/forcing/adcirc_met.py @@ -641,7 +641,7 @@ def sfcf_valid_times(cycle: datetime, segment: str, start: Optional[datetime] = def _interp_bracket(cycle: datetime, valid: datetime) -> Optional[Tuple[int, int, float]]: - """(lo, hi, weight) forecast hours when ``valid`` falls between two 3-hourly GFS records, else None.""" + """(lo, hi, weight) forecast hours when ``valid`` falls between two 3-hourly GFS records, else None. MJ (10/11/26)""" fh = int((valid - cycle).total_seconds() // 3600) if fh <= 120 or fh % 3 == 0: return None @@ -677,7 +677,7 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix up to ``wait_s`` seconds in total. Output files appear under their final names only when all three are complete, 221 last (the ops skip test is on 221). Records are copied one at a time. In real time (within 12 h of ``cycle``) the cycle-time record first waits for f000 so getges does not fall back to f006 of the - previous cycle. + previous cycle. MJ (10/11/26) """ import netCDF4 as nc diff --git a/tests/test_adcirc_met.py b/tests/test_adcirc_met.py index 348f128..93f23bc 100644 --- a/tests/test_adcirc_met.py +++ b/tests/test_adcirc_met.py @@ -497,9 +497,9 @@ def test_fcst_is_hourly_with_interpolation_after_f120(tmp_path): p = _read(f["221"], "pressfc") assert p.shape[0] == 181 s = lambda h, k: _valid_seed(CYC, h) * 10 + k - for h in (0, 7, 120, 123, 126, 180): # native records are exact copies + for h in (0, 7, 120, 123, 126, 180): # native records are exact copies. MJ (10/11/26) assert p[h][0, 0] == s(h, 2) - # 121 and 122 are 1/3 and 2/3 of the way from f120 to f123; the seeds step 10 per hour + # 121 and 122 are 1/3 and 2/3 of the way from f120 to f123; the seeds step 10 per hour. MJ (10/11/26) assert p[121][0, 0] == pytest.approx(s(120, 2) + 10.0) assert p[122][0, 0] == pytest.approx(s(120, 2) + 20.0) q = _read(f["222"], "ugrd10m") @@ -520,7 +520,7 @@ def test_fcst_interpolation_keeps_fill_values(tmp_path): def test_forecast_waits_then_fails_clearly(tmp_path): tank = tmp_path / "g" - _fcst_tank(tank, last=177) # f180 is missing + _fcst_tank(tank, last=177) # f180 is missing. MJ (10/11/26) clock = {"t": 0.0} naps = [] @@ -544,7 +544,7 @@ def sleep_and_land(s): def test_current_cycle_f000_is_not_waited_for(tmp_path): - # fcst record 0 (valid == cycle) goes through getges, so it never waits + # fcst record 0 (valid == cycle) goes through getges, so it never waits. MJ (10/11/26) _fill_tank(tmp_path / "g", [(CYC - timedelta(hours=6), 6)] + [(CYC, h) for h in range(1, 121)] + [(CYC, h) for h in range(123, 181, 3)]) files = am.build_sfcf_forcing(tmp_path / "g", CYC, "fcst", tmp_path / "o", wait_s=0, From 7628ffe2451ec919e8f30cc545c0412987305fae Mon Sep 17 00:00:00 2001 From: mansurjisan Date: Sun, 11 Oct 2026 20:32:13 -0400 Subject: [PATCH 3/4] adcirc_met: check the later bracket file's grid; manifest lists both bracket files of interpolated records --- nos_utils/forcing/adcirc_met.py | 4 +++- tests/test_adcirc_met.py | 3 +++ 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/nos_utils/forcing/adcirc_met.py b/nos_utils/forcing/adcirc_met.py index 6a4103b..1dfeee1 100644 --- a/nos_utils/forcing/adcirc_met.py +++ b/nos_utils/forcing/adcirc_met.py @@ -707,7 +707,7 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix src = sfcf_path(comin_gfs, cycle, int((valid - cycle).total_seconds() // 3600)) if not _wait_for(src, deadline, poll_s, sleep, clock): raise FileNotAvailableError("GFS file did not appear before the wait expired: {}".format(src)) - used.append((valid, src)) + used.extend((valid, p) for p in (srcs if br else [src])) # interpolated records list both brackets. MJ (10/11/26) log.info("record %d valid %s <- %s", rec, valid, src) with nc.Dataset(str(src)) as ds: ds.set_auto_maskandscale(False) @@ -719,6 +719,8 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix try: if hi is not None: hi.set_auto_maskandscale(False) + if hi.variables["lat"].shape != lat.shape: + raise MetForcingError("{}: grid differs from {}".format(srcs[1], src)) for n, names in SFC_OUTPUTS: o = outs[n] if lat.shape != (len(o.dimensions["grid_yt"]), len(o.dimensions["grid_xt"])): # ncecat needs equal shapes. MJ (10/06/26) diff --git a/tests/test_adcirc_met.py b/tests/test_adcirc_met.py index 93f23bc..d04698e 100644 --- a/tests/test_adcirc_met.py +++ b/tests/test_adcirc_met.py @@ -505,6 +505,9 @@ def test_fcst_is_hourly_with_interpolation_after_f120(tmp_path): q = _read(f["222"], "ugrd10m") assert q[121][1, 2] == pytest.approx(s(120, 0) + 10.0 + 6) assert _read(f["225"], "icec").shape[0] == 181 + man = json.loads((tmp_path / "o" / "stofs_2d_glo_fcst.sfcf_manifest.json").read_text())["files"] + assert len(man) == 181 + 40 and [e["source_location"] for e in man[121:123]] == [ + str(am.sfcf_path(tank, CYC, 120)), str(am.sfcf_path(tank, CYC, 123))] def test_fcst_interpolation_keeps_fill_values(tmp_path): From 460ee0056d1fcecdb2b16d11b4637c0a7e6e8b21 Mon Sep 17 00:00:00 2001 From: mansurjisan Date: Sun, 11 Oct 2026 21:22:01 -0400 Subject: [PATCH 4/4] adcirc_met: bracket grid check compares lat/lon values, not only shape --- nos_utils/forcing/adcirc_met.py | 2 +- tests/test_adcirc_met.py | 11 +++++++++++ 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/nos_utils/forcing/adcirc_met.py b/nos_utils/forcing/adcirc_met.py index 1dfeee1..825cd23 100644 --- a/nos_utils/forcing/adcirc_met.py +++ b/nos_utils/forcing/adcirc_met.py @@ -719,7 +719,7 @@ def build_sfcf_forcing(comin_gfs, cycle: datetime, segment: str, out_dir, prefix try: if hi is not None: hi.set_auto_maskandscale(False) - if hi.variables["lat"].shape != lat.shape: + if not (np.array_equal(hi.variables["lat"][:], lat) and np.array_equal(hi.variables["lon"][:], lon)): raise MetForcingError("{}: grid differs from {}".format(srcs[1], src)) for n, names in SFC_OUTPUTS: o = outs[n] diff --git a/tests/test_adcirc_met.py b/tests/test_adcirc_met.py index d04698e..022595a 100644 --- a/tests/test_adcirc_met.py +++ b/tests/test_adcirc_met.py @@ -521,6 +521,17 @@ def test_fcst_interpolation_keeps_fill_values(tmp_path): assert i[121][0, 0] == np.float32(9.99e20) and i[120][0, 0] < 1e6 and i[121][0, 1] < 1e6 +def test_fcst_rejects_a_bracket_on_another_grid(tmp_path): + import netCDF4 as nc + tank = tmp_path / "g" + _fcst_tank(tank) + with nc.Dataset(str(am.sfcf_path(tank, CYC, 123)), "a") as d: + d.variables["lon"][:] = d.variables["lon"][:] + 0.5 + with pytest.raises(am.MetForcingError, match="grid differs"): + am.build_sfcf_forcing(tank, CYC, "fcst", tmp_path / "o") + assert list((tmp_path / "o").glob("*.nc*")) == [] + + def test_forecast_waits_then_fails_clearly(tmp_path): tank = tmp_path / "g" _fcst_tank(tank, last=177) # f180 is missing. MJ (10/11/26)