Quality Control Checks - xOPR Documentation

This notebook demonstrates the xopr.qc module, which provides standardized quality control checks for radar datasets. Each check produces a per-trace boolean mask that can be used to filter or visualize data quality.

The built-in checks are:

ice_thickness_threshold — flags traces where ice thickness is below a minimum
snr_bed_pick — flags traces where the bed return is weak relative to the noise floor
heading_change — flags traces with rapid aircraft heading changes (if the season has no Heading, one is reconstructed from the GPS track)
minimum_agl — flags traces where the platform is too close to the surface

We’ll load a flight line, run all checks, and visualize which checks pass or fail along the radargram.

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

import xopr
from xopr.qc import run_qc

opr = xopr.OPRConnection(cache_dir="radar_cache")

Load an example flight¶

We’ll start by querying the first few frames of an example flight.

collection, segment = "2019_Antarctica_GV", "20191105_01"

stac_items = opr.query_frames(collections=collection, segment_paths=[segment], max_items=5)
stac_items

We’ll load all frames from the results, merge them into one dataset, and resample to uniform 2-second spacing.

Heading is an angle, so a plain mean is wrong wherever it wraps at ±π. We drop it before resampling and let the QC reconstruct it from the resampled GPS track instead.

frames = opr.load_frames(stac_items)
flight_line = xopr.merge_frames(frames)
flight_line = flight_line.drop_vars("Heading").resample(slow_time='2s').mean()
flight_line.xopr

Run QC checks¶

run_qc() automatically loads layer picks (standard:surface and standard:bottom) when they are not already in the dataset — just pass the opr connection. It also fills in Heading from the GPS positions when a season did not record one (see below); pass "heading_change": {"source": "gps"} to force the reconstruction even when a measured heading exists.

Passing checks=None (the default) runs all registered checks with their default parameters. You can override parameters for specific checks by passing a dict:

# Run all default checks (with a custom ice thickness threshold)
qc_ds = run_qc(
    flight_line,
    opr=opr,
    checks={
        "ice_thickness_threshold": {"min_thickness_m": 500},
        "snr_bed_pick": {"min_snr_db": 5.0},
        "heading_change": {"max_deg_per_km": 5.0},
        "minimum_agl": {"min_agl_m": 100.0},
    },
)

# Summary
qc_vars = [v for v in qc_ds.data_vars if v.startswith("qc_")]
for v in qc_vars:
    n = int(qc_ds[v].sum())
    print(f"  {v}: {n}/{qc_ds.sizes['slow_time']} passed")

n_pass = int(qc_ds["qc"].sum())
n_total = qc_ds.sizes["slow_time"]
print(f"\nCombined: {n_pass}/{n_total} traces passed all checks ({100*n_pass/n_total:.1f}%)")

  qc_ice_thickness_threshold: 481/969 passed
  qc_snr_bed_pick: 944/969 passed
  qc_heading_change: 938/969 passed
  qc_minimum_agl: 969/969 passed

Combined: 430/969 traces passed all checks (44.4%)

Visualize the radargram with per-check QC flags¶

The top panel shows the radargram with layer picks overlaid and red shading on traces that fail the combined QC. The bottom panel shows a stacked plot of each individual check. Each check’s trace is high when the check passes and low when it fails.

fig, (ax_rg, ax_qc) = plt.subplots(
    2, 1, figsize=(15, 7),
    gridspec_kw={"height_ratios": [2, 1], "hspace": 0.05},
    sharex=True,
)

# --- Top panel: radargram with red shading ---
radargram = 10 * np.log10(np.abs(qc_ds["Data"]))
radargram.plot.imshow(x="slow_time", cmap="gray", ax=ax_rg, add_colorbar=False)
ax_rg.invert_yaxis()

if "standard:surface" in qc_ds:
    qc_ds["standard:surface"].plot(ax=ax_rg, x="slow_time", color="cyan", linewidth=0.5, label="Surface")
if "standard:bottom" in qc_ds:
    qc_ds["standard:bottom"].plot(ax=ax_rg, x="slow_time", color="yellow", linewidth=0.5, label="Bottom")

ax_rg.set_ylabel("TWTT [s]")

# Shade failing traces
fail_mask = ~qc_ds["qc"].values
if fail_mask.any():
    slow_times = qc_ds.slow_time.values
    half_dt = (slow_times[1] - slow_times[0]) / 2 if len(slow_times) > 1 else np.timedelta64(1, "s")
    diff = np.diff(fail_mask.astype(int))
    starts = np.where(np.concatenate(([fail_mask[0]], diff == 1)))[0]
    ends = np.where(np.concatenate((diff == -1, [fail_mask[-1]])))[0]
    for s, e in zip(starts, ends):
        ax_rg.axvspan(slow_times[s] - half_dt, slow_times[e] + half_dt, color="tab:red", alpha=0.25)

ax_rg.legend(loc="lower right")
ax_rg.set_title(f"{flight_line.attrs['collection']} — {flight_line.attrs['segment_path']}")
ax_rg.set_xlabel("")

# --- Bottom panel: per-check pass/fail ---
qc_vars = sorted(v for v in qc_ds.data_vars if v.startswith("qc_")) + ["qc"]
for idx, var in enumerate(qc_vars):
    (qc_ds[var]+(1.5*idx)).plot(ax=ax_qc, x="slow_time", linewidth=0.5, label=var)

ax_qc.set_yticks(np.arange(len(qc_vars))*1.5)
ax_qc.set_yticklabels(qc_vars, fontsize=8)

plt.show()

Flight track map¶

The flight track projected onto EPSG:3031 with coastline context. Traces that pass all QC checks are shown in green; traces filtered out by any check are shown in red.

import cartopy.crs as ccrs
import cartopy.feature as cfeature

proj = xopr.geometry.project_dataset(qc_ds, target_crs="EPSG:3031")
passed = proj["qc"].values

fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={"projection": ccrs.SouthPolarStereo()})
ax.add_feature(cfeature.COASTLINE, linewidth=0.5)
ax.add_feature(cfeature.OCEAN)

ax.scatter(proj["x"].values[~passed], proj["y"].values[~passed],
           c="tab:red", s=4, label="Failed QC", transform=ccrs.epsg(3031))
ax.scatter(proj["x"].values[passed], proj["y"].values[passed],
           c="tab:green", s=4, label="Passed QC", transform=ccrs.epsg(3031))

ax.legend(loc="upper right")
ax.set_title(f"{collection} — {segment}")
plt.show()

/home/runner/work/xopr/xopr/.venv/lib/python3.12/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_ocean.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

/home/runner/work/xopr/xopr/.venv/lib/python3.12/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Running a subset of checks¶

You can run only specific checks by passing a subset of the registry keys. You can also mix in custom check functions as callable keys:

# Run only ice thickness and SNR checks with stricter thresholds
strict_ds = run_qc(
    flight_line,
    opr=opr,
    checks={
        "ice_thickness_threshold": {"min_thickness_m": 1000},
        "snr_bed_pick": {"min_snr_db": 10.0},
    },
)

n_pass_strict = int(strict_ds["qc"].sum())
print(f"Strict thresholds: {n_pass_strict}/{n_total} traces passed ({100*n_pass_strict/n_total:.1f}%)")
print(f"Default thresholds: {n_pass}/{n_total} traces passed ({100*n_pass/n_total:.1f}%)")

Strict thresholds: 288/969 traces passed (29.7%)
Default thresholds: 430/969 traces passed (44.4%)

Seasons without a heading¶

Many older seasons (for example every Greenland P3 season before 2002, and the 2009–2012 DC8 and TO seasons) have no Heading, Roll, or Pitch in their CSARP files. For heading_change, run_qc reconstructs a heading from the GPS positions using xopr.add_heading. The dataset records where the heading came from in attrs["heading_source"].

The reconstruction is the course over ground, not the true platform heading: it omits the crab angle from crosswind. That offset is typically a few degrees and changes slowly, so it cancels in the rate of change that heading_change uses. Traces where the GPS track is unreliable (stationary platform, position dropouts, time gaps) get a NaN heading and fail the check.

old_items = opr.query_frames(collections="2012_Greenland_P3", max_items=2)
old_line = xopr.merge_frames(opr.load_frames(old_items))
print("Heading in file:", "Heading" in old_line)

old_qc = run_qc(old_line, checks={"heading_change": {"max_deg_per_km": 5.0}})
print("heading_source:", old_qc.attrs["heading_source"])
print(f"qc_heading_change: {int(old_qc['qc_heading_change'].sum())}/{old_qc.sizes['slow_time']} passed")

Heading in file: False
heading_source: gps
qc_heading_change: 3338/3910 passed

/home/runner/work/xopr/xopr/src/xopr/opr_tools.py:57: FutureWarning: In a future version of xarray the default value for data_vars will change from data_vars='all' to data_vars=None. This is likely to lead to different results when multiple datasets have matching variables with overlapping values. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set data_vars explicitly.
  merged_segment = xr.concat(segment_frames, dim='slow_time', combine_attrs=merge_dicts_no_conflicts).sortby('slow_time')

On a season that does have a measured heading we can compare the two. The difference is the crab angle: small, slowly varying, and irrelevant to the heading rate.

line = xopr.merge_frames(frames)
gps = xopr.add_heading(line, overwrite=True)
crab = np.degrees((gps["Heading"] - line["Heading"] + np.pi) % (2 * np.pi) - np.pi)

fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(12, 5), sharex=True)
ax0.plot(line.slow_time, np.degrees(line["Heading"]), lw=0.8, label="measured (INS)")
ax0.plot(gps.slow_time, np.degrees(gps["Heading"]), lw=0.8, label="reconstructed from GPS")
ax0.set_ylabel("heading [deg]"); ax0.legend()
ax1.plot(crab.slow_time, crab, lw=0.8)
ax1.set_ylabel("GPS − INS [deg]"); ax1.set_xlabel("slow_time")
plt.show()

Custom check functions¶

You can also pass a callable as a check key. The function must accept a dataset as its first argument and return a dataset with QC mask variables added (via apply_qc_mask).

Here we define a simple roll-angle check that flags traces where the platform roll exceeds 5 degrees:

Note that Roll (like Heading) is only present for seasons with an inertial navigation system; unlike heading, it cannot be reconstructed from GPS, so this check raises on older seasons.

from xopr.qc import apply_qc_mask

def roll_check(ds, max_roll_deg=5.0):
    """Flag traces where abs(Roll) exceeds a threshold."""
    roll_deg = np.degrees(ds["Roll"].values)
    mask = xr.DataArray(np.abs(roll_deg) <= max_roll_deg, dims="slow_time")
    return apply_qc_mask(ds, mask, "roll")

custom_ds = run_qc(
    flight_line,
    opr=opr,
    checks={
        "ice_thickness_threshold": {"min_thickness_m": 500},
        roll_check: {"max_roll_deg": 5.0},
    },
)

for v in sorted(v for v in custom_ds.data_vars if v.startswith("qc")):
    n = int(custom_ds[v].sum())
    print(f"  {v}: {n}/{custom_ds.sizes['slow_time']} passed")

  qc: 438/969 passed
  qc_ice_thickness_threshold: 481/969 passed
  qc_roll: 917/969 passed