Retrieval

Retrievals are done using a series of interchangable “blocks” which carry out the many functions required for a DOAS-style retrieval. The hierarchy of classes is shown below.

The following example shows how to simulate radiances using SASKTRAN, package the radiances in a Radiance object that skdoas.retrieval can use, and then perform a total column NO2 retrieval.

[1]:
import numpy as np
import xarray as xr
from scipy.integrate import trapz

import sasktran as sk
import skdoas.retrieval as skr
[2]:
### SIMULATE RADIANCES ###
mjd = 54000.6
lats, lons = [52.1], [-106.7]  # Saskatoon

# create lines of sight
observer = sk.Geodetic()
observer.from_lat_lon_alt(0., -100., 35786e3)  # approximate TEMPO location
geometry = sk.NadirGeometry()
geometry.from_lat_lon(mjd=mjd, observer=observer, lats=lats, lons=lons)
geometry.reference_point = lats + lons + [0.0, mjd]

# create atmosphere
atmosphere = sk.Atmosphere()
atmosphere.atmospheric_state = sk.MSIS90()
atmosphere['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())
atmosphere['ozone'] = sk.Species(sk.O3DBM(), sk.Labow())
atmosphere['no2'] = sk.Species(sk.NO2Vandaele1998(), sk.Pratmo())

wl = np.arange(405., 465., 0.2)
engine = sk.EngineHR(geometry, atmosphere, wl)
rad = engine.calculate_radiance('xarray')

measurement_l1 = skr.Radiance_SASKTRAN(sk_result_xarray=rad)
[3]:
### PERFORM RETRIEVAL ON SIMULATED RADIANCES ###

print('setting up retrieval')
# create SCD fitting object
instrument = skr.DeltaFunction()
cross_sections = skr.CrossSection_Default(instrument=instrument, wavelength=np.arange(405., 465.1, 0.1))
cross_sections.add_cross_section(sk.NO2Vandaele1998(), 220., 1e5)
cross_sections.add_cross_section(sk.O3DBM(), 220., 1e5)
cross_sections.add_cross_section(sk.HITRANChemical('H2O'), 220., 1e5)
cross_sections.release_weights()
solar_spectrum = skr.SolarSpectrum_Default(instrument=instrument, wavelength=np.arange(405., 465.1, 0.1))
ring_spectrum = skr.RingSpectrum_Default(solar_spectrum=solar_spectrum, wavelength=np.arange(405., 465.1, 0.1))
scd_object = skr.SlantColumn_OMNO2(cross_sections, ring_spectrum)

# create AMF object
amf_atmosphere = skr.Atmosphere_Default(atmosphere=atmosphere)
amf_table = skr.BoxAirMassFactor_440nm50kmClearTable()
amf_object = skr.AirMassFactor_OMNO2(amf_table, amf_atmosphere)

retrieval = skr.Retrieval_OMNO2_TotalColumnOnly(scd_object, amf_object)

print('performing retrieval')
vcdr, vcdr_error = retrieval.retrieve(measurement_l1)

alts = amf_table.altitude_edges()
print('\nretrieval results:')
for i in range(measurement_l1.num_lines_of_sight()):
    profile = sk.Pratmo().get_parameter('SKCLIMATOLOGY_NO2_CM3', lats[i], lons[i], alts, mjd)
    vcd = trapz(profile, 1e2 * alts)
    print(f'{vcdr[i]:10.4e} | {vcd:10.4e} | {1e2*(vcdr[i]/vcd - 1):6.2f}%')

setting up retrieval
performing retrieval

retrieval results:
2.2799e+15 | 2.2922e+15 |  -0.54%