IASI

The Infrared Atmospheric Sounding Interferometer (IASI) is an instrument on the METOP satellites. It is a hyper-spectral infrared sensor measuring upwelling radiation from a sun-synchronous orbit in 8461 spectral channels (645.0 - 2760.0 cm-1). These measurements are used to derive a plethora of geophysical variables (e.g. temperature and humidity profiles). This makes IASI a key data source for numerical weather prediction (NWP) and applications in atmospheric chemistry and monitoring of essential climate variables. See IASI Level 1 product page, IASI Level 1 Principal Component Scores product page and IASI Level 2 product page for more information.

Static plot of L1C spectra

This example is made using the following packages.

[13f3f980] CairoMakie v0.12.16
[db073c08] GeoMakie v0.7.9

The key variable is "gs1cspect" which contains the radiance spectra measured by IASI. The spectra from a full orbit are almost 2 GB of data. In this example we will just load one data record with observation locations and plot the spectra of two observations. The spectra are converted from radiances to brightness temperature since this is often convenient when interpreting the IASI spectra. The cloud cover in each individual field of view (from AVHRR) is also read from the file and included as legends on the plot of the spectra.

using MetopDatasets
using CairoMakie, GeoMakie

ds = MetopDataset("IASI_xxx_1C_M01_20240925202059Z_20240925220258Z_N_O_20240925211316Z.nat", maskingvalue=NaN);

# Select a single data record
data_record_index = 670

# Select 2 points to plot the full spectrum 
selected_points = CartesianIndex.([(2,15), (1,27)])
selected_colors = [:red, :black]

# Read the geolocation of the data record
longitude, latitude = let 
    lon_lat = ds["ggeosondloc"][:,:,:,data_record_index]
    lon_lat[1,:,:], lon_lat[2,:,:]
end

# read spectrum as brightness temperature
wavenumber = ds["spectra_wavenumber"][:, data_record_index]
selected_T_b = let 
    spectra = ds["gs1cspect"][:,:,:,data_record_index]
    # convert select spectra to brightness temperature
    [brightness_temperature.(spectra[:,i] , wavenumber) 
        for i in selected_points]
end

# Read cloud cover for the selected points
selected_cloud_cover = ds["geumavhrr1bcldfrac"][:,:,data_record_index][selected_points]

# plot the data
fig = let

    fig = Figure()

    # axis to plot geolocation
    ax = GeoAxis(fig[1, 1],
        title = "Observations",
        xlabel = "longitude",
        ylabel = "latitude",
        limits = (extrema(longitude), extrema(latitude)))

    # plot all observations from data record in gray
    scatter!(ax, longitude[:], latitude[:], color=:gray)
    # plot selected observations in color
    scatter!(ax, longitude[selected_points], latitude[selected_points], 
        color=selected_colors, marker=:xcross, markersize = 15)

    # Add coastlines
    lines!(ax, GeoMakie.coastlines()) 

    # Plot the selected spectra
    ax2 = Axis(fig[2, 1],
        title = "Selected spectra",
        xlabel = "Wavenumber (cm-1)",
        ylabel = "Brightness Temperature (K)",
        limits=((500,3500), (150,300)))

    for i in eachindex(selected_points)
        cloud_cover_i = Int(selected_cloud_cover[i])
        lines!(ax2, wavenumber./100, selected_T_b[i], color = selected_colors[i],
            label = "$(cloud_cover_i)% cloud cover")
    end

    # Add legends
    axislegend(ax2, position = :rb)

    fig
end

The top plot shows a row of observations west of the French coast. Two observation have been marked with a colored X. The spectrum of these observations are shown below.

Advanced interactive plots of L1C spectra

This example shows how the IASI spectra can be shown interactively on a responsive map. The example relies on Tyler.jl and GLMakie.jl. It is recommend to first try some examples from the Tyler documentation if you are new to this package.

[e9467ef8] GLMakie v0.10.16
[e170d443] Tyler v0.2.0

This example will use Tyler to create an interactive background map similar to google maps in the top panel of the figure. Then all the locations of the observation will be plotted on the map colored based on the cloud cover. The bottom panel will display the spectrum of a single observation. Clicking an observation in the top panel will load that spectrum and plot it in the bottom panel. The loading of the spectrum will be done lazily to conserve memory.

using Tyler, GLMakie, MetopDatasets

ds = MetopDataset("IASI_xxx_1C_M01_20240925202059Z_20240925220258Z_N_O_20240925211316Z.nat", maskingvalue=NaN);

# read geolocation points and data shape
pts, pts_original_size = let
    lon_lat = ds["ggeosondloc"][:,:,:,:]

    # make sure longitude is from -180 to 180
    lon_lat[1, :,:,:][ lon_lat[1,:,:,:] .>180] .-= 360

    # convert the points to tuples
    lon_lat = tuple.(lon_lat[1, :,:,:],lon_lat[2, :,:,:])

    # store the original shape of the points
    pts_original_size = size(lon_lat)

    # Flatten the points and convert them to web_mercator (the coordinate system used by Tyler)
    MT = Tyler.MapTiles
    pts = [Point2f(MT.project(lon_lat[i], MT.wgs84, MT.web_mercator)) 
            for i in eachindex(lon_lat)]
    
    pts, pts_original_size
end;

# read cloud fraction
cloud_fraction = Float32.(ds["geumavhrr1bcldfrac"][:]);

# helper function to read the spectrum for a single point
function read_spectrum_pts(ds, index::CartesianIndex)
    # read spectrum and wavenumber
    spectrum = ds["gs1cspect"][:,Tuple(index)...]
    wavenumber = ds["spectra_wavenumber"][:, Tuple(index)[end]]

    # covert to brightness temperature
    T_B = brightness_temperature.(spectrum, wavenumber)
    wavenumber_cm = wavenumber./100

    # join brightness temperature and wavenumber to points
    spectrum_pts = Point2f.(tuple.(wavenumber_cm, T_B))
    return spectrum_pts
end

# read an initial spectrum
spectrum_pts = read_spectrum_pts(ds, CartesianIndex(1,1,1));

# create the inter active plot.
fig = let
    fig = Figure()
    
    # select background map and initial zoom
    provider = Tyler.TileProviders.Esri(:WorldImagery);
    extent = Tyler.Extent(X = (-10, 10), Y = (-10, 10));

    # create background map
    ax1 = Axis(fig[1,1])
    m = Tyler.Map(extent; provider, figure=fig, axis=ax1);
    wait(m);

    # Plot observation points with cloud cover
    objscatter = scatter!(ax1, pts, color = cloud_fraction, 
        colorrange = (0,100), colormap=:grays, markersize=15)
    # hack from https://github.com/MakieOrg/Tyler.jl/issues/109
    translate!(objscatter, 0, 0, 10) 

    # Plot a red cross on top of the selected point
    selected_point = Observable(pts[1:1])
    selected_scatter = scatter!(ax1, selected_point, 
        color = :red, markersize=10, marker =:xcross)
    # hack from https://github.com/MakieOrg/Tyler.jl/issues/109
    translate!(selected_scatter, 0, 0, 11)

    # Add colorbar
    Colorbar(fig[1, 2], limits = (0.0,100.0), colormap =:grays, label = "Cloud Fraction")
    hidedecorations!(ax1)

    # Create the second plot for the spectrum
    ax2 = Axis(fig[2, 1],
        title = "Observation spectrum",
        xlabel = "Wavenumber (cm-1)",
        ylabel = "Brightness Temperature (K)")

    # plot the spectrum
    spectrum_observable = Observable(spectrum_pts)
    lines!(ax2, spectrum_observable)
    
    # Add event handler to update the plot when the user click on a new observation 
    obs_func = on(events(m.axis).mousebutton) do event
        if event.button == Mouse.left && event.action == Mouse.press
            plt, i = pick(m.axis)
            if plt == objscatter # check if an observation was clicked
                # get the CartesianIndex
                cartesian_i = CartesianIndices(pts_original_size)[i]
                # load the selected spectrum and update the spectrum plot
                spectrum_observable[] = read_spectrum_pts(ds, cartesian_i)
                # update the red x
                selected_point[] = pts[i:i]
            end
        end
    end

    fig
end

It is now possible to interactively explore the nearly 100 000 observations from an obit of IASI with the background map giving important context.

Using IASI Principal Component Scores

IASI radiance spectra are also available as Principal Component Scores (PCS). A single orbit of PCS data is approximately 100 MB, compared to an uncompressed IASI Level 1 orbit, which is around 2 GB.

MetopDataset can be used to read IASI PCS products and provides methods to reconstruct the full IASI spectra from the scores. Note that the IASI principal component vectors are required for spectral reconstruction. These vectors are included in this package as a 200 MB lazy artifact, which is downloaded automatically the first time reconstruct_iasi_spectrum is called. This approach provides a smooth user experience for IASI PCS users while avoiding unnecessary disk usage for users who do not work with IASI data.

The example below compares a reconstructed IASI spectrum with the original spectrum.

using MetopDatasets
using CairoMakie

ds      = MetopDataset("IASI_xxx_1C_M01_20260416213255Z_20260416231454Z_N_O_20260416222234Z.nat", maskingvalue=NaN)
ds_pcs  = MetopDataset("IASI_PCS_1C_M01_20260416213255Z_20260416231454Z_N_O_20260416222424Z.nat", maskingvalue=NaN)

# chose an observation to test.
data_record_index = 670
sounder_index = 3
xtrack_index = 8

# check the observation time and location match.
lon_lat_original = ds["ggeosondloc"][:, sounder_index, xtrack_index, data_record_index]
lon_lat_pcs = ds_pcs["ggeosondloc"][:, sounder_index, xtrack_index, data_record_index]
if !isapprox(lon_lat_original, lon_lat_pcs)
    error("Location not matching")
end

timestamp_original = ds["record_start_time"][data_record_index]
timestamp_pcs = ds_pcs["record_start_time"][data_record_index]
if timestamp_original != timestamp_pcs
    error("Timestamp not matching")
end

# Create variable for the reconstruct spectrum.
# Note: This downloads the 200 MB principal component vectors the first time.
reconstruct_spectra_variable = reconstruct_iasi_spectrum(ds_pcs)

# Indexing the variable reads the PCS and reconstructs the spectra.
observation_reconstructed = reconstruct_spectra_variable[:, sounder_index, xtrack_index, data_record_index]
# Reading the original spectrum.
observation_original = ds["gs1cspect"][:, sounder_index, xtrack_index, data_record_index]

# get wave numbers
wavenumber_cm_original = ds["spectra_wavenumber"][:, data_record_index]./100
wavenumber_cm_pcs = ds_pcs["spectra_wavenumber"][:]./100

# plot figure
let

    fig = Figure()
    ax1 = Axis(fig[1, 1],
        limits= (nothing,nothing),
        xlabel = "Wavenumber (cm-1)",
        ylabel = "Radiance (W/m2/sr/m-1)",
        )

    lines!(ax1, wavenumber_cm_original, observation_original, label = "original")
    lines!(ax1, wavenumber_cm_pcs, observation_reconstructed, label = "reconstruct")
    lines!(ax1, wavenumber_cm_pcs, 
        observation_reconstructed .- observation_original[1:length(observation_reconstructed)], label = "difference")
    axislegend(ax1, position = :rt)

    # add a second plot with zoom-in.
    ax2 = Axis(fig[2, 1],
        limits= ((1585,1600),(-5e-6,5e-5)),
        xlabel = "Wavenumber (cm-1)",
        ylabel = "Radiance (W/m2/sr/m-1)",
        )

    lines!(ax2, wavenumber_cm_original, observation_original, label = "original")
    lines!(ax2, wavenumber_cm_pcs, observation_reconstructed, label = "reconstruct")
    lines!(ax2, wavenumber_cm_pcs, 
        observation_reconstructed .- observation_original[1:length(observation_reconstructed)], label = "difference")

    fig
end

The top plot compares an IASI PCS reconstructed spectrum with the original Level 1 spectrum. The agreement is so good that the reconstructed curve completely overlaps the original. The bottom plot shows a zoomed-in view, where small differences between the original and reconstructed spectra become visible.

Level 2 Combined Sounding Products

The IASI level 2 products contains derived atmospheric profiles of temperature, water vapour, ozone and trace gases. The availability of these profiles depend on cloud cover and therefore the number of profiles will vary through out the product. These variables of changing size are padded with fill values (default to missing) to generate an array that fits with the MetopDataset interface. This example plots the temperature and water vapour profiles next to a map showing the location of the observation. The example uses maskingvalue = NaN for selected variables to avoid missing values. We read the "first guess" variables (fg_atmospheric_water_vapour, fg_atmospheric_temperature), because they contain the data from the statistical all-sky retrieval. The 1DVar algorithm that generates the "non-first guess" profiles will be phased out in favor of the statistical retrieval in the future.

using MetopDatasets
using CairoMakie, GeoMakie

ds = MetopDataset("IASI_SND_02_M03_20250120105357Z_20250120123253Z_N_O_20250120123416Z.nat");

# Select color to plot temperature and humidity profile
selected_color = :red

# Select a single data record
data_record_index = 105

# Read temperature and humidity profiles;
# we read the profiles from the statistical retrieval with associated quality indicators
temperature, humidity, selected_point  = let
    temp_var = cfvariable(ds, "fg_atmospheric_temperature", maskingvalue = NaN)
    qi_temp_var = cfvariable(ds, "fg_qi_atmospheric_temperature", maskingvalue = NaN)
    humidity_var = cfvariable(ds, "fg_atmospheric_water_vapour", maskingvalue = NaN)
    qi_humidity_var = cfvariable(ds, "fg_qi_atmospheric_water_vapour", maskingvalue = NaN)
    # keep only soundings where the quality indicators are low
    # thresholds can be relaxed to increase data yield
    # only look at selected data_record_index to save time
    good_retrievals = findall(qi_temp_var[:,data_record_index] .< 2 
                               .&& qi_humidity_var[:,data_record_index] .< 4)
    # we select the 60th retrieval for plotting
    selected_point = good_retrievals[60]
    
    temp_var[:, selected_point, data_record_index], 
    humidity_var[:, selected_point, data_record_index], selected_point
end

# Read the geolocation of the data record
longitude, latitude = let 
    lat_lon_var = cfvariable(ds, "earth_location", maskingvalue = NaN)
    lat_lon = lat_lon_var[:,:,data_record_index]
    lat_lon[2,:], lat_lon[1,:]
end

# Read the pressure levels 
temp_pressure_levels = cfvariable(ds, "pressure_levels_temp", maskingvalue = NaN)[:]
hum_pressure_levels  = cfvariable(ds, "pressure_levels_humidity", maskingvalue = NaN)[:]

# Plot figure
fig = let

    fig = Figure()

    # axis to plot geolocation
    ax = GeoAxis(fig[1, 1:2],
        title = "Observations",
        xlabel = "longitude",
        ylabel = "latitude",
        limits = (extrema(longitude), extrema(latitude)))

    # plot all observations from data record in gray
    scatter!(ax, longitude[:], latitude[:], color=:gray)
    # plot selected observations in color
    scatter!(ax, longitude[[selected_point]], latitude[[selected_point]], 
        color=[selected_color], marker=:xcross, markersize = 15)

    # Add coastlines
    lines!(ax, GeoMakie.coastlines()) 

    # Plot the temperature profile 
    y_limits = (0.001,1900.0)

    ax2 = Axis(fig[2, 1],
        title = "Temperature profile",
        ylabel = "Pressure (hPa)",
        xlabel = "(K)", yreversed = true,
        limits = (nothing, y_limits),
	yscale= log10)

    lines!(ax2,  temperature, temp_pressure_levels/100, color = selected_color)
    
    # Plot the temperature profile 
    ax3 = Axis(fig[2, 2],
        title = "Water vapour  profile",
        ylabel = "Pressure (hPa)",
        xlabel = "(kg/kg)", yreversed = true,
        limits = (nothing, y_limits),
	yscale = log10)

    lines!(ax3,  humidity, hum_pressure_levels/100, color = selected_color)
    hideydecorations!(ax3, grid=false)

    fig
end

The top plot shows a row of observations west of the west African coast. One observation has been marked with a colored X. The corresponding temperature and water vapour profile of the observation are shown below.

Level 2 Carbon monoxide

This example shows how to plot retrieved carbon monoxide for the IASI L2 products.

using MetopDatasets
using CairoMakie, GeoMakie, Statistics

iasi_file = "IASI_SND_02_M03_20250120105357Z_20250120123253Z_N_O_20250120123416Z.nat"
ds = MetopDataset(iasi_file , maskingvalue=NaN);

# select the layer for the carbon monoxide
layer_index = 15
height_in_km = Int.(ds["forli_layer_heights_co"][layer_index])/1000

# get locations
longitude, latitude = let 
    lat_lon = ds["earth_location"][:,:,:]
    lat_lon[2,:,:], lat_lon[1,:,:]
end

# read carbon monoxide data
co_integrated = ds["integrated_co"][:,:].*1000 # convert to g/m2
co_at_height = ds["co_cp_co_a"][layer_index,:,:] .* ds["co_x_co"][layer_index,:,:]

# get color range for plotting
get_quantiles(a) = quantile(a,0.02), quantile(a,0.98)
co_range = get_quantiles(co_at_height[.!isnan.(co_at_height)])
co_integrated_range = get_quantiles(co_integrated[.!isnan.(co_integrated)])

# plot the data
fig = let

    # Create figure and axis
    fig = Figure()
    ax1 = GeoAxis(fig[1, 1],
        title = "Carbon monoxide at $height_in_km km",
        xlabel = "longitude",
        ylabel = "latitude")
    ax2 = GeoAxis(fig[2, 1],
        title = "Carbon monoxide integrated",
        xlabel = "longitude",
        ylabel = "latitude")

    # Add coastlines
    lines!(ax1, GeoMakie.coastlines(), color=:black)
    lines!(ax2, GeoMakie.coastlines(), color=:black) 

    # plot observations
    scatter!(ax1, longitude[:], latitude[:],
        color = co_at_height[:], colorrange = co_range, markersize = 2)
    scatter!(ax2, longitude[:], latitude[:],
        color = co_integrated[:], colorrange = co_integrated_range, markersize = 2)

    # Add colorbar
    c1 = Colorbar(fig[1,2],colorrange = co_range .*10^-16, label="molecules/cm2 x 10^16")
    c2 = Colorbar(fig[2,2],colorrange = co_integrated_range, label="g/m2")

    # hide ticks
    hidedecorations!(ax1, grid = false)
    hidedecorations!(ax2, grid = false)

    fig
end