Creating Synthetic Observations with perfect_model_obs

Sample a model state as if it were the real ocean - the first step in an Observing System Simulation Experiment (OSSE)!

A typical synthetic-observation workflow consists of five main steps:

Understand why (and when) synthetic observations are useful.
Build DART’s MOM6 executables.
Design an observing network — the locations from Tutorial 1, plus some random ones.
Run perfect_model_obs through dartobsgen.
Compare synthetic and real observations.

This is Part 2 of the DART tutorial series: 1. Creating Observations · 2. Synthetic Observations · 3. Cycling DART–CESM

Missed Tutorial 1?

You can run this notebook standalone:

Point REAL_OBS_DIR at the staged workshop copy, <CROC_DART_OBS>/real, in Step 3.1, or
Skip the harvest entirely — Step 3.1 falls back to a regular-grid network if it finds no Tutorial 1 output.

SECTION 1: Why synthetic observations?¶

In an OSSE you treat one model run as the “truth”, sample it where (and how accurately) your instruments would sample the real ocean, and assimilate those samples into a different ensemble. Because you know the truth exactly, you can measure precisely how much of it the assimilation recovers.

OSSEs answer questions like:

Is my DA plumbing working at all? (If you can’t recover a known truth, real obs won’t help.)
Where should new instruments go to constrain the circulation feature I care about?
How does observation accuracy trade off against observation count?

perfect_model_obs is DART’s tool for the sampling step: it reads a model state, applies the same forward operators the filter would use, and (optionally) perturbs each value with noise drawn from the observation error variance you assign.

SECTION 2: Build DART for MOM6¶

Step 2.1: Compile perfect_model_obs¶

You cloned DART in Tutorial 1 (“Before you start”). Now build the MOM6 executables on Derecho — load the same compiler and netCDF modules you use for CESM, then:

cd <DART_SRC>/models/MOM6/work
./quickbuild.sh

This compiles perfect_model_obs, filter, and the supporting utilities into the work directory. See the DART getting-started docs for Derecho-specific build settings (mkmf.template).

Step 2.2: What perfect_model_obs needs¶

dartobsgen’s PerfectModelSource drives perfect_model_obs for you, one assimilation window at a time. It needs three things in the DART work directory:

the compiled perfect_model_obs executable,
an input.nml containing a &perfect_model_obs_nml block (the one shipped in models/MOM6/work is a good starting point), and
a MOM6 state file to sample — the “truth”. Edit input.nml so the model reads the staged panama restart, <CROC_DART_OBS>/panama1_restart.nc (copy it into the work directory first).

For each window, PerfectModelSource writes a template obs_seq.in holding your network (locations, types, times, error variances — with placeholder values of 0.0), patches the namelist, runs the executable in an isolated subdirectory under <DART_SRC>/models/MOM6/work/windows/, and collects the resulting obs_seq.out.

Step 2.3: Parameters and sanity check¶

The shared series parameters, then a quick check that the pieces from Steps 2.1–2.2 are in place.

# --- CROCODILE DART tutorial series parameters (same cell in all 3 notebooks) ---
from pathlib import Path
import datetime

WORKDIR = Path("<DART_WORKDIR>")        # your scratch working directory on Derecho

START = datetime.datetime(2013, 4, 1)   # must match RUN_STARTDATE in Tutorial 3
END   = datetime.datetime(2013, 4, 4)   # 3 days -> 12 six-hour assimilation windows
FREQ  = datetime.timedelta(hours=6)

# Bounding box: the panama1 domain (lon 278-281E, lat 7-10N) padded by ~2 degrees
LAT_MIN, LAT_MAX = 5.0, 12.0
LON_MIN, LON_MAX = -84.0, -77.0

OBS_TYPES = ["ARGO_TEMPERATURE", "ARGO_SALINITY"]

REAL_OBS_DIR      = WORKDIR / "obs" / "real"       # Tutorial 1 output
SYNTHETIC_OBS_DIR = WORKDIR / "obs" / "synthetic"  # Tutorial 2 output

DART_WORK_DIR = Path("<DART_SRC>") / "models" / "MOM6" / "work"

for required in ["perfect_model_obs", "input.nml"]:
    f = DART_WORK_DIR / required
    status = "OK     " if f.exists() else "MISSING"
    print(f"{status} {f}")

SECTION 3: Design the observation network¶

Step 3.1: Harvest the Tutorial 1 locations¶

A realistic OSSE samples the truth where instruments actually were. We read the first Tutorial 1 window with pyDARTdiags and turn its unique locations into ObsNetworkEntry objects — one per location per observation type.

One thing changes compared to Tutorial 1: the error variance is now yours to choose. Real converters carry the instrument error with the data; in an OSSE the “instrument” is imaginary, so you decide how good it is. We use (0.2 °C)² for temperature and (0.1 PSU)² for salinity — typical Argo-like values.

import numpy as np
import pydartdiags.obs_sequence.obs_sequence as obsq
from dartobsgen import ObsNetworkEntry

OBS_ERR_VAR = {
    "ARGO_TEMPERATURE": 0.04,  # (0.2 degC)^2
    "ARGO_SALINITY":    0.01,  # (0.1 PSU)^2
}
PROFILE_DEPTHS = [10.0, 50.0, 100.0, 200.0, 500.0, 1000.0]  # Argo-like profile depths (m)

nb1_files = sorted(REAL_OBS_DIR.glob("obs_seq.*.out"))
# Missed Tutorial 1? Use the staged workshop copy instead:
# nb1_files = sorted(Path("<CROC_DART_OBS>/real").glob("obs_seq.*.out"))

network = []
if nb1_files:
    real = obsq.ObsSequence(str(nb1_files[0]))
    locations = (real.df[["longitude", "latitude", "vertical"]]
                 .drop_duplicates()
                 .reset_index(drop=True))
    print(f"Harvested {len(locations)} unique locations from {nb1_files[0].name}")
    for _, row in locations.iterrows():
        lon = row.longitude if row.longitude <= 180 else row.longitude - 360
        for obs_type in OBS_TYPES:
            network.append(ObsNetworkEntry(
                obs_type=obs_type,
                lat=float(row.latitude),
                lon=float(lon),
                vertical=float(row.vertical),
                vert_unit="height (m)",
                obs_err_var=OBS_ERR_VAR[obs_type],
            ))
else:
    # Fallback: no Tutorial 1 output found -- build a coarse regular-grid network.
    print("No Tutorial 1 output found; building a 1-degree grid network instead.")
    for lat in np.arange(LAT_MIN + 0.5, LAT_MAX, 1.0):
        for lon in np.arange(LON_MIN + 0.5, LON_MAX, 1.0):
            for depth in PROFILE_DEPTHS:
                for obs_type in OBS_TYPES:
                    network.append(ObsNetworkEntry(
                        obs_type=obs_type,
                        lat=float(lat), lon=float(lon),
                        vertical=depth, vert_unit="height (m)",
                        obs_err_var=OBS_ERR_VAR[obs_type],
                    ))

n_harvested = len(network)
print(f"Network so far: {n_harvested} observations")

Step 3.2: Add random locations¶

Now the part you can’t do with real data: invent instruments. We add 20 random profile locations inside the domain, each sampling the standard depths. The seeded random generator means everyone in the workshop gets the same “random” network — change the seed and your network is yours alone.

rng = np.random.default_rng(seed=42)   # fixed seed: reproducible "random" network
N_RANDOM = 20

random_lats = rng.uniform(LAT_MIN, LAT_MAX, N_RANDOM)
random_lons = rng.uniform(LON_MIN, LON_MAX, N_RANDOM)

for lat, lon in zip(random_lats, random_lons):
    for depth in PROFILE_DEPTHS:
        for obs_type in OBS_TYPES:
            network.append(ObsNetworkEntry(
                obs_type=obs_type,
                lat=float(lat), lon=float(lon),
                vertical=depth, vert_unit="height (m)",
                obs_err_var=OBS_ERR_VAR[obs_type],
            ))

print(f"Added {len(network) - n_harvested} obs at {N_RANDOM} random profile locations "
      f"({len(network)} total)")

Step 3.3: Map and save the network¶

Plot the two populations, and save the network to CSV — a compact, shareable record of your observing-system design.

import pandas as pd
import matplotlib.pyplot as plt

net_df = pd.DataFrame([{
    "obs_type": e.obs_type, "lat": e.lat, "lon": e.lon,
    "vertical": e.vertical, "obs_err_var": e.obs_err_var,
} for e in network])

WORKDIR.mkdir(parents=True, exist_ok=True)
net_df.to_csv(WORKDIR / "obs_network.csv", index=False)

harvested = net_df.iloc[:n_harvested]
random_part = net_df.iloc[n_harvested:]

fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(harvested["lon"], harvested["lat"], s=18, color="tab:blue",
           label=f"Tutorial 1 locations ({len(harvested)})")
ax.scatter(random_part["lon"], random_part["lat"], s=30, color="tab:orange",
           marker="^", label=f"random additions ({len(random_part)})")
ax.set_xlim(LON_MIN, LON_MAX)
ax.set_ylim(LAT_MIN, LAT_MAX)
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")
ax.set_title("Synthetic observing network")
ax.legend()
plt.show()

Hint

random_lats = rng.uniform(8.0, 9.0, N_RANDOM)
random_lons = rng.uniform(-80.0, -79.0, N_RANDOM)

Dense local sampling constrains the feature you target, but the rest of the domain is corrected only through ensemble covariances — and in Tutorial 3, localization deliberately limits how far one observation can reach.

SECTION 4: Generate the synthetic observations¶

Step 4.1: Run perfect_model_obs through dartobsgen¶

Same ObsGenConfig pattern as Tutorial 1 — only the source changes. That symmetry is the point: to the rest of the pipeline, an observation source is an observation source.

from dartobsgen import ObsGenConfig, PerfectModelSource, generate_obs_sequences

config = ObsGenConfig(
    start=START, end=END,
    lat_min=LAT_MIN, lat_max=LAT_MAX,
    lon_min=LON_MIN, lon_max=LON_MAX,
    obs_types=OBS_TYPES,
    assimilation_frequency=FREQ,
    output_dir=SYNTHETIC_OBS_DIR,
)

source = PerfectModelSource(
    dart_work_dir=str(DART_WORK_DIR),
    obs_network=network,
)

SYNTHETIC_OBS_DIR.mkdir(parents=True, exist_ok=True)

# max_workers=1 runs windows sequentially; set to None to run them in parallel.
written = generate_obs_sequences(config, source, max_workers=1)

print(f"{len(written)} obs_seq file(s) written to {SYNTHETIC_OBS_DIR.name}/")
for p in written:
    print("  ", Path(p).name)

Step 4.2: Look inside¶

The template values (0.0) have been replaced by the model sampled at each location, plus noise drawn from your obs_err_var.

syn = obsq.ObsSequence(str(Path(written[0])))

print(f"{Path(written[0]).name}: {len(syn.df)} observations")
syn.df[["type", "longitude", "latitude", "vertical", "observation", "obs_err_var"]].head(8)

SECTION 5: Real vs. synthetic, side by side¶

For the locations harvested from Tutorial 1 we now have two sets of values: what the real ocean said (CrocoLake) and what the model + noise says (perfect_model_obs). Comparing them makes “perfect model plus noise” concrete — and any systematic offset you see is a preview of the model bias the assimilation in Tutorial 3 will fight.

real_T = real.df[real.df["type"] == "ARGO_TEMPERATURE"]["observation"]
syn_T  = syn.df[syn.df["type"] == "ARGO_TEMPERATURE"]["observation"]

fig, ax = plt.subplots(figsize=(6, 4))
ax.hist(real_T, bins=30, alpha=0.6, label=f"real, CrocoLake ({len(real_T)})")
ax.hist(syn_T, bins=30, alpha=0.6, label=f"synthetic, perfect_model_obs ({len(syn_T)})")
ax.set_xlabel("Temperature (degC)")
ax.set_ylabel("Count")
ax.set_title("Real vs. synthetic temperature observations, first window")
ax.legend()
plt.show()

Recap¶

Your takeaway artifacts: obs/synthetic/ and obs_network.csv under your WORKDIR.

Where to go from here?¶

Tutorial 3: Cycling DART–CESM — assimilate the real observations, then swap in obs/synthetic/ to run it as an OSSE.
The script version of this workflow: mom6_perfect_model.py in the dartobsgen repository.
DART perfect_model_obs documentation.