Aerosol model analysis with ALG

Suspended particles (aerosols) are one of the main contributors of light scattering in the atmosphere radiative transfer. Understanding the radiometric effect of aerosol optical properties in top-of-atmosphere (TOA) radiance is therefore necessary to properly compensate for it and improve the quality of atmospheric correction.

In this tutorial, partly based on (Vicent et al, 2017), we will learn how to generate atmospheric LUTs with different options for the nodes distribution. These LUTs will be based on the coupling between the OPAC aerosol database and MODTRAN. In addition, we will also show how to extract and analyse the OPAC-related data produced by ALG.

Note: On this page, we will use the MATLAB® notation (and commands) for the operations with the LUT data (e.g., reading files, plotting). However, similar commands can be found in other programing languages (e.g., Python, C/C++, Fortran).

Creating the atmospheric LUTs

In this tutorial, we will use two different LUTs: a reference LUT, based on the coupling of OPAC and MODTRAN, and a retrieval LUT, based on parametrization of aerosol optical properties. Both of them are simulated at a solar zenith angle (SZA) of 45 deg and an aerosol optical thickness (AOT) of 0.15. The simulated data will cover the spectral range between 400 nm and 1000 nm at a spectral sampling of 15 cm-1 (approx. 0.7 nm at 700 nm).

Let's start by creating the reference LUT pressing the New LUT configuration button . On the General Configuration step, introduce the file name on the Saved LUT file(s) text box (e.g., C:\ALG\LUTs_folder\reference). Secondly, set the user level as advanced to have access to additional LUT configuration parameters. We can then select MODTRAN5 in the RTM pop-up menu, and the Transfer Functions as operation mode. Finally, we choose the Manual (gridded) sampling as our LUT sampling strategy. We can finally click on the Next button at the bottom-right of the screen to proceed introducing the key atmospheric variables.

We continue by defining the LUT key input variables and their values. We start by pressing on the Geometric group and selecting the Solar zenith angle parameter. On the distribution sub-panel, we introduce the 45 both for the minimum and maximum values and 2 for the number of samples. This will set a single static value of 45 deg for the SZA. Click on the Add parameter button (or press ) to add these value on the table:

We now press on the Aerosols group and select the Aerosol model parameter. We select the Continental (average) Urban (OPAC), Desert (OPAC), Maritime (clean) and User defined options on the list by Crtl.+click them.

After pressing , the user-defined OPAC aerosols configuration window will pop-up. Here, we will randomly define two new aerosol mixtures between continental and desert and other two with a mixture between continental and maritime (see Table 4 in Hess et al, 1998 for further information).

  • For the continental-desert mixture, select the following components and Nmin-Nmax values: Insoluble (0.3-0.6), Water-soluble (2000-7000), and the three minerals (nuclei, acc. mode and coa. mode) with the default Nmin-Nmax values. Then introduce the number 2 on the Number of aerosol mixtures (N) text box and click on the Mix particles (N aerosols) button.

  • For the continental-maritime mixture, repeat the previous steps but now selecting the following components and Nmin-Nmax values: Insoluble (0.15-0.6), Water-soluble (500-15000), soot (2500-8000) and the two sea-salt (acc. mode and coa. mode) with their default Nmin-Nmax values. Add 2 new aerosol mixtures on the table.

Notice that particle density values of the new four aerosol mixtures might differ since they are randomly generated betwween the introduced Nmin-Nmax values. However, you can manually edit these values on the table.

Once defined the configuration of these user-defined aerosol mixtures, click on the Accept. This will add these values on the key input variables table, under the aerosol model variable IHAZE. Notice that these user-defined OPAC aerosols take the ID numbers ≤ -11.

On the Sensor Spectral Configuration step, we will create a single continuous spectral band from 400 nm to 1000 nm at 15 cm-1 spectral resolution. On the Advanced configuration we will not modify the default advanced parameter and we will simply click on the Next button.

ALG will automatically create the input and output OPAC and MODTRAN files (.tp5 and .tp7 respectively for MODTRAN). In the case of MODTRAN, these input/output files will have as name combNrhop5 and combNrhop15, where N is an index between identifying the LUT node combination. A progress bar will show the writing of the MODTRAN input files. After its completion, OPAC will run with the configuration of the selected default and user-defined aerosols. This will create an OPAC LUT file named reference_OPAC.h5 that contains all the relevant aerosol optical properties for the selected OPAC aerosols. MODTRAN will then run the simulations for these RTM input files, including the aerosol optical properties extracted from this OPAC-related file.

With respect to the configuration of the retrieval LUT, on the General Configuration step, introduce the file name on the Saved LUT file(s) text box (e.g., C:\ALG\LUTs_folder\retrieval). As with the reference LUT, set the user level as advanced, MODTRAN5 as the selected RTM, and Transfer Functions as the operation mode. Finally, we choose the Latin hypercube as our LUT sampling strategy. We can finally click on the Next button at the bottom-right of the screen to proceed introducing the key atmospheric variables.

We continue by defining the LUT key input variables and their values. We start by introducing a SZA of 45 deg as we did for the reference LUT. We then press on the Aerosols group and start by selecting the Henyey-Greestein asymmetry factor, setting the values of 0.6 and 0.99 respectively on the min/max text boxes and clicking on the Add parameter button to introduce these values on the variables table. We repeat this process with the following aerosol parameters:

  • Angstrom coefficient. Min: 0.5; Max: 2.
  • Single Scattering Albedo. Min: 0.85; Max: 0.99.

Since the distribution is Latin hypercube we need to define the total # of grid points to be included in the LUT. We therefore introduce a value of 500 points in the corresponding text box.

For the spectral configuration and advanced parameters, we follow the same steps as in the reference LUT.

Once both (reference and retrieval) RTM simulations have finished, we can generate the atmospheric LUTs from the output RTM files. On the menu bar, press the New LUT generation button and select the LUT configuration reference.xml file. ALG will automatically load the LUT configuration, generating the LUT key input parameters into a header file (reference.hdr) and processing/storing the RTM output files into the LUT output file (reference.h5). Repeat the process for the retrieval.xml file generating the retrieval.hdr file*.

*Note: From version v.1.2, ALG stores the LUT header information in the same .h5 file as the LUT data. All steps related with reading the LUT header are valid replacing the extension .hdr by .h5.

Analysis of aerosol optical properties

In this second section, we will show how to read the OPAC output data stored in the reference_OPAC.h5 file. This file contains information about the phase function, the extinction coefficient (in km-1), the absorption coeffiecient (in km-1) and the asymmetry parameter. These optical properties are given in a wavelength VARSPC (in μm) and relative humidity RH (in %) grid. The phase function is also given at an scattering angle ANGFgrid (in deg). We can load this information by reading the corresponding datasets from the .h5 file:

VARSPC  = h5read(filename,'/VARSPC')*1000;
ANGF    = h5read(filename,'/ANGF');
RH      = h5read(filename,'/RH');
SPEC    = h5read(filename,'/SPEC'); 
auxF    = h5read(filename,'/F');

The aerosol phase function (F) is a 4-dimensional variable (ANGF×VARSPC×RH×aerosol number). We can normalize its values by applying the following equation (1):

$$F' = {F}/{2π∫_0^π{F·sin(θ)·dθ}}$$

The SPEC dataset contains the asymmetry parameter (ASYM), the extinction coefficient (EXTC) and the absorption coefficient (ABSC). Therefore, it is a 4-dimensional variable (VARSPC×RH×optical property×aerosol number). We can extract each individual optical property with the following command lines:

auxEXTC = squeeze(SPEC(:,:,1,:));
auxABSC = squeeze(SPEC(:,:,2,:));
auxASYM = squeeze(SPEC(:,:,3,:));

We can interpolate these optical properties to the RH of our atmosphere (mid-latitude summer by default; RH=74.4%). In the following figures, we can plot these optical properties as function of wavelength or scattering angle.

Note: In the following figures, the color code identifies the aerosol types: Continental + Maritime (#2), Continental + Maritime (#1), Continental + Desert (#2), Continental + Desert (#1), Maritime (clean), Desert, Urban, Continental (average).

In the case of the extinction coefficient values (see circles in figure below), we can fit (dashed line) an Ångström law \({ε(λ)}/{ε(550 nm)}={({550 nm}/{λ})}^α\). We can observe a good fitting of this parametric approximation except for the Continental (average) curve.

In the case of the phase function values (see solid line in figure below), we can fit (dashed line) a parametric Henyey-Greenstein function \({1}/{4π}{1-g^2}/({1+g^2-2g·cosθ})^{3/2}\). We can observe an important discrepancy between the real phase function and the parametric approximation, accentuated by the use of a log-log scale in the plot:

With respect to the asymmetry parameter, we notice that aerosols affected by minerals (desert) or by soot (Urban and Continental (average)) deviates from a nearly spectrally constant value:

With respect to the single scattering albedo (i.e., (EXTC-ABSC)/EXTC), we notice again that aerosols affected by minerals (desert) or by soot (Urban and Continental (average)) deviates from a nearly spectrally constant value.

Error analysis in atmospheric correction

Let us start by reading the LUT header and output files for the reference and retrieval LUTs:

hdr1 = readLUThdr([fileref,'.hdr']);
hdr2 = readLUThdr([fileret,'.hdr']);
Y1   = h5read([fileref,'.h5'],'/LUTdata');
Y2   = h5read([fileret,'.h5'],'/LUTdata');

In the files above, fileref and fileret represent the path and file name of our reference and retrieval LUT files (e.g, C:\ALG\LUTs_folder\reference).

As it can be read in the LUT header file, both variables Y1 and Y2 contain the following atmospheric transfer functions:

  • Atmospheric path radiance (L0) in [mW·m-2·sr-1·nm-1].
  • At-surface direct and diffuse solar irradiance (Edir and Edif) in [mW·m-2·nm-1].
  • Target-to-sensor direct and diffuse transmittance (Tdir and Tdif) (unit-less).
  • Spherical albedo (S) (unit-less).

We can therefore extract these atmospheric transfer functions taking into account that each of these functions have the same number of wavelengths:

L0   = Y1(1:numel(hdr1.wvl),:);
Edir = Y1(numel(hdr1.wvl)+1:2*numel(hdr1.wvl),:);
Edif = Y1(2*numel(hdr1.wvl)+1:3*numel(hdr1.wvl),:);
S    = Y1(3*numel(hdr1.wvl)+1:4*numel(hdr1.wvl),:);
Tdir = Y1(4*numel(hdr1.wvl)+1:5*numel(hdr1.wvl),:);
Tdif = Y1(5*numel(hdr1.wvl)+1:6*numel(hdr1.wvl),:);

We can propagate the top-of-canopy Lambertian reflectance ρ to construct TOA radiance (L_{ref}) with the following equation (2):

$$L_{ref} = L_0 + {(E_{dir}\cosθ_{il}+E_{dif})·(T_{dir}+T_{dif})·ρ}/{π(1-Sρ)}=$$ $$= L_0 + {E_{tot}·T_{tot}·ρ}/{π(1-Sρ)}$$

where θil=45 deg is the solar zenith angle.

Let us consider the following surface reflectance, which has been previously resampled to the wavelengths in hdr1.wvl =hdr2.wvl:

For this surface reflectance, let us now construct the TOA radiance for all combinations in our LUT:

LTOA_ref = L0 + (1/pi)*(bsxfun(@times,(Edir*cosd(hdr1.SZA)+Edif).*(Tdir+Tdif),rho))./(1-bsxfun(@times,S,rho));

We can repeat these steps to extract the atmospheric transfer functions and calculate the TOA radiance for the retrieval LUT.

The following figure shows the resulting TOA radiances overlapped for both LUTs:

We can observe that the use of parametric approximations of aerosol optical properties in the retrieval LUT (in red) covers the radiance range obtained by the complex OPAC aerosol properties introduced in the reference LUT. We can therefore find the combinations in the retrieval LUT that best matches each spectral in the reference LUT. In this way, we can analyse the error propagation in the atmospherically corrected surface reflectance.

So for each TOA radiance spectra (i) in the reference LUT, we will find the TOA radiance spectra (j) in the retrieval LUT that minimizes the following equation (3):

$$j_{min} = min_j∑↙{λ}{[L_{ref,i}(λ) - L_{ret,j}(λ)]^2}$$

In the equation above, we will filter out the wavelengths inside the absorption bands of the O2-B (686-708 nm), O2-A (758-772 nm) and the H2O (713-753 nm, 780-850 nm and 891-986 nm).

Lets now perform an atmospheric correction of the reference TOA radiances using the bets matching combinations (jmin) from the atmospheric transfer functions in the retrieval LUT. The correction is achieved by inverting the equation (2), which leads to the following equation (4):

$$ρ = {π(L_{ref} - L_0)}/{E_{tot}·T_{tot} + πS·(L_{ref} - L_0)}$$

In Matlab:

A 	 = pi*(LTOA - L0(:,filter2));
B 	 = Tgas.*Etot.*Ttot; B = B(:,filter2);
rho_corr = A./(B+A.*S(:,filter2));

We can therefore calculate the relative error in reflectance after doing the atmospheric correction:

We observe that the higher relative errors are obtained for the Continental + Desert (#2), Continental + Desert (#1), Maritime (clean), Desert aerosols, i.e., .those with higher spectral variability of the asymmetry parameter and single scattering albedo. Also, we can observe the higher errors in absorption bands for the Urban and Continental (average) aerosols, particularly in the H2O band at 950 nm. This is due to the higher absorption of these aerosols.

References