Consistency of ozone mixing ratio from limb, nadir, and in situ measurements

logo

Please note that this repository is used for development and review, so quality assessments should be considered work in progress until they are merged into the main branch

1.2.1. Consistency of ozone mixing ratio from limb, nadir, and in situ measurements#

Production date: 20-07-2025

Produced by: Chiara Volta (ENEA, Italy) and Sandro Calmanti (ENEA, Italy)

🌍 Use case: Support policy in detecting large ozone depletion events.#

❓ Quality assessment questions#

  • How well do limb and nadir satellite datasets agree with ozonesondes in monitoring stratospheric ozone mixing ratio in high-latitude regions?

  • Can differences between limb and nadir measurements affect the interpretation of ozone depletion events for policymaking?

Stratospheric ozone plays a vital role in protecting life on Earth from harmful solar ultraviolet radiation, and assessing its return to historical levels and its full recovery from ozone-depleting substances (ODS) is a critical objective [1]. Since the adoption of the 1987 Montreal Protocol, the international scientific community has been tasked with providing governments with regular updates on the state of the ozone layer and guidance for strengthening the Protocol’s provisions. In this context, the availability of high-quality ozone observations is essential for monitoring trends, evaluating the effectiveness of policy measures, and supporting robust scientific assessments.

The goal of the present assessment is to evaluate the consistency between ozone observations from limb- and nadir-viewing satellite sensors and in situ ozonesonde measurements of stratospheric ozone. Additionally, it aims to assess whether differences between limb and nadir data could influence the interpretation of ozone depletion events, with potential implication for the scientific advice that informs policy decisions.

A full description of ozone data from limb and nadir satellite sensors and derived datasets is provided in [2]. Ozonesonde data are freely available from the World Ozone and Ultraviolet Radiation Centre (WOUDC).

📢 Quality assessment statements#

These are the key outcomes of this assessment

  • Both limb- and nadir-viewing satellite datasets effectively capture the seasonal variability of stratospheric ozone mixing ratios at mid- and high-latitudes.

  • Satellite datasets demonstrate significantly higher agreement with ozonesondes at mid-to-high Northern latitudes than in the Antarctic region.

  • Despite the extended temporal coverage provided by successive satellite missions, inconsistencies between different versions of the same product and between products from similar sensor types (i.e., limb or nadir) limit their suitability for detecting long-term ozone trends in the stratospheric ozone layer.

  • The consistency between satellite and in situ measurements depends on both the ozone concentration in the stratosphere and latitude.

  • The tendency for satellites to overestimate ozone during depletion events at high Southern latitudes could lead to an underestimation of the severity of ozone depletion events.

📋 Methodology#

The analysis performed in this notebook focuses on comparing ozone mixing ratios measured by limb and nadir satellite sensors with in situ ozonesonde observations from two stations: Edmonton (EDT - Canada) and Marambio (MBI - Antarctica). These stations were selected due to the availability of long-term time series and can be considered as representative of mid- to high-latitude regions in the Northern and Southern Hemispheres, respectively. To evaluate consistency across datasets, time series from satellite and ozonesonde measurements were compared, and statistical indicators such as Root Mean Square Error (RMSE), Mean Bias Error (MBE), Mean Absolute Percentage Error (MAPE) and Standard Error (SE) were calculated.
The analysis excludes merged satellite datasets (e.g., MSR and AMZM) to avoid potential inconsistencies related to cross-sensor averaging and algorithmic harmonization.

1. Set up the code and choose the satellite data to use

  • Import libraries

  • Define data requests and other parameters

2. Satellite data caching functions and retrieval

  • Define required functions

  • Data retrieval

3. Download external ozonesonde data

  • Define function to extract ozonesonde data

  • Download ozonesonde data by station

4. Compute statistics

  • Define statistical function

5. Display and discuss results

  • Display results

  • Discussion

📈 Analysis and results#

1. Set up the code and choose the satellite data to use#

Import libraries#

Besides the standard libraries used to manage and analyse multidimensional arrays, the C3S EQC custom function c3s_eqc_automatic_quality_control is imported to download satellite data and calculate their spatially weighted means. The requests module is imported to download ozonesonde data.

Hide code cell source

 
import os
import cdsapi
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
import numpy as np
import pandas as pd
import requests
import xarray as xr
from c3s_eqc_automatic_quality_control import diagnostics, download

plt.style.use("seaborn-v0_8-notebook")

Define data requests and other parameters#

The analysis performed in this notebook focuses on ozone data from both limb and nadir satellite sensors [2]. Specifically, it uses limb and nadir data extracted over latitude bands centred on the ozonesonde station locations (see ‘Download ozonesonde data by station’ section below), at an atmospheric pressure level of 50 hPa. This pressure level, corresponding to an altitude of about 20 km, lies within the lower stratosphere and is commonly used to evaluate ozone profiles due to the high ozone concentrations typically found at this altitude [3].

Hide code cell source

 
# Client configuration
os.environ["CDSAPI_RC"] = os.path.expanduser("~/calmanti_sandro/.cdsapirc")

# Latitudes
lat_slices = [slice(54.5, 55.5), slice(-65.5, -64.5)]

# Pressure level in hPa
pressure = 50

# Set satellite data requests
collection_id = "satellite-ozone-v1"
request_dict = {
    "limb": {
        "processing_level": "level_3",
        "variable": "mole_concentration_of_ozone_in_air",
        "vertical_aggregation": "vertical_profiles_from_limb_sensors",
    },
    "nadir": {
        "processing_level": "level_3",
        "variable": "mole_content_of_ozone_in_atmosphere_layer",
        "vertical_aggregation": "vertical_profiles_from_nadir_sensors",
    },
}

2. Satellite data caching functions and retrieval#

Define required functions#

The satellite datasets require the dedicated preprocess function to properly handle their time dimension. While the nadir datasets directly provide ozone volume mixing ratio (i.e., mole fraction in ppv), the limb datasets require it to be calculated from ozone concentration (mol m\(^{-3}\)), along with the corresponding temperature (in K) and atmospheric pressure (in hPa). In this case, an ideal gas constant (rconst) of 8.314 J mol\(^{-1}\) K\(^{-1}\) is used. Moreover, while the limb data are already provided as zonal means, the nadir data must be spatially averaged. Appropriate unit conversions are applied to express ozone mixing ratio in parts per million by volume (ppmv).

Hide code cell source

 
def preprocess(ds):
    if "time" in ds.dims:
        return ds
    time_str = ds.attrs["time_coverage_start"][:6]
    time = pd.to_datetime(time_str, format="%Y%m")
    return ds.expand_dims(time=[time])

def compute_ozone_mixing_ratio(ds, lat_slices, pressure, rconst=8.314):
    # Select pressure
    if "pressure" in ds.dims:
        ds = ds.sel(pressure=pressure)
    else:
        ds = ds.where(ds["pressure"].notnull())
        diff = np.abs(ds["pressure"] - pressure).fillna(9999)
        ds = ds.isel(altitude=diff.argmin("altitude").compute())

    dataarrays = []
    for lat_slice in lat_slices:
        ds_lat = ds.sel(latitude=lat_slice)

        # Spatially weighted mean
        ds_lat = ds_lat.mean(
            set(ds_lat.dims) - {"time", "latitude", "longitude"},
            keep_attrs=True,
        )
        if ds_lat.sizes["latitude"] != 1:
            latitude = (
                ds_lat["latitude"].coarsen(latitude=ds_lat.sizes["latitude"]).mean()
            )
            ds_lat = diagnostics.spatial_weighted_mean(ds_lat)
            ds_lat = ds_lat.expand_dims(latitude=latitude)
            ds_lat["latitude"].attrs = latitude.attrs

        # Compute mixing ratio
        if "ozone_mixing_ratio" in ds_lat:
            da = ds_lat["ozone_mixing_ratio"] * 1.0e6 # convert volume mixing ratio (ppv) to ppmv
        else:
            da = (
                (ds_lat["ozone_concentration"] * rconst * ds_lat["temperature"])
                / ds_lat["pressure"]
            ) * 1e4 # compute and convert mole fraction to ppmv
        dataarrays.append(da.compute())

    da = xr.concat(dataarrays, "latitude")
    da.attrs = {"long_name": "O₃ mixing ratio", "units": "ppm"}
    return da.to_dataset(name="ozone_mixing_ratio")

def download_combined_dataset(collection_id, request):
    collection = cdsapi.Client(quiet=True).client.get_collection(collection_id)

    datasets = []
    for sensor in collection.apply_constraints(request)["sensor"]:
        if sensor in ["cllg", "cmzm", "merged_np"]:
            continue  # no merged products
        sensor_r = request | {"sensor": sensor}
        for algo in collection.apply_constraints(sensor_r)["algorithm"] or [None]:
            algo_r = sensor_r | {"algorithm": algo or []}
            for version in collection.apply_constraints(algo_r)["version"]:
                version_r = algo_r | {"version": version}
                request_list = []
                for year in collection.apply_constraints(version_r)["year"]:
                    year_r = version_r | {"year": year}
                    months = collection.apply_constraints(year_r)["month"]
                    request_list.append(year_r | {"month": months})

                product = "-".join([sensor, version] + ([algo] if algo else []))
                if product == "omps-v0002-usask":
                    continue  # time_coverage_start attribute is missing

                ds = download.download_and_transform(
                    collection_id,
                    request_list,
                    preprocess=preprocess,
                    transform_func=compute_ozone_mixing_ratio,
                    transform_func_kwargs={
                        "lat_slices": lat_slices,
                        "pressure": pressure,
                    },
                    quiet=True,
                )
                datasets.append(ds.expand_dims(product=[product]))
    return xr.concat(datasets, "product")

Data retrieval#

Limb and nadir datasets are downloaded and combined into two separate datasets. Overall, 17 limb and 19 nadir datasets are downloaded. The limb dataset OMPS-v0002-usask is excluded from the analysis, as most of its files (106 out of 130) lack the required time attribute.

Hide code cell source

 
datasets = {}
for sensor, request in request_dict.items():
    datasets[sensor] = download_combined_dataset(collection_id, request).compute()

2025-06-26 08:35:07,870 WARNING [2025-06-23T00:00:00] Scheduled System Session affecting Service reliability - 30 June 2025. Please follow status [here](https://status.ecmwf.int/) or in our [forum](https://forum.ecmwf.int/t/scheduled-maintenance-of-the-cloud-infrastructure-on-30-june-2025/13598)
2025-06-26 08:36:06,558 WARNING [2025-06-23T00:00:00] Scheduled System Session affecting Service reliability - 30 June 2025. Please follow status [here](https://status.ecmwf.int/) or in our [forum](https://forum.ecmwf.int/t/scheduled-maintenance-of-the-cloud-infrastructure-on-30-june-2025/13598)

3. Download external ozonesonde data#

Ozonesonde data are retrieved from the World Ozone and Ultraviolet Radiation Centre (WOUDC).

Define function to extract ozonesonde data#

The get_station function is implemented to retrieve in situ ozone partial pressure (mPa) and atmospheric pressure (hPa) data from the WOUDC database. Only data at a pressure level of about 50 hPa are extracted. Ozonesonde data are combined to calculate the ozone mixing ratio (ppmv), which is then averaged monthly.

Hide code cell source

 
def get_station(station_gaw_id, pressure):
    response = requests.get(
        url="https://api.woudc.org/collections/ozonesonde/items",
        params={"station_gaw_id": station_gaw_id, "f": "json", "limit": 100000},
    )
    response.raise_for_status()
    data = response.json()
    dataarrays = []
    for feature in data["features"]:
        properties = feature["properties"]
        da = xr.DataArray(
            [
                o3pp if o3pp is None else float(o3pp)
                for o3pp in properties["o3partialpressure"]
            ],
            dims="pressure",
            coords={"pressure": list(map(float, properties["pressure"]))},
        )
        if pressure not in da["pressure"]:
            continue

        da = da.sel(pressure=pressure, drop=True)
        if "pressure" in da.dims:
            da = da.mean("pressure")  # Mean of duplicates
        time = pd.to_datetime(properties["timestamp_date"]).tz_localize(None)
        da = da.expand_dims(time=[time])
        dataarrays.append(da)

    da = xr.concat(dataarrays, "time").sortby("time").resample(time="1MS").mean()
    da = (da / pressure) * 10 # convert mole fraction to ppmv
    da.attrs = {"long_name": "O₃ mixing ratio", "units": "ppm"}
    return da.rename("ozone_mixing_ratio")

Download ozonesonde data by station#

Ozonesonde data are downloaded as single array from two stations: Edmonton (EDT) and Marambio (MBI). The EDT station is located in Canada (lat 53.55°N, lon 114.09°W) and represents mid- to subpolar latitudes in the Northern Hemisphere, while MBI station is in Antarctica (lat 64.24°S, lon 56.62°W) and represents high southern latitudes. Among the mid- and high-latitude stations in the WOUDC database, these two offer the longest time series available. More information can be found on the WOUDC website.

Hide code cell source

 
# Select ozonesonde stations
stations = {55: "EDT", -65: "MBI"}

# Download ozonesonde data
da_station = xr.concat(
    [
        get_station(station, pressure=pressure).expand_dims(station=[station])
        for station in stations.values()
    ],
    "station",
)

/data/common/miniforge3/envs/wp5/lib/python3.11/site-packages/xarray/computation/nanops.py:118: RuntimeWarning: invalid value encountered in scalar divide
  data = data / (valid_count - ddof)
/data/common/miniforge3/envs/wp5/lib/python3.11/site-packages/xarray/computation/nanops.py:118: RuntimeWarning: invalid value encountered in scalar divide
  data = data / (valid_count - ddof)
/data/common/miniforge3/envs/wp5/lib/python3.11/site-packages/xarray/computation/nanops.py:118: RuntimeWarning: invalid value encountered in scalar divide
  data = data / (valid_count - ddof)

4. Compute statistics#

Define statistical function#

The compute_statistics function is implemented to calculate the Root Mean Square Error (RMSE), Mean Bias Error (MBE), Mean Absolute Percentage Error (MAPE) and Standard Error (SE) between ozone monthly mixing ratios from limb and nadir datasets and in situ ozonesonde observations.

Hide code cell source

 
def compute_statistics(reference, comparison):
    mean_dims = ["product", "time"]
    diff = comparison - reference
    rmse = np.sqrt((diff**2).groupby(ds["sensor"]).mean(mean_dims)) #Root Mean Squared Error
    mbe = diff.groupby(ds["sensor"]).mean(mean_dims) #Mean bias error
    mape = np.abs(diff / reference).groupby(ds["sensor"]).mean(mean_dims) * 100 #Mean Absolute Percentage Error
    se = (diff.groupby(ds["sensor"]).std(mean_dims))/(np.sqrt(diff.groupby(ds["sensor"]).count(mean_dims))) #Standard Error
    dataarrays = [rmse.rename("RMSE"), mbe.rename("MBE"), mape.rename("MAPE"), se.rename("SE")]
    return xr.merge(dataarrays).to_pandas()

5. Display and discuss results#

Display results#

Monthly time series are displayed to compare satellite-derived ozone mixing ratios against ozonesonde measurements at two representative latitudes corresponding to the EDT (55°N) and MBI (65°S) stations. The satellite datasets are grouped by geometry: limb-viewing sensors in the upper panels and nadir-viewing in the lower panels. Each colored lines represents a satellite product, while black line corresponds to ozonesonde data. For each station, a scatter plot comparing monthly limb and nadir ozone mixing ratios (in red and green, respectively) with ozonesonde data is also presented, with statistical metrics indicated in the legend. The statistical analysis is performed by excluding monthly data before year 2002 due to the high discontinuity in both satellite and ozonesonde measurements. Dashed black 1:1 line in the scatter plots indicates perfect agreement between satellite and ozonesonde ozone mixing ratios.

Hide code cell source

 
# Display ozonesonde and satellite data time series
for sensor, ds in datasets.items():
    (da,) = ds.data_vars.values()
    da = da.dropna("product", how="all")
    da_station_cutout = da_station.sel(
        time=slice(*da["time"].dt.strftime("%Y-%m").values[[0, -1]])
    )
    cmap = plt.get_cmap("tab10", da.sizes["product"])
    colors = [cmap(i) for i in range(da.sizes["product"])]
    with plt.rc_context({"axes.prop_cycle": plt.cycler(color=colors)}):
        facet = da.plot(row="latitude", hue="product", figsize=(10, 8))
    for ax, sel_dict in zip(facet.axs.flatten(), facet.name_dicts.flatten()):
        station = stations[sel_dict["latitude"]]
        ax.plot(
            da_station_cutout["time"],
            da_station_cutout.sel(station=station),
            color="k",
            lw=1.5,
            label=station,
            zorder=1,
        )
        ax.legend()
        ax.grid()
        ax.set_title(f"Latitude: {sel_dict['latitude']}°", x=0.1,  fontsize=12)
        ax.set_ylabel("O$_3$ mixing ratio (ppmv)")
    facet.fig.suptitle(f"{sensor.capitalize()}-Ozonesonde datasets comparison", y=1.01, x=0.4, fontsize=15)
    plt.show()

# Display scatter plots and statistics
for latitude, station in stations.items():
    fig, ax = plt.subplots(1, 1)
    ds = xr.concat(
        [
            ds.sel(latitude=latitude)
            .reset_coords(drop=True)
            .assign_coords(sensor=("product", [sensor] * ds.sizes["product"]))
            for sensor, ds in datasets.items()
        ],
        dim="product",
    )
    ds = ds.assign_coords(ozonesonde=da_station.sel(station=station))
    ds = ds.sel(time=slice("2002-01-01","2024-12-01"))
    for var, prefix in zip(["ozone_mixing_ratio", "ozonesonde"], ["Satellite", "Ozonesonde"]):
        ds[var].attrs["long_name"] = f"{prefix} {ds[var].long_name}"
    sc = ds.plot.scatter(
        x="ozonesonde",
        y="ozone_mixing_ratio",
        hue="sensor",
        s=10,
        ax=ax,
        alpha=0.8,
        cmap="prism", 
        add_legend=False 
    )
    if sc.colorbar:
        sc.colorbar.remove()
    sensors = np.unique(ds.sensor.values)
    colors = plt.get_cmap("prism")(np.linspace(0, 1, len(sensors)))
    handles = [mpatches.Patch(color=colors[i], label=sensor) for i, sensor in enumerate(sensors)]
    handles.append(mlines.Line2D([], [], color='k', linestyle='--', label='1:1 line'))
    limits = ax.get_xlim() + ax.get_ylim()
    coords = [min(limits), max(limits)]
    ax.plot(coords, coords, color="k", ls="--", lw=2, zorder=1)
    ax.set_ylabel("Satellite O$_3$ mixing ratio (ppmv)")
    ax.set_xlabel("Ozonesonde O$_3$ mixing ratio (ppmv)")
    ax.legend(handles=handles, loc="best")
    ax.grid()
    # Compute and display statistics
    df = compute_statistics(ds["ozonesonde"], ds["ozone_mixing_ratio"])
    plt.figtext(
        1,
        0.5,
        str(df.T),
        ha="left",
        va="center",
        fontfamily="monospace",
    )
    plt.tight_layout()
    plt.show()
../../_images/2e11f9d3e2673c9b9ed9f3c60c4384d3ca1f392e9ad1ddf4d08ff3076a69f3d9.png ../../_images/63b341ee6601797ec4d2179bf903be7d02a08ad9758ebcd8de73e8435e14bcdc.png ../../_images/b3fec69c783d65bc9f8a1e227c41fbb20c1c0e5a934a94e44bab2f878883c4b5.png ../../_images/b703c0ea4e0ec79311bff53d17bb11210b6ee78c667fa8a6d9d15a30c6ee2f3b.png

Discussion#

Monthly time series

Time series figures show the seasonal cycle of ozone mixing ratio at two latitudes (55°N and 65°S), approximately corresponding to the EDT and MBI ozonesonde stations, respectively. The seasonal cycle is more pronounced at 65°S, and the results demonstrate the consistency of different satellite products in capturing this variability. These figures also provide a long-term perspective by displaying overlapping datasets from multiple satellite missions, including both limb- and nadir-viewing sensors, thereby increasing confidence in trend analyses. However, some inconsistencies are evident between datasets, even when derived from sensors with the same viewing geometry. For example, notable differences are observed between the ACE and OSIRIS limb products at both stations due to different retrieval techniques, as well as temporal, spatial and spectral resolution [4]. A similar analysis reveals discrepancies between successive versions of the same product (e.g., the GOME nadir dataset), which arise from algorithm improvements and extended temporal coverage [2] [4]. These results highlight the improved performance and stability of updated dataset versions , emphasising the critical need for using the latest reprocessed records to ensure data stability for accurate long-term ozone trend assessments (see Product’s Overview [here]). Additionally, the figures indicate a clear discontinuity in both satellite and ozonesonde records before 2002 at 65°S, likely due to the limited satellite coverage and the sporadic nature of in situ campaigns during the austral winter.

Scatter plots and statistics

Although limb-viewing sensors are generally regarded as more accurate than nadir-viewing sensors for stratospheric ozone retrievals, the scatter plots and associated statistical metrics reveal only minor differences between the two observation geometries in their agreement with ozonesonde measurements at the two high-latitude stations.
However, the overall agreement is substantially poorer at the MBI station (65°S) than at the EDT station (55°N). This difference is particularly evident in the RMSE, which increases from approximately 0.34 ppmv at EDT to more than 0.52 ppmv at MBI for both sensor types. This pattern indicates that satellite products encounter greater difficulty under the extreme atmospheric conditions of the Antarctic region, likely related to the polar vortex, and strong temperature and ozone variability, than at Northern mid-to-high latitudes. Across all datasets, the mean bias error (MBE) at EDT remains relatively small, with a slight underestimation by nadir products (−0.085 ppmv). At MBI, nadir products instead show a small positive bias (0.077 ppmv). Although this value is numerically smaller than at EDT, the substantially higher RMSE and MAPE (exceeding 22%) at MBI indicate a larger overall uncertainty. In addition, limb products at MBI exhibit a considerably larger positive bias (0.283 ppmv), nearly four times higher than that of the nadir products. While this may appear counterintuitive, the smaller standard error (SE) associated with the limb products (0.012 vs 0.017 ppmv for nadir) suggests a more precise and consistent representation of the in situ observations. This likely reflects a more stable and homogeneous set of satellite–sonde matchups.
The scatter plots further show that the agreement between satellite and ozonesonde ozone mixing ratios depends on the absolute ozone concentration. At the EDT station, nadir products exhibit a relatively larger overall negative bias, yet they agree better with ozonesonde data when ozone mixing ratios are below 3.0 ppmv, whereas at higher concentrations, all satellite products tend to underestimate ozone. In contrast, at the MBI station, both limb and nadir products generally overestimate ozone mixing ratios, particularly for values below 2.5 ppmv.
These concentration-dependent biases have important implications for interpreting ozone recovery. As satellites tend to overestimate ozone during low-value events at MBI, they may underestimate the severity of Antarctic ozone depletion episodes. This could result in an overly optimistic assessment of ozone layer recovery in policy contexts. In addition, the discontinuities and version-to-version differences (e.g., GOME-v0004 vs. v0006) evident in the time series highlight the risk of attributing apparent recovery trends to physical changes when they may instead reflect algorithm updates. The consistent use of the most recent, fully reprocessed, and harmonized data records is therefore essential for robust long-term trend analyses and for supporting evidence-based decisions under the Montreal Protocol.

ℹ️ If you want to know more#

Key resources#

Code libraries used:

Further readings:

References#

[1] Bernhard, G.H., et al. (2023). Stratospheric ozone, UV radiation, and climate interactions. Photochemical and Photobiological Science, 16(10), 1519-1523.

[2] Van Roozendael, M. (2024). Product User Guide and Specification (PUGS): Level 3 and Level 4 ozone columns and profiles. E.U. Copernicus Climate Change Service. Document ref. C3S2_312a_Lot2_D-WP2-FDDP-PUGS_202311_O3_v3.5.

[3] Bekki, S., & Lefevre, F. (2009). Stratospheric ozone: History and concepts and interactions with climate. European Physical Journal Conferences, 1, 113-136.

[4] Van Roozendael, M. (2024). Ozone Algorithm Theoretical Basis Document (ATBD), Version: Multiple Ozone products and related versions. E.U. Copernicus Climate Change Service. Document ref. C3S2_312a_Lot2_D-WP2-FDDP-ATBD_202311_O3_v3.6.