Go to notebook file

[1]:

%reload_ext eradiate.notebook.tutorials

Last updated: 2022-10-05 13:34 (eradiate v0.22.4.post102+g84ef86b.d20221005)

Particle layer basics#

Overview

In this tutorial, we introduce Eradiate’s particle layer simulation features. We also introduce the general heterogeneous atmosphere container.

Prerequisites

What you will learn

  • How to add one or several aerosol layers to a molecular atmosphere.

The particle layer interface#

We start by activating the IPython extension and importing and aliasing a few useful components. We also select the CKD mode.

[2]:

%load_ext eradiate
import eradiate
import numpy as np
import matplotlib.pyplot as plt
from eradiate.notebook.tutorials import plot_sigma_t

eradiate.set_mode("ckd")

We alias the eradiate.scenes module and the unit registry for convenience:

[3]:

import eradiate.scenes as ertsc
from eradiate import unit_registry as ureg

The ParticleLayer class is used to configure particle layers. It is parametrised by:

  • an extent, defined by its bottom and top parameters;

  • a number of layers (n_layers) the particle layer is discretised;

  • a particle spatial density distribution (e.g. uniform or gaussian)

  • an optical thickness value (tau_ref) at a reference wavelength set by its w_ref parameter (by default 550 nm).

We start by creating a particle layer with default parameters:

[4]:

particle_layer_default = ertsc.atmosphere.ParticleLayer()
particle_layer_default

ParticleLayer(
    id='atmosphere',
    geometry=None,
    scale=None,
    _bottom=0.0 km,
    _top=1.0 km,
    distribution=UniformParticleDistribution(bounds=array([0., 1.])),
    w_ref=550 nm,
    tau_ref=0.2,
    n_layers=16,
    dataset=<xarray.Dataset>
Dimensions:  (w: 20, mu: 181, i: 1, j: 1)
Coordinates:
  * w        (w) float64 350.0 400.0 412.0 443.0 ... 1.95e+03 2.25e+03 3.75e+03
  * i        (i) int64 0
  * j        (j) int64 0
  * mu       (mu) float64 1.0 1.0 1.0 1.0 1.0 ... -0.9983 -0.9993 -0.9998 -1.0
Data variables:
    sigma_t  (w) float64 0.05645 0.04158 0.03882 ... 0.00578 0.005448 0.005256
    albedo   (w) float64 0.9808 0.9798 0.9795 0.9784 ... 0.9129 0.9011 0.8525
    phase    (w, mu, i, j) float64 178.4 173.6 161.7 ... 0.2595 0.2662 0.2686
Attributes:
    convention:  CF-1.8
    title:       Particles radiative properties
    history:     2021-11-03 15:29:43 - data set creation - kata, version 0.0.2,
    _phase=TabulatedPhaseFunction(
        id='phase_atmosphere',
        data=<xarray.DataArray 'phase' (w: 20, mu: 181, i: 1, j: 1)>
array([[[[2.1040000e-01]],

        [[2.1451000e-01]],

        [[2.2864000e-01]],

        ...,

        [[1.6167650e+02]],

        [[1.7361663e+02]],

        [[1.7842933e+02]]],


       [[[2.1991000e-01]],

        [[2.1822000e-01]],

        [[2.2853000e-01]],
...
        [[6.7030030e+01]],

        [[6.7152350e+01]],

        [[6.7193330e+01]]],


       [[[2.6860000e-01]],

        [[2.6621000e-01]],

        [[2.5953000e-01]],

        ...,

        [[3.0493370e+01]],

        [[3.0512870e+01]],

        [[3.0519370e+01]]]])
Coordinates:
  * w        (w) float64 350.0 400.0 412.0 443.0 ... 1.95e+03 2.25e+03 3.75e+03
  * i        (i) int64 0
  * j        (j) int64 0
  * mu       (mu) float64 -1.0 -0.9998 -0.9993 -0.9983 ... 1.0 1.0 1.0 1.0
Attributes:
    standard_name:  scattering_phase_matrix
    long_name:      scattering phase matrix
    units:          str^-1
    )
)

We see that our particle_layer_default variable represents an aerosol layer with:

  • a vertical extent from 0 to 1 km;

  • a uniform particle density distribution within that extent;

  • a division into 16 uniform cells;

  • an optical thickness at 550 nm equal to 0.2.

We can now visualise the extinction coefficient using the following convenience function:

[5]:

plot_sigma_t(particle_layer_default)
../../_images/tutorials_getting_started_particle_layer_9_0.png

Now, let’s define a particle layer with a different optical thickness:

[6]:

particle_layer_thick = ertsc.atmosphere.ParticleLayer(tau_ref=1.0)
plot_sigma_t(
    particle_layer_default,
    particle_layer_thick,
    labels=["default", "thick"],
)
../../_images/tutorials_getting_started_particle_layer_11_0.png

We can also modify the extent of our particle using the bottom and top parameters:

[7]:

particle_layer_high = ertsc.atmosphere.ParticleLayer(tau_ref=0.5, bottom=1 * ureg.km, top=2 * ureg.km)
plot_sigma_t(
    particle_layer_default,
    particle_layer_thick,
    particle_layer_high,
    labels=["default", "thick", "high"],
)
../../_images/tutorials_getting_started_particle_layer_13_0.png

We can vary the particle density distribution within the defined bounds:

[8]:

particle_layer_gaussian = ertsc.atmosphere.ParticleLayer(
    tau_ref=0.5,
    bottom=0.5 * ureg.km,
    top=1.5 * ureg.km,
    distribution="gaussian",
    n_layers=31,
)

plot_sigma_t(
    particle_layer_default,
    particle_layer_thick,
    particle_layer_high,
    particle_layer_gaussian,
    labels=["default", "thick", "high", "gaussian"],
)
../../_images/tutorials_getting_started_particle_layer_15_0.png

Adding a particle layer to a molecular atmosphere#

Although Eradiate can simulate radiative transfer in a single particle layer, realistic use cases rather involve adding a particle layer to a molecular atmosphere. For this, we use the HeterogeneousAtmosphere class. This class acts as a container which can mix a MolecularAtmosphere and an arbitrary number of ParticleLayers.

[9]:

exponential = ertsc.atmosphere.ParticleLayer(
    tau_ref=0.05,
    bottom=0 * ureg.km,
    top=6 * ureg.km,
    distribution="exponential",
)

gaussian = ertsc.atmosphere.ParticleLayer(
    tau_ref=0.02,
    bottom=1 * ureg.km,
    top=2 * ureg.km,
    distribution="gaussian",
)

atmosphere = ertsc.atmosphere.HeterogeneousAtmosphere(
    molecular_atmosphere={"type": "molecular"},
    particle_layers=[exponential, gaussian]
)

We can display the extinction coefficient for the full atmospheric profile and all its components as follows (note that we restrict the view range to the [0, 10] km extent for clarity

[10]:

plot_sigma_t(
    atmosphere,
    atmosphere.molecular_atmosphere,
    exponential,
    gaussian,
    labels=["atmosphere", "molecular", "exponential", "gaussian"],
    altitude_extent=(0, 10),
)
../../_images/tutorials_getting_started_particle_layer_19_0.png

We see on this plot that the various components of the HeterogeneousAtmosphere container are stacked as expected. We also see that the components are all assumed piecewise-constant on their own grid, then regridded on a fine mesh prior to merge.

Running the simulation#

We can now use this atmosphere definition to run a 1D simulation and compute the top-of-atmosphere BRF:

[11]:

# Show only spectral loop progress
eradiate.config.progress = "spectral_loop"

exp = eradiate.experiments.AtmosphereExperiment(
    surface={"type": "lambertian", "reflectance": 1.0},
    atmosphere=atmosphere,
    illumination={"type": "directional", "zenith": 30.0, "azimuth": 0.0},
    measures={
        "type": "mdistant",
        "construct": "from_viewing_angles",
        "zeniths": np.arange(-75, 76, 5),
        "azimuths": 0.0,
        "spectral_cfg": {"bins": ["550"]},  # Run the simulation for the 550 nm bin
        "spp": 10000,
    },
)
result = eradiate.run(exp)
result

[11]:

<xarray.Dataset>
Dimensions:         (sza: 1, saa: 1, y_index: 1, x_index: 31, w: 1, srf_w: 1)
Coordinates: (12/13)
  * sza             (sza) float64 30.0
  * saa             (saa) float64 0.0
  * y_index         (y_index) int64 0
  * x_index         (x_index) int64 0 1 2 3 4 5 6 7 ... 23 24 25 26 27 28 29 30
    bin             (w) object '550'
    y               (y_index) float64 0.0
    ...              ...
  * w               (w) float64 550.0
    bin_wmin        (w) float64 545.0
    bin_wmax        (w) float64 555.0
    vza             (x_index, y_index) float64 -75.0 -70.0 -65.0 ... 70.0 75.0
    vaa             (x_index, y_index) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
  * srf_w           (srf_w) float64 550.0
Data variables:
    radiance        (sza, saa, w, y_index, x_index) float64 0.4107 ... 0.415
    spp             (sza, saa, w) float64 1e+04
    irradiance      (sza, saa, w) float64 1.621
    srf             (srf_w) float64 1.0
    brdf            (sza, saa, w, y_index, x_index) float64 0.2534 ... 0.2561
    brf             (sza, saa, w, y_index, x_index) float64 0.7961 ... 0.8044
    radiance_srf    (sza, saa, y_index, x_index) float64 0.4107 0.433 ... 0.415
    irradiance_srf  (sza, saa) float64 1.621
    brdf_srf        (sza, saa, y_index, x_index) float64 0.2534 ... 0.2561
    brf_srf         (sza, saa, y_index, x_index) float64 0.7961 ... 0.8044
Attributes:
    convention:  CF-1.8
    source:      eradiate, version 0.22.4.post102+g84ef86b.d20221005
    history:     2022-10-05 13:34:35 - data creation - AtmosphereExperiment.p...
    references:
    title:       Top-of-atmosphere simulation results

And we can plot the computed TOA BRF:

[12]:

with plt.rc_context({"lines.marker": ".", "lines.linestyle": ":"}):
    result.brf.squeeze(drop=True).plot(x="vza");
../../_images/tutorials_getting_started_particle_layer_25_0.png

Final words#

This quick introduction provides just an overview of how particle layers are parametrised. In particular, particle radiative property datasets (the dataset parameter of the ParticleLayer constructor) are not (yet) covered by this tutorial. Until this is done, you can have a look in the resources/data/spectra/particles and explore those datasets.

Further reading#