MicaSense Image Processing Tutorial 3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 3 - DLS Sensor Basic Usage\n",
    "\n",
    "\n",
    "This tutorial will walk through how to convert RedEdge data from raw images to radiace and then use the DLS information in the image metadata to convert that to reflectance.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Update for DLS2\n",
    "\n",
    "This tutorial was originally written as a guide to basic processing of DLS1 data for RedEdge-3 and RedeEdge-M.  Since October 2018, Altum and RedEdge-MX have been shipped standard with DLS2.  DLS1 is a small, red square with a single diffuser on top.  DLS2 is a larger, thinner black rectangle with two small diffusers on top and 8 small diffusers on different up-looking faces.\n",
    "\n",
    "For DLS2 data, we recommend using the `Capture` class to access the `Capture.dls_irradiance()` method. This will provide a compensated horizontal irradiance useful for radiometrically correcting imagery. We leave the below intact for legacy users and as a tutorial on remote sensing calibrations and terms.  \n",
    "\n",
    "For RedEdge and Altum with up-to-date firmware, the `EXIF:HorizontalIrradiance` tag and the appropriate scale factor will provide the necessary information for reflectance compensation. If necessary, this can be further field-calibrated by using a pane image to tie the DLS2 and camera calibrations together in-situ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remote Sensing Theory\n",
    "\n",
    "In [Tutorial 1](./MicaSense%20Image%20Processing%20Tutorial%201.html) we covered conversion of the image radiance to reflecance using the average radiance from an area of a specially made lambertian reflectance panel.  As we did not get very far into remote sensing theory, we will cover that some here as we move on to the more complex problem of using the data from the RedEdge Downwelling Light Sensor (DLS).  To get started, we will review some of the fundamental measurements in remote sensing.\n",
    "\n",
    "## Definition of Terms\n",
    "\n",
    "\n",
    "The DLS measures the energy coming from the sky, while the camera measures the energy coming from the fields below.  The property that we are generally looking to estimate in agricultural remote sensing is the surface reflectance of the plants and ground below the sensor. \n",
    "\n",
    "The MicaSense RedEdge measures spectral radiance, a measure of the power incident on the CMOS image sensor at each pixel, and has units power per unit area per unit bandwidth per unit solid angle (W/m^2/nm/sr). We will call the radiance image L.\n",
    "\n",
    "The MicaSense Downwelling Light Sensor (DLS) measures incoming spectral irradiance on the top surface of the device for the same 5 spectral bands as the camera. Irradiance has the units of power per unit area per unit bandwidth, or Watts per meter squared per nanometer (W/m^2/nm). We will call this downwelling quantity I.\n",
    "\n",
    "If we look closely at the units, they differ by the last term - the unit over solid angle. This is because each pixel of the camera sensor only sees a specific solid angle, while the dls detector has a specific area.  \n",
    "\n",
    "\n",
    "## Understanding the measurements\n",
    "\n",
    "In order to use downwelling irradiance and upwelling radiance to determine reflectance, we simply need to divide pi times the irradiance by the radiance.  Our units cancel out, and we are left with the unitless quantity of reflectance.  However, there is one caveat.  In most systems, the DLS is not held level to the ground at all times, so the amount of light that shines on it's surface changes as the aircraft moves about.  In oder for that simple math to work, it's necessary that the irradiance and radiance be compensated for a very important property: the illumination angle.  \n",
    "\n",
    "Both the radiance and irradiance are sensitive to changes in the direction of the lighting. Both quantites decrease with the cosine of the angle from the surface normal, eventually reaching zero when the light source is parallel to the surface. So when measuring irradiance, this effect needs to be corrected in order to estimate the actual irradiance.\n",
    "\n",
    "## Different types of irradiance\n",
    "\n",
    "To complicate matters further, the light reaching the DLS and the ground is not coming from a single source, even under sunny conditions.  The direct sunlight provides what we will call the direct source, resulting in a direct irradiance `I_direct`.  The rest of the (usually blue) sky provides what we term the indirect, or diffuse, irradiance which we will call `I_diffuse`.  Under even the sunniest conditions with the summer sun directly overhead, only about 85% of the light is direct, while the other 15% is diffuse, which we will express as the ratio 6:1.  Under thin cloud conditions, this ratio drops to close 1:1, meaning half of the measured irradiance is coming from parts of the sky other than the direction of the sun. The ratio goes even lower under fully cloudy skies, becoming fully diffuse.\n",
    "\n",
    "As we ultimately want to estimate the time-varying irradiance, we will need to estimate both the direct and diffuse irradiance terms as well as correct for the cosine effect of the DLS sensor itself.\n",
    "\n",
    "## Sensor-specific challenges\n",
    "\n",
    "There's one final piece of theory required that pertains to irradiance sensing, and that is the Fresnel Effect.  The Fresnel Effect describes a physical property of the diffuser that gathers light on top of the DLS.  This effect results from the fact that much of the light striking the diffuser is reflected, and the amount reflected changes with the angle of the incident light.  So, instead of our diffuser behaving as a perfect cosine detector, we need to also compensate for the material imperfections.\n",
    "\n",
    "## DLS pose estimation\n",
    "\n",
    "An onboard orientation sensor provides a fairly rough version of the earth-fixed orientation of the DLS.  The quality of this orientation varies based on the quality of the sensor installation and calibration, but in even the best cases it has a few degrees of error. If the DLS and the camera are mounted fixed to one another, the photogrammetric workflow can be used to augment the orientation of the irradiance sensor.  Our examples here will use the sensor measurements directly, but significant improvements can be made using the very precise angles photorammetry can provide.\n",
    "\n",
    "## Radiometric calibration options\n",
    "\n",
    "As we've discussed in the first tutorial, the method of using a calibrated reflectance panel is the scientific standard for remote sensing measurements. Spectroradiometers which cost tens of thousands of dollars are calibrated in this method in the field.  Panels images may also be used both before and after each flight to track slowly changing solar irradiance. In many cases of direct sunlight with no clouds, panel images can provide the best results when taken before and after each flight. As the diffuse lighting increases, the utility of the DLS increases as well.\n",
    "\n",
    "In cases of more varied irradiance, the DLS can be used for irradiance estimation, using only the internal pose estimation.  The DLS irradiance can be be used in conjunction with panels to ensure that the DLS and panel measurements match at the time of panel capture.  In order for this to be effective, it's important that the DLS and panel are exposed to the same irradiance.  This is best accomplished by ensuring panel images are taken with the operator's back to the sun and the operator shadow next to the panel. The aircraft and operator should block as little of the sky as possible.\n",
    "\n",
    "Finally, in fixed mount systems, the photogrammetric angles can be used, along with the offset between the camera and DLS sensors, to most accurately estimate irradiance. This can again be coupled with one or more panel measurements for best radiometric results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Putting it all Together\n",
    "\n",
    "Below, we will use the `micasense.Capture` class to load a set of images and use tilt-compensated DLS irradiance data to roughly translate those images to reflectance. We will then compare these to the known reflectances of a lambertian panel, and develop a site-specific compensation factor to tie together the DLS irradiance and camera radiance measurements through the panel.\n",
    "\n",
    "## Accessing and processing DLS data\n",
    "\n",
    "Below, we show how to access DLS data and compute the minimum required corrections to the DLS data for use.  The code in this cell is present in the `micasense.dls` library, and is used by the `Image` and `Capture` class to expose corrected DLS data.  If your intent is to simply use this code to correct your data, you can skip this section and access pose-corrected DLS data through the `capture.dls_irradiance()` method.\n",
    "\n",
    "For the purpose of this tutorial, we assume a constant irradiance (over time) and a ratio of direct to diffuse irradiance of 6:1.  These assumptions only hold approximately, and only for clear sky. Cloud cover can lead to dramatically different results. \n",
    "\n",
    "## Loading Captures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import micasense.dls as dls\n",
    "\n",
    "import os, glob\n",
    "import micasense.capture as capture\n",
    "\n",
    "images_path = os.path.join('.','data','0000SET','000')\n",
    "image_names = glob.glob(os.path.join(images_path,'IMG_0000_*.tif'))\n",
    "cap = capture.Capture.from_filelist(image_names)\n",
    "# set panel corners manually if zbar is not installed\n",
    "panelCorners = [[[809,613],[648,615],[646,454],[808,452]],\n",
    "                [[772,623],[613,625],[610,464],[770,462]],\n",
    "                [[771,651],[611,653],[610,492],[770,490]],\n",
    "                [[829,658],[668,659],[668,496],[829,496]],\n",
    "                [[807,632],[648,634],[645,473],[805,471]]]\n",
    "\n",
    "cap.set_panelCorners(panelCorners)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing Solar orientation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define DLS sensor orientation vector relative to dls pose frame\n",
    "dls_orientation_vector = np.array([0,0,-1])\n",
    "# compute sun orientation and sun-sensor angles\n",
    "(\n",
    "    sun_vector_ned,    # Solar vector in North-East-Down coordinates\n",
    "    sensor_vector_ned, # DLS vector in North-East-Down coordinates\n",
    "    sun_sensor_angle,  # Angle between DLS vector and sun vector\n",
    "    solar_elevation,   # Elevation of the sun above the horizon\n",
    "    solar_azimuth,     # Azimuth (heading) of the sun\n",
    ") = dls.compute_sun_angle(cap.location(),\n",
    "                      cap.dls_pose(),\n",
    "                      cap.utc_time(),\n",
    "                      dls_orientation_vector)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correcting DLS readings for orientations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since the diffuser reflects more light at shallow angles than at steep angles,\n",
    "# we compute a correction for this\n",
    "fresnel_correction = dls.fresnel(sun_sensor_angle)\n",
    "\n",
    "# Now we can correct the raw DLS readings and compute the irradiance on level ground\n",
    "dls_irradiances = []\n",
    "center_wavelengths = []\n",
    "for img in cap.images:\n",
    "    dir_dif_ratio = 6.0\n",
    "    percent_diffuse = 1.0/dir_dif_ratio\n",
    "    # measured Irradiance / fresnelCorrection\n",
    "    sensor_irradiance = img.spectral_irradiance / fresnel_correction\n",
    "    untilted_direct_irr = sensor_irradiance / (percent_diffuse + np.cos(sun_sensor_angle))\n",
    "    # compute irradiance on the ground using the solar altitude angle\n",
    "    dls_irr = untilted_direct_irr * (percent_diffuse + np.sin(solar_elevation))\n",
    "    dls_irradiances.append(dls_irr)\n",
    "    center_wavelengths.append(img.center_wavelength)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "plt.scatter(center_wavelengths,dls_irradiances)\n",
    "plt.xlabel('Wavelength (nm)')\n",
    "plt.ylabel('Irradiance ($W/m^2/nm$)')\n",
    "plt.show();\n",
    "\n",
    "cap.plot_undistorted_reflectance(dls_irradiances)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tying the camera and DLS together via Panels\n",
    "\n",
    "Finally, if we want to use the panel as the ultimate reference tying the panel and the DLS readings together with a known lambertian reflector, we can calculate a DLS factor which we can then apply to DLS derived irradiances over the whole flight.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "panel_reflectance_by_band = [0.67, 0.69, 0.68, 0.61, 0.67] #RedEdge band_index order\n",
    "\n",
    "panel_radiances = np.array(cap.panel_radiance())\n",
    "irr_from_panel = math.pi * panel_radiances / panel_reflectance_by_band\n",
    "dls_correction = irr_from_panel/dls_irradiances\n",
    "cap.plot_undistorted_reflectance(dls_irradiances*dls_correction)\n",
    "\n",
    "plt.scatter(cap.center_wavelengths(), cap.panel_reflectance())\n",
    "plt.title(\"Panel Reflectances\")\n",
    "plt.xlabel(\"Wavelength (nm)\")\n",
    "plt.ylabel(\"Reflectance\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Future improvements\n",
    "\n",
    "While the need to compensate the DLS for movement increases the complexity of the problem, it also provides us more information than might be available for other diffuser types. This is because as the DLS moves about on the aircraft, it provides many measurements per second that should all have consistent effects assuming the light isn't changing over those very fast measurements. This movement can then be used to develop an irradiance series which uses all of the available information from the on board sensors.  A future tutorial will provide an overview of some methods for improving irradiance sensing further."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "Copyright (c) 2017-2019 MicaSense, Inc.  For licensing information see the [project git repository](https://github.com/micasense/imageprocessing)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.1"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {
    "height": "1002px",
    "left": "0px",
    "right": "1117px",
    "top": "107px",
    "width": "212px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}