The observation database

The high-fidelity observations \(\mathbf{y}_k\) used by the analysis phase (see Data Assimilation) are not passed to CONES directly: they are read once at start-up from a netCDF file, whose path is set by obsFile in The conesDict reference. This page describes the file’s expected structure and shows two ways to build one.

Expected structure

The file is read by conesStaticObservation and the netCDF helpers in conesToolBox.conesFunctions (cones_count_netcdf_obs(), cones_get_netcdf_obs_coord()), which impose the following schema:

Dimension	Meaning
obs	Number of observation points. Required, and must be named exactly obs — cones_count_netcdf_obs() reads it directly as ds.sizes["obs"].
time	Number of time samples available for each observed variable.

Variable	Meaning
x, y, z	Coordinates of each observation point, shaped (obs,). Static in time: an observation is a fixed point in space.
one per obsVar letter (e.g. u, v, w)	The observed value of that variable at each point and time, shaped (time, obs). The variable names must match obsVar in The conesDict reference exactly.
the _std companion of each of the above (e.g. u_std)	The standard deviation \(\sigma_{i,k}\) of that observation, same shape. This is what fills the observation covariance \(\mathbf{R}_k = \mathrm{diag}(\sigma_{i,k}^2)\) (see The conesDict reference) and the perturbation drawn by setObs(). Every observed variable needs one — get_vars_from_netcdf() looks it up unconditionally and raises a KeyError if it is missing.

Only the variables named in obsVar (plus their _std companion) are actually read at run time; extra data variables in the file are ignored.

A CF-decoding gotcha

If you set a dtype entry directly in a variable’s .attrs (rather than its .encoding) — for instance to force time to be interpreted as timedelta64[ns] — some xarray versions refuse to CF-decode the file on open. Both examples below set ds.time.attrs['dtype'] for exactly this reason; cones_open_netcdf_dataset() is CONES’s read-side workaround, moving any such dtype attribute into encoding before decoding. Reuse it if you write your own reader.

Example 1: synthetic observations

utils/generate_observation_file.py builds a database of 5 point velocity observations, constant in time — this is what the cavity tutorial uses. Run it with the relative standard deviation to apply to each value (5% here):

python3 utils/generate_observation_file.py 0.05

import xarray as xr
import numpy as np
import sys

if len(sys.argv) != 2:
    print("Usage: python3 generate_observation_file.py <relative_std>")
    print("  relative_std: fraction of each observed value used as its standard deviation, e.g. 0.05")
    sys.exit(1)
relative_std = float(sys.argv[1])

# Define number of observations and coordinates
num_obs = 5
n_time = 101
time = np.linspace(0, 100, n_time)

obs_ids = np.arange(num_obs)

x = np.array([0.0425, 0.0475, 0.0525, 0.0575, 0.0675])
y = np.ones(num_obs)*0.0975
z = np.ones(num_obs)*0.005

# Define observations values and std

u_static = np.array([0.753300, 0.799275, 0.824227, 0.851823, 0.874961])
v_static = np.array([0.0455249, 0.0403204, 0.0357992, 0.0322402, 0.0263673])
w_static = np.zeros(num_obs)

u_std_static = u_static*relative_std
v_std_static = v_static*relative_std
w_std_static = w_static*relative_std

# Broadcasting
u = np.zeros((n_time, num_obs))
v = np.zeros((n_time, num_obs))
w = np.zeros((n_time, num_obs))
u_std = np.zeros((n_time, num_obs))
v_std = np.zeros((n_time, num_obs))
w_std = np.zeros((n_time, num_obs))


for t_idx, t in enumerate(time):
    u[t_idx, :] = u_static
    v[t_idx, :] = v_static
    w[t_idx, :] = w_static
    u_std[t_idx, :] = u_std_static
    v_std[t_idx, :] = v_std_static
    w_std[t_idx, :] = w_std_static


# Assemble dataset

ds = xr.Dataset(
    data_vars={
        "u": (["time", "obs"], u),
        "v": (["time", "obs"], v),
        "w": (["time", "obs"], w),
        "u_std": (["time", "obs"], u_std),
        "v_std": (["time", "obs"], v_std),
        "w_std": (["time", "obs"], w_std),
        "x": (["obs"], x),
        "y": (["obs"], y),
        "z": (["obs"], z),
        # "obs_type":(["obs"], obs_types)
    },
    coords={
        "time": time,
        "obs": obs_ids,
        }
)

# ds = ds_static.expand_dims(time=t)
ds.time.attrs['units'] = 'nanoseconds'
ds.x.attrs = {'axis': 'X', 'units': 'm', 'standard_name': 'projection_x_coordinate'}
ds.y.attrs = {'axis': 'Y', 'units': 'm', 'standard_name': 'projection_y_coordinate'}
ds.z.attrs = {'axis': 'Z', 'units': 'm', 'standard_name': 'projection_z_coordinate'}

ds.attrs = {'feature_type': 'point', 'Conventions': 'CF-1.6'}
ds.time.attrs['dtype'] = 'timedelta64[ns]'
# Add Metadata
# ds.u.attrs = {"units": "m/s", "standard_name": "u"}
# ds.u_std.attrs = {"units": "m/s", "descriptions": "Standard deviation of u"}
# ds.v.attrs = {"units": "m/s", "standard_name": "v"}
# ds.v_std.attrs = {"units": "m/s", "descriptions": "Standard deviation of v"}
# ds.w.attrs = {"units": "m/s", "standard_name": "w"}
# ds.w_std.attrs = {"units": "m/s", "descriptions": "Standard deviation of w"}
# ds.x.attrs = {"units": "m"}
# ds.y.attrs = {"units": "m"}
# ds.z.attrs = {"units": "m"}
# ds.attrs = {"description": "Database of point observations over time"}
ds.to_netcdf("observation_database.nc")
print("File saved!")
print(ds)

Adapt the coordinates and observed values to your own case; keep every data variable shaped (time, obs) (or (obs,) for x/y/z) and make sure every observed variable has its _std companion.

Example 2: observations from OpenFOAM probes

If you already have a reference/high-fidelity run, utils/probes_to_observation_netcdf.py converts the output of OpenFOAM’s probes function object (as configured in system/probes) into an observation database, reading the probe coordinates from the file’s own header comments:

python3 utils/probes_to_observation_netcdf.py postProcessing/probes/0/U

It currently hardcodes num_obs = 5 and a flat 5% relative standard deviation — open the script and adjust both to match your system/probes setup before relying on it.