Extinction Retrieval

A simple extinction retrieval using a single wavelength channel and two polarization states.

from skretrieval.retrieval.rodgers import Rodgers
from skretrieval.retrieval.statevector import StateVector
from ali.retrieval.aerosol import AerosolRetrieval
from ali.retrieval.statevector import StateVectorAerosolProfile
from ali.retrieval.measvec import (MeasurementVector, MeasurementVectorElement)
from ali.retrieval.measvec.transformer import AltitudeNormalization, LogRadiance, LinearCombination
from ali.instrument.simulator import ImagingSimulator
from ali.instrument.sensor import ALISensor
from ali.util.geometry import optical_axis_from_geometry
import numpy as np
import sasktran as sk
import matplotlib.pyplot as plt


def basic_extinction_retrieval(image_center_alts=(22.5,), mjd=54372):

    altitudes = np.arange(500.0, 45001.0, 200)
    sensor_wavel = np.array([750.0])
    integration_time = 0.1

    # ---------------------------------------------------------------------------------------- #
    # Setup measurement geometry
    # ---------------------------------------------------------------------------------------- #
    geometry = sk.VerticalImage()
    geometry.from_sza_ssa(sza=60, ssa=90, lat=45.0, lon=0,
                          tanalts_km=image_center_alts, mjd=mjd, locallook=0, satalt_km=500)
    optical_geometry = optical_axis_from_geometry(geometry)

    # ---------------------------------------------------------------------------------------- #
    # Create simulated level 1 data
    # ---------------------------------------------------------------------------------------- #
    atmo_sim = sk.Atmosphere()
    atmo_sim['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())
    atmo_sim['ozone'] = sk.Species(sk.O3DBM(), sk.Labow())
    atmo_sim['aerosol'] = sk.SpeciesAerosolGloSSAC()

    sensor = ALISensor(wavelength_nm=sensor_wavel,
                       image_horiz_fov=5, image_vert_fov=1.1,
                       num_columns=1, num_rows=512)
    sensor._simulate_pixel_averaging = 50  # simulate the binning of 50 detector columns
    sensor.ccd = 'teledyneccd97'

    simulator = ImagingSimulator(sensors=[sensor], optical_axis=[optical_geometry],
                                 atmosphere=atmo_sim, options={'numordersofscatter': 1,
                                                               'polarization': True})
    simulator.sun = geometry.sun
    simulator.grid_sampling = True  # model every line of sight
    simulator.dual_polarization = True
    simulator.integration_time = integration_time

    measurement_l1 = simulator.calculate_radiance()

    # ---------------------------------------------------------------------------------------- #
    # Setup forward model
    # ---------------------------------------------------------------------------------------- #
    atmo_ret = sk.Atmosphere()
    atmo_ret['air'] = sk.Species(sk.Rayleigh(), sk.MSIS90())
    atmo_ret['ozone'] = sk.Species(sk.O3DBM(), sk.Labow())

    ret_sensor = ALISensor(wavelength_nm=sensor_wavel,
                           image_horiz_fov=5 / 512, image_vert_fov=1.1,
                           num_columns=1, num_rows=512)
    ret_sensor._simulate_pixel_averaging = 50
    ret_sensor.ccd = 'teledyneccd97'
    ret_sensor.add_adc = False
    ret_sensor.add_noise = False
    ret_sensor.add_dark_current = False

    # radiative transfer options
    options = {'numordersofscatter': 1, 'polarization': True,
               'calcwf': 2, 'wfheights': altitudes,
               'wfwidths': np.concatenate([np.diff(altitudes), [altitudes[-1] - altitudes[-2]]])}
    forward_model = ImagingSimulator(sensors=[ret_sensor], optical_axis=[optical_geometry],
                                     atmosphere=atmo_ret, options=options)
    forward_model.grid_sampling = True
    forward_model.store_radiance = False  # recompute radiance every call
    forward_model.sun = geometry.sun
    forward_model.dual_polarization = True
    forward_model.integration_time = integration_time

    # ---------------------------------------------------------------------------------------- #
    # Setup retrieval
    # ---------------------------------------------------------------------------------------- #
    # Setup measurement vector
    aerosol_mv = MeasurementVectorElement()
    aerosol_mv.add_transform(LinearCombination({0: 1, 1: 1}))
    aerosol_mv.add_transform(AltitudeNormalization(norm_alts=(35000.0, 40000.0)))
    aerosol_mv.add_transform(LogRadiance())
    meas_vec = MeasurementVector([aerosol_mv])

    # setup state vector
    ps_clim = {'SKCLIMATOLOGY_LOGNORMAL_MODERADIUS_MICRONS': np.ones_like(altitudes) * 0.08,
               'SKCLIMATOLOGY_LOGNORMAL_MODEWIDTH': np.ones_like(altitudes) * 1.6}
    ps_clim = sk.ClimatologyUserDefined(altitudes, ps_clim)
    aerosol_opt_prop = sk.MieAerosol(particlesize_climatology=ps_clim, species='H2SO4')
    apriori_profile = np.loadtxt('aerosol_apriori.txt')
    apriori_profile = np.interp(altitudes, apriori_profile[:, 0], apriori_profile[:, 1])
    aerosol_state_element = StateVectorAerosolProfile(altitudes_m=altitudes,
                                                      values=np.log(apriori_profile),
                                                      species_name='aerosol',
                                                      optical_property=aerosol_opt_prop,
                                                      lowerbound=5500,
                                                      upperbound=40000,
                                                      second_order_tikhonov_factor=np.array([20, 10, 30]),
                                                      second_order_tikhonov_altitude=np.array([5000.0, 25000.0, 40000.0]))
    aerosol_ret = AerosolRetrieval(StateVector([aerosol_state_element]), meas_vec,
                                   cloud_vector=None, retrieval_altitudes=altitudes)
    aerosol_state_element.add_to_atmosphere(atmo_ret)  # add aerosol and weighting functions to atmosphere

    # Do the retrieval
    rodgers = Rodgers(max_iter=5, lm_damping=0.01)
    aerosol_ret.configure_from_model(forward_model, measurement_l1)
    output = rodgers.retrieve(measurement_l1, forward_model, aerosol_ret)

    # ---------------------------------------------------------------------------------------- #
    # Plot retrieval results
    # ---------------------------------------------------------------------------------------- #
    tanalts = measurement_l1[0].tangent_locations().altitude / 1000
    true_state = atmo_sim['aerosol'].climatology.get_parameter('SKCLIMATOLOGY_AEROSOL_EXTINCTIONPERKM', 45.0, 0.0,
                                                               altitudes, mjd)
    xsec = atmo_sim['aerosol'].optical_property.calculate_cross_sections(sk.MSIS90(), 45.0, 0.0,
                                                                         1000.0, 54372, 525.0).total * 1e5
    fig, ax = plt.subplots(1, 2, sharey=True, dpi=200)
    ax[0].plot(np.exp(output['xs'][-1]), altitudes / 1000)
    ax[0].plot(true_state / xsec, altitudes / 1000)
    ax[1].plot(output['y_meas'], tanalts)
    ax[1].plot(output['ys'][-1], tanalts)
    ax[0].set_ylabel('Altitude [km]')
    ax[0].set_xlabel('Number Density [cm$^-3$]')
    ax[1].set_xlabel('Measurement Vector')


if __name__ == '__main__':
    basic_extinction_retrieval()
../../_images/example_ext_ret.png