05-DA_wofost.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Introduction\" data-toc-modified-id=\"Introduction-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Introduction</a></span></li><li><span><a href=\"#DA-to-the-rescue:-EnKF\" data-toc-modified-id=\"DA-to-the-rescue:-EnKF-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>DA to the rescue: EnKF</a></span></li><li><span><a href=\"#EnKF-for-WOFOST\" data-toc-modified-id=\"EnKF-for-WOFOST-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>EnKF for WOFOST</a></span><ul class=\"toc-item\"><li><span><a href=\"#The-experiment\" data-toc-modified-id=\"The-experiment-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>The experiment</a></span></li></ul></li><li><span><a href=\"#Exercises\" data-toc-modified-id=\"Exercises-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Exercises</a></span><ul class=\"toc-item\"><li><span><a href=\"#No-observations\" data-toc-modified-id=\"No-observations-4.1\"><span class=\"toc-item-num\">4.1&nbsp;&nbsp;</span>No observations</a></span><ul class=\"toc-item\"><li><span><a href=\"#Uncertainty-in-the-model\" data-toc-modified-id=\"Uncertainty-in-the-model-4.1.1\"><span class=\"toc-item-num\">4.1.1&nbsp;&nbsp;</span>Uncertainty in the model</a></span></li></ul></li><li><span><a href=\"#The-type-of-observations\" data-toc-modified-id=\"The-type-of-observations-4.2\"><span class=\"toc-item-num\">4.2&nbsp;&nbsp;</span>The type of observations</a></span></li><li><span><a href=\"#The-properties-of-the-observations\" data-toc-modified-id=\"The-properties-of-the-observations-4.3\"><span class=\"toc-item-num\">4.3&nbsp;&nbsp;</span>The properties of the observations</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ucfajlg/.local/lib/python3.6/site-packages/tqdm/autonotebook/__init__.py:14: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)\n",
      "  \" (e.g. in jupyter console)\", TqdmExperimentalWarning)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import sys, os\n",
    "import copy\n",
    "import datetime as dt\n",
    "import ipywidgets as widgets\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "from wofost_da import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"float:right\">\n",
    "<table>\n",
    "<tr>\n",
    "    <td> <img src=\"figs/nceo_logo.png\" alt=\"Drawing\" style=\"width:200px;height:40px;\"/> </td>\n",
    "    <td> <img src=\"figs/gssti_logo.png\" alt=\"Drawing\" style=\"width:200px;height:40px;\"/> </td>\n",
    "    <td> <img src=\"figs/multiply_logo.png\" alt=\"Drawing\" style=\"width:40px;height:40px;\"/> </td>\n",
    "</tr>\n",
    "    </table>\n",
    "</div>\n",
    "\n",
    "# Data assimilation with the WOFOST model: EnKF filter\n",
    "\n",
    "**Author:** Jose Gómez-Dans (NCEO & UCL)  `j.gomez-dans@ucl.ac.uk`\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "So far, we've had a look at EO data of croplands, and we had a brief play with crop models. Although these topics are major research areas on their own right, it is possible to do a first sum up of what we've seen.\n",
    "\n",
    "**EO data**\n",
    "- Retrieving agronomically meaningful parameters from e.g. reflectance or radiance data is a complex task. A lot of assumptions go into this task. \n",
    "- Perhaps the most interesting parameter is LAI (and is closely related cousing, fAPAR). One reason for its popularity is that LAI is also a parameter used and/or predicted by mechanistic growth models and a major indicator of crop development\n",
    "- Optical EO data is hampered by clouds, which can result in not much data being available (e.g. large data gaps in agronomically meaningful periods). Other properties of the sensor that affect the retrievals are its spatial resolution, its spectral properties\n",
    "- There are other observables from space. These include thermal data (e.g. land surface temperature), passive microwaves (often used for soil moisture). SAR is also used for understanding both vegetation and soil properties. New measurements, such as SIF (solar induced fluorescence) are also becoming practical.\n",
    "- A number of *products* that translate the \"raw\" EO signal into useful parameters have been developed. It makes sense to use them!\n",
    "- ...but always remember to evaluate them to see how they perform. Also be aware that they make make assumptions that are incompatible with other bits of your processing!\n",
    "- EO data only gives a partial window into crops: we infer LAI, maybe soil moisture, but we are interested in biomass or grain yield!\n",
    "\n",
    "**Crop models**\n",
    "- There are many crop models available. We can go from statistical models (that e.g. link anomalies in temperature, precipitation and a vegetation index and yield), to production efficiency models to more complex mechanistic models.\n",
    "- PEM are fun to play around with, but are critically dependent on calibration and having sufficient EO data.\n",
    "- Mechanistic models tend to be preferred, as they encompass a lot of our understanding of crop physiology, soil hydraulics. When we refer to crop models, we'll refer to this type of crop models.\n",
    "- Generally speaking, a crop model takes meteorology, soil and crop type/cultivar parameterisations and calculates the evolution of a crop.\n",
    "- All this requires defining many processes, their interfaces and the parameters that govern such processes and how they are affected by meteorology, soil properties, management decisions, etc.\n",
    "- Calibration of models to local conditions is critical\n",
    "- However, the vast variation of farming practices, local microclimate and soil properties, as well as the many shortcomings of models themselves, of driver data and soil maps, often result in a model is that \"approximately OK\" for some area, but completely wrong for any individual field!\n",
    "\n",
    "DA aims to merge these two independent views of the world, acknowledging the limitations of each. So if we have frequent observations of LAI over a field, how can we inform our crop model of that evidence, in a way that acknowledges the fact that LAI measurements may be uncertain, and that the model is also a bit ropey?\n",
    "\n",
    "## DA to the rescue: EnKF\n",
    "\n",
    "You've already encountered the EnKF. To introduce a common mathematical notation, let $y$ be our observations of e.g. LAI and/or soil moisture. Let our state vector be the stacking of state parameters of the crop model (e.g. LAI, SM, TAGP, DVS, ...). We will assume that our observations are corrupted by zero mean Gaussian noise with some uncertainty, and that we can encode the uncertainty in our state by using a multivariate Gaussian. If we call our state mean vector at a given timestep $\\mu$ and its associated covariance matrix $\\mathbf{P}$, then we can draw a number of realisations of the state using a random number generator. Each of these can be propagated by the crop model (run possibly sampling different soil properties, or sampling maybe uncertainties in meteorological drivers or crop parameters), to provide a population of states after some time. We can then *assimilate* the LAI observations to refine our guess of the state at that point, and go back to running the crop model.\n",
    "\n",
    "<img src=\"http://robertsy.github.io/ensemble-da/assets/img/EnKF.png\" width=\"400px\" />\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EnKF for WOFOST\n",
    "\n",
    "### The experiment\n",
    "\n",
    "Rather than using real data (who needs that, amirite!!!?), we will use a synthetic WOFOST run from where we will derive some observations. This synthetic experiment allows us to compare the result of assimilating the data to the \"truth\", and also allows us to check the effect of the imperfections in the observations.\n",
    "\n",
    "We will consider observations of LAI and SM. As a start, we will define an *ensemble* of models where each member will be run with different parameters and initial conditions. In effect, the ensemble will be made up of our uncertainty in the initial state of the system, as well as the uncertainty in the model parameters. Other factors we could have included would be different meteorology to account for meteo uncertainty. This basically means that each ensemble will run forward in a different way, but the ensemble should capture the typical variability of the crop trajectories we want to observe. This again underlines the necessity of having a relisitic, uncertainty-quantified calibration of the model: DA cannot make miracles!\n",
    "\n",
    "In essence, each of these members will be run forward until an observation is found. At that point, since both LAI and SM are in the observations and in the model, the standard Kalman update formula for the identity observation operator will be used to update the LAI and/or SM value at that timestep of every ensemble member, which will then be run forward using WOFOST, until a new observation is available.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "29bf9046077d41c1a77372deacbb490f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=50, description='Number of ensemble members', max=1000, min=2), Selectio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "run_enkf_widget();\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "The previous framework allows you to do some experiments to test the different parts of a DA system. Next, we suggest a set of experiments. The widget above saves PDFs of the plots that you can get through the [file navigator](./da_plots/) and download for further use. The names are based on the parameter choices.\n",
    "\n",
    "\n",
    "### No observations\n",
    "\n",
    "Perhaps the simplest thing to do is to run an *open loop* simulation. Set a reasonable number of ensembles (>50) to get an idea of the dispersion in the predictions. \n",
    "\n",
    "#### Uncertainty in the model\n",
    "\n",
    "Any model is an approximation to reality. In DA, this is acknowledged directly. Here, our uncertainty comes from the parameters and initial conditions of the simulation. For a number of parameters, we draw random values given a particular distribution. While you cannot change the shape of the distribution (it is Gaussian), we can change the standard deviation (the *dispersion* bit). Have a play with that, with no assimilation and see what the ensemble looks like.\n",
    "\n",
    "### The type of observations\n",
    "\n",
    "It'd be dissappointing if after all this talk of data and so on, we didn't give you observations!!! You can assimilate LAI, SM or both simultaneously. See what effect assimilating either or both observational streams has in the results.\n",
    "\n",
    "\n",
    "\n",
    "### The properties of the observations\n",
    "\n",
    "The number and timing of the observations is crucial: clearly if your observations fall outside the period of interest, then you have a problem! You can see reduce the number of observations (they're selected randomly) to see the effect of gappiness in the data. Remember that e.g. for S2, the gaps in the rainy season are substantial.\n",
    "\n",
    "Another important characteristic of observations is that they are inferences, often fairly far removed from the actual measurement. As such, they are *uncertain*. You can change the relative unertainty of each observatinal stream, making the data more or less certain.\n",
    "\n",
    "Consider runs using many observations but very uncertain, and runs using few observations but very precise.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\">Creative Commons Attribution-NonCommercial 4.0 International License</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}