Weighting Functions

The weighting function classes are used in the same way as the sasktran.Engine classes: specify a sasktran.Geometry object, a sasktran.Atmosphere object, wavelengths, and engine options (passed as keyword arguments rather than as a dictionary), and they will produce a numpy array or an xarray dataset.

When calculating weighting functions using sasktran.Engine, the atmosphere.wf_species property must be set. Here, this property should be left blank (it will be ignored) and the species should be specified through the wf_species property of the weighting function object.

When using built-in weighting function methods these classes provide little advantage over simply using a sasktran.Engine object. The main advantage of these classes is that they provide a consistent interface for both built-in methods and for finite difference weighting functions.

import skdoas
import sasktran as sk
import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline
plt.rcParams['figure.figsize'] = [15, 9]
plt.rcParams['font.size'] = 20

slat, slon, mjd = 52.1, -106.7, 54000.6  # Saskatoon
mlat, mlon, mjd = 25.8, -80.2, 54000.6  # Miami

# create line 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=[slat, mlat], lons=[slon, mlon])

# 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 = [440.]

heights = np.arange(25e2, 6e4, 5e3)
widths = 25e2 * np.ones(12)

hr = skdoas.WeightingFunctionFiniteDifferenceHR(
    geometry=geometry, atmosphere=atmosphere, wf_species='no2', wfheights=heights, wfwidths=widths)
do = skdoas.WeightingFunctionDO(
    geometry=geometry, atmosphere=atmosphere, wf_species='no2', wfaltitudes=heights, wfwidths=widths)

wfhr = hr(wl).unstack('perturbation').wf_no2
wfdo = do(wl).unstack('perturbation').wf_no2

plt.plot(wfhr.isel(wavelength=0, los=0), wfhr.altitude, label='HR, Saskatoon')
plt.plot(wfdo.isel(wavelength=0, los=0), wfdo.altitude, label='DO, Saskatoon')
plt.plot(wfhr.isel(wavelength=0, los=1), wfhr.altitude, label='HR, Miami')
plt.plot(wfdo.isel(wavelength=0, los=1), wfdo.altitude, label='DO, Miami')
plt.title('Weighting functions from TEMPO line of sight')
plt.xlabel('Weighting function dI/dn (photons.s$^{-1}$.nm$^{-1}$.cm$^{-2}$.sr$^{-1}$/cm$^{-3}$)')
plt.ylabel('Altitude (km)')
plt.legend()
plt.show()