2-canopy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Canopy simulation tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains sample cases for the Eradiate preview session of Dec 7th 2021. In this tutorial, we introduce the canopy simulation capabilities of Eradiate. We will explore the features and interfaces of different objects used to construct a scene, simulate radiative transfer and visualise results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We load the Rich notebook extension for improved object inspection\n",
    "%load_ext rich\n",
    "\n",
    "# We set the Matplotlib inline backend to SVG for better figure quality\n",
    "from matplotlib_inline import backend_inline\n",
    "backend_inline.set_matplotlib_formats(\"svg\")\n",
    "\n",
    "# We import general processing and plotting libraries\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import pandas as pd\n",
    "\n",
    "# We import the top-level Eradiate module\n",
    "import eradiate\n",
    "\n",
    "# For convenience, we alias a few components\n",
    "from eradiate import unit_registry as ureg\n",
    "import eradiate.scenes as ertsc\n",
    "import eradiate.experiments as ertxp\n",
    "\n",
    "# We'll use the monochromatic mode in double precision\n",
    "eradiate.set_mode(\"mono_double\")\n",
    "\n",
    "# We activate a few convenient presets\n",
    "eradiate.notebook.install()\n",
    "eradiate.plot.set_style(rc={\"figure.dpi\": 96})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with the definition of a couple of utility functions which will simplify the visualisation of results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_camera(\n",
    "    exp,\n",
    "    measure_id,\n",
    "    robust=True,\n",
    "    add_colorbar=False,\n",
    "    vmin=None,\n",
    "    vmax=None\n",
    "):\n",
    "    \"\"\"\n",
    "    Display the output of a monochromatic camera measure.\n",
    "    \"\"\"\n",
    "    _, ax = plt.subplots(1,1)\n",
    "    exp.results[measure_id][\"radiance\"].squeeze(drop=True).plot.imshow(\n",
    "        ax=ax,\n",
    "        origin=\"upper\",\n",
    "        cmap=\"Greys_r\",\n",
    "        vmin=vmin,\n",
    "        vmax=vmax,\n",
    "        robust=robust,\n",
    "        add_colorbar=add_colorbar,\n",
    "    )\n",
    "    ax.set_aspect(1)\n",
    "    plt.show()\n",
    "    plt.close()\n",
    "\n",
    "def show_brf(exp, measure_id):\n",
    "    \"\"\"\n",
    "    Display the BRF output of a distant radiance measure.\n",
    "    \"\"\"\n",
    "    _, ax = plt.subplots(1,1)\n",
    "    with plt.style.context({\"lines.linestyle\": \":\", \"lines.marker\": \".\"}):\n",
    "        exp.results[measure_id][\"brf\"].squeeze(drop=True).plot(ax=ax, x=\"vza\")\n",
    "    plt.show()\n",
    "    plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rendering a surface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "First, let's create a very simple scene consisting of a default Lambertian surface observed by a camera. We start by defining the surface and set its width to 10 m (the default is much larger)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lambertian_surface=ertsc.surface.LambertianSurface(\n",
    "    width=10 * ureg.m, reflectance=0.5\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's make a monochromatic render of this scene. We'll produce a 4:3 ratio image with a fairly low resolution sample count—all we want is to visualise the scene.\n",
    "\n",
    "**Note**: Using the correlated *k* band mode to render scenes is more computationally intensive. Therefore, we will stick to Eradiate's monochromatic mode for renders in this tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "camera_oblique = ertsc.measure.PerspectiveCameraMeasure(\n",
    "    id=\"camera_oblique\",\n",
    "    origin=[15, 15, 15] * ureg.m,\n",
    "    target=[0, 0, 0] * ureg.m,\n",
    "    up=[0, 0, 1],\n",
    "    film_resolution=(320, 240),\n",
    "    spp=512,\n",
    ")\n",
    "exp = ertxp.RamiExperiment(\n",
    "    surface=lambertian_surface, measures=camera_oblique\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp.run()\n",
    "show_camera(exp, \"camera_oblique\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might notice rendering artifacts. This is because Eradiate's camera settings are currently not optimised for rendering. This will be improved in future releases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rendering a surface and canopy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Let's add a canopy above this ground patch. We'll configure a homogeneous canopy composed of randomly oriented floating disks. This is done using the [`DiscreteCanopy` class](https://eradiate.readthedocs.io/en/latest/rst/reference/generated/autosummary/eradiate.scenes.biosphere.DiscreteCanopy.html), or, rather, its [`homogeneous()`](https://eradiate.readthedocs.io/en/latest/rst/reference/generated/autosummary/eradiate.scenes.biosphere.DiscreteCanopy.html) class constructor. We will use abstract parameters for the sake of demonstrating this feature:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "homogeneous_canopy = ertsc.biosphere.DiscreteCanopy.homogeneous(\n",
    "    l_vertical=1.0 * ureg.m,\n",
    "    l_horizontal=10.0 * ureg.m,\n",
    "    lai=2.0,\n",
    "    leaf_radius=10 * ureg.cm,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We add this canopy to a scene with a surface using the [`RamiExperiment` class](https://eradiate.readthedocs.io/en/latest/rst/reference/generated/autosummary/eradiate.experiments.RamiExperiment.html). Note that this time, we don't set the surface size: it is automatically adjusted to match the extents we set for the canopy object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp = ertxp.RamiExperiment(\n",
    "    surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "    canopy=homogeneous_canopy,\n",
    "    measures=camera_oblique,\n",
    ")\n",
    "exp.run()\n",
    "show_camera(exp, \"camera_oblique\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rendering a surface and padded canopy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In order to simulate the fact that our canopy is periodic, we can pad it with clones of itself. Let's start with just of few of them.\n",
    "\n",
    "**Tip**: If the render time is too long, you can reduce the number of radiance samples per image pixel by adjusting the measure's `spp` parameter.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "camera_oblique = ertsc.measure.PerspectiveCameraMeasure(\n",
    "    id=\"camera_oblique\",\n",
    "    origin=[15, 15, 15] * ureg.m,\n",
    "    target=[0, 0, 0],\n",
    "    up=[0, 0, 1],\n",
    "    film_resolution=(320, 240),\n",
    "    spp=64,  # Lower sample count value ensures quicker rendering\n",
    ")\n",
    "\n",
    "exp = ertxp.RamiExperiment(\n",
    "    surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "    canopy=homogeneous_canopy,\n",
    "    padding=1,\n",
    "    measures=camera_oblique,\n",
    ")\n",
    "exp.run()\n",
    "show_camera(exp, \"camera_oblique\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's relocate our camera to the top of canopy. We will observe it at a grazing angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "camera_toc = ertsc.measure.PerspectiveCameraMeasure(\n",
    "    id=\"camera_toc\",\n",
    "    origin=[15, 15, 2] * ureg.m,\n",
    "    target=[0, 0, 0],\n",
    "    up=[0, 0, 1],\n",
    "    film_resolution=(320, 240),\n",
    "    spp=64,  # Lower sample count value ensures quicker rendering\n",
    ")\n",
    "exp = ertxp.RamiExperiment(\n",
    "    surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "    canopy=homogeneous_canopy,\n",
    "    padding=2,\n",
    "    measures=[camera_toc],\n",
    ")\n",
    "exp.run()\n",
    "show_camera(exp, \"camera_toc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing the BRF in the principal plane\n",
    "\n",
    "Now that we know how to build our scene, we are going to compute the BRF in the principal plane for a similar scene. This is similar to the test cases for which we can submit results to the RAMI online checker (ROMC.)\n",
    "\n",
    "In order to observe something meaningful, let's set the Sun zenith angle to 30°."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "illumination = ertsc.illumination.DirectionalIllumination(\n",
    "    zenith=30.0, azimuth=0.0\n",
    ")\n",
    "\n",
    "brf_pplane = ertsc.measure.MultiDistantMeasure.from_viewing_angles(\n",
    "    id=\"brf_pplane\",\n",
    "    azimuths=0.0,\n",
    "    zeniths=np.arange(-75, 76, 5),\n",
    "    spp=50000,\n",
    ")\n",
    "\n",
    "exp = ertxp.RamiExperiment(\n",
    "    illumination=illumination,\n",
    "    surface=ertsc.surface.LambertianSurface(),\n",
    "    canopy=homogeneous_canopy,\n",
    "    padding=10,\n",
    "    measures=[brf_pplane],\n",
    ")\n",
    "exp.run()\n",
    "show_brf(exp, \"brf_pplane\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There still is some Monte Carlo noise, but we're keeping the number of samples low in order so that the simulation remains short.\n",
    "\n",
    "**Exercises for the interested reader**\n",
    "\n",
    "1. Change the viewing zenith angle and verify that the hot spot moves accordingly.\n",
    "2. Change the viewing azimuth angle and adjust the measure so as to keep it aligned with the principal plane.\n",
    "3. Compute the BRF in a plane orthgonal to the principal plane.\n",
    "4. Gradually increase padding and determine from which value the quasi-periodic approximation can be considered valid."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing the BRF in the principal plane with an atmosphere"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can combine our canopy with an atmosphere. This uses a different experiment class ([`Rami4ATMExperiment`](https://eradiate.readthedocs.io/en/latest/rst/reference/generated/autosummary/eradiate.experiments.Rami4ATMExperiment.html)), which merges the capabilities of the `RamiExperiment` and `OneDimExperiment` classes. We will only use a purely scattering atmosphere and perform monochromatic simulations so as to limit the runtime of our examples.\n",
    "\n",
    "**Note**: This section demonstrates experimental features. We are aware of several issues, notably intersection inaccuracies. This will be fixed in the coming weeks!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "eradiate.set_mode(\"mono_double\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will first add a default atmosphere and reuse the top-of-canopy view from our earlier example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "camera_toc = ertsc.measure.PerspectiveCameraMeasure(\n",
    "    id=\"camera_toc\",\n",
    "    origin=[15, 15, 2] * ureg.m,\n",
    "    target=[0, 0, 0],\n",
    "    up=[0, 0, 1],\n",
    "    film_resolution=(320, 240),\n",
    "    spp=64,  # Lower sample count value ensures quicker rendering\n",
    ")\n",
    "exp = ertxp.Rami4ATMExperiment(\n",
    "    surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "    atmosphere=ertsc.atmosphere.MolecularAtmosphere(has_absorption=False),\n",
    "    canopy=homogeneous_canopy,\n",
    "    padding=2,\n",
    "    measures=[camera_toc],\n",
    ")\n",
    "exp.run()\n",
    "show_camera(exp, \"camera_toc\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our render is very noisy, but we can clearly see that the addition of the atmosphere to our previous setup added a haze. We can now reposition our camera to the ToA and select a wavelength in the UV, where Rayleigh scattering is strong:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for atmosphere in [\n",
    "    None,  # Atmosphere-free reference\n",
    "    ertsc.atmosphere.MolecularAtmosphere(has_absorption=False)  # Rayleigh scattering only (no gaseous absorption, no aerosols)\n",
    "]:\n",
    "    camera_toa = ertsc.measure.PerspectiveCameraMeasure(\n",
    "        id=\"camera_toa\",\n",
    "        origin=[0, 0, 101] * ureg.km,\n",
    "        target=[0, 0, 0],\n",
    "        up=[0, 1, 0],\n",
    "        film_resolution=(128, 128),\n",
    "        fov=2e-2,  # The camera is far away, reducing the field of view is like zooming\n",
    "        spp=512,  # Lower sample count value ensures quicker rendering\n",
    "        spectral_cfg={\"wavelengths\": [300.0] * ureg.nm},\n",
    "    )\n",
    "    exp = ertxp.Rami4ATMExperiment(\n",
    "        surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "        atmosphere=atmosphere,\n",
    "        canopy=homogeneous_canopy,\n",
    "        padding=2,\n",
    "        measures=[camera_toa],\n",
    "    )\n",
    "    exp.run()\n",
    "    show_camera(exp, \"camera_toa\", add_colorbar=True, vmin=0.07, vmax=0.15)\n",
    "    # Note: Normalise v-range so that we can compare both images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Despite the noise, we can see that Rayleigh scattering dims the signal at this wavelength."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercises for the interested reader**\n",
    "\n",
    "1. With a non-zenith illumination, set up a TOA BRF measurement for a few wavelengths ranging from the UV to NIR regions and see how Rayleigh scattering broadens the hot spot reflective lobe.\n",
    "1. Increase padding value and see how adjacency effects impact simulation results differently depending on the wavelength.\n",
    "2. Switch to the CKD mode and compute the TOA BRF for Sentinel-2 bands 1 and 8A. You will have to use the AFGL 1986 atmosphere model. You will also have to adust simulation parameters to keep the computational time reasonable while setting up your simulation (decrease the sample count while preparing it, then increase it to compute the final results)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution to exercise 1\n",
    "\n",
    "toa_brf = ertsc.measure.MultiDistantMeasure.from_viewing_angles(\n",
    "    id=\"toa_brf\",\n",
    "    zeniths=np.arange(-75, 76, 5),\n",
    "    azimuths=0.0,\n",
    "    spp=10000,\n",
    "    spectral_cfg={\"wavelengths\": [300, 325, 350, 400, 450, 500]}\n",
    ")\n",
    "\n",
    "exp = ertxp.Rami4ATMExperiment(\n",
    "    illumination=ertsc.illumination.DirectionalIllumination(zenith=30.0),\n",
    "    surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "    atmosphere=ertsc.atmosphere.MolecularAtmosphere(has_absorption=False),\n",
    "    canopy=homogeneous_canopy,\n",
    "    padding=2,\n",
    "    measures=[toa_brf],\n",
    ")\n",
    "exp.run()\n",
    "ds = exp.results[\"toa_brf\"]\n",
    "\n",
    "with plt.style.context({\"lines.linestyle\": \":\", \"lines.marker\": \".\"}):\n",
    "    ds.brf.squeeze(drop=True).plot(hue=\"w\", x=\"vza\")\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution to exercise 3\n",
    "\n",
    "eradiate.set_mode(\"ckd_double\")\n",
    "eradiate.config.progress = 1  # Do not show Mitsuba-level progress bar\n",
    "\n",
    "toa_brf = ertsc.measure.MultiDistantMeasure.from_viewing_angles(\n",
    "    id=\"toa_brf\",\n",
    "    zeniths=np.arange(-75, 76, 5),\n",
    "    azimuths=0.0,\n",
    "    spp=1000,\n",
    "    spectral_cfg={\"srf\": \"sentinel_2a-msi-8a\"}\n",
    ")\n",
    "\n",
    "exp = ertxp.Rami4ATMExperiment(\n",
    "    illumination=ertsc.illumination.DirectionalIllumination(zenith=30.0),\n",
    "    surface=ertsc.surface.LambertianSurface(reflectance=0.5),\n",
    "    atmosphere=ertsc.atmosphere.MolecularAtmosphere.afgl_1986(),\n",
    "    canopy=homogeneous_canopy,\n",
    "    padding=2,\n",
    "    measures=[toa_brf],\n",
    ")\n",
    "exp.run()\n",
    "ds = exp.results[\"toa_brf\"]\n",
    "\n",
    "with plt.style.context({\"lines.linestyle\": \":\", \"lines.marker\": \".\"}):\n",
    "    ds.brf_srf.squeeze(drop=True).plot(x=\"vza\")\n",
    "plt.show()\n",
    "plt.close()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "5eafca8c7657d92f106998161825f52f7ae6a82058c49ab253eccd48b77a16ba"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}