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.1.1. Satellite Earth’s outgoing longwave and shortwave radiation budget intercomparison for climate monitoring

Production date: 11-06-2025

Produced by: CNR-ISMAR, Andrea Storto, Vincenzo de Toma

🌍 Use case: Consistency assessment of Earth Radiation Budget products

❓ Quality assessment question(s)

• How consistent are the Earth Radiation Budget products provided by C3S?

The Earth Radiation Budget (ERB) Outgoing Longwave Radiation (OLR) and Outgoing Shortwave Radiation (OSR) products allow quantifying outgoing longwave and shortwave radiation, respectively, which are in turn important for several climate applications. In particular, a non-exhaustive list of possible applications include: i) evaluation and potential improvement of climate models, through routinely validation of CMIP individual simulations and ensembles against OLR and OSR data, allowing to assess the net radiative forcing and possibly leading to improved configurations through re-tuning of e.g. aerosol-cloud interactions, radiative transfer schemes, etc.; ii) tracking the Earth’s Energy Imbalance (EEI), which consists of long-term mean imbalance of incoming and outgoing energy at the top-of-the-atmosphere (TOA). This imbalance is mostly accumulated into the oceans, and its analysis is crucial to monitor recent global warming; iii) diagnosing El-Niño Southern Oscillation (ENSO) and other climate variability modes, which are intimately linked to anomalies of OLR (e.g. enhanced convection); iv) analyse cloud radiative effects (e.g. in combination with clear-sky radiative data) and their impact on climatic changes; and many more. However, their temporal coverage and spatial detail varies across the ERB products. Over the overlapping periods, the products are in good agreement, except for the OLR at high latitudes (e.g. north of 60*N / south of 60°S), while discrepancies in OSR arise at both mid- (generally less than 10 W m-2) and high-latitudes.

📢 Quality assessment statement

These are the key outcomes of this assessment

• The OLR/OSR products capture the spatial and temporal variability of the outgoing longwave radiation / outgoing shortwave radiation (at the top of the atmosphere) for use in climate analyses, climate model evaluations, and radiative forcing studies.

• The ERB OLR/OSR products have different temporal coverage and spatial resolution, namely their combined use is not straightforward and it is recommended for advanced users only, or for selected periods.

• The consistency between the products concerning spatially averaged values is generally acceptable. Climatological maps and zonally averaged values indicate that the products are sufficiently consistent in most areas, although differences exist and tend to increase near the polar regions.

The estimation of outgoing longwave radiation (OLR) at the top of the atmosphere is crucial for understanding the Earth’s radiation budget (ERB), which plays a significant role in climate dynamics and energy balance. Despite its importance, accurately assessing OLR-ERB variability using remote sensing data and satellite missions is fraught with challenges due to calibration inconsistencies, data gaps, and the complex interactions of radiation with atmospheric constituents ([1]; [2]; [3]; [4]).

Fig. 1.1.1.1 Monthly mean time series of OLR (first row), OSR (second row) Ensemble mean Net Top Of the Atmosphere Radiation (NTOAR) (third row) for different geographical regions (columns).

📋 Methodology

In this notebook we inter-compare the following Earth Radiation Budget Products, available through the Climate Data Store (CDS) by C3S:

1. CERES Energy Balanced and Filled (EBAF) TOA Monthly means data in netCDF Edition4.1. NASA Langley Atmospheric Science Data Center DAAC (NASA CERES EBAF)

2. NOAA CDR Program, (2018): NOAA Climate Data Record (CDR) of Monthly Outgoing Longwave Radiation (OLR), Version 2.7 (NOAA/NCEI HIRS).

3. Earth’s radiation budget from 1979 to present derived from satellite observations. CCI ICDR product version 3.1. Copernicus Climate Change Service Climate Data Store (C3S CCI).

4. ESA Cloud Climate Change Initiative (ESA Cloud_cci) data: Cloud_cci ATSR2-AATSR L3C/L3U CLD_PRODUCTS v3.0 (ESA Cloud CCI).

CERES is the only product which gives also the incoming shortwave radiation: thus, in order to evaluate the Net Top Of the Atmosphere Radiation, the ensemble mean budget was built using outgoing longwave and outgoing shortwave from each product, and the incoming shortwave from CERES, and then averaging across the product dimension.

The analysis and results are organised in the following steps, which are detailed in the sections below:

1. Choose the data to use and setup code

2. Download and transform

3. Plot spatial weighted time series

4. Plot time weighted means

5. Plot spatial weighted zonal means

📈 Analysis and results

1. Choose the data to use and setup code

In this section, we import the required packages and set up the dataset names for further use in the following sections. Processing functions are also defined. We select only complete years available for each dataset in order to have a fair comparison between different products.

Show code cell source

Hide code cell source

import os
os.environ["CDSAPI_RC"] = os.path.expanduser("~/detoma_vincenzo/.cdsapirc")
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
import datetime
from cartopy import crs as ccrs, feature as cfeature
import cartopy.mpl.ticker as cticker
from c3s_eqc_automatic_quality_control import download, plot, utils, diagnostics
import warnings
warnings.filterwarnings("ignore")

# Variable to analyse
variables = ["outgoing_longwave", "outgoing_shortwave", "incoming_shortwave"]
request_dicts = []

for variable in variables:

    # Region for timeseries
    region_slices = {
        "global": {"lat_slice": slice(-90, 90), "lon_slice": slice(0, 360)},
        "northern hemisphere": {"lat_slice": slice(0, 90), "lon_slice": slice(0, 360)},
        "southern hemisphere": {"lat_slice": slice(-90, 0), "lon_slice": slice(0, 360)},
    }

    collection_id = "satellite-earth-radiation-budget"
    chunks = {"year": 1}
    varnames = (
        {"olr", "toa_lw_all_mon"} if variable == "outgoing_longwave" else {"rsf", "toa_sw_all_mon"}
    )

    collection_id = "satellite-earth-radiation-budget"
    chunks = {"year": 1}
    varnames = (
        {"olr", "toa_lw_all_mon", "LW_flux"}
        if variable == "outgoint_longwave"
        else {"rsf", "toa_sw_all_mon", "SW_flux"}
    )

    request_dict = {
        "CERES": {
            "start": "2001-01",
            "stop": "2023-12",
            "product_family": "ceres_ebaf",
            "climate_data_record_type": "thematic_climate_data_record",
            "time_aggregation": "monthly_mean",
            "format": "zip",
            "origin": "nasa",
            "variable": f"{variable}_radiation",
        },
        "Sentinel 3A": {
            "start": "2017-01",
            "stop": "2021-12",
            "format": "zip",
            "origin": "c3s",
            "sensor_on_satellite": "slstr_on_sentinel_3a",
            "variable": "all_variables",
            "product_family": "cci",
            "time_aggregation": "monthly_mean",
            "climate_data_record_type": "interim_climate_data_record",
        },
        "Sentinel 3B": {
            "start": "2019-01",
            "stop": "2021-12",
            "format": "zip",
            "origin": "c3s",
            "sensor_on_satellite": "slstr_on_sentinel_3b",
            "variable": "all_variables",
            "product_family": "cci",
            "time_aggregation": "monthly_mean",
            "climate_data_record_type": "interim_climate_data_record",
        },
        "Sentinel 3A_3B": {
            "start": "2019-01",
            "stop": "2021-12",
            "format": "zip",
            "origin": "c3s",
            "sensor_on_satellite": "slstr_on_sentinel_3a_3b",
            "variable": "all_variables",
            "product_family": "cci",
            "time_aggregation": "monthly_mean",
            "climate_data_record_type": "interim_climate_data_record",
        },
        "ESA ENVISAT": {
            "start": "2003-01",
            "stop": "2011-12",
            "format": "zip",
            "origin": "esa",
            "product_family": "cci",
            "climate_data_record_type": "thematic_climate_data_record",
            "time_aggregation": "monthly_mean",
            "sensor_on_satellite": "aatsr",
            "variable": "all_variables",
        },
        "CLARA_A3": {
            "start": "1979-01",
            "stop": "2020-12",
            "product_family": "clara_a3",
            "origin": "eumetsat",
            "variable": f"{variable}_radiation",
            "climate_data_record_type": "thematic_climate_data_record",
            "time_aggregation": "monthly_mean",
        },
        "HIRS": {
            "start": "1979-01",
            "stop": "2023-12",
            "format": "zip",
            "origin": "noaa_ncei",
            "product_family": "hirs",
            "climate_data_record_type": "thematic_climate_data_record",
            "time_aggregation": "monthly_mean",
            "version": "2_7_reprocessed",
            "variable": f"{variable}_radiation",
        },
        "CERES_noMask": {
            "start": "2001-01",
            "stop": "2023-12",
            "product_family": "ceres_ebaf",
            "climate_data_record_type": "thematic_climate_data_record",
            "time_aggregation": "monthly_mean",
            "format": "zip",
            "origin": "nasa",
            "variable": f"{variable}_radiation",
        },
    }
    request_dicts.append(request_dict)

def preprocess_time(ds):
    if "time" in ds and "units" in ds["time"].attrs:
        # Could not decode
        ds = ds.squeeze("time", drop=True)
    if "time" not in ds:
        time = pd.to_datetime(ds.attrs["time_coverage_start"].replace("Z", ""))
        ds = ds.expand_dims(time=[time])
    return ds


def spatial_weighted_mean(ds, lon_slice, lat_slice, regrid, ds_obj, SKIP, mask=None):
    if mask is None:
        if regrid==True:
            ds = diagnostics.regrid(ds, ds_obj, method='bilinear', skipna=SKIP)
        ds = utils.regionalise(ds, lon_slice=lon_slice, lat_slice=lat_slice)
        return diagnostics.spatial_weighted_mean(ds, weights=True, skipna=SKIP)
    else:
        if regrid==True:
            ds = diagnostics.regrid(ds, ds_obj, method='bilinear', skipna=SKIP)
        ds = utils.regionalise(ds, lon_slice=lon_slice, lat_slice=lat_slice)
        mask = utils.regionalise(mask, lon_slice=lon_slice, lat_slice=lat_slice)
        return diagnostics.spatial_weighted_mean(ds.where(~np.isnan(mask)), weights=True, skipna=SKIP)

xarray_kwargs = {
    "drop_variables": ["time_bounds", "record_status"],
    "preprocess": preprocess_time,
}

2. Download and transform

The code below will download the products. Notice the definition of common periods for the different products, and their inclusion in the mean map dictionary. When the variable to download is the shortwave component, we skip the HIRS dataset, since only longwave is available within it.

Show code cell source

Hide code cell source

weights=True
skipna=True
p = []
start_var = {}
stop_var = {}
DS = {}
for variable, request_dict in zip(variables, request_dicts):
    Ds = {}
    start_prod = {}
    stop_prod = {}
    for product, request in request_dict.items():
        if product == "HIRS" and variable != "outgoing_longwave":
            continue
        if "CERES" not in product and variable=="incoming_shortwave":
            continue
        p.append(product)
    
        start = request.pop("start")
        stop = request.pop("stop")
        start_prod[product] = start
        stop_prod[product] = stop
        
        requests = download.update_request_date(
            request, start=start, stop=stop, stringify_dates=True
        )
    
        ds = download.download_and_transform(
            collection_id,
            requests,
            #transform_func=diagnostics.time_weighted_mean,
            #transform_func_kwargs={'weights':True, 'skipna':False},
            chunks=chunks,
            transform_chunks=False,
            **xarray_kwargs,
            quiet=True,
        )
        Ds[product] = ds
    start_var[variable] = start_prod
    stop_var[variable] = stop_prod
    DS[variable] = Ds

Show code cell source

Hide code cell source

clara_a3_start, clara_a3_end ="1980-01", "2020-12"
ceres_start, ceres_end = "2001-01", "2023-12"
esa_envisat_start, esa_envisat_end = "2003-01", "2011-12" 
sentinel_start, sentinel_end = "2019-01", "2021-12"

def mean_map_period(ds, start, stop, W, S):
    selection = ds.sel(time=slice(start, stop))
    time_coord =  np.array(pd.date_range(start, stop, freq='MS'), dtype=np.datetime64)
    ds = selection.assign_coords({'time':('time', time_coord, selection.time.attrs)})
    ds_map = diagnostics.time_weighted_mean(ds, weights=W, skipna=S)
    return ds_map

Show code cell source

Hide code cell source

def get_varname(ds, variable):
    if variable=="outgoing_longwave":
        varnames = (
            {"olr", "toa_lw_all_mon", "LW_flux"}
        )
    elif variable=="outgoing_shortwave":
        varnames = (
            {"rsf", "toa_sw_all_mon", "SW_flux"}
        )
    else:
        varnames = (
            {"solar_mon"}
        )
        
    (varname,) = set(ds.data_vars) & varnames
    return varname

ds_maps = {}
ds_maps_p1 = {}
ds_maps_p2 = {}
ds_maps_p3 = {}
ds_maps_p4 = {}
p = []
for variable, request_dict in zip(variables, request_dicts):
    da_maps = {}
    da_maps_p1 = {}
    da_maps_p2 = {}
    da_maps_p3 = {}
    da_maps_p4 = {}
    for product, request in request_dict.items():
        if product == "HIRS" and variable != "outgoing_longwave":
            continue
        if "CERES" not in product and variable=="incoming_shortwave":
            continue
        p.append(product)
    
        start = start_var[variable][product]
        stop = stop_var[variable][product]

        ds = DS[variable][product]
        
        varname=get_varname(ds, variable)
        
        da = ds[varname]
        if product=="CERES":
            regrid=False
            da_obj = ds[varname]
        else:
            regrid=True
            da = diagnostics.regrid(da, da_obj, method='bilinear')
        da.attrs.update({"start": start, "stop": stop})
        da_maps[product] = diagnostics.time_weighted_mean(da, weights=True, skipna=False)

        curr_time = pd.date_range(start, stop)
        ceres_time = pd.date_range(ceres_start, ceres_end)
        clara_a3_time = pd.date_range(clara_a3_start, clara_a3_end)
        esa_envisat_time = pd.date_range(esa_envisat_start, esa_envisat_end)
        sentinel_time = pd.date_range(sentinel_start, sentinel_end)
        if clara_a3_time[0] in curr_time and clara_a3_time[-1] in curr_time:
            da_maps_p1[product] = mean_map_period(da, clara_a3_start, clara_a3_end, W=weights, S=skipna)
            da_maps_p1[product] = diagnostics.regrid(da_maps_p1[product].rename(variable), da_obj, method='bilinear')
        if ceres_time[0] in curr_time and ceres_time[-1] in curr_time:
            da_maps_p2[product] = mean_map_period(da, ceres_start, ceres_end, W=weights, S=skipna)
            da_maps_p2[product] = diagnostics.regrid(da_maps_p2[product].rename(variable), da_obj, method='bilinear')
        if esa_envisat_time[0] in curr_time and esa_envisat_time[-1] in curr_time:
            da_maps_p3[product] = mean_map_period(da, esa_envisat_start, esa_envisat_end, W=weights, S=skipna) 
            da_maps_p3[product] = diagnostics.regrid(da_maps_p3[product].rename(variable), da_obj, method='bilinear')
        if sentinel_time[0] in curr_time and sentinel_time[-1] in curr_time:
            da_maps_p4[product] = mean_map_period(da, sentinel_start, sentinel_end, W=weights, S=skipna) 
            da_maps_p4[product] = diagnostics.regrid(da_maps_p4[product].rename(variable), da_obj, method='bilinear')

    ds_maps_p1[variable] = da_maps_p1
    ds_maps_p2[variable] = da_maps_p2
    ds_maps_p3[variable] = da_maps_p3
    ds_maps_p4[variable] = da_maps_p4
    ds_maps[variable] = da_maps

Show code cell source

Hide code cell source

listmasks = [ds_maps[k] for k in ds_maps.keys()]
lw_mask = xr.concat([xr.where(~np.isnan(listmasks[0][l]), 1, np.nan) for l in listmasks[0].keys()], dim='product').mean(['product'], skipna=False)
sw_mask = xr.concat([xr.where(~np.isnan(listmasks[1][l]), 1, np.nan) for l in listmasks[1].keys()], dim='product').mean(['product'], skipna=False)
mask = xr.concat([lw_mask, sw_mask], dim='variable').assign_coords({'variable':variables[:2]})

Show code cell source

Hide code cell source

ds_timeseries = {}
p = []
for variable, request_dict in zip(variables, request_dicts):
    da_timeseries = {}
    for product, request in request_dict.items():
        if product == "HIRS" and variable != "outgoing_longwave":
            continue
        if "CERES" not in product and variable=="incoming_shortwave":
            continue
        p.append(product)
        
        start = start_var[variable][product]
        stop = stop_var[variable][product]

        ds = DS[variable][product]
        
        varname=get_varname(ds, variable)
        
        da = ds[varname]
        if product=="CERES":
            regrid=False
            da_obj = ds[varname]
        else:
            regrid=True
            da = diagnostics.regrid(da, da_obj, method='bilinear')
        da.attrs.update({"start": start, "stop": stop})
        
        regrid_dict={'regrid':regrid, 'ds_obj':da_obj}

        SkipNa={'SKIP': True}

        if product=="CERES_noMask":
            Mask = {'mask': None}
        elif variable=='incoming_shortwave':
            Mask = {'mask': None}  
        else:
            Mask = {'mask': mask.sel(variable=variable)}
        
        dataarrays = []
        for region, slices in region_slices.items():
            dall = {}
            dall.update(slices)
            dall.update(regrid_dict)
            dall.update(SkipNa)
            dall.update(Mask)
            da_ts = spatial_weighted_mean(da, lon_slice=dall['lon_slice'], lat_slice=dall['lat_slice'], regrid=dall['regrid'], ds_obj=dall['ds_obj'], SKIP=dall['SKIP'], mask=dall['mask'])
            dataarrays.append(da_ts.expand_dims(region=[region]))
        
        da = xr.concat(dataarrays, "region")
        da.attrs.update({"start": start, "stop": stop})
        da_timeseries[product] = da.rename(variable)
    ds_timeseries[variable] = da_timeseries

Show code cell source

Hide code cell source

ds_periods_list = [ds_maps_p1, ds_maps_p2, ds_maps_p3, ds_maps_p4]
periods = [clara_a3_time, ceres_time, esa_envisat_time, sentinel_time]

3. Plot spatial weighted time series

Below, we calculate and plot spatially weighted means for the different earth radiation budget products, in terms of montly climatologies, although the sampling period of each dataset is, in general, different from the others (i.e. CERES: 2001-2023, Sentinel 3A: 2017-2021, Sentinel 3B: 2019-2021, Sentinel 3A_3B: 2019-2021, ESA_Envisat: 2003-2011, CLARA_A3: 1979-2020, HIRS: 1979-2023, CERES_noMask: 2001-2023). Despite each dataset has a different length, we want here to inspect how the different products reproduce the seasonal cycle, both for the global mean, and for each hemisphere. Notice we also used the non-masked version of CERES (CERES_noMask), for two reasons: fist, we want to have an idea of what is the main effect of the masking (i.e. if it affects mainly the northern and southern hemisphere), and second to have the possibility of giving the best possible estimate of the Net Top of the Atmosphere Radiation (i.e. as much as possible close to a sort of Earth Energy Imbalance).

The pronounced seasonality in OLR is mostly linked with the larger land area in the Northern Hemisphere, while that one in OSR linked with the solar cycle.

Show code cell source

Hide code cell source

lw_ts = xr.concat([ds_timeseries['outgoing_longwave'][k] for k in ds_timeseries['outgoing_longwave'].keys()], dim='product', coords='minimal').assign_coords({'product': [k for k in ds_timeseries['outgoing_longwave'].keys()]})
sw_ts = xr.concat([ds_timeseries['outgoing_shortwave'][k] for k in ds_timeseries['outgoing_shortwave'].keys()], dim='product', coords='minimal').assign_coords({'product': [k for k in ds_timeseries['outgoing_shortwave'].keys()]})

t1 = lw_ts.groupby('time.month').mean('time').plot(col='region', hue='product', figsize=(12,4), marker='o', markeredgecolor='k') #.resample(time='1M').mean()
for ax in t1.axs.flatten():
    ax.set_xticks([1,2,3,4,5,6,7,8,9,10,11,12])
    ax.set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
    ax.grid()
    

Show code cell source

Hide code cell source

t2 = sw_ts.groupby('time.month').mean('time').plot(col='region', hue='product', figsize=(12,4), marker='o', markeredgecolor='k') #.resample(time='1M').mean()
for ax in t2.axs.flatten():
    ax.set_xticks([1,2,3,4,5,6,7,8,9,10,11,12])
    ax.set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
    ax.grid()

The offsets are greater in OLR than in OSR, with the ESA product (ENVISAT) being positively biased compared to the other products.

Show code cell source

Hide code cell source

p_nohirs = sw_ts.product.values
res_var = {}

def fix_lon(ds):
    if ds.longitude.max().values>350.:
        ds['longitude'] = ((ds['longitude'] + 180) % 360) - 180
        # Step 2: Sort by the new longitude values to maintain increasing order
        ds = ds.sortby('longitude')
    else:
        ds=ds
    return ds
    
for request_dict in request_dicts:
    prod = {}
    for product, request in request_dict.items():
        if product == "HIRS" and variable != "outgoing_longwave":
            continue
        p.append(product)

        lw_ds = DS['outgoing_longwave'][product].groupby('time.month').mean('time')
        sw_ds = DS['outgoing_shortwave'][product].groupby('time.month').mean('time')
        sw_ds_in = DS['incoming_shortwave']['CERES'].groupby('time.month').mean('time')
        sw_ds_in = fix_lon(sw_ds_in)
        lw_ds = fix_lon(lw_ds)
        sw_ds = fix_lon(sw_ds)

        shortname_out = get_varname(sw_ds, 'outgoing_shortwave')
        longname = get_varname(lw_ds, 'outgoing_longwave')
        shortname_in = get_varname(sw_ds_in, 'incoming_shortwave')
        da = np.abs(sw_ds[shortname_out]) - np.abs(sw_ds_in[shortname_in]) + np.abs(lw_ds[longname])
        dataarrays = []
        for region, slices in region_slices.items():
            dall = {}
            dall.update(slices)
            dall.update(regrid_dict)
            dall.update(SkipNa)
            dall.update(Mask)
            da_sel = da.sel(longitude=dall['lon_slice'], latitude=dall['lat_slice'])
            da_ts = da_sel.weighted(np.cos(da_sel.latitude*np.pi/180.)).mean(['longitude', 'latitude'])
            #da_ts = spatial_weighted_mean(da, lon_slice=dall['lon_slice'], lat_slice=dall['lat_slice'], regrid=dall['regrid'], ds_obj=dall['ds_obj'], SKIP=dall['SKIP'], mask=dall['mask'])
            dataarrays.append(da_ts.expand_dims(region=[region]))
        
        da = xr.concat(dataarrays, "region")
        prod[product] = da
    res_var['EEI'] = prod

Most satellite datasets available from the Climate Data Store (CDS), such as those from the CLARA or HIRS record, are not gap-filled and do not always provide full global coverage. Therefore, they cannot be directly used to compute global EEI without additional assumptions, spatial extrapolation, or fusion with other datasets. Here, CERES-EBAF is referenced as a benchmark for EEI estimation, but the datasets from the CDS are not used to compute EEI directly, since they do not currently support a comparable global radiative balance analysis. Still, it is possible to look at the Imbalance annual cycles for the Net Top Of the Atmosphere Radiation, which is shown below. Specifically, NTOAR is calculated for each of the datasets in this notebook as the difference between the net shortwave and net longwave component, but calculating the net shortwave using always the incoming component from CERES-EBAF. In the end, the mean across product is extracted, constituting a sort of Ensemble mean.

Show code cell source

Hide code cell source

res = xr.concat([res_var['EEI'][p] for p in p_nohirs], dim='product').assign_coords({'product': p_nohirs}).mean('product')
t3 = res.plot(col='region', figsize=(12,4), marker='o', markeredgecolor='k', label='NTOAR EnsMean') #.resample(time='1M').mean()
for a, ax in enumerate(t3.axs.flatten()):
    if a==0:
        ylabel = (
        "Net\n"
        "Top Of the Atmosphere\n"
        r"Radiation = $\langle (sw_{\mathrm{out}} - sw_{\mathrm{in}}^{[\mathrm{CERES}]}) - lw \rangle$"
        "\n[W m$^{-2}$]"
        )
        ax.set_ylabel(ylabel)
    ax.set_xticks([1,2,3,4,5,6,7,8,9,10,11,12])
    ax.set_xticklabels(['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec'])
    ax.legend(loc='best')
    ax.grid()

Next, we show the Net Top Of the Atmosphere Radiation annual cycle for the whole available domain (with masking and regridding), and for the Northern and Southern Hemisphere. The EEI is the small difference between the sunlight the Earth absorbs and the thermal (longwave) plus reflected (shortwave) radiation it emits, and represents a fundamental “engine” of climate change. Even a tiny positive imbalance (on the order of 0.5-1 W m-2, according to recent estimates [13]) sustained year after year is what drives ocean heat uptake, ice melt, sea-level rise and ultimately the global temperature increase we observe.

The global average of Net Top Of the Atmosphere Radiation indicates a small seasonal cycle (about 15 W m-2 peak-to-peak) compared to hemispheric averages, with a slight maximum Net Top Of the Atmosphere Radiation in June-July. That means mid-year the Earth is absorbing just a bit more energy than at other times, largely because Northern Hemisphere land is warm, and snow-free.

4. Plot time weighted means

Below, we calculate and plot time-weighted means for the different radiation budget products. Please note that different periods were considered, in order to make fair comparison between products wich have a different temporal extent.

Show code cell source

Hide code cell source

listmasks = [ds_maps[k] for k in ds_maps.keys()]
lw_mask = xr.concat([xr.where(~np.isnan(listmasks[0][l]), 1, np.nan) for l in listmasks[0].keys()], dim='product').mean(['product'], skipna=False)
sw_mask = xr.concat([xr.where(~np.isnan(listmasks[1][l]), 1, np.nan) for l in listmasks[1].keys()], dim='product').mean(['product'], skipna=False)

mask = xr.concat([lw_mask, sw_mask], dim='variable').assign_coords({'variable':variables[:2]})

Show code cell source

Hide code cell source

for ds_maps, times in zip(ds_periods_list, periods):
    p_to_show = []
    for k in ds_maps['outgoing_longwave'].keys():
        if k=='CERES_noMask':
            continue
        else:
            p_to_show.append(k)
    
    p_to_show_whirs = p_to_show.copy()
    p_to_show.remove("HIRS")
    lw_maps = xr.concat([ds_maps['outgoing_longwave'][p] for p in p_to_show_whirs], dim='product').assign_coords({'product': p_to_show_whirs})
    lw_maps_masked = lw_maps.where(~np.isnan(mask.sel(variable='outgoing_longwave')))
    p1 = lw_maps_masked.plot(col='product', 
                             col_wrap=2,
                             subplot_kws={'projection': ccrs.PlateCarree()},
                             levels=range(150, 315, 5),
                             extend="both",
                             cmap="nipy_spectral",
                             add_colorbar=False,
                             transform=ccrs.PlateCarree())
    for ax in p1.axs.flat:  # loop through the map axes
        subplotspec = ax.get_subplotspec()
        if subplotspec.is_last_row():
            ax.xaxis.set_visible(True)
        if subplotspec.is_first_col():
            ax.yaxis.set_visible(True)
        ax.coastlines('50m')
        ax.gridlines()
    p1.fig.set_layout_engine("compressed")
    cax = p1.fig.add_axes([1.05, 0.1, 0.025, 0.75])
    p1.add_colorbar(cax=cax, orientation='vertical', label=r'longwave [W m-2]')
    p1.fig.suptitle(f'Mean map on the period {times[0].year}-{times[-1].year}')
    if len(p_to_show)>0:
        sw_maps = xr.concat([ds_maps['outgoing_shortwave'][p] for p in p_to_show], dim='product').assign_coords({'product': p_to_show})
        sw_maps_masked = sw_maps.where(~np.isnan(mask.sel(variable='outgoing_shortwave')))
        p2 = sw_maps_masked.plot(col='product',
                                 col_wrap=2, 
                                 subplot_kws={'projection': ccrs.PlateCarree()},
                                 levels= range(80, 210, 5),
                                 extend="both",
                                 cmap="gist_ncar",
                                 add_colorbar=False,
                                 transform=ccrs.PlateCarree())
        for ax in p2.axs.flat:  # loop through the map axes
            subplotspec = ax.get_subplotspec()
            if subplotspec.is_last_row():
                ax.xaxis.set_visible(True)
            if subplotspec.is_first_col():
                ax.yaxis.set_visible(True)
            ax.coastlines('50m')
            ax.gridlines()
        p2.fig.set_layout_engine("compressed")
        cax = p2.fig.add_axes([1.05, 0.1, 0.025, 0.75])
        p2.add_colorbar(cax=cax, orientation='vertical', label=r'net shortwave [W m-2]')
        p2.fig.suptitle(f'Mean map on the period {times[0].year}-{times[-1].year}')

Time mean spatial maps show consistent spatial features across the products, although with large offsets and difference (e.g. near the Equator, at high latitudes). Differences are also due to varying spatial resolution, and different temporal and spatial sampling across the products.

All products show highest OLR over the subtropical deserts (e.g. Sahara, Arabian Peninsula, central Australia) and clear‐skies over the eastern subtropical ocean gyres, and lowest OLR where deep convection occurs (ITCZ, tropical rainforests) and at high latitudes under clouds. All OLR products agree on the major “warm” desert maxima and “cold” convective minima. Differences are very small (<5 W m-2) in the well‐sampled tropical belt; HIRS can show a bit more structure in the Intertropical Convergence Zone (likely due to higher spectral resolution), whereas the other products tend to smooth the strongest minima and maxima slightly. Concerning the OSR, all products feature bright subtropical desert maxima (200–220 W m-2), low OSR (80–100 W m-2) in persistently cloudy regions (ITCZ, midlatitude storm tracks, equatorial Pacific), and moderate values (120–160 W m-2) over clear‐sky oceans. The spatial patterns are remarkably consistent across the products, with differences generally of the order of less than 4 W m-2.

5. Plot spatial weighted zonal means

The code below will calculate and plot weighted zonal means for the radiation budget products. For the sake of comparison, different periods are shown consistently with different temporal overlapping windows.

Show code cell source

Hide code cell source

i=0
for ds_maps, times in zip(ds_periods_list, periods):
    fig=plt.figure(i, figsize=(12,5))
    ax1 = fig.add_subplot(121)
    ax2 = fig.add_subplot(122)
    fig.suptitle(f'Zonal Mean on the period {times[0].year}-{times[-1].year}')
    for variable, da_maps in ds_maps.items():
        if variable=='incoming_shortwave':
            continue
        ax=ax1 if variable=='outgoing_longwave' else ax2
        for product, da in da_maps.items():
            if product=='CERES_noMask':
                continue
            else:
                da = diagnostics.spatial_weighted_mean(da, dim="longitude")
                da = da.where(da!=0, np.nan)
                da.plot(ax=ax, y="latitude", label=product)
        ax.legend(loc='best')
        if variable=='outgoing_shortwave':
            ax.set_ylabel('')
        ax.grid()
        ax.set_xlabel(f"{variable} [W m-2]")
    i=i+1

Zonal averages indicate the good consistency among the products, except at thig latitudes (e.g., south of 50°S and north of 50°N). In particular, CERES and ESA ENVISAT show good agreement in the OLR, with CLARA_A3 and HIRS possibly deviating slightly, especially at higher latitudes. Also, OSR shows the largest discrepancies going poleward, which may arise from differences in cloud detection or albedo algorithms.

Discussion and applications

The differences identified across the OLR and OSR datasets primarily reflect known characteristics of the satellite instruments and retrieval algorithms from which the various climate data records are constructed. For OLR, the spread between the CERES EBAF record and the HIRS and CLARA-A3 products is consistent with the long-recognized sensitivity of longwave flux retrievals to spectral sampling, cloud detection, and scene-dependent angular corrections ([1], [4], [5]). CERES derives broadband fluxes using detailed angular-distribution models (ADMs) developed specifically for top-of-atmosphere applications [3], which typically results in smaller biases in cloudy and high-latitude regions. In contrast, HIRS and the AVHRR-based CLARA-A3 rely on narrower channels and empirical regressions to estimate broadband longwave radiation, a known source of latitudinal and surface-type–dependent discrepancies documented in their validation reports (HIRS PVR, CLARA-A3 PVR). These differences help explain the stronger gradients observed in the CLARA-A3 OLR fields and the slight warm-scene bias in HIRS over subtropical dry regions.

The Sentinel-3 SLSTR (3A, 3B, and the merged 3A+3B) interim CDRs show overall good consistency with the CERES reference but display enhanced variability in the tropics and over land surfaces. This behaviour is in line with their validation documentation, which highlights sensitivity to cloud screening and to the limited temporal extent of the interim SLSTR record relative to multi-decadal products (CERES Quality Summary, Sentinel Quality Assurance Document). The relatively short data period (2017–2021) also means that interannual variability associated with ENSO events—known to modulate tropical OLR [7]—can disproportionately influence the mean differences when compared to long records such as CERES or HIRS.

For OSR, discrepancies among products similarly reflect differences in sensor spectral coverage, calibration stability in the visible and near-infrared, and the treatment of anisotropy in the shortwave domain. CERES shortwave fluxes benefit from scene-dependent ADMs and long-term calibration stability [1], [4], [9], which contributes to the coherent latitudinal reflection patterns typically used as a benchmark. The CLARA-A3 product, derived from AVHRR, shows slightly higher reflected shortwave radiation in bright subtropical and high-latitude regions, consistent with known limitations in AVHRR-based angular corrections and cloud-phase retrievals noted in its product validation reports. Similarly, OSR from ESA’s ENVISAT AATSR and from SLSTR shows small but systematic differences over marine stratocumulus regions—areas where shortwave anisotropy and cloud optical properties strongly influence flux retrievals. These patterns are consistent with prior evaluations of shortwave absorption and model–observation differences [11], [12].

Overall, the magnitude and spatial structure of the inter-product differences observed here align with previously documented behaviour of broadband radiation datasets and with theoretical expectations derived from radiative-transfer and climate-feedback studies [6], [8], [10], [13]. The consistency of many features across independent retrievals—particularly the subtropical shortwave reflection maxima and the tropical longwave minima—supports the robustness of these products for climatological analyses, while the remaining discrepancies can be traced to well-understood algorithmic and instrumental differences identified in the respective validation reports.

ℹ️ If you want to know more

Key resources

Code libraries used:

• C3S EQC custom functions, c3s_eqc_automatic_quality_control, prepared by B-Open

References

[1] Loeb, N. G., Wielicki, B. A., Doelling, D. R., Smith, G. L., Keyes, D. F., Kato, S., … & Wong, T. (2009). Toward optimal closure of the Earth’s top-of-atmosphere radiation budget. Journal of Climate, 22(3), 748-766

[2] Harries, J. E., Brindley, H. E., Sagoo, P. J., & Bantges, R. J. (2001). Increases in greenhouse forcing inferred from the outgoing longwave radiation spectra of the Earth in 1970 and 1997. Nature, 410(6826), 355-357.

[3] Wielicki, B. A., Barkstrom, B. R., Harrison, E. F., Lee III, R. B., Smith, G. L., & Cooper, J. E. (1996). Clouds and the Earth’s Radiant Energy System (CERES): An Earth observing system experiment. Bulletin of the American Meteorological Society, 77(5), 853-868.

[4] Loeb, N. G., Wielicki, B. A., Su, W., Loukachine, K., Sun, W., Wong, T., … & Lin, B. (2012). Advances in understanding top-of-atmosphere radiation variability from satellite observations. Surveys in Geophysics, 33(3-4), 359-385

[5] Susskind, J., G. Molnar, L. Iredell, et al., 2012: Interannual variability of outgoing longwave radiation as observed by AIRS and CERES. J. Geophys. Res. Atmos., 117, D23107, doi: 10.1029/2012JD017997.

[6] Donohoe, A., Armour, K. C., Pendergrass, A. G., & Battisti, D. S. (2014). Shortwave and longwave radiative contributions to global warming under increasing CO2. Proceedings of the National Academy of Sciences, 111(47), 16700-16705. https://doi.org/10.1073/pnas.1412190111

[7] Su, W., N. G. Loeb, L. Liang, N. Liu, and C. Liu (2017), The El Niño–Southern Oscillation effect on tropical outgoing longwave radiation: A daytime versus nighttime perspective, J. Geophys. Res. Atmos., 122, 7820–7833, https://doi.org/10.1002/2017JD027002

[8] Soden, B. J., and I. M. Held, 2006: An Assessment of Climate Feedbacks in Coupled Ocean–Atmosphere Models. J. Climate, 19, 3354–3360, https://doi.org/10.1175/JCLI3799.1

[9] Loeb, N. G., J. M. Lyman, G. C. Johnson, R. P. Allan, D. R. Doelling, T. Wong, B. J. Soden, and G. L. Stephens, 2012: Observed changes in top-of-the-atmosphere radiation and upper-ocean heating consistent within uncertainty. Nat. Geosci., 5, 110–113, https://doi.org/10.1038/ngeo1375

[10] von Schuckmann, K., and Coauthors, 2016: An imperative to monitor Earth’s energy balance. Nat. Climate Change, 6, 138–144, https://doi.org/10.1038/nclimate2876

[11] Schwarz, M., D. Folini, S. Yang, and M. Wild, 2019: The Annual Cycle of Fractional Atmospheric Shortwave Absorption in Observations and Models: Spatial Structure, Magnitude, and Timing. J. Climate, 32, 6729–6748, https://doi.org/10.1175/JCLI-D-19-0212.1

[12] Yang, S., Zhang, X., Guan, S., Zhao, W., Duan, Y., Yao, Y., … Jiang, B. (2023). A review and comparison of surface incident shortwave radiation from multiple data sources: satellite retrievals, reanalysis data and GCM simulations. International Journal of Digital Earth, 16(1), 1332–1357. https://doi.org/10.1080/17538947.2023.2198262

[13] von Schuckmann, K. and Coauthors, 2023: Heat stored in the Earth system 1960–2020: where does the energy go? Earth System Science Data, 15(4), 1675-1709, https://doi.org/10.5194/essd-15-1675-2023