{ "cells": [ { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "9363186e-7292-46fa-b9b9-ab7233a0d2ee" }, "slideshow": { "slide_type": "slide" } }, "source": [ "# Assimilating EO data into an ecosystem model with a particle filter " ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "e9923300-1d60-4b9a-98f5-441363bc8216" }, "slideshow": { "slide_type": "slide" } }, "source": [ "* The aim of today's practical is to explore the blending of **observations** and **models**\n", "* You will use the **DALEC ecosystem model**, that tracks the fate of C through photosynthesis, biomass allocation and decomposition.\n", "* We will use **observations of LAI** from space-borne sensors to keep the model in track\n", "* The observations will be combined with the model using a type of **particle filter** (a **sequential Metropolis-Hastings filter**)\n", "* We will focus on the Metolius Young Pine FLUXNET site, and will compare assimilation results to flux tower data.\n" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "dfa52bce-5f45-4cd7-97b1-6a9f335a846a" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## The DALEC model\n", "The DALEC model is a simple box model of C cycling.\n", "![DALEC](http://ars.els-cdn.com/content/image/1-s2.0-S0034425707003276-gr2.jpg)" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "749f2688-434c-4dba-926a-0a75b7c2fa9e" }, "slideshow": { "slide_type": "slide" } }, "source": [ "1. Carbon is acquired through the GPP box, which uses the very simple ACM model to model photosynthesis as a function of $LAI$, and temperature, incoming radiation, $CO_2$ concentration and VPD.\n", "2. Part of the GPP is lost as **autototrophic respiration**, to get NPP.\n", "3. The NPP is allocated to the **foliar**, **woody** and **roots** pools.\n", "4. The foliar and root pools loose C to the **litter** pool.\n", "5. The woody biomass pool looses C to the **SOM/Woody debris** pool.\n", "6. The litter pool decomposes into the **SOM** pool\n", "7. Both **litter** and **SOM** pools release C back to the atmopshere through microbial activity and **heterotrophic respiration**" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "e7484ffc-4770-4bad-af58-23a2206602b6" }, "slideshow": { "slide_type": "slide" } }, "source": [ "* The DALEC model has been calibrated for the Metolius site\n", "* We have a fairly rough understanding of the different C pools at around 2000.\n", "* However, the DALEC model is **very simplistic**!!! \n", "* The idea is that this simple model can be made to track observations, and hence combine heterogeneous observational streams into a consistent story of what's happening to C." ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "f1459bcd-1fac-4d48-99e8-c92973cf78bf" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## The MODIS observations" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "88ebb92f-326a-4e88-9e16-5c9bf3bc2d39" }, "slideshow": { "slide_type": "slide" } }, "source": [ "* We will use the MODIS LAI product, which provides both an estimate of LAI and of uncertainty in LAI.\n", "* This product takes the **red** and **nir** reflectance, and maps them to LAI.\n", "* The uncertainty of the LAI product is controlled by things like cloudiness, snow, etc." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2018-03-23T10:48:13.067257Z", "start_time": "2018-03-23T10:48:11.532151Z" }, "nbpresent": { "id": "b6c95724-8651-4e73-bddb-ea78264de730" } }, "outputs": [ { "data": { "text/plain": [ "Text(0,0.5,'LAI $m^2m^{-2}$')" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAJiCAYAAABHDf8zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xu4HXV97/HP7IRA3QQbglTpichF\nv8Z6KUrPYyTV9GgsB602VVK8pdqq2CvFemhr9MCxpbejLW29FNGqecSmoKa2iKdNtdEGtrR4AuJJ\n+FWIIlAftAFryGMMyV7nj1k7WXvvmbXm8pvLb+b9eh4ediZrzXzXb+Y3+5tZsz4rGgwGAgAAQDim\nmi4AAAAA+dDAAQAABIYGDgAAIDA0cAAAAIGhgQMAAAgMDRwAAEBgaOAAAAACQwMHAAAQGBo4AACA\nwNDAAQAABIYGDgAAIDA0cAAAAIGhgQMAAAgMDRwAAEBgaOAAAAACQwMHAAAQmKVNF+CDmQ2cc02X\nAQAAkEVUdgVcgQMAAAgMDRwAAEBgaOAAAAACQwMHAAAQGBo4AACAwNDAAQAABIYGDgAAIDA0cAAA\nAIGhgQMAAAgMDRwAAEBgaOAAAAACQwMHAAAQGBo4AACAwNDAAQAABIYGDgAAIDA0cAAAAIGhgQMA\nAAgMDRwAAEBgaOAAAAACQwMHAAAQGBo4AACAwNDAAQAABIYGDgAAIDA0cAAAAIGhgQMAAJ0XRZGi\nKGq6DG9o4AAAAAJDAwcAABAYGjgAAIDA0MABAAAEhgYOAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo\n4AAAAAJDAwcAABAYGjgAAIDALG1io2b2Rkm/LekRSb/unLuxiToAAABCVHsDZ2ZnSPpVSU+TdKqk\nv5dEAwcAAJBRE1fgXiLpw865hyU9bGYXNlADAABAsJpo4M6WNGtmt0g6XtJvNlADAABAsJpo4H5A\n0mmSnidplaTPm9npzrlHFj7QzNZJWldrdQBQQhRFkqTBYNBwJQC6rIkG7iFJX3LOHZT0VTO7V9Jj\nJd278IHOuR2SdkxaoZld7rlGAACA1moiRuRzks43syVmdpqkR0v69wbqAAAACFLtDZxz7jOSbpP0\n/yR9WtLFzrkjddcBAAAQqqgL92mY2cA513QZAMA9cEBLtWxuRmVXwDcxAAAABIYGDgAAIDA0cAAA\nAIGhgQOAFomi6Oi9OgCQhgYOAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo4AAAAAJDAwcAABAYGjgA\nAIDA0MABAAAEhgYOAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo4AAAAAJDAwcAABAYGjgAAIDA0MAB\nAAAEhgYOAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo4AAAAAJDAwcAABAYGjgAAIDA0MABAAAEhgYO\nAAB4FUWRoihquoxOo4EDAAAIDA0cAABAYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAw\nNHAAAACBoYEDAAAIDA0cAABAYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACB\noYEDAAAIDA0cAABAYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACBoYEDAAAI\nDA0cAABAYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACBoYEDAAAIDA0cAABA\nYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACBoYEDAAAIDA0cAABAYGjgAAAA\nArO0iY2a2a2SThz+8Tbn3EVN1AEAABCi2hs4M1si6bvOuXPr3jYAAOif2T0zunvzRq1aMa0j11ym\naO0GTa1e03RZpTRxBe40Sd9sYLsAAKBnZvfMaLB9i04/efjG3/59Gmzfolkp6CauiXvgTpf0DDO7\n3cxuMbPzGqgBAAD0wGDnNunwofkLDx+KlwesiStwByR9QNJ7JD1V0t+Y2dnOuUcWPtDM1klaV2t1\nAACgO/bvy7c8EE00cHdKusM5d1jSLjN7QNIPSbpv4QOdczsk7Zi0QjO73HONAACgC5avTG7Wlq+s\nvxaPmngL9RJJfyhJZnaWpEdL+vcG6gAAAB0Xrd0gLV02f+HSZfHygDXRwL1b0hPN7C5J10n6Befc\nbAN1AACAjptavUbR+k2658GHNTsYSMtXKlq/KegPMEhSNBgMmq6hNDMbOOeaLgMAFEWRJKnoubXs\n84E2aONx3LKaorIr4JsYAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACBoYEDAAAI\nDA0cAABAYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACBoYEDAAAIDA0cAABA\nYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAFARaIoUhRFTZcBoINo4AAAAAJDAwcAABAYGjgAAIDA0MAB\nAAAEhgYOAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo4AAAAAJDAwcAABAYGjgAAIDA0MABAAAEhgYO\nAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo4IAWiaJIURQ1XQYAtF7fz5c0cAAAAIGhgQMAAAgMDRwA\nAEBgaOAAAAACQwMHAAAQGBo4AACAwNDAAQAABIYGDuiIvmciAU1i/qFuNHAAAACBoYEDAAAIDA0c\nAABAYGjgAAAAAkMDBwCezO6Z0d2bN+rQO1+nI9dcpovOObPpkgB01NKmCwCALpjdM6PB9i06/eQT\n4wX79+nqjec1WxSAzuIKHAB4MNi5TTp8aN6y6WXH6coLzm2oIgBdRgMH1Iy8qI7avy9x8aoV0zUX\nAqAPaOAAwIflKxMX3/vQgZoLAdAHNHAA4EG0doO0dNm8ZQcOPaLNN97aUEUAuowGDgA8mFq9RtH6\nTbrnwYc1OxhIy1fq4utu0tZde5suDUAHRYPBoOkaSjOzgXOu6TKATObuf0uae+P+rsx6UZ/R/VBk\nn7Afw8R+m6+O8ci7jZbto9I3QnMFDgAAIDA0cAAAAIGhgQMAAAgMDVwB5HgBAIAm0cABAAAEhgYO\nAAAgMDRwAAAAgaGBAwAACAwNHAAAQGBo4AAAAAJDAwcAABAYGjgAANBJXc5tpYEDAAAIDA0cAABA\nYGjgAAAAAtNYA2dmU2b2RTM7v6kaAAAAQtTkFbhLJD2xwe0DAAAEqZEGzszOkLRe0g1NbB8AACBk\nTV2B+3NJl0oaNLR9AACAYC2te4Nm9vOSbnHOOTOb9Nh1ktbVUBYAAGjIXFbbYMB1nayiugfLzLZI\nerakWUmPk/RdSa91zn22xDoHzjlPFU7GgYYyxh0/ZY4tjst2GN0PRfYJ+zFM7Lf58o5HVXOl7Hys\nUOl04dqvwDnnNs39bGYflrS1TPMGAADQN+TAAQAABKb2K3CjnHOvbXL7AAAAIWq0gQMAHDO7Z0Z3\nb96oVSumdeSayxSt3aCp1WuaLgtAC9HAAUALzO6Z0WD7Fp1+8onxgv37NNi+RbMSTRyARbgHDgBa\nYLBzm3T40PyFhw/FywFgARo4oIOiKDr6kXkEYv++fMsBFNaFcyQNHAC0wfKV+ZYD6DUaOABogWjt\nBmnpsvkLly6LlwPAAjRwANACU6vXKFq/Sfc8+LBmBwNp+UpF6zfxAQYAifgUKgC0xNTqNTrryudI\nas3X/QBoKa7AAQAABIYGDgAAIDA0cAAAAIHhHjgAABC0uUy3cfeOdu2+Uq7AAQAABIYGDgAAIDA0\ncAAAAIGhgQMAAAgMDRwAAEBgaOAAAAACQwMHAAAQGBo4AACAwNDAAQAABIYGDgAAIDA0cAAAAIGh\ngQMAAAgMDRwAAEBgljZdAAAA6I7ZPTO6e/NGrVoxrSPXXKZo7QZNrV7TdFmdQwMHAAC8mN0zo8H2\nLTr95BPjBfv3abB9i2YlmjjPeAsVAAB4Mdi5TTp8aP7Cw4fi5fCKBg4IWBRFiqKo6TIAILZ/X77l\nKIwGDgAA+LF8Zb7lKIwGDgAAeBGt3SAtXTZ/4dJl8XJ4RQMHAAC8mFq9RtH6TbrnwYc1OxhIy1cq\nWr+JDzBUgE+hAgAAb6ZWr9FZVz5HkjQYDBqupru4AgcAABCYsVfgzGyJpJdKmpV0g3Pu8HD5hc65\n62uoDwAAAAtMugK3RdIzJf2opJ1mdvZw+S9WWhUAAABSTboH7jTn3Kskycw+IulDZnZF5VUBQEDm\nsvi43wdAXSZdgTvezI6XJOfc1yT9lKS3SHpq1YUBAAAg2aQG7s2SVsz9wTm3X/E9cW+usigAAACk\nG/sWqnPui6N/NrNTnXPfkvTRSqsCAABAqrwxIh+vpAoAAABklreB41uzAQAAGpa3geMjVgAAAA3j\nmxgAAAACw1uoDYqi6Gh+FAAgLJzD0aS8DdzLK6kCAAAAmeVq4JxzD1RVCAAAALKZ9FVai5jZT0t6\n2/C5d0q6TdLtkm5zzn3Tb3kAAABYKHcDJ+n9kn5X0lckPUXS0yVtkPQjkk70VxoAAACSFGngvifp\nPc65I5I+N7fQzLiTEwAAoAZFYkTeJemXFi50zpERBwAAUIMiV+C2S/qMma2XdIPie+C+7Jw76LUy\nAAAAJCpyBW6bpFsVN24vlPQxSfvNbLfPwgAAQL/VnbUXUrZfkStwp0q6cPQtUzNbrvjDDAAAAKhY\nkQZuq6TnSdoxt8A5t1/STZ5qAgAAwBhF3kI9U9L1ZnapmT3Jd0EAAAAYr8gVuOslOUkvlfR2MztO\n0h2Kg3wXfToVAAAAfuVu4JxzHxz9s5k9XtIzhv8BAACgYkWuwM3jnPuGpG9I+rvy5QAAAGCSsQ2c\nmS1R/FbprKQbnHOHh8svdM5dX0N9AAAAWGDShxi2SHqmpB+VtNPMzh4u/8VKq0LrhJSNg3Tsx35i\nvwPdM+kt1NOcc6+SJDP7iKQPmdkVlVcFAACAVJOuwB1vZsdLknPua5J+StJbJD216sIAAACQbFID\n92ZJK+b+MAzsfelwOQAAABoQDQaDyY9qOTMbOOdq297cvSRlx87XeuoQUq1tN24s847z6OPTfi6z\nfmSTdz9kXVcV9cGfsvu6y8qcy7Kup8w5rwXnyNI3peaOETGzn5b0tuFz71T8pfa3Kw7y/WbZggAA\nADBekRy490v6XUlfkfQUxV9iv0HSj0g60V9pAAAASFKkgfuepPc4545I+tzcQjPjM+oAAAA1KPJl\n9u+StOg7T51zvPEPAAWQ0xYO9hXaosgVuO2SPmNm6yXdoPgeuC875w56rQwAAACJilyB2ybpVsWN\n2wslfUzSfjPb7bMwAAAAJCtyBe5USReOvmVqZssVf5gBAAAAFSvSwG2V9DxJO+YWDAN+b8ry5GGz\n9zFJT5J0QNIbnHNfKlAHAABALxV5C/VMSdeb2aVm9qQCz3+zpC8650xxntz/KrAOAGi12T0zunvz\nRh165+t05JrLNLtnpumSgM5YOL8uOufMpkuqXZErcNdLcoq/UuvtZnacpDsUB/ku+nRqgn+QtHf4\n8w9K+s8CNQBAa110zpkabN+i008eRmPu36fB9i2alTS1ek2jtQGhm90zs2h+Xb3xvGaLakDuBs45\n98HRP5vZ4yU9Y/hflufPDJ+3a/ic8/PWAABtduUF50qHD81fePiQBju3STRwQCmDndsWza/pZcfF\n865HilyBm8c59w1J35D0dzmfd46ZPU/SX0o6K+kxZrZO0rqSJQJArVatmE7+i/37JPHdpEApw3m0\nUOq8myDU+Vi6gcvLzP5E0rucc/c55z5vZieYWZQUBOyc26GRD0uMWefl/isFgGLufejAsbd3Ri1f\nWX8xQNcsX5nYxN370AH16U64Ih9iKGuJpI2SZGZrJH2db3EA0CWbb7xVWrps/sKlyxSt3dBMQUCH\nRGs3LJpfBw49Es+7HmmigbtS0ovM7N8Ufy3XGxuoAQAqs3XXXkXrN+meBx/W7GAgLV+paP0mPsAA\neDC1es2i+XXxdTdp6669k5/cIVHW93zN7LmTHuOc+0Lpigows4Fzrrbt+Xq/PKT33UOqte3GjWXe\ncR59fNrPZdaPbMruhyzP8VUfysm7r/uqzLks63rKnPNasB9Lf6FunnvgPjLh7wdSr95+BgAAaETm\nBs45d0aVhQAAACCb0p9CNbPjJf2UpFc4515WviQAAACMU6iBM7Mlkl4o6RWKv5HhQeXMgQMAAEAx\nuRq44QcZXinpZZK+K+k0Sec75z5fQW0AAABIkDlGxMzuk/RRSQ9LusA5d9bw569WVBsAAAAS5MmB\n+5qkEyQdr/gTpxr5PwAAAGqSuYFzzv24pHMlfVPSFjO7U9JyER0CAABQq8xBvguZ2bMkbZL0s5K+\nI+kTzrnNHmvLUwtBvhULqda2I8i3ewjy7Y+kfXVk9836+rVXadWKaU2ddIqitRt6/60bBPlOVDrI\nt/BXaTnnvuScu0TSD0v6DUmXlC0GAICQXHTOmRps36LTTz5RU1Ek7d+nwfYtmt0z03Rp6LjS34Xq\nnDsi6TZJP1C+HAAAwnHlBedKhw/NX3j4kAY7tzVTEHqjiS+zR8WiKDp6GRjNYT+gKhxb7bFqxXTy\nX+zfV28h6B0aOAAACrr3oQPJf7F8Zb2FoHdo4AAAKGjzjbdKS5fNX7h0maK1G5opCL2R+ZsYzOx7\nSs9941o+AKB3tu7aq2uv3cSnUFG7PF+l9eTKqgAAIFBTq9forCufI4moFtQncwPnnLunykIAFDe7\nZ0Z3b96oVSumdeSay3TROWdq6669TZcFAKhIri+zB9A+ozlUkqT9+3T1xvOaLQoAUCk+xAAELimH\nanrZcfFyAEAn0cAhNzKo2iUthyo1nwqNYN4A8IkGDghcWg5Vaj4VACB4NHBA4JJyqA4ceiReDgDo\nJBo4IHBbd+1VtH6T7nnwYc0OBtLylbr4upv4FCoAdFjUhcwaMxs452rb3tx9LGXHztd66ljv6Dqr\nqrtr0sZp3PjlHdu0/VJk2yguy34Yt3/yzq8yxwnKyTvn+srXMZr1fJn3+S3Yj6VviOUKHAAAQGC4\nAlcAV+D4l2YWef9FOLtnJvfX8XAFrh24AtcfLbhyE4Q2XIFLO6cuXP7q931CW3ftDe4KHEG+QAvM\n7plZFMY72L5FsxLfqQgAOaWdU4/cf5e0++ZOBJ/zFirQAoOd2xaF8erwoXg5ACCXtHOq7vhCZ4LP\nuQKHUnjbwJP9+/ItR2sxF4AWSDt3DmYTF4cYfM4VOKANlq/MtxwAkC7t3Bkltz0hBp/TwAEtEK3d\nsCiMV0uXxcsBALmknVP1tOd2JvicBg5oganVaxaF8UbrN/EBBgAoIO2cuuQFr+lM8DkxIgUQI8JH\n57MgyLc/yoaFEiMSDs6F2bQhRiTrekMN8uVDDDnN7pnR3Zs3atWKaR255rJMWV3AHI4fAIAPNHA5\nkNWFMjh+AAC+cA9cDmR1oQyOHwCALzRweQSY1RVF0dH39Lu4vaAEePyge5ijqArHVr1o4PIgqwtl\ncPwAADyhgcuBrC6UwfEDAPCFBi4HsrpQBscPAMAXcuAKCCkHropas+ZZ9R05cP0RUg4cx0A5LcgP\na60y40EOXIEVdOFgo4FLXy8NXHOqbuBm98zo69depVUrpjV10il69fs+oa279tLANYAGrj9a8Iu/\ntYqOx8Jz2Wg+Zpbz5ZHdN+d+fgv2I0G+QB8lZcpdvfG8ZosCgJzK5mNedM6Zvc3X5B44IEBJmXLT\ny47TlRec21BFAJBf2XzMKy84t7f5mjRwgUnL2SF/p2dSsuNWrZiuuZD2YA6gCI6bY7KOxejjSo9f\nyXzM1HNeD/I1aeCAEKVkx9370IGaCwGAEkrmY6ae83qQr0kDBwQoKVPuwKFHtPnGWxuqCADyK5uP\nufnGW3ubr0kDBwQoKVPu4utu0tZde5suDQAyK5uPuXXX3t7maxIjUkCTMSJNfSSaGJH8yIGrV5Ov\njRiRcFU950KS9TVkGQNy4CYiRgTAfLN7ZnT35o1atWJaR665bF4mEgCgG2jggA7pcyYSAPQJ98AB\nHdLnTCQA6BMauICRX4SF+pyJhP4Yd+7zmlHWA1nHEu1DAwd0SJ8zkQCgT2jggA7pcyYSAPQJDRzQ\nIX3ORAKAPuFTqEDHTK1eo7OufI6ksDOpAADpaOAAAEDvhZahSQMHAAB6LcQMTe6BAwAAvRZihiYN\nXM2K5OpUkcWTts6+5v6Ued19HbO2Y7+EqY37rS01NXnebssYVCXEDE0aOAAA0GshZmjSwAEAgF4L\nMUOTBg4AAPRaiBmafAoVAAD0XmgZmjRwDQktbwZAsrJzeeHzLzrnTG3dtbfCiuED+w1No4FrQIh5\nMwAWm90zU2ouJz3/6o3nVVcwvEg6h7PfUDfugWtAiHkzABYb7NxWai4nPX962XHxOQKtlXQOZ7+h\nbjRwFRmXmdNU3kzeHJ+25/70Ncuu668vKGlzNutcTnlc6jkihyaPkzq3nXVbPmtK2z9l9hvzGnnR\nwDUgxLwZAAnS5mzWuZzyuNRzBFohbf+w31AnGrgGhJg3A2CxaO2GUnM56fkHDj0SnyPQWknncPYb\n6kYD14AQ82YALDa1ek2puZz0/Iuvu4lPM7Zc0jmc/Ya6RXVnnZjZlKQPSnq+pIck/ZJz7qaS6xw4\n53yUl8ncfQrjxi7tMaPLs6xn3HN8/exrW1nHxpcsr6HsuqqoY9y28h5bZfZv17T9WMw6b8qcX4ps\nuwpVzKe0x2QdyzJzosycy8rX75Uy2y0ylr7Go+z5ssx5uO75MVdC2RU0ESOyQdLJkk6XtFrSVklP\nb6AOZETeUTeQPdgOSfuhzu0xfxGKOs5ZIc+PJhq4x0r6iHNuIGm3mT3GzKLhn9Ey5B11Q9m8MviR\nlgFZ1S8N5i9CVcc5K/T5Ufs9cM659zjnPilJZvYGSV+leWsv8o66oWxeGfxIy4Csaj4xfxGqOs5Z\noc+PRr6JwcymJV0laY2kl4x53DpJ6+qpqpgq3i9v031KVeQdhapN+yW3snllLVLH/SpVrbfu+RT6\n/A16zjWkjWNW5H64Q+98naaScvE8nrNCnx+1N3Bm9ihJ/yzp85LOdc4dTHusc26HpB0Z1nm5r/ow\n370PHTh2eXnB8jMbqAcFLV+ZfOIje7BW4+ZT3dtj/qLN0o5dn+es0OdHEzEib5K00zl36bjmDe1A\n3lE3lM0rgx9pGZBVzSfmL0JVR15q6POjiQbuWZI2mNmdc/81UAMyIu+oG8rmlcGPtAzIquYT8xeh\nqiMvNfT5UXsOXBWazIHzlT0z7rnkwCWrOwfOV94bOXDl1HEsVpFhRQ5c/ucWOY/2PQeuzLFbx++e\nvHWQAzdmBSGfyOfQwB37+cjum/X1a6/SqhXTmjrhRH173z6tnD5eUyedole/7xPaumsvDdwCs3tm\njo3ZSaekZg3RwLXDxON+uA+XPOU5Rx9XZhtlHl9FA1fkeO1jA7fwmBg9/+WtNcs6u9bALTzO5l5r\n2hiUHRsauPwa+RQqqrEo0+b7B/SYE0+Ifw4s36Yu5KOFr+5stSZxvGZTRb5X6JlheSQdZ1dvPE9r\nnnBqruWoFt+F2iGJGVMjQsq3qQv5aOGrO1utSRyv2VSR7xV6ZlgeScfZ9LLjdPGaJ+da3sWxaROu\nwGVU9WVVH+vPkl0TSr5NbTzno3XhLcoiqnjdWdeZN8upobdLMhtbS4B5fk2McRX5Xk1lhjVyjKYc\nT0umkt/1S1s+aWyqem1tnNdV4Apch2TJkqoqbypYaZlC5KMFI+2Y7uSxzvGaSRXHBMeZdGQ2uSFK\nW97JsWkRGrgOSczNGRFSvk1dyEcLX93Zak3ieM2minyv0DPD8kg6zg4cekRXz9yZa3kXx6ZNaOA6\nZFGmzfHT+vbDB4PMt6kL+WjhqztbrUkcr9lUke8VemZYHknH2cXX3aRLtn0x1/Iujk2bECOSka+4\nhrLRHFXn8vQxRiTrc3x9HL2LMSJ13gNX9lisYpyqqjXv9opsuwpVxFpkfUwVc6LMnMvK1/Hqc7zb\n9vsma+2+tl2x0jEiNHAZ+TgppOXqZD25l83fmVR3kfpo4LLnks1dJaGBK7dOGrjuNnDjzkFZ11/m\nnFcmUy/0Bm5SlmKW81zW32nj6s6a3zfudfelgeNTqDVJy9XJKi2DyFf+Ttn6cExaLhlZXUC6us9B\nZOodUyRLkay95nEPXE3ScnWy5uSkZRD5yt8pWx+OScslI6sLSFf3OYhMvWOKZCmStdc8Gri6pOTq\nZM0QSntc0fydUVEUafa7/1F6PSGIoujoJfLRn7M8PstyKX3M0sY477brVLaGMmPss44siuxrHJN3\nPy56fMlzZG4TMvXavt/z1lfknDVu7Ivm4vmuo89o4OqSkquTNScn7XG+8nd6lXFUMcYSKKDkOdLX\n9vqYqVfknEXWXvNo4GqSlquTNScnLYPIV/5OnzKOqsZYAvmVPUf62F5fM/WKZCmStdc8GriapOXq\nZM3JScsg8pW/06eMo6oxlkB+Zc+RPrbX10y9IlmKZO01jxiRjHzFNRT52HPe5zT187jXXYUiERBV\nfNS87Jhl2V6ZmrKOTZKsjy+zvSrGvsr1+n4N456fd3tFtu2Lr3Et+5imztVlx7XqOqoay7p/x5QZ\nm7zPrVjpGyxp4DKqcrKMy71Jy0Uqk/2W5bmZMudOOFHf3rdPK6ePT8wNqlresc+ac5e2rqJjSQNX\nTfOTZ96k5fDlXW/R1zBpnXU1cOPyCcsqMz8mrWf0MfP27YJzkM54mu7ZcUPh7SUtL5t1VmT88tSR\n5XhNmw+Tcvd8nOfG/c6o6nzZlwaOHLiGjcu9SctFKpP9liVnJ3Pm3PcP6DEnnnD0MZNyg5pURcYU\nmUXNyTtvsuZ7dTnbqu58wipe96J9u+AcpC/v6M2+y/q7IG0+HLn/Lmn3zRNfW5kxWPTcBfuL82U5\n3APXsHG5N2m5SGWy37Lk7OTJnJtnQm5Qk6rImCKzqDl5503WfK8uZ1vVnU9YxetO3LdjdHnfZf1d\nkDYfdMcXMr22MmOQeMwVWA+S0cB5FEX5c4PG5t6k5BSVyX7LkrOTN3Mubw15FRnXRSrImPKRWZT1\ntXkZgwBkfZ1F5k3q8qzrXcBl8Zb7AAAfA0lEQVRLrR6UrSPLuKRtb9y2fb/uKErPqxynzHysYt9l\n2V9Z68j0uyBt/w5mJz93zLbL/L7J+xgko4Fr2Njcm5Q8ojLZb1lydvJmzuWtoREVZEyRWdScIvMm\nS75Xl7OtUrdXUe5ZnWPp+zmTntuWfZfpd0Ha/o2Sf/0v3FaZMfD1GCSjgWvYuNybtFykMtlvWXJ2\n8mTOzTMhN6hJVWRMkVnUnLzzJmu+V5ezrdKyvqrKPatrLMfp8r7L+rsgbT7oac/N9NrKjMGk/cX5\nshwauIaNy71Jy0Uqk/2WJWcnc+bc8dP69sMHM+cGNamKjCkyi5qTd95kzffqcrZVWtZXVblntYzl\ngnOQnr6uN/su6++CtPmw5AWvyfTayozBpP3F+bIcYkQyqiKLJ+tz25C/U6RWX8pEtYxbVxX7y+f4\n+VxvlrFJUneMSB3HYtVzOe+2so5TVdvOq85jtKo5VOfrGTd+ddZR5Dhuw++VsufzuudHDsk3MObA\nFTgAAIDAcAUuo6QOPWuwYpkQ0YXhuEX/9VIkBDLvtsaFp/oa+zlZgykX1jBpPNL+lVYmjDJPmGnR\n/djGK3BZjwcfx/G4sTyy++Y4SmH/vvgtnLl9l2E/FpnLeV9bnnGetN60sc8SaJtmds/MsfFbvnLy\nOFU0rmWuwOXdnu+A7lFl9leReTBuu+PGu+y+4ArcRKWvwBHkW1CeYMUyIaI+wnHLhEBmVSY81de2\nkoIpk2rIG0xZJozSZxBoW0JFs6jieCjy+ueeczSLanTfTdiPdQff+uIztHpuXUfHL8s4tWxc826v\n7nmWd3+FdB5AtXgLtaBcYbclQkR9hOOWCYHMqkx4qq9tJQVTJtWQN5iyTBilzyDQtoSKZlHF8VDk\n9U/ad+PWVXfwrS8+Q6uzrCvv/GhDoPC47dU9z/Lur5DOA6gWDVxBucNuPYeI5lEqBDKrEuGp3raV\nEkw5+vjBYHDsX64LpC33FZCcd51VrsuHKEoOJI2i9MDVueVpzx2nyOvPOzZZAq2LhMnWymdodYZ1\n5Z0fdY9r3gDj2ueZpzom1TcYDOp4W3Di9rLUUXetoaKBKyh32G1DIaLjnl8mEHiREuGp3raVEky5\n6PE5a/UVkJx3nVWuq2ptCcTNu70sgdZtHO95fIZWZ1hX3vlR97jmDTCufb97qqP1xyW8o4ErKFfY\nbYkQUR/huGVCILMqE57qa1tJwZRJNeSttUwYpc8g0LaEimbRlkDcPMGvWQOt2zjeo3yGVmdZV975\nUfe45g0wrru+vPsr1OMS/tHAFZQnWLFMiKiPcNwyIZBZlQlP9bWtpGDKpBry1lomjNJnEGhbQkWz\naEsg7txzjl7lGN13E/ZjSOM9ymdo9dy6jo5flnFq2bjmDTCuu768+yvU4xL+ESOSURVhimU+Ft/G\nnxfW6kuZj8UXWW/aY+oevyo+eu8zRqQNx1yRY7HOqIOmYkTybDuvNuzfsnOo6teTdfyarGPSdsvW\nVOY1Z3lukfryPrdi+W4CTloBDdxkeTN6aOD62cDlzY4al/lXJnfOVwPnI1Ov6Dj5yOHKMwd91jRp\nWzRw1cynsg2cz0zMrOM3sY6SmZNZsjibauDG/V7NMl5p9eXN+Ay5gSMHbgKfmUrortLZciOZf5IK\n5875kjXHzcfrzpITVuVrbmNNfdfGY6AOvjMn25pbWNXv1bbsx7rQwE3gM1MJ3eUlW240869g7pwv\nY3PcRn4Z+HjdWXLCqnzNbayp79p4DNTBd+Zk0pxtg6p+r7ZlP9aFDzFM4iFTKYomZ14NBuTeZJFl\nLJvYhq/MplUrpkvlziXJ+3qiaHKO26Q6yiyvO4erjTU1oaq5VeTc5mu807a98LU2tX+z1jEqb61J\ncznrvq7s91LB36sL615YXx/m6SgauEl8Ziqhs3xlNt370IFSuXO+ZK3b5+su+tyy2lhT37XxGKhD\nWzInK1fR79WgxsADGrgJfGYqobvyZjONy/wrkzvnS9bX4+N1tzEnrOma+q6Nx0Ad2pI5WbWqfq+G\nNAY+0MBN4DNTCd2VN5tpXOZfmdy5ul+Pj9fdxpywpmvquzYeA3VoS+Zk1ar6vRrSGPhAjEhGdURw\nJG2rjm2HFCPiMyqjTPxCHeNX5zGX9/U3OU5lj8Wm5rLPeeNz20mPL7JtXzWFut+rGr8mj8sydac9\nN8s68z6myGsqux89iEqvoIYiK9fmBq5sdlSZbbfpl+Yks3tm4k8m7d8XX40a5hdlzeAbN35pkh43\nb3sl89e63MCl7Ze04z3v8iJzqM5f5L6zpkJs4PKe28rmaXatgfOZO5flGC1zTi0yfmnPzbLOvI/J\nWl+Zx1egdANHjEiF+pZJU9RcJtDRj38P84uO3H+XtPvm2sZvUTZRA/lrIUjLcFrzhFO9LB/VxjnU\nxprqlncMyNOcr+pjKC0Tru5zKqpFA1ehvmXSFJWWOaY7viANZuctrnL8EuuoadshSctwunjNk70s\nb3v+WhtrqlveMSBPc76qj6G0TLi6z6moFh9iqFBa9szRf/20xGBwLEtn9Oc6RFF65tjCE80cX5k+\ni15rSjZR0rbzjlPd41qplHFaMpX8jsDSJcmnmbTlefLXqjp209Y7GAxS529Xs6aS5D63pRwzSY9v\ny1ypsg4feWXj6ktdj6dzaq/Pfy1CA1eh1OyZlAycvkodpyj58Kws0yfDfulqnlAuKeN0ZDblBJ2y\nH7Ps31bmOpENmf/clnd5x1V9XLfmnIpK0cBVKC3rK1q7oZmCWiptnPS059aa6ZOUTVTXtkOSluF0\n9cydmfdj1v3bxlwnsiHzn9sS51aPz4VVH9dtOaeiWjRwFUrL+mrbFws3bW6cjv5rfDhOS17wmloz\nfRZlEzWQvxaCtAynS7Z9MfF4T9qPWfdvG3OdyIbMf25LGrM+nwurPq7T9k/d51RUixiRjKr4eHlb\nYkR8veZxrynLNsrUkXedWWuqOq4h1BiRsq877/bKxDVUFR/g89gIMUbEx7ktz/K6Y0R8PXehqn+X\nVBWhUTbmI89YFFl/H2NEaOAymrSzfeS9pS1v8y9yGrh2N3BZ88qKvB7fDVzRvLcsr6PI65ukDQ1c\nWtZXkdeZpYZxGWJZ1pW0fOE6k15DnvNrmxs4n1l4eesr8/iFz8myrixo4MohRsQDcqHQRiEdlyHV\n2hZpWV+zUiVvTVaR5Za0zoWvoSvHBll48I174DwgFwptFNJxGVKtbZGW9TXYua2S7VWR5ZaWATn6\nGrpybJCFB9+4AudBlqyqppTZdtm6y1yKrvEy9iJ1bjPrtorU5CNrqi511jo6llXv6yrnUOrYZMgz\nLCRlvaX2UVqtI8vrPL9WejzkzMLzIW09Vb3OkG7JCqnWNFyB86CVWVXovZCOy5BqbYvacyaryL/L\nkA/XmWODLDx4RgPnQRuzqoCQjsuQam2LunMmq8i/y5IP15Vjgyw8+EYD50Ebs6qAkI7LkGpti7pz\nJqvIv8uSD9eVY4MsPPhGjEhGVcd0jFvua9tVfNw7a61lxqDuGJEsdeStNWt9VUQaVBEf4/N1V/Ea\nxm2viseXGddx2y4zl6uaN772dZnX6TPuwtecKHtcVrGvs/A5rnk1GSNSZBueESOCcM3umdHdmzdq\n1YppHbnmsnn5T1kef9E5Zwb3r/BQJe0rIEScR8KU9/dFH9DAoRFZ8p9GdSULKkRpeWP84kNoOI+E\nqe7Mw1BwDxwakSX/aVRXsqBClJY3xtgjNJxHwlR35mEoaOBKGgwG3t8rX7jO0T+nbS9vHVXUnUuG\n/KdRbcg08zX2ZbdXx7ZH15+UUyXF+VV5xsNXrVnXE9ycaJEoio7eC9QlaeeLtGN8Thvm4rht1DG3\n6j7/jao98zAQNHBoRs5MpM5kQYWI/Cp0RO3ZefCC/ZaMBg6NyJuJ1JUsqBCRX4WuqDs7D36w35LR\nwKEReTORupIFFSLyq9AVdWfnwQ/2WzJy4DKqIounbN6Mr6ylvNvqaw5c3seXzYHLW1Pd+YR5686y\nvbIZW3Xqcg6cz6xCX8uryoHLq+x6qj5vV/14n8iBK4cYEQBHkbXUT2SjAeGhgQMgKX82H7ohab+T\njQa0H/fAAZCUP5sP3ZC038lGA9qPK3AZ5X0/vI77CdpYUxXG5RLlWV6HoPd7zmy+QtvIoe3Hq8/6\nfK1r9P6dsvs9S8ZiU3Nw3PrbeF7Ioo1jWbWqtj263i7/nuQKHIAYeW/9lDN7EUA70MABkETeW18l\n7XcyFoH2o4EDIIm8t75K2u9kLALt12gOnJm9VNIa59xvlVxP5Tlwo5rK68mqivU2mQOXZXlVfOVF\nZV1n1flUVeXAlVF3tl8VqsqaqjPjrWzOWhZV5N0V0cYcOF/aPlfStPG4rFjpHLhGGjgziyT9saRX\nSvoQDVz3Grgju2/W16+9SqtWTGvqpFNS88SSap3dM5P6XBq4bNujgatXFQ1c1nngu4GbNHdp4Gjg\nqkADl1+Tn0L9bIPbRoUuOufMwnliZJEBzc2DMnMXQL0auQfOOTdwzt0g6fYmto9qXXnBuYXzxMgi\nA5qbB2XmLoB6tToHzszWSVrXcBnBafrSeWp+VELe1KJaS2SRVaHpscwjb60Nff/fIlm2GdJ+yGLi\n65kwD6oajzxzt4i8WWeh7ve2z6cuKvO6Qx2zVjdwzrkdknZMepyZXV55Mcjs3ocOHHsLZlSWPLHl\nK5N/WZBFhj5paB6UmrsAakWMCLzbfOOthfPEyCIDmpsHZeYugHrRwMG7rbv2Fs4TI4sMaG4elJm7\nAOrVaA6cL8SIVK+qHLgs28uyvI3qjhHJu+0yURRVCWn/jioTz+JzvW3PgSsjpBiRNo5f29UdJdMC\nQceIoIdm98zEn2jbvy/+131KPhyAfGb3zOjuzRu1asW0jlxzGXML6DgaONRmLtvqaEwBGVOAF+S3\nAf3DPXCoDRlvQDXIbwP6hwbOo8Fg0JX35qvhIeMt7xhHUXT0vgmUU8fxzRzKbzAYJEd/SKnZi2Uz\ns5rKN/O13SaPM47x4hi7+WjgUJ+0LCkypoBymFtA79DAoTZkvAHVYG4B/UMDh9rMZVsdvSpAxhTg\nBfmJQP+QA1dA1VldbdSWHLim1uNz26HmwCFdkzlwWR7vK08O5XPgkKyHY1n65mwauAK60FzkRQPn\nb9ttaOCO7L5ZX7/2Kq1aMa2pE07Ut/ft08rp4zV10il69fs+oa279gZxXLYFDVx/TBqn2T0zx+bW\nSaeQx5cRDVx+5MABPbMoM+z7B/SYE0+If96/T1dvPK+54oCAzWVdkseHOnAPHNAziZlhI6aXHRc/\nBkAuZF2iTjRwBfjKogkp0yat1ra8hjbmvbVlbEaNzQwbsWrFdA3VoC5px2LavGnjsdtGi8bJQ9Zl\nX3HM5UcDB/RNhmywex86UEMhQMeQx4ca0cABPZOYGTbiwKFHtPnGW2usCOgG8vhQJxo4oGcWZYYd\nP61vP3zwaH7YxdfdpK279jZdJhAc8vhQJ2JEkFvWGIw6Y0RC+gh6kzEiaY8hQqKctseIZHk8+z2b\nqvY1knV4LMmBk2jg6ta2Bi4td6mteUw0cN1DA9cfNHD16vBYkgOHfkvLXTpy/13S7pvJYwIAdBL3\nwCFoablLuuML5DEBADqLK3DIbdyl7NG/K3PJO/Nz0/KVBrP5Hl+jrONX1TaqfG5fZRmzIuOa9zns\n9+oxTvVivNNxBQ5hS8tXilIObfKYAAAdQAOHoKXlLulpzyWPCQDQWTRwCFpa7tKSF7yGPCYAQGcR\nI4JOCCkHri2IkAD7vRqMKzIoHSPCFTgAAIDA8ClUAOih2T0zunvzRq1aMa0j11zWmqBrANnQwAFA\nz6QFYBN0DYSDt1ABoGfSArAJugbCwRU41KLqm3q5WTg/X6HLqIfXOZQWaN2CoOsuYD6hDlyBA4C+\nSQu0JugaCAYNHAD0TFoANkHXQDho4ACgZ9ICsPkAAxAOgnxRi6aCLQnURFdUcSwzP4DGlA7y5UMM\nANByZLYBWIgGDgBajMw2AEm4Bw4AWozMNgBJuAKHTuPeHgSPzLagcZ8hqsIVOABoMzLbACSggQOA\nFiOzDUASGjgAaDEy2wAkIQcOteA+EKAccuDCxBgjBTlwaD8yrAAA8IsGDpUiwwoAAP+4Bw6VIsMK\nAAD/uAKHapFhBXhRxT1U3JcFhIsrcKgWGVYAAHhHA4dKkWEFAIB/NHCoFBlWAAD4Rw4cakEWEoA+\n4tyHFOTAAQDQRmRgoko0cAAAeEYGJqrGPXAAAHhGBiaqRgMHAIBvZGCiYjRwAAD4RgYmKkYDBwCA\nZ2Rgomo0cAAAeEYGJqpGDhxqQRYSgD7i3IcUpXPguAIHAAAQGBo4AACAwNDAAQAABIZvYkAtuP8D\nAAB/uAIHAAAQGBo4AACAwNDAAQAABIYGDgAAIDA0cAAAAIGhgQMAAAgMDRwAAEBgyIEDAKAiZGCi\nKlyBAwAACAwNHAAAQGBo4AAAAAJT+z1wZhZJeq+kCyR9U9LLnXP31V0HAABAqJq4AvcSSY+R9ARJ\nfyLpygZqAAAACFYTDdwFkj7inBtI2iZpXQM1AAAABKuJBu7xku6TJOfcIUlLzIx78QAAADJqIgdu\nIOnwyJ8PO+dmkx5oZuvEFToAAIB5mmjg7pd0mqQ7zOw4SQfTHuic2yFpx6QVmtnlvooDAABouybe\nuvy0pFcOf36lpO0N1AAAABCsJq7AfUrSi81sr6R7Jb2sgRoAAACCVXsDN/z06evr3i4AAEBX8OlP\nAACAwNDAAQAABIYGDgAAIDA0cAAAAIGhgQMAAAgMDRwAAEBgaOAAAAACQwMHAAAQGBo4AACAwNDA\nAQAABIYGDgAAIDA0cAAAAIGp/cvsq2JmTZcAAACQxcA5F5VZQTQYDHwV03lmdoVz7oqm6wgJY5Yf\nY5YP45UfY5YfY5YfY1Yt3kIFAAAIDA0cAABAYGjgAAAAAkMDBwAAEBgaOAAAgMDQwAEAAASGBg4A\nACAwNHAAAACBoYHLZ0fTBQRoR9MFBGhH0wUEZkfTBQRoR9MFBGhH0wUEaEfTBXQZ38QAAAAQGK7A\nAQAABIYGDgAAIDA0cAAAAIGhgQMAAAjM0qYLaDsziyS9V9IFkr4p6eXOufuarap9zGxK0gclPV/S\nQ5J+SdLpkn5X0qHhw17vnNvZTIXtZGa3Sjpx+MfbJL1V0sclnSLp4865NzdVWxuZ2e9IunBk0aMl\nXSTpI5IODpdd7pz767praysze6mkNc653zKzH1M8Vo+S9OfOuXdxjltswZj9uKQPSIokfUHSxc65\nIwvnrnPuoobKbYUFY/ZKJZz7zeztkt4g6T8lvco59+WGyu0EGrjJXiLpMZKeIGmjpCsl/VyTBbXU\nBkknK27aVkvaKukTiifu55osrK3MbImk7zrnzh1Ztk3S5ZJukHSDmf2Ec+6fmqqxbZxzb5f0dkky\ns9WS/kTx3HyHc+4vGyytdYaN2R9LeqWkDw0Xv0/SKyTtlvQvw+PtaeIcJyl1zN4j6VWSviTpWkk/\nY2af1IK521cpY/ZELTj3m9kzJL1I0pmSniXpzyStq7XYjuEt1MkukPQR59xA0jZxwKV5rIbj5Jzb\nrfgXwuMl3dtsWa12muIrHpKONnTPknTD8Hi7XtILG6otBH8q6TfFcTbOZyV9TJLM7HGSIufc7c65\nRyR9StJ6cY5baHTMIknfcs7dOhyfmyU9WQvmLo6N2VDSnLxA0rXOucPOuVskrTKzR9VVYBfRwE32\neEn3SZJz7pCkJcO3CzHCOfce59wnJcnM3iDpq4qvxr3XzPaY2Z+b2XGNFtk+p0t6hpndbma3SFor\n6aHhLwpJ+nfFjTEWMLO1kg46525XPI6bh8fZR83spIbLa4XhP6ZukHT7cNHRc9nQ3PHFOW5o4ZgN\n//wCSTKzkyW9UdJNWjB3zey8xopuWMJxJiWf+xcefw8o/oc+CurlJM1pIOnwyJ8PO+dmmyqmzcxs\n2syukXSJpNcqTuG+RPFbND8o6ZcbK66dDii+t+ZcSW+S9H80/1gbSDrSQF0huFjxFThJ+r+S3ibp\nKZK+IekdTRXVcgvPZXPHF+e4Cczs+YqPs78evi24cO5+jH+gzrNDi8/9accfCqKBm+x+xZfLNZyg\nB8c/vJ+Gl8L/WdLDik9qX5N0lXNut3PusKS/kvTUBktsozslvds594hzbpekr0h6xsjfn6a4IcEI\nM1su6Xk69jU9H3XO7RxeufyoOM7SHD2XDc0dX5zjxjCzV0u6WtLPOueuHC5eOHcfkPRDTdXYJsO3\nnZPO/QuPv5MlfauBEjuDBm6yTyu+OVPD/29vsJY2e5Oknc65S51zByUtkfRvZjY3YV8i6ZbGqmun\nSyT9oSSZ2VmSTpJ0o5n9xPB+uFcrPv4w3xrFx9rcv95vGd4gLUkvFcdZIufc/ZKWmtmTzGxa8Zzc\nLs5xqcxsqaTfk7RueN/WnIVz99GK35JG+rn/05J+1symzOx5kr46fMseBfEp1Mk+JenFZrZX8U2Z\nL2u4nrZ6lqTnmtnoTfeXSNphZocVX537cBOFtdi7Fb/1cpfij9X/guJ/yW+VtFLSh51zX2qwvrZa\no/gTgXPeJOnaYdP7FcXjiGS/KumTkqYl/Z5z7ptmxjku3RmK79P6RzObW/ZuJcxd3naOOecOm9mi\nc79z7hEz+2dJd0n6juJPPKMEvsweAAAgMLyFCgAAEBgaOAAAgMDQwAEAAASGBg4AACAwNHAAAACB\noYEDUAkzi8zsEjPbbWYHzexbZvZXZnbmyGO+bmbPbrLOUWZ2hZn9xfDnTWbmLYfPzHaY2UXDnz9o\nZpd5Wu9rzWx2GGmx8O/+p5kdMbMdPrYFoD1o4ABU5fcl/Y/hfysV57dNSbrJzCr9DsRhJlwpzrkt\nzrkX+agnYd2/4Jz7I4+r3OKcOzthO++Q9HyP2wHQEgT5AvBumMJ+qaSfdM7tGC6+28xeqfhLry+V\n9Nbh8heb2bWSHiXpXc65dw7X8WuSflvSMknbJP2Kc+6gmf2opL9Q/PU8M5Je75y7x8yukPQESWdL\n+qqZvULSmuFXHcnM/kbSTkl/LOkqSa9SfA78R0mvk3S+pP85fOwRSf8q6SLn3PlmdqKkP1OcKr9f\n0h84564ePnag+Lsef2P4Gq6Y+7uR8XiP4q//+nEzm5V0gaQ7nXN/MHz+b0l6i6TvSfo1Sb+i+Cvp\nbpS0aRiO+gRJ10h6tqTdkt7onBv9AnEAPcIVOABVeL6k/xhp3iRJw6+/+riktSOLnyvpOZL+u6S3\nmtl5ZnaGpMslPVPSWZKeJOnnh18B9WnFX2/0WMVN1odH1vUySZdJ+nlJfy/pRZJkZidIeoGk64bb\nWSfpyZJOV5y0/3POueskvUPS+51zv7zg9bxL8bcXnK24iXubma0f+fv/Jukcxd8K8admdvyC1/3L\nkj4v6VXD7Sz0dMWp/58Y/vdHw9qeKel8M5uS9LfDv3uspI8obmoB9BQNHIAq/BfFX8uU5AFJp4z8\n+Urn3APOudskXa+4wYoUX826YPiY8xU3ai+WdJtz7m+dcw9LukLSs81sxfBxn3bO3Tz8YvuPa9jA\nSVov6cvOuW9I+pfher6r+AvIDyp+izfR8Mu5XyPpt51z33HO3aH4CuDoV05d4Zz7rqQbFF/Vy/sW\n8f8evp5/lXS7c+4fnHP/qfirwU6T9GOSljjn/sI5d8A5915JUyPfAQugZ3gLFUAVvqP0Jubxkr4+\n8ufRLwG/X9LjnHN7zexCxW8lXqW46foVxVelftLMDo48Z4niRmxuu3P+VtLVZnaK4qtmfz1cfpKk\nD0n6YUl7JK3QeKdKOl7S10aWfVvx9//O+Q8pvsI4/P7HvOfWB4b/n1XcWI6aUvy6n5zwuh87+sDh\nl4R/dvjHzzvnuP8N6CiuwAGowucknWFmzxxdOHwr8ELF93bNOW3k58dLutfMzlb8Fuz5ipszJ+l3\nJH1L0ieccyc4506Q9AOK3479t4UFDK9gfU7xVbgLFF/dk+K3SW92zp3lnHuxpLsnvJYHJQ0kPW5k\n2RM1vwkta9KXUn9L0q1zr3v42v+rpC+MPsg593nn3NLhfzRvQIdxBQ6Ad845Z2YfkHS9mb1e0i2S\nTlTchB1UfDP+nLeY2S2K7wH7acU36T9B0jVm9kJJ9ym+AvZtSZ+R9PvD6JHbJP26pAudc88ys6RS\nrpe0WdJdzrm5K31LJR0/cl/cTyr+gEUk6bCklcOf517LI2b2KUlXmtklks6UtEnx27p5HNb8t47z\n+KKkU83sxZK2S3q54vvkzii4PgCB4wocgKr8oqQPSHq/pIckfVnxlaZ1zrnvjzxuRvEVtL+X9Fbn\n3F3OuX9U3HztVPz24qmSLnfOPSDp1cP1Pqj4frmLxtTwKcXN4OgHB35H8T1x+yS9UvFbs5cq/iDB\nP0l6oeK3bUe9SfFbrfcrvs/t8rlPt+Zwo6R3mdnP5HyenHMHJW2Q9HbFY/kbkjY45w7lXReAbogG\ng0lX7gEAbWVmr1XcFL825e/XKf6Qxbr6qgJQNa7AAQAABIYGDgDCtyntq7R07FOpADqEt1ABAAAC\nwxU4AACAwNDAAQAABIYGDgAAIDA0cAAAAIGhgQMAAAgMDRwAAEBg/j8/9s2kPAGvLAAAAABJRU5E\nrkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "d=np.loadtxt(\"data/Metolius_MOD15_LAI.txt\")\n", "plt.figure(figsize=(10,10))\n", "x = np.arange( 166 )\n", "plt.plot ( x, d[:, 2], 'o' )\n", "plt.vlines ( x, d[:, 2] - d[:, 3], d[:, 2] + d[:, 3] )\n", "plt.xlabel(\"Observation time[-]\")\n", "plt.ylabel(r'LAI $m^2m^{-2}$')" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "2600177b-24b4-429e-96fb-b173c84aaa6d" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## Data assimilation as Bayesian update\n", "\n", "* Carbon pools are parameterised by a PDF\n", "* Spread of the PDF defines our belief where the true value of the pool lies\n", " * This is inherently Bayesian approach!\n", "* We start with a rough estimate of the pool sizes...\n", "* We *update* the pool sizes conditioned on observations of LAI\n", "* This new estimate is now used as **the prior** for the next step\n", "* The DALEC model is used to **propagate** the state from one time step to the next\n", " * We assume the DALEC model is wrong, so propagation implies noise inflation" ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "4f3d587f-41e3-4bc6-a751-20b002d6529a" }, "slideshow": { "slide_type": "slide" } }, "source": [ "## The particle filter\n", "\n", "* We will use the Dowd 2007 **sequential MH PF**\n", "* The model of our dynamic system is given by \n", "\n", "\\begin{split}\n", " \\mathbf{x}_{k+1} &= \\mathcal{M}(\\mathbf{x}_{k},I) + \\nu\\\\\n", " \\mathbf{y}_{k+1} &=\\mathcal{H}(\\mathbf{x}_{k+1},S) + \\eta,\n", "\\end{split}\n", "\n", "* The first equation encodes the forward propagation of the state in time using the model\n", "* The second couples the state with the observations\n", "* We will couple the foliar C pool $C_f$ with observed LAI from MODIS, noting that they are related by the SLA (assumed known)." ] }, { "cell_type": "markdown", "metadata": { "nbpresent": { "id": "a832bf6f-1524-4f84-ab5e-e53b68cab9d1" }, "slideshow": { "slide_type": "slide" } }, "source": [ "\n", "* The filter propagates a number of *particles* that determine the state using the model, and addss a random forcing\n", "* **If** an observation is available, then particles which are able to better explain the observation (i.e., where the foliar C pool scaled by SLA is close to the observed LAI) will tend to dominate, whereas particles that are far from the LAI observation will tend to be rejected.\n", "* After a few iterations, the trajectory ought to start tracking the observations.\n", "\n", "![pf trajectory](http://jgomezdans.github.io/dalec_pf/pf_trajectories.png)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true, "nbpresent": { "id": "7a19668d-5e0e-4838-bac7-4e42dd017c46" } }, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "Ctrl-E", "itemize": "Ctrl-I" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "nbpresent": { "slides": { "023a064c-30fe-487e-a964-49193a36edc5": { "id": "023a064c-30fe-487e-a964-49193a36edc5", "prev": null, "regions": { "4fc23a2e-7301-400b-b048-d79b08fd4932": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "9363186e-7292-46fa-b9b9-ab7233a0d2ee", "part": "whole" }, "id": "4fc23a2e-7301-400b-b048-d79b08fd4932" } } }, "070591dc-eb39-470b-aee1-722e83352bf2": { "id": "070591dc-eb39-470b-aee1-722e83352bf2", "prev": "3bb3eff1-8081-4347-b490-14da7e168882", "regions": { "7c94a893-83b2-4ce4-bd36-0fe221d46836": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "dfa52bce-5f45-4cd7-97b1-6a9f335a846a", "part": "whole" }, "id": "7c94a893-83b2-4ce4-bd36-0fe221d46836" } } }, "0863736b-15d4-487e-a2b7-2512d4b9f82c": { "id": "0863736b-15d4-487e-a2b7-2512d4b9f82c", "prev": "88f4b6ec-121a-4810-a6b3-c8773635398f", "regions": { "6c8dc0ae-b785-4068-8ef6-acdd4441e4ac": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "749f2688-434c-4dba-926a-0a75b7c2fa9e", "part": "whole" }, "id": "6c8dc0ae-b785-4068-8ef6-acdd4441e4ac" } } }, "2042f4ce-8198-4be0-be71-4a4308641ba3": { "id": "2042f4ce-8198-4be0-be71-4a4308641ba3", "prev": "cf0e8db1-d407-49d9-a857-7738f957a9a9", "regions": { "7c590b5b-c851-484f-826a-2a22c37dde87": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "4f3d587f-41e3-4bc6-a751-20b002d6529a", "part": "whole" }, "id": "7c590b5b-c851-484f-826a-2a22c37dde87" } } }, "23c96bc3-a828-4777-871c-6d0a991e1c45": { "id": "23c96bc3-a828-4777-871c-6d0a991e1c45", "prev": "070591dc-eb39-470b-aee1-722e83352bf2", "regions": { "2159a7a9-40ee-4027-9f43-b531fbf2d323": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "749f2688-434c-4dba-926a-0a75b7c2fa9e", "part": "whole" }, "id": "2159a7a9-40ee-4027-9f43-b531fbf2d323" } } }, "37c6c2df-049f-4b8f-92a2-10c7d98db1fd": { "id": "37c6c2df-049f-4b8f-92a2-10c7d98db1fd", "prev": "23c96bc3-a828-4777-871c-6d0a991e1c45", "regions": { "d02ff7a6-76e8-4783-88af-06d4b0bcb21f": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "e7484ffc-4770-4bad-af58-23a2206602b6", "part": "whole" }, "id": "d02ff7a6-76e8-4783-88af-06d4b0bcb21f" } } }, "3bb3eff1-8081-4347-b490-14da7e168882": { "id": "3bb3eff1-8081-4347-b490-14da7e168882", "prev": "414bfd9a-aa9a-4705-a8dc-f2f503a79dc6", "regions": { "8db3e0df-5ff0-405b-a614-a9ace8546719": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "e9923300-1d60-4b9a-98f5-441363bc8216", "part": "whole" }, "id": "8db3e0df-5ff0-405b-a614-a9ace8546719" } } }, "414bfd9a-aa9a-4705-a8dc-f2f503a79dc6": { "id": "414bfd9a-aa9a-4705-a8dc-f2f503a79dc6", "prev": null, "regions": { "d1258b1b-90d9-45db-b7de-4d28d1fe27ef": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "9363186e-7292-46fa-b9b9-ab7233a0d2ee", "part": "whole" }, "id": "d1258b1b-90d9-45db-b7de-4d28d1fe27ef" } } }, "571e27ed-a38d-49f5-8850-21e1c2469dbe": { "id": "571e27ed-a38d-49f5-8850-21e1c2469dbe", "prev": "6bbe74ca-4eab-44a3-8fe7-a314641d12db", "regions": { "798b35a9-f171-4987-9080-a56d38b825a9": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "88ebb92f-326a-4e88-9e16-5c9bf3bc2d39", "part": "whole" }, "id": "798b35a9-f171-4987-9080-a56d38b825a9" } } }, "60f50e5c-e8a7-48df-8482-6067a05f5935": { "id": "60f50e5c-e8a7-48df-8482-6067a05f5935", "prev": "c88353ff-39b6-496f-9919-9bdd0e8f280f", "regions": { "c81f12bb-b382-46cc-bfbe-338f7f69c410": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "a832bf6f-1524-4f84-ab5e-e53b68cab9d1", "part": "whole" }, "id": "c81f12bb-b382-46cc-bfbe-338f7f69c410" } } }, "6b51198c-e50e-473a-b932-a0457a130061": { "id": "6b51198c-e50e-473a-b932-a0457a130061", "prev": "7d5ab888-dd42-448b-a5f1-656a818c36f5", "regions": { "3490ebd0-9501-4590-a103-8664be2fa8fa": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "88ebb92f-326a-4e88-9e16-5c9bf3bc2d39", "part": "whole" }, "id": "3490ebd0-9501-4590-a103-8664be2fa8fa" } } }, "6bbe74ca-4eab-44a3-8fe7-a314641d12db": { "id": "6bbe74ca-4eab-44a3-8fe7-a314641d12db", "prev": "37c6c2df-049f-4b8f-92a2-10c7d98db1fd", "regions": { "390541cd-755d-40b6-8f51-da6835c80a32": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "f1459bcd-1fac-4d48-99e8-c92973cf78bf", "part": "whole" }, "id": "390541cd-755d-40b6-8f51-da6835c80a32" } } }, "7d5ab888-dd42-448b-a5f1-656a818c36f5": { "id": "7d5ab888-dd42-448b-a5f1-656a818c36f5", "prev": "bd9ad0ca-5bb0-49bb-8989-54106a549964", "regions": { "4954b36d-d76b-4168-b8a4-3a049fe23ad1": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "f1459bcd-1fac-4d48-99e8-c92973cf78bf", "part": "whole" }, "id": "4954b36d-d76b-4168-b8a4-3a049fe23ad1" } } }, "86cb0672-54cb-4470-be7d-318b947e48d2": { "id": "86cb0672-54cb-4470-be7d-318b947e48d2", "prev": "a107cbbe-7e5f-4233-993a-3136269c78aa", "regions": { "09ad547b-3eb6-4c00-87bd-13c19e30ced7": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "2600177b-24b4-429e-96fb-b173c84aaa6d", "part": "whole" }, "id": "09ad547b-3eb6-4c00-87bd-13c19e30ced7" } } }, "88f4b6ec-121a-4810-a6b3-c8773635398f": { "id": "88f4b6ec-121a-4810-a6b3-c8773635398f", "prev": "d2c0e02b-bba3-471c-b12d-7dab65a89d90", "regions": { "785a331e-999e-4428-a097-86e70fc4ef4c": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "dfa52bce-5f45-4cd7-97b1-6a9f335a846a", "part": "whole" }, "id": "785a331e-999e-4428-a097-86e70fc4ef4c" } } }, "a107cbbe-7e5f-4233-993a-3136269c78aa": { "id": "a107cbbe-7e5f-4233-993a-3136269c78aa", "prev": "571e27ed-a38d-49f5-8850-21e1c2469dbe", "regions": { "93b56c2e-6db9-49af-b933-ef8e11872ef0": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "b6c95724-8651-4e73-bddb-ea78264de730", "part": "whole" }, "id": "93b56c2e-6db9-49af-b933-ef8e11872ef0" } } }, "bd9ad0ca-5bb0-49bb-8989-54106a549964": { "id": "bd9ad0ca-5bb0-49bb-8989-54106a549964", "prev": "0863736b-15d4-487e-a2b7-2512d4b9f82c", "regions": { "cf64a6c7-c852-4587-8ec5-ff335d19576a": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "e7484ffc-4770-4bad-af58-23a2206602b6", "part": "whole" }, "id": "cf64a6c7-c852-4587-8ec5-ff335d19576a" } } }, "c04783f0-1718-431c-9a8c-66d001281e5d": { "id": "c04783f0-1718-431c-9a8c-66d001281e5d", "prev": "cf80b4fa-056e-4516-a3e4-8befed2ecf56", "regions": { "df9fda1c-3b4e-46b5-9be6-0c2d2d04b89e": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "7a19668d-5e0e-4838-bac7-4e42dd017c46", "part": "whole" }, "id": "df9fda1c-3b4e-46b5-9be6-0c2d2d04b89e" } } }, "c88353ff-39b6-496f-9919-9bdd0e8f280f": { "id": "c88353ff-39b6-496f-9919-9bdd0e8f280f", "prev": "86cb0672-54cb-4470-be7d-318b947e48d2", "regions": { "ceb550c6-7a3a-4163-a295-1e345f968d74": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "4f3d587f-41e3-4bc6-a751-20b002d6529a", "part": "whole" }, "id": "ceb550c6-7a3a-4163-a295-1e345f968d74" } } }, "cf01f1e4-0d21-4780-b769-69356c4026c8": { "id": "cf01f1e4-0d21-4780-b769-69356c4026c8", "prev": "60f50e5c-e8a7-48df-8482-6067a05f5935", "regions": { "26b8e775-3bf2-4f72-8cfc-56f956606d95": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "7a19668d-5e0e-4838-bac7-4e42dd017c46", "part": "whole" }, "id": "26b8e775-3bf2-4f72-8cfc-56f956606d95" } } }, "cf0e8db1-d407-49d9-a857-7738f957a9a9": { "id": "cf0e8db1-d407-49d9-a857-7738f957a9a9", "prev": "d9b3bb7a-72f8-4370-9c67-50002f71f6c3", "regions": { "47c4e162-761b-4a97-b42e-10757a354aed": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "2600177b-24b4-429e-96fb-b173c84aaa6d", "part": "whole" }, "id": "47c4e162-761b-4a97-b42e-10757a354aed" } } }, "cf80b4fa-056e-4516-a3e4-8befed2ecf56": { "id": "cf80b4fa-056e-4516-a3e4-8befed2ecf56", "prev": "2042f4ce-8198-4be0-be71-4a4308641ba3", "regions": { "8176f27d-d6d6-4c79-accf-b6fc31dacb56": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "a832bf6f-1524-4f84-ab5e-e53b68cab9d1", "part": "whole" }, "id": "8176f27d-d6d6-4c79-accf-b6fc31dacb56" } } }, "d2c0e02b-bba3-471c-b12d-7dab65a89d90": { "id": "d2c0e02b-bba3-471c-b12d-7dab65a89d90", "prev": "023a064c-30fe-487e-a964-49193a36edc5", "regions": { "31aef8b8-6889-4c0e-9971-1f8c822fd133": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "e9923300-1d60-4b9a-98f5-441363bc8216", "part": "whole" }, "id": "31aef8b8-6889-4c0e-9971-1f8c822fd133" } } }, "d9b3bb7a-72f8-4370-9c67-50002f71f6c3": { "id": "d9b3bb7a-72f8-4370-9c67-50002f71f6c3", "prev": "6b51198c-e50e-473a-b932-a0457a130061", "regions": { "97ae3cd1-a726-467b-a3ff-b3a844be8a65": { "attrs": { "height": 0.8, "width": 0.8, "x": 0.1, "y": 0.1 }, "content": { "cell": "b6c95724-8651-4e73-bddb-ea78264de730", "part": "whole" }, "id": "97ae3cd1-a726-467b-a3ff-b3a844be8a65" } } } }, "themes": {} } }, "nbformat": 4, "nbformat_minor": 1 }