Read EXF binary files in review and animate utils instead of ERA5 NetCDF

  Adds load_exf_binaries() to common.py which memory-maps the written .bin
  files and reconstructs coordinates from NetCDF metadata only. Updates
  review_exf_conditions.py and animate_exf_conditions.py to use this so QC
  checks and animations reflect the actual data that MITgcm will consume.

Author: fluidnumericsJoe

spectre_utils/animate_exf_conditions.py

@@ -7,8 +7,6 @@
 import matplotlib.pyplot as plt
 from matplotlib.animation import FFMpegWriter
 from datetime import datetime
-from metpy.calc import specific_humidity_from_dewpoint
-from metpy.units import units
 
 from spectre_utils import common
 
@@ -106,15 +104,9 @@ def main():
     t1 = datetime.strptime(config["domain"]["time"]["start"], "%Y-%m-%d")
     t2 = datetime.strptime(config["domain"]["time"]["end"], "%Y-%m-%d")
 
-    print("Loading dataset...")
-    ds = common.load_atm_dataset(working_directory, prefix, years, atm_vars, t1, t2)
-
-    if "sp" in ds and "d2m" in ds:
-        d2m_celsius = ds["d2m"] - 273.15
-        ds["aqh"] = specific_humidity_from_dewpoint(ds["sp"] * units.Pa, d2m_celsius * units.degC)
-
     animations_dir = os.path.join(simulation_directory, "animations", "atmosphere")
     os.makedirs(animations_dir, exist_ok=True)
+    simulation_input_dir = os.path.join(simulation_directory, "input")
 
     # Collect all variable names to animate (configured + computed), deduplicated
     seen = set()
@@ -130,6 +122,11 @@ def main():
             seen.add(n)
             to_animate.append(n)
 
+    print("Loading EXF binary files...")
+    ds = common.load_exf_binaries(
+        simulation_input_dir, to_animate, working_directory, prefix, years, atm_vars, t1, t2
+    )
+
     for name in to_animate:
         if name not in ds:
             print(f"  Skipping {name} (not in dataset)")

spectre_utils/common.py

@@ -68,6 +68,54 @@ def from_copernicus(dataset_id, variables, start_date, end_date, min_long, max_l
     return ds_full, grid
 
 
+def load_exf_binaries(simulation_input_dir, var_names, working_directory, prefix, years, atm_vars, t1, t2):
+    """Load EXF binary files into a labelled xarray Dataset.
+
+    Coordinates (time, latitude, longitude) are extracted from the ERA5 NetCDF
+    source files (metadata only — no data are read from them).  The binary data
+    are memory-mapped so only the portions actually accessed are read into RAM.
+    """
+    import xarray as xr
+    import numpy as np
+    import os
+
+    # Obtain grid coordinates from the first available NetCDF variable.
+    coord_ds = None
+    for var in atm_vars:
+        mitgcm_name = var["mitgcm_name"]
+        files = [f"{working_directory}/{prefix}_{mitgcm_name}_{year}.nc" for year in years]
+        if all(os.path.exists(f) for f in files):
+            coord_ds = xr.open_mfdataset(files, combine="by_coords").sel(
+                valid_time=slice(t1, t2)
+            )
+            break
+    if coord_ds is None:
+        raise FileNotFoundError("No ERA5 NetCDF files found for coordinate extraction.")
+
+    times = coord_ds["valid_time"].values
+    lats = coord_ds["latitude"].values
+    lons = coord_ds["longitude"].values
+    coord_ds.close()
+
+    nt, ny, nx = len(times), len(lats), len(lons)
+
+    data_vars = {}
+    for name in var_names:
+        bin_path = os.path.join(simulation_input_dir, f"{name}.bin")
+        if not os.path.exists(bin_path):
+            print(f"Warning: binary file not found, skipping: {bin_path}", file=sys.stderr)
+            continue
+        arr = np.memmap(bin_path, dtype=">f4", mode="r", shape=(nt, ny, nx))
+        data_vars[name] = xr.DataArray(
+            arr,
+            dims=["valid_time", "latitude", "longitude"],
+            coords={"valid_time": times, "latitude": lats, "longitude": lons},
+            name=name,
+        )
+
+    return xr.Dataset(data_vars)
+
+
 def load_atm_dataset(working_directory, prefix, years, atm_vars, t1, t2):
     """Load ERA5 atmospheric variables per MITgcm name, rename, and apply optional scale factors.
 

spectre_utils/review_exf_conditions.py

@@ -1,10 +1,9 @@
 """
 QC review of EXF atmospheric forcing fields.
 
-Loads the ERA5 xarray dataset (same config-driven pipeline as
-mk_exf_conditions.py), computes derived variables, runs physical
-range / consistency / temporal checks, verifies binary file sizes, and
-writes a Markdown report plus three diagnostic figures to:
+Loads the written binary forcing files, runs physical range / consistency /
+temporal checks, verifies binary file sizes, and writes a Markdown report
+plus three diagnostic figures to:
 
     <simulation_directory>/review/atmosphere/
         report.md
@@ -23,9 +22,6 @@
 import matplotlib.pyplot as plt
 from datetime import datetime, timezone
 
-from metpy.calc import specific_humidity_from_dewpoint
-from metpy.units import units as metpy_units
-
 from spectre_utils import common
 
 # ---------------------------------------------------------------------------
@@ -216,22 +212,8 @@ def main():
     os.makedirs(review_dir, exist_ok=True)
 
     # ------------------------------------------------------------------
-    # Load dataset
+    # Build ordered list of all mitgcm names (configured + computed)
     # ------------------------------------------------------------------
-    print("Loading ERA5 dataset...")
-    ds = common.load_atm_dataset(working_directory, prefix, years, atm_vars, t1, t2)
-
-    if "sp" in ds and "d2m" in ds:
-        d2m_celsius = ds["d2m"] - 273.15
-        ds["aqh"] = specific_humidity_from_dewpoint(
-            ds["sp"] * metpy_units.Pa, d2m_celsius * metpy_units.degC
-        )
-
-    nt = ds.sizes["valid_time"]
-    ny = ds.sizes.get("latitude")
-    nx = ds.sizes.get("longitude")
-
-    # Ordered list of all mitgcm names (configured + computed), deduplicated
     seen: set[str] = set()
     all_names: list[str] = []
     for var in atm_vars:
@@ -245,6 +227,18 @@ def main():
             seen.add(n)
             all_names.append(n)
 
+    # ------------------------------------------------------------------
+    # Load binary forcing files
+    # ------------------------------------------------------------------
+    print("Loading EXF binary files...")
+    ds = common.load_exf_binaries(
+        simulation_input_dir, all_names, working_directory, prefix, years, atm_vars, t1, t2
+    )
+
+    nt = ds.sizes["valid_time"]
+    ny = ds.sizes.get("latitude")
+    nx = ds.sizes.get("longitude")
+
     # ------------------------------------------------------------------
     # Diagnostic figures
     # ------------------------------------------------------------------