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 workflow for generating synthetic observations consists of the following steps:

Design an observing network: we’ll use the locations from Tutorial 1, plus some random ones.
Run a “truth” or nature run of your model: MOM6 for this tutorial
Use dartobsgen with perfect_model_obs create the synthetic observations
Examine the synthetic observations you have created.

This is Part 2 of the DART tutorial series:
1. Working with Real Observations · 2. Creating 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 a model run as the “truth”. You sample your model run where your instruments would sample the real ocean, and add realistic errors such as instrument noise and representativeness error. Because you know the truth exactly, you can use OSSEs to estimate how new instruments and observing networks will impact your data assimilation.

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 assimilation would use, and perturbs each value with noise drawn from the observation error variance you assign.

SECTION 2: Set Up Your Experiment Parameters¶

dartobsgen’s PerfectModelSource drives perfect_model_obs for you, one assimilation window at a time. It needs three things in the DA_PROJECT_DIR 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 to sample, the “truth”.

For this tutorial, we’ll use a MOM6 run of the Hawaii domain as model truth, and the perfect_model_obs from /glade/work/hkershaw/croc26/DART/.

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

DA_PROJECT_DIR = Path("<DART_DA_PROJECT_DIR>")        # your DA project working directory on Derecho

START = datetime.datetime(2023, 6, 15)   # must match RUN_STARTDATE in Tutorial 3
END   = datetime.datetime(2023, 6, 18)   # 3 days -> 3 one-day assimilation windows, centered on midnight
FREQ  = datetime.timedelta(hours=24)

# Bounding box:
LAT_MIN, LAT_MAX = 20.0, 25.0
LON_MIN, LON_MAX = -160.0, -155.0

OBS_TYPES = ["ARGO_TEMPERATURE", "ARGO_SALINITY"]

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

PMO_RUN_DIR = Path("<DART_DA_PROJECT_DIR>") / "pmo_runs" / "windows"

For each window, PerfectModelSource writes a template obs_seq.in holding your network (locations, types, times, error variances), sets the DART namelist for each window, runs the executable in an isolated subdirectory under <DART_DA_PROJECT_DIR>/pmo_runs/windows/, and collects the resulting obs_seq.out.

SECTION 3: Design the observation network¶

Step 3.1: Harvest the Tutorial 1 locations¶

Often in an OSSE we want to sample the truth where instruments actually were. Let’s use the locations of the real observations from Tutorial 1. 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
}
# Argo-like profile depths (m), used as a fallback if no Tutorial 1 obs
PROFILE_DEPTHS = [10.0, 50.0, 100.0, 200.0, 500.0, 1000.0]  

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. You can change the seed to generate a different network.

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 the network¶

Plot the two networks, the locations harvested from Tutorial 1 and the random additions, together.

import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

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])

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

fig, ax = plt.subplots(figsize=(6, 5), subplot_kw={"projection": ccrs.PlateCarree()})
ax.scatter(harvested["lon"], harvested["lat"], s=18, color="tab:blue",
           label=f"Tutorial 1 locations ({len(harvested)})", transform=ccrs.PlateCarree())
ax.scatter(random_part["lon"], random_part["lat"], s=30, color="tab:orange",
           marker="^", label=f"random additions ({len(random_part)})", transform=ccrs.PlateCarree())

ax.coastlines(resolution="10m")
ax.add_feature(cfeature.LAND, facecolor="lightgray")
ax.set_extent([LON_MIN, LON_MAX, LAT_MIN, LAT_MAX], crs=ccrs.PlateCarree())
ax.gridlines(draw_labels=True)

ax.set_title("Synthetic observing network")
ax.legend()
plt.show()

SECTION 4: Generate The Synthetic Observations¶

Step 4.1: Configure `dartobsgen` For Our New Observation Network¶

ObsGenConfig is almost the same as Tutorial 1. We want the same time period, same bounding box, same observation types, same assimilation frequency, but we want our output written to a different directory: SYNTHETIC_OBS_DIR. The source for “observations” is now a MOM6 run rather than CrocoLake, so we use PerfectModelSource rather than CrocoLake.

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)

Step 4.3: Map every window¶

Loop over your synthetic observation sequences and map the ARGO_TEMPERATURE observations to see the observation value sampled from MOM6 over the time window.

temp_dfs = []
for p in written:
    obs_seq = obsq.ObsSequence(str(Path(p)))
    temp_dfs.append(obs_seq.df[obs_seq.df["type"] == "ARGO_TEMPERATURE"])

vmin = min(df["observation"].min() for df in temp_dfs)
vmax = max(df["observation"].max() for df in temp_dfs)

ncols = 4
nrows = -(-len(written) // ncols)  # ceil division
fig, axes = plt.subplots(nrows, ncols, figsize=(4 * ncols, 3.2 * nrows),
                          sharex=True, sharey=True)
axes = np.atleast_1d(axes).ravel()

for ax, p, df in zip(axes, written, temp_dfs):
    lon = np.where(df["longitude"] > 180, df["longitude"] - 360, df["longitude"])
    sc = ax.scatter(lon, df["latitude"], c=df["observation"],
                     cmap="viridis", vmin=vmin, vmax=vmax, s=14)
    ax.set_xlim(LON_MIN, LON_MAX)
    ax.set_ylim(LAT_MIN, LAT_MAX)
    ax.set_title(Path(p).name.removeprefix("obs_seq.").removesuffix(".out"), fontsize=8)

for ax in axes[len(written):]:
    ax.axis("off")

fig.colorbar(sc, ax=axes[:len(written)].tolist(), shrink=0.6,
             label="ARGO_TEMPERATURE (degC)")
fig.suptitle("perfect_model_obs temperature, every assimilation window")
plt.show()

Recap¶

Your takeaway artifact: A directory obs/synthetic/ocn_obs_seq containing synthetic observations under your DA_PROJECT_DIR.

Where to go from here?¶

Tutorial 3: Cycling DART–CESM: assimilate the real observations, then swap in obs/synthetic/ocn_obs_seq to run the experiment as an OSSE.
Take a look at the script version of this workflow: mom6_perfect_model.py in the dartobsgen repository.
Dive deeper into the DART perfect_model_obs documentation.