Introduction and Guide¶

This is a notebook to demonstrate the RetrievalToolbox software and how to create a trace gas retrieval based on measurements from NASA's EMIT instrument. Users can SHIFT-ENTER through each of the cells below. This demo sets up all the required inputs and then lets the user manually perform a few iterations of the inverse problem to adjust a state vector which best fits the measurement.

In order to run this demo, users have to download a few additional files (spectroscopy, solar model, EMIT L1B). Instructions to download those files are available on the Github page.

In [1]:
using Pkg; Pkg.activate();

using CSV
using DataFrames
using Dates
using HDF5
using Interpolations
using LinearAlgebra
using LoopVectorization
using NCDatasets
using Plots
default(
    titlefontsize = 10,
    legendfontsize = 8,
    guidefontsize = 8,
    tickfontsize = 8,
    
    )
using Polynomials
using Printf
using ProgressMeter
using Statistics
using StatsBase
using Unitful

using RetrievalToolbox; const RE = RetrievalToolbox

include("helpers.jl")
include("forward_model.jl")

  Activating project at `~/.julia/environments/v1.11`
┌ Warning: XRTM_PATH environment variable was not found! Not loading XRTM!
└ @ RetrievalToolbox ~/Work/ghgc/RetrievalToolbox.jl/src/RetrievalToolbox.jl:82
┌ Warning: Could not find XRTM module!
└ @ RetrievalToolbox ~/Work/ghgc/RetrievalToolbox.jl/src/RetrievalToolbox.jl:117
┌ Warning: Calls to XRTM library functions will crash the session!
└ @ RetrievalToolbox ~/Work/ghgc/RetrievalToolbox.jl/src/RetrievalToolbox.jl:118
Out[1]:
forward_model! (generic function with 1 method)

Preparation steps¶

Read inputs¶

In [2]:
# Read the noise coefficients from TXT file
noise_csv = CSV.File(
    "data/emit_noise.txt", skipto=2,
    header=["wl", "a", "b", "c", "rmse"]
);

USERS MUST EDIT THIS LINE IN CASE THEY DOWNLOAD A DIFFERENT EMIT L1B FILE

In [3]:
# Open the L1B file
nc_l1b = NCDataset(joinpath("L1B", "EMIT_L1B_RAD_demo.nc"));
In [4]:
# Read wavelength
wavelengths = nc_l1b.group["sensor_band_parameters"]["wavelengths"].var[:];

# Produce linear noise interpolation for each coeff a,b,c
noise_itps = Dict(
    "a" => linear_interpolation(noise_csv.wl, noise_csv.a, extrapolation_bc=Line()),
    "b" => linear_interpolation(noise_csv.wl, noise_csv.b, extrapolation_bc=Line()),
    "c" => linear_interpolation(noise_csv.wl, noise_csv.c, extrapolation_bc=Line()),
);
    
# Number of bands
Npix = length(wavelengths)
Out[4]:
285
In [5]:
# Read ISRF FWHMs (needed later for creation of ISRFs)
fwhms = nc_l1b.group["sensor_band_parameters"]["fwhm"].var[:];
In [6]:
# This is the start time of the EMIT L1B
time_coverage_start = nc_l1b.attrib["time_coverage_start"]
Out[6]:
"2023-06-12T16:21:03+0000"
In [7]:
# Read all Lat/Lon/Alt
nc_lon = replace(nc_l1b.group["location"]["lon"][:,:], missing => NaN);
nc_lat = replace(nc_l1b.group["location"]["lat"][:,:], missing => NaN);
nc_alt = replace(nc_l1b.group["location"]["elev"][:,:], missing => NaN);

# Read radiances!
# needs .var to ignore missing
nc_rad = replace(nc_l1b["radiance"][:,:,:], missing => NaN);
In [8]:
# Create the solar model
solar_model = RE.TSISSolarModel(
    "data/hybrid_reference_spectrum_p005nm_resolution_c2022-11-30_with_unc.nc",
    spectral_unit=:Wavelength
);
In [9]:
# Create the spectroscopy objects

ABSCO_CH4 = RE.load_ABSCOAER_spectroscopy(
    "data/CH4_03900-05250_v0.0_init.nc";
    spectral_unit=:Wavelength, distributed=true
);

ABSCO_CO2 = RE.load_ABSCOAER_spectroscopy(
    "data/CO2_03900-05250_v0.0_init.nc";
    spectral_unit=:Wavelength, distributed=true
);

ABSCO_H2O = RE.load_ABSCOAER_spectroscopy(
    "data/H2O_03900-05250_v0.0_init.nc";
    spectral_unit=:Wavelength, distributed=true
);

Prepare required objects¶

In [10]:
#=
    Retrieval windows
    =================
=#

window_dict = Dict{String, RE.SpectralWindow}()
Out[10]:
Dict{String, SpectralWindow}()

Users can make a choice regarding the so-called retrieval windows: to retrieve only the methane-affected portions of the spectrum, execute the next cell which will add a spectral window named CH4 which covers the relevant part of the spectrum. The cell following that contains a spectral window that contains CO$_2$ absorption, named CO2.

Users may choose to retrieve either CH$_4$, CO$_2$ or both; the retrieval application will adjust dynamically. Changing the retrieval windows after doing a first run is generally possible, however all subsequent cells must be executed before doing another retrieval.

In [11]:
window_dict["CH4"] = RE.spectralwindow_from_ABSCO(
    "CH4",
    2225.0, # min
    2350.0, # max
    2300.0, # reference
    35.0, # buffer (for ISRF convolution)
    ABSCO_CH4, # Spectroscopy from which to get the wavelength grid from
    u"nm"
);
In [12]:
window_dict["CO2"] = RE.spectralwindow_from_ABSCO(
    "CO2",
    1925.0,
    2125.0,
    2000.0,
    35.0,
    ABSCO_CH4,
    u"nm"
);
In [13]:
#=
    Dispersion objects
    ==================
=#

pixels = collect(1:Npix)
# Make a polynomial fit to the wavelengths from the L1B file

disp_fit = Polynomials.fit(pixels, wavelengths, 2) # order >2 does not seem to work nicely!
disp_coeffs = disp_fit.coeffs .|> Float64 # Get the coeffs

dispersion_dict = Dict{RE.SpectralWindow, RE.SimplePolynomialDispersion}()
for (wname, swin) in window_dict

    dispersion_dict[swin] = RE.SimplePolynomialDispersion(
        copy(disp_coeffs) * u"nm",
        1:Npix,
        swin # Reference this spectral window
    )

end
In [14]:
for (swin, disp) in dispersion_dict
    println("Dispersion coefficients for `$(swin)`:")
    display(disp.coefficients)
end

Dispersion coefficients for `SpectralWindow: CO2`:

3-element Vector{Float64}:
 372.4044494628906
   7.4781622886657715
  -0.00012527835497166961
Dispersion coefficients for `SpectralWindow: CH4`:

3-element Vector{Float64}:
 372.4044494628906
   7.4781622886657715
  -0.00012527835497166961
In [15]:
#=
    Instrument Spectral Response Function
    =====================================

    Build the ISRF table - we need a table since the ISRF changes with band/pixel

=#


Ndelta = 200
wl_delta_unit = u"nm"
wl_delta = zeros(Ndelta, Npix)
rr = zeros(Ndelta, Npix)

for i in 1:Npix

    σ = FWHM_to_sigma(fwhms[i]) # this is in nm

    wl_delta[:,i] = collect(LinRange(-5*σ, 5*σ, Ndelta))
    rr[:,i] = @. 1 / (σ * sqrt(2*pi)) * exp(-0.5 * ( - wl_delta[:,i])^2 / σ^2)
end

isrf = RE.TableISRF(
    wl_delta, # Δλ
    wl_delta_unit, # unit of Δλ
    rr # Relative response
)

# Create the dispersion => ISRF dictionary
isrf_dict = Dict(disp => isrf for (swin, disp) in dispersion_dict)
Out[15]:
Dict{SimplePolynomialDispersion{SpectralWindow{Float64}, Float64, UnitRange{Int64}}, TableISRF{Float64}} with 2 entries:
  SimplePolynomialDispersion{SpectralWindow{Float64}, … => TableISRF: (200, 285)
  SimplePolynomialDispersion{SpectralWindow{Float64}, … => TableISRF: (200, 285)
In [16]:
p = plot()
for i in [1,100,285]
    plot!(p,
        isrf.ww_delta[:,i],
        isrf.relative_response[:,i],
        #yscale=:log10,
        label="Pos. $(i)"
    )
end
xlabel!("Δλ [$(isrf.ww_delta_unit)]")
ylabel!("Relative response")
Out[16]:
No description has been provided for this image
In [17]:
#=
    Gases
    =====

    Take representative profiles from the shipped atmosphere file (H2O, CH4, CO2 only)
=#

# Some user defined function to generate a pressure grid:
function generate_plevels(psurf)
    return vcat(
        collect(LinRange(0.001u"hPa", 200.0u"hPa", 4)),
        collect(LinRange(300.0u"hPa", psurf, 6))
    )
end

plevels = generate_plevels(1000.0u"hPa")
N_RT_lev = length(plevels)

gases = RE.GasAbsorber[]

gas_ch4 = RE.create_example_gas_profile("./data/EMIT-example.csv", "CH4",
    ABSCO_CH4, plevels; is_example=false
)
ch4_prior = copy(gas_ch4.vmr_levels)

gas_h2o = RE.create_example_gas_profile("./data/EMIT-example.csv", "H2O",
    ABSCO_H2O, plevels; is_example=false
)
h2o_prior = copy(gas_h2o.vmr_levels)

gas_co2 = RE.create_example_gas_profile("./data/EMIT-example.csv", "CO2",
    ABSCO_CO2, plevels;  is_example=false)
co2_prior = copy(gas_co2.vmr_levels)

push!(gases, gas_ch4);
push!(gases, gas_h2o);
push!(gases, gas_co2);
In [18]:
#=
    Atmosphere
    ==========
=#

# Use `atm_orig` as the original reference atmosphere
atm_orig = RE.create_example_atmosphere("./data/EMIT-example.csv", N_RT_lev;
    T=Float64, is_example=false
);

# Use `atm` as a working copy
atm = deepcopy(atm_orig);
N_MET_lev = atm.N_met_level;

# Ingest the retrieval grid..
RE.ingest!(atm, :pressure_levels, plevels);

# Make sure we have layer quantites everywhere
# (they are calculated from level quantities)
RE.calculate_layers!(atm);

# Add gases to atmosphere
push!(atm.atm_elements, gases...);
In [19]:
# Plot the gases in the atmosphere
p = []
for gas in gases
    _plot = plot(
        gas.vmr_levels,
        atm.pressure_levels,
        label=nothing,
        #label="$(gas.gas_name)",
        marker=:o,
        #legend=:bottomleft,
    )
    yflip!(_plot)
    ylabel!(_plot, "Pressure level [$(atm.pressure_unit)]")
    xlabel!(_plot, "$(gas.gas_name) [$(gas.vmr_unit)]")
    push!(p, _plot)
end
plot(p...)
Out[19]:
No description has been provided for this image

State Vector set-up¶

In [20]:
sv_ch4_scaler = RE.GasLevelScalingFactorSVE(
    1,
    N_RT_lev,
    gas_ch4,
    Unitful.NoUnits,
    1.0,
    1.0,
    1.0
)

#Uncomment for CO2
sv_co2_scaler = RE.GasLevelScalingFactorSVE(
    1,
    N_RT_lev,
    gas_co2,
    Unitful.NoUnits,
    1.0,
    1.0,
    1.0
)


sv_h2o_scaler = RE.GasLevelScalingFactorSVE(
    1,
    N_RT_lev,
    gas_h2o,
    Unitful.NoUnits,
    1.0,
    1.0,
    1.0e-2
)

# Retrieve a polynomial for the Lambertian surface albedo
sv_surf = RE.SurfaceAlbedoPolynomialSVE[]
surf_order = 2

for (win_name, swin) in window_dict
    for o in 0:surf_order

        o == 0 ? fg = 0.25 : fg = 0.0

        push!(sv_surf,
            RE.SurfaceAlbedoPolynomialSVE(
                swin,
                o,
                u"nm",
                fg,
                fg,
                1.0^(-2*o)
            )
        )

    end
end

# Construct the state vector
state_vector = RE.RetrievalStateVector([
    sv_ch4_scaler,
    sv_h2o_scaler,
    sv_co2_scaler,
    sv_surf..., # Expand list
])
Out[20]:
                  State Vector (current)
┌────┬──────────────────────────────────────┬──────┬──────┐
│ #1 │ GasLevelScalingFactorSVE CH4 (1, 10) │  1.0 │      │
│ #2 │ GasLevelScalingFactorSVE H2O (1, 10) │  1.0 │      │
│ #3 │ GasLevelScalingFactorSVE CO2 (1, 10) │  1.0 │      │
│ #4 │ SurfaceAlbedoPolynomialSVE (0) [CH4] │ 0.25 │      │
│ #5 │ SurfaceAlbedoPolynomialSVE (1) [CH4] │  0.0 │ nm⁻¹ │
│ #6 │ SurfaceAlbedoPolynomialSVE (2) [CH4] │  0.0 │ nm⁻² │
│ #7 │ SurfaceAlbedoPolynomialSVE (0) [CO2] │ 0.25 │      │
│ #8 │ SurfaceAlbedoPolynomialSVE (1) [CO2] │  0.0 │ nm⁻¹ │
│ #9 │ SurfaceAlbedoPolynomialSVE (2) [CO2] │  0.0 │ nm⁻² │
└────┴──────────────────────────────────────┴──────┴──────┘
In [21]:
#=
    Buffer
    ======
=#

N1 = 50_000 # Number of spectral points needed for monochromatic calculations, e.g. convolution
N2 = Npix # Number of spectral points at instrument level needed

my_type = Float64

# Will contain outputs of ISRF application
inst_buf = RE.InstrumentBuffer(
    zeros(my_type, N1),
    zeros(my_type, N1),
    zeros(my_type, N2),
)

# Buffer needed for the monochromatic radiance calculations
rt_buf = RE.ScalarRTBuffer(
    dispersion_dict, # Already a SpectralWindows -> Dispersion dictionary
    RE.ScalarRadiance(my_type, N2), # Hold the radiance - we use ScalarRadiance because we don't need polarization
    Dict(sve => RE.ScalarRadiance(my_type, N2) for sve in state_vector.state_vector_elements),
    Dict(swin => zeros(Int, 0) for swin in values(window_dict)), # Hold the detector indices
    u"mW * cm^-2 * sr^-1 * nm^-1" # Radiance units for the forward model
)

# Create the EarthAtmospherBuffer using this helper function rather than doing it manually
buf = RE.EarthAtmosphereBuffer(
    state_vector, # The state vector
    values(window_dict) |> collect, # The spectral window (or a list of them)
    [(:Lambert, 7) for x in window_dict], # Surfaces
    atm.atm_elements, # All atmospheric elements
    Dict(swin => solar_model for swin in values(window_dict)), # Solar model dictionary (spectral window -> solar model)
    [:BeerLambert for swin in window_dict], # Use the speedy Beer-Lambert RT model
    RE.ScalarRadiance, # Use ScalarRadiance for high-res radiance calculations
    rt_buf,
    inst_buf,
    N_RT_lev, # The number of retrieval or RT pressure levels
    N_MET_lev, # The number of meteorological pressure levels, as given by the atmospheric inputs
    my_type # The chosen Float data type (e.g. Float16, Float32, Float64)
);
In [22]:
# Supplement the scene object with sensible values

# Put our atmosphere in here
buf.scene.atmosphere = atm

# Set the time (roughly)
buf.scene.time = DateTime(split(time_coverage_start, "+")[1]);

# Set satellite observer
# (the EMIT L1B files do not have satellite angles in them, so we just pretend it's
#  all at nadir)
buf.scene.observer = RE.SatelliteObserver(
    0.0, # zenith angle (let's assume it's close to nadir)
    0.0, # azimuth angle (does not matter here)
    [0., 0., 0.], # Satellite position
    [0., 0., 0.] # Satellite velocity
);

Forward model set-up¶

In [23]:
#=
    Forward model definitions
=#

# Forward model keyword arguments (to be supplied with the forward model function)
fm_kwargs = (
    buf=buf,
    inst_buf=inst_buf,
    rt_buf=rt_buf,
    dispersions=dispersion_dict,
    isrf_dict=isrf_dict,
    solar_doppler_factor=nothing, # Let RetrievalToolbox calculate the solar Doppler shift
    solar_distance=1.0
    );
In [24]:
#=
    Noise model
=#

# EMIT noise model
function noise_model!(noise, rad, a, b, c)
    @turbo for i in eachindex(noise)
        noise[i] = abs(a[i] * sqrt(rad[i] + b[i]) + c[i])
    end

    # Bump up negative noise to this (as per reference algorithm)
    noise[noise .<= 0] .= 1e-5
end

# Pre-allocate a vector to contain the noise, so we don't have to do it over and over
this_noise = zeros(size(nc_rad, 1))

# Pre-calculate the a,b,c as function of wavelength!
noise_a = noise_itps["a"].(wavelengths);
noise_b = noise_itps["b"].(wavelengths);
noise_c = noise_itps["c"].(wavelengths);
In [25]:
# Function to calculate pressure from elevation (valid only for this model atmosphere)
# (we want this to obtain the surface pressure given the L1B elevation)
itp_logp_from_alt = linear_interpolation(
    reverse(buf.scene.atmosphere.altitude_levels),
    log10.(reverse(buf.scene.atmosphere.met_pressure_levels)),
    extrapolation_bc=Line()
);

# Function to calculate altitude from pressure
# (we want this to produce the column-integrated CH4)
# Note that we don't reverse here since the knots need to be in ascending order
# and pressure levels increase when going towards the surface.
itp_alt_from_logp = linear_interpolation(
    log10.(buf.scene.atmosphere.met_pressure_levels),
    buf.scene.atmosphere.altitude_levels,
    extrapolation_bc=Line()
);

# Before going into the scene loop, let's establish the solar strength for each
# spectral window, so we can calculate a good initial guess for the surface albedo.
solar_strength_guess = Dict{RE.SpectralWindow, Float64}()

for (swin, rt) in buf.rt
    solar_idx = searchsortedfirst.(Ref(rt.solar_model.ww), swin.ww_grid[:] / 1000.0)
    solar_strength_guess[swin] = maximum(rt.solar_model.irradiance[solar_idx])
end

Perform the retrieval¶

In [26]:
# Pick a scene index here

idx = CartesianIndex(237, 1649)
Out[26]:
CartesianIndex(237, 1649)
In [27]:
# We must first take all the per-scene data from the various shared arrays.
this_meas = @view nc_rad[:,idx] # Radiance / measurement
this_lon = nc_lon[idx]
this_lat = nc_lat[idx]
this_alt = nc_alt[idx]

# Adjust surface pressure for this scene to follow the terrain..
this_psurf = 10^(itp_logp_from_alt(this_alt));
In [28]:
# Move new retrieval grid into buffer
new_plevels = generate_plevels(this_psurf * buf.scene.atmosphere.pressure_unit)
RE.ingest!(buf.scene.atmosphere, :pressure_levels, new_plevels)

# Calculate layer quantities from level quantities
RE.calculate_layers!(buf.scene.atmosphere)

# Set the scene location
loc = RE.EarthLocation(
    this_lon,
    this_lat,
    this_alt,
    u"m"
)
buf.scene.location = loc
Out[28]:
EarthLocation
Longitude: -111.09999412186187
Latitude:  41.24026978970571
Altitude:  2063.561698153218 m
In [29]:
# Calculate solar angles from the location and time
RE.update_solar_angles!(buf.scene)
println("Solar zenith  = $(buf.scene.solar_zenith)")
println("Solar azimuth = $(buf.scene.solar_azimuth)")

Solar zenith  = 47.61081360746263
Solar azimuth = 100.57536348382311
In [30]:
# Set the gas profiles back to their original prior state
gas_co2.vmr_levels[:] .= co2_prior;
gas_ch4.vmr_levels[:] .= ch4_prior;
gas_h2o.vmr_levels[:] .= h2o_prior;
In [31]:
# Calculate the noise (in-place)!
noise_model!(this_noise, this_meas,
    noise_itps["a"].(wavelengths),
    noise_itps["b"].(wavelengths),
    noise_itps["c"].(wavelengths),
);
In [32]:
# Create solver
solver = RE.IMAPSolver(
    forward_model!,
    state_vector,
    Diagonal(RE.get_prior_covariance(state_vector)), # Prior covariance matrix - just use diagonal
    20, # number of iterations
    1.0, # dsigma scale
    dispersion_dict,
    rt_buf.indices,
    rt_buf.radiance,
    rt_buf.jacobians,
    Dict(disp => this_meas for disp in values(dispersion_dict)), # the measurement (full spectrometer)
    Dict(disp => this_noise for disp in values(dispersion_dict)) # Noise (full spectrometer)
);
In [33]:
# Adjust surface reflectance first guess for every retrieval window
for (swin, rt) in buf.rt

    # Find out any unit conversion factor between measured radiance (rt_buf.radiance_unit)
    # and the radiance we are using internally.
    rad_unit_fac = 1.0 * rt_buf.radiance_unit / rt.radiance_unit |> upreferred

    # Calculate apparent albedo from the measured radiances
    signal = percentile(RE.get_measured(solver, swin; view=true), 99)
    albedo_prior = pi * signal / (
        solar_strength_guess[swin] * cosd(buf.scene.solar_zenith)) * rad_unit_fac
    
    for (sve_idx, sve) in RE.StateVectorIterator( # loop through all albedo SVEs
        state_vector, RE.SurfaceAlbedoPolynomialSVE)

        # Skip non-matching state vector element
        sve.swin != swin && continue
        
        if sve.coefficient_order == 0
            sve.first_guess = albedo_prior
            sve.prior_value = albedo_prior
            sve.iterations[1] = albedo_prior

            println("Setting albedo prior for $(swin) to $(albedo_prior)")
        end
    end
end

Setting albedo prior for SpectralWindow: CO2 to 0.11092285670307869
Setting albedo prior for SpectralWindow: CH4 to 0.16081745179442844
In [34]:
# Re-set the state vector to first guess values (empty the iterations)
for sve in state_vector.state_vector_elements
    empty!(sve.iterations)
    push!(sve.iterations, sve.first_guess)
end
state_vector
Out[34]:
                    State Vector (current)
┌────┬──────────────────────────────────────┬──────────┬──────┐
│ #1 │ GasLevelScalingFactorSVE CH4 (1, 10) │      1.0 │      │
│ #2 │ GasLevelScalingFactorSVE H2O (1, 10) │      1.0 │      │
│ #3 │ GasLevelScalingFactorSVE CO2 (1, 10) │      1.0 │      │
│ #4 │ SurfaceAlbedoPolynomialSVE (0) [CH4] │ 0.160817 │      │
│ #5 │ SurfaceAlbedoPolynomialSVE (1) [CH4] │      0.0 │ nm⁻¹ │
│ #6 │ SurfaceAlbedoPolynomialSVE (2) [CH4] │      0.0 │ nm⁻² │
│ #7 │ SurfaceAlbedoPolynomialSVE (0) [CO2] │ 0.110923 │      │
│ #8 │ SurfaceAlbedoPolynomialSVE (1) [CO2] │      0.0 │ nm⁻¹ │
│ #9 │ SurfaceAlbedoPolynomialSVE (2) [CO2] │      0.0 │ nm⁻² │
└────┴──────────────────────────────────────┴──────────┴──────┘

Manually iterate!

Below cell will execute perform an iteration - meaning the forward model will be evaluated and a state vector update is computed from the forward model result (through the difference to the measurement). The new state vector values are added to each state vector element, such that the current state vector value represents that updated state.

Note that the first time this cell is executed in a fresh session will cause a lot of underlying code to be compiled, and hence will take ~30s. Subsequent calls to the forward model will generally only take ~50ms.

Convergence should be achieved in ~3-5 iterations. To re-set the state vectors, execute the cell above.

In [38]:
# Repeat this step until a desirable fit is achieved..
@time RE.next_iteration!(solver; fm_kwargs)


# Restrict all gas scalers > 0.01
# (this is mostly to keep the H2O scale factor > 0)
for (sve_idx, sve) in RE.StateVectorIterator(
    state_vector, RE.GasLevelScalingFactorSVE)

    current = RE.get_current_value(sve)

    if (current < 0.01)
        sve.iterations[end] = 0.01
    end
end


forward_model!(state_vector; fm_kwargs...)

@info "Converged? $(RE.check_convergence(solver))"

p = plot()

chi2 = RE.calculate_chi2(solver);

for (cnt, (swin, disp)) in enumerate(dispersion_dict)
    
    meas = RE.get_measured(solver, swin);
    conv = RE.get_modeled(solver, swin);

    plot!(disp.detector_samples[disp.index], meas, linewidth=5,
            label="Measurement $(cnt)", color=:black)
    plot!(disp.detector_samples[disp.index], conv, linewidth=2,
            marker=:o, label="Fit $(cnt)", color=:red, legend=:bottom)

end
xlabel!("Spectral sample / channel");
ylabel!("Radiance\n[$(buf.rt_buf.radiance_unit)]")
tstring = ""
for (k,v) in chi2
    tstring *= @sprintf "%s, χ\$^2\$ = %0.4f\n" k v
end
title!(tstring, fontsize=6);

display(p)

# Display the state vector
display(state_vector)

  0.124465 seconds (170.90 k allocations: 36.160 MiB)

[ Info: Converged? true
No description has been provided for this image 
                      State Vector (current)
┌────┬──────────────────────────────────────┬──────────────┬──────┐
│ #1 │ GasLevelScalingFactorSVE CH4 (1, 10) │     0.770976 │      │
│ #2 │ GasLevelScalingFactorSVE H2O (1, 10) │     0.838197 │      │
│ #3 │ GasLevelScalingFactorSVE CO2 (1, 10) │       1.0148 │      │
│ #4 │ SurfaceAlbedoPolynomialSVE (0) [CH4] │     0.143683 │      │
│ #5 │ SurfaceAlbedoPolynomialSVE (1) [CH4] │ -0.000346305 │ nm⁻¹ │
│ #6 │ SurfaceAlbedoPolynomialSVE (2) [CH4] │   2.28398e-6 │ nm⁻² │
│ #7 │ SurfaceAlbedoPolynomialSVE (0) [CO2] │     0.146061 │      │
│ #8 │ SurfaceAlbedoPolynomialSVE (1) [CO2] │  0.000119099 │ nm⁻¹ │
│ #9 │ SurfaceAlbedoPolynomialSVE (2) [CO2] │   4.56715e-7 │ nm⁻² │
└────┴──────────────────────────────────────┴──────────────┴──────┘
In [39]:
# Show some posterior analysis
q = RE.calculate_OE_quantities(solver);
RE.print_posterior(q)

                                               Posterior state vector
┌──────────────────────────────────────┬───────┬──────────────┬─────────────┬─────────────┬─────────────┬──────────┐
│                                 Name │ Units │        Value │ Uncertainty │ Uncertainty │ Uncertainty │       AK │
│                                      │       │              │     (total) │ (smoothing) │     (noise) │          │
├──────────────────────────────────────┼───────┼──────────────┼─────────────┼─────────────┼─────────────┼──────────┤
│ GasLevelScalingFactorSVE CH4 (1, 10) │       │     0.770976 │   0.0687225 │  0.00495206 │   0.0685438 │ 0.995277 │
│ GasLevelScalingFactorSVE H2O (1, 10) │       │     0.838197 │   0.0207111 │  0.00429214 │   0.0202615 │ 0.957105 │
│ GasLevelScalingFactorSVE CO2 (1, 10) │       │       1.0148 │  0.00990371 │ 0.000199172 │   0.0099017 │ 0.999902 │
│ SurfaceAlbedoPolynomialSVE (0) [CH4] │       │     0.143683 │  0.00130415 │  8.71199e-5 │  0.00130124 │ 0.999998 │
│ SurfaceAlbedoPolynomialSVE (1) [CH4] │  nm⁻¹ │ -0.000346305 │  8.13553e-6 │   7.1073e-7 │  8.10443e-6 │      1.0 │
│ SurfaceAlbedoPolynomialSVE (2) [CH4] │  nm⁻² │   2.28398e-6 │  1.76512e-7 │  1.24388e-8 │  1.76073e-7 │      1.0 │
│ SurfaceAlbedoPolynomialSVE (0) [CO2] │       │     0.146061 │ 0.000922116 │ 0.000170098 │ 0.000906292 │ 0.999999 │
│ SurfaceAlbedoPolynomialSVE (1) [CO2] │  nm⁻¹ │  0.000119099 │  1.91717e-5 │  3.41837e-6 │  1.88645e-5 │      1.0 │
│ SurfaceAlbedoPolynomialSVE (2) [CO2] │  nm⁻² │   4.56715e-7 │  1.23661e-7 │  1.85625e-8 │  1.22259e-7 │      1.0 │
└──────────────────────────────────────┴───────┴──────────────┴─────────────┴─────────────┴─────────────┴──────────┘
In [40]:
# Posterior covariance
q.Shat
Out[40]:
9×9 Matrix{Float64}:
  0.00472278   -0.000148685   3.40917e-6   …   1.11838e-7   -5.5759e-10
 -0.000148685   0.00042895   -1.733e-5        -3.4163e-7     1.85484e-9
  3.40917e-6   -1.733e-5      9.80834e-5       3.85463e-8   -4.55803e-10
  8.67709e-5   -7.60191e-7   -1.72739e-8       4.83614e-10  -1.70711e-12
  2.72265e-7    6.56482e-8   -2.8225e-9       -5.27015e-11   2.89287e-13
 -6.03576e-9    1.08753e-9   -4.0722e-11   …  -8.58152e-13   4.59896e-15
 -5.97719e-6    1.69977e-5    2.14735e-6      -1.41683e-8    6.65657e-11
  1.11838e-7   -3.4163e-7     3.85463e-8       3.67555e-10  -2.23936e-12
 -5.5759e-10    1.85484e-9   -4.55803e-10     -2.23936e-12   1.52919e-14

Correlation matrix

The cell below produces a correlation matrix derived from the posterior covariance. It reveals how "confounding" state vector element pairs are with respect to others, and includes the effect of instrument noise and prior covarince. The correlation matrix is symmetric, and for easier reading, we only produce the upper triangle here.

A value closer to -1 or 1 suggests that the retrieval "cannot distinguish" between adjusting one over the other as a means to minimize the cost function. Highly correlated state vector elements will thus result in systematic biases. In our case, observe the correlation between the zero-order surface albedo SurfaceAlbedoPolynomialSVE (0) and any of the gas profile scalers, such as GasLevelScalingFactorSVE CH4. This analysis shows that the inverse method can easily distinguish between making the surface slightly brighter, or increasing the gas concentration.

Note however, that a high correlation does not imply that the effect is large, only that it is likely systematic. Due to the presense of instrument noise, it might require a larger set of retrievals to show the true impact of this correlation.

In [41]:
# Produce a correlation matrix from Shat!
C = similar(q.Shat)
for i in 1:size(C,1), j in 1:size(C,2)
    if i < j 
        C[i,j] = q.Shat[i,j] / sqrt(q.Shat[i,i] * q.Shat[j,j])
    else
        C[i,j] = NaN
    end
    
end
xl = [split("$(sv)",":")[1] for sv in state_vector.state_vector_elements]
heatmap(xl, xl, C, c=:PuOr, clim=(-1,1), xrotation=70, size=(650, 600), aspect=1)
yflip!()

for i in 1:size(C,1), j in 1:size(C,2)
    i == j && continue
    if abs(C[i,j]) < 0.5
        tcolor = :black
    else
        tcolor = :white
    end
    annotate!(xl[j], xl[i], text((@sprintf "%.2f" C[i,j]), tcolor, 7))
end

title!("Posterior correlation matrix")
Out[41]:
No description has been provided for this image
In [42]:
# Plot the high-res model along with the at-instrument radiance
p_list = []
for (swin, rt) in buf.rt

    conv = RE.get_modeled(solver, swin);
    
    p = plot()
    
    plot!(p, swin.ww_grid, rt.hires_radiance, label="High-res model")
    
    #plot!(p, swin.ww_grid, rt.hires_solar, label="Solar model")
    plot!(p,
        RE.get_wavelength(solver, swin),
        conv * ((1.0 * buf.rt_buf.radiance_unit / rt.radiance_unit) |> upreferred),
        linewidth=5, label="After instrument model"
    )
    title!("$(swin)")
    xlabel!("Wavelength [$(swin.ww_unit)]")
    ylabel!("Radiance\n[$(buf.rt_buf.radiance_unit)]")
    push!(p_list, p)
end

plot(p_list..., size=(800, 500))
Out[42]:
No description has been provided for this image