{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Performing NCGR at a gridpoint\n", "This notebook walks you through how to use the `NCGR` package to perform non-homogenous Gaussian regression (NCGR) on an ice-free date or freeze-up date ensemble forecast at a single location using the `NCGR.ncgr_gridpoint` module. This option is perhaps preferred if you don't want to rely on the requirements involved in the `ncgr_fullfield()` module, if memory limitations prevent you from using `ncgr_fullfield()`, or if you would like to modify or have more control over different steps of the NCGR algorithm. \n", "\n", "Outline\n", "-----------\n", "* The first section of this noetbook, **Main Recipe**, provides the relevant function calls to perform NCGR and produce relevant forecast variables.\n", "\n", "* The second section of this notebook, **Working Script**, offers a copy and pasteable section of code that loads some sample hindcasts and observations at a single gridpoint and performs NCGR on a forecast for a specified year.\n", "\n", "* The third section of the notebook, **Detailed explanation & plotting**, goes through step-by-step and provides insights on each piece of code, the data, and creates plots showing the results. \n", "\n", "*Hint: You can hold \"shift+tab\" inside the parenthesis of any function in the notebook to see its help/docstring.*\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Main Recipe\n", "------------------\n", "```python \n", "from NCGR import ncgr\n", "\n", "# instatiate the ncgr_gridpoint and fcst_vs_clim modules\n", "ngp = ncgr.ncgr_gridpoint(a, b)\n", "fvc = ncgr.fcst_vs_clim(a,b) \n", "\n", "# build model\n", "predictors_tau, predictors_t, coeffs0 = ngp.build_model(X_t, X_tau, Y_tau, tau, t)\n", "\n", "# optimize/minimize the CRPS to estimate regression coefficients\n", "coeffs = ngp.optimizer(Y_tau, predictors_tau, coeffs0)\n", "\n", "# Compute calibrated forecast distribution parameters for year t\n", "mu_t, sigma_t = ngp.forecast_mode(predictors_t, coeffs)\n", "\n", "# get relevant forecast information based on the forecast distribution\n", "result = fvc.event_probs(mu_t, sigma_t, Y_clim)\n", "fcst_probs, clim_terc, clim_params = result.probs, result.terciles, result.params\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bare-bones code" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Calibrated forecast distribution parameters (mu,sigma): [360.06167043316805, 15.537755970468954]\n", "Probilities for early, normal, late freeze-up: [0.08214745 0.24499739 0.67285516]\n", "Climatology terciles: [338.45221123 353.10369408]\n", "Climatology distribution parameters (mu,sigma): [345.77795265 17.00784101]\n" ] } ], "source": [ "from NCGR import ncgr, sitdates\n", "import numpy as np\n", "\n", "# load data\n", "X_all = np.load('Data/fud_X.npy') \n", "Y_all = np.load('Data/fud_Y.npy') \n", "\n", "# define time variables\n", "event='fud' \n", "years_all = np.arange(1979, 2017+1)\n", "im = 12 # initialization month\n", "\n", "# instantiate the sitdates class\n", "si_time = sitdates.sitdates(event=event)\n", "a = si_time.get_min(im) # minimum date possible\n", "b = si_time.get_max(im) # maximum date possible\n", "\n", "\n", "# instatiate the ncgr_gridpoint module\n", "ngp = ncgr.ncgr_gridpoint(a, b)\n", "# instantiate the fcst_vs_clim module\n", "fvc = ncgr.fcst_vs_clim(a,b) \n", "\n", "\n", "# the forecast year to calibrate\n", "t=2016 \n", "\n", "# set up relevant time variables for calibrating the forecast\n", "# using data only from previous\n", "t_idx = np.where(years_all==t)[0][0]\n", "tau = years_all[years_all" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAGDCAYAAAA8mveiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd5xU5fX48c+Z2dnCUpayIG1ZqhCwIVEMCgJRUVGxYIklRo0p+ovEb0yiKdYUk9iixljQ2JKooEYRLMkKCoK6oKBU6bC0BXZZtk47vz/m7rp1dmZ379bzfr3m5cxzy3PGxDl7n/vc54iqYowxxsTK09IBGGOMaVsscRhjjImLJQ5jjDFxscRhjDEmLpY4jDHGxMUShzHGmLhY4jAdmogsFJHrGnhshogUioi3qeMypjWzxGFMjERkq4h8u/yzqm5X1c6qGmrGGHqIyGsiUiQi20TkO/XsP1ZEPnAS3F4RuanStrtF5AsRCYrIHdWO6ysib4jILhFREcl05QuZNskSh2kXRCQhlrZ24FHAD/QBLgceE5HRte0oIr2At4HHgZ7AMODdSrtsBH4OvFXL4WHn2AubLHLTbljiMK2CiAwUkVdFJFdEDojIIyLiEZFfO39Z7xOR50Skm7N/pvOX8LUish3Iqq3N2Xe8iHwkIvkislJETq0jhqEikuX0v19EXhSRNGfb80AG8Kbz1/vPK/WX4OzTz/kr/aCIbBSR71c69x0i8rLzHQ6LyGoRGRfnv6NUIj/kv1HVQlVdDLwBXFnHITcD76jqi6papqqHVXVt+UZVfVZVFwCHqx+oqntV9W/Ap/HEaDoGSxymxTn3COYB24BMoD/wb+Bq5zUZGAJ0Bh6pdvgkYBRwRm1tItKfyF/U9wA9gJ8Bc0UkvbZQgD8A/ZzjBwJ3AKjqlcB24BxneOpPtRz/L2Cnc/xFwO9FZGql7ec63yuNyA9+xXcRkXlOYqvtNc/ZbQQQUtUNlc65Eqj1igMYDxx0kuY+EXlTRDLq2NeYmFniMK3BCUR+bG9R1SJVLXX+mr4cuF9VN6tqIXArcGm1Iag7nGNK6mi7ApivqvNVNayq7wHZwFnVg1DVjar6nvPXeS5wP5EkVC8RGQicDPzCif9z4CmqXg0sduIIAc8Dx1Tqe7qqptXxmu7s1hk4VK3rQ0CXOsIaAHwXuInI1dIWIsnNmEZpj2PApu0ZCGxT1WC19n5ErkLKbSPy/9k+ldp21HK+ym2DgJkick6lNh/wfvWDRKQ38FfgFCI/xh4gL8bv0A84qKqVh322AZWHo/ZUel8MJItIQi3fuy6FQNdqbV2pZajJUQK8pqqfAojIncB+EemmqtUTkDExsysO0xrsADJquZm9i8gPf7kMIAjsrdRW2/LOldt2AM9X+ws+VVX/WMtxf3COPVpVuxK5WpF6+qocaw8RqfzXfwaQE+WYCiKywLl3UttrgbPbBiBBRIZXOvQYYHUdp11VLeby91LLvsbEzBKHaQ0+AXYDfxSRVBFJFpEJRIZVfioig0WkM/B74KU4/kIHeAE4R0TOEBGvc+5TRWRALft2IfJXfb5zb+SWatv3ErnXUoOq7gA+Av7g9HE0cC3wYixBquqZzr2T2l5nOvsUAa8Cdzn/niYA5xEZ9qrNM8D5InKsiPiA3xAZLssHEBGfiCQT+R1IcOKueCbF2ZbkfExyPhtjicO0PGfM/xwi00W3E7nBfAnwNJEfxQ+IjM+XAv8vznPvIPLjehuQS+QK5BZq///+ncBYIvcN3iLyI13ZH4BfOzesf1bL8ZcRubm/C3gNuN25p9KUfgykAPuIJNYfqepqABE5RUQKy3dU1Swi3/stZ/9hQOXnPp4kMpx1GfAr533lezIlRBIpwDrnszGIFXIyxhgTD7viMMYYExfXE4czrvxZ+Vx0ifidiGwQkbUi8pNK7X91HpxaJSJj3Y7NGGNM/JpjOu5NwFq+nkZ4NZHplyNVNexMgQQ4ExjuvE4EHnP+aYwxphVx9YrDmblyNpEHocr9CLhLVcMAqrrPaT8PeE4jlgFpItLXzfiMMcbEz+2hqgeJLKIWrtQ2FLhERLKduevlc9L7U/XBrZ1OmzHGmFbEtaEqEZkO7FPV5dUWlUsCSlV1nIhcQGTK5SnU/lBSjSlfInI9cD1Aamrq8SNHjmzy2I0xpj1bvnz5flWtbb22mLh5j2MCcK6InAUkA11F5AUiVxJznX1eI/KQEk77wErHDyAyH74KVX0CeAJg3Lhxmp2d7U70xhjTTonItvr3qptrQ1WqequqDlDVTOBSIEtVrwBeB6Y4u00isowCRFYLvcqZXTUeOKSqu92KzxhjTMO0xCKHfwReFJGfEnkqtbxs53wiK5ZuJLIA3PdaIDZjjDH1aJbEoaoLgYXO+3wiM62q76PADc0RjzHGmIZrd8uql5SUsHnzZkKhZisD3Sp4vV6GDBlCSkpKS4dijGnn2l3i2Lx5M7169SI9PR2Pp2OsqBIOh8nNzWXz5s2MHl1XMThjjGka7e6XNRQKdaikAeDxeEhPT+9wV1nGmJbRLn9dO1LSKNcRv7MxpmXYr40xxpi4WOIwxhgTF0scLrn44ovp0aMHw4cPr9J+zz33MHz4cIYNG8bdd99d0X7XXXcxbNgwhg8fzjnnnENxcXHFtrlz5zJ48GAyMjK47bbbmu07GGNMbTp84nhrdhaXDf0JZ3S+isuG/oS3Zmc1yXmvueYa3nzzzSpt2dnZPPvss6xYsYK1a9eyYMECvvzyS7Zs2cLjjz/OypUr+eqrrwiFQsyePRuAYDDIrFmzmD9/Phs2bGDu3LmsWLGiSWI0xpiG6NCJ463ZWTx2y4sc3JMPCgf35PPYLS82SfKYNm0avXr1qtL2xRdfMHbsWLp06YLP52PChAm89NJLQGQ2WFFREYFAgJKSEgYMGADAokWLyMzMZNSoUSQnJ3PhhRcyZ86cRsdnjDEN1aETxwu/f51AWaBKW6AswAu/f92V/o499lg+/vhj9u7dy+HDh3nvvffYsWMHgwcP5sYbbyQzM5PevXvTtWtXzj//fAB27NhBv379Ks4xcOBAcnJyXInPGGNi0aETx8G9+XG1N9Zxxx3HrFmzmDx5MpMnT2b06NEkJCSQm5vLvHnz2LhxI3v27KG4uJjHHnsMgMhKLFWJ1LYCvTHGNI8OnTh69EmLq70pzJo1izVr1pCdnU337t0ZPnw48+bNY9CgQfTr14+kpCRmzJjBRx99BEBGRga7dn29unz1KxBjjGluHTpxXHHbDHxJviptviQfV9w2w7U+y4eZvvrqK9566y2uueYaMjMzWb58OYcPHyYcDpOVlUV5gaqJEyeyZcsW1q1bR2lpKXPnzuXCCy90LT5jjKlPu1urKh5nXxspC/LC71/n4N58evRJ44rbZlS0N8Y555zDsmXLyMvLo0+fPtx6663MmjWL8847j7y8PBISEnjooYdIT09n8uTJnHPOORx99NEkJCQwZswYbr75ZgB8Ph8PPPAA06ZNIxQKcfnll3P88cc3Oj5jjGkoqW0Mva2orQLgqlWrOProo1soopbVkb+7MSZ2IrJcVcc19PgOPVRljDEmfpY4jDHGxMUShzHGmLhY4jDGGBMXSxzGGGPiYonDGGNMXCxxGGOMiYslDmOMMXGxxOGSpizk1L9/f0aMGMHIkSMZM2ZMs30HY4ypTYdPHPNey+bSs+/n9PF3cenZ9zPvtez6D4pBUxVyKrdo0SLWrVvHl19+2STxGWNMQ3XoxDHvtWwee+BdDh4oBODggUIee+DdJkkeTVXIyRhjWpsOnTheeOoDAv5glbaAP8gLT33gSn8NKeRUburUqYwePZr77rvPldiMMSZWHTpxlF9pxNreWA0p5ASwZMkS1qxZw7vvvssTTzzB22+/7Up8xhgTiw6dOHr07BxXe1OIt5ATQGZmJhC5ST59+nSWLl3qWnzGGFOfDp04rrhuIr7EqiVJfIkJXHHdRNf6jLeQU0FBAfn5+RXvs7KybOl0Y0yL6tCFnKafH1mO/oWnPuDggUJ69OzMFddNrGhvjKYq5JSTk8OMGZGKhKFQiIsuusgqABpjWpQVcmpHOvJ3N8bEzgo5GWOMaVaWOIwxxsTFEocxxpi4WOIwxhgTF0scxhhj4mKJwxhjTFxcTxwi4hWRz0RkXrX2h0WksNLnJBF5SUQ2isjHIpLpdmzGGGPi1xxXHDcBays3iMg4IK3aftcCeao6DHgAuLcZYjPGGBMnVxOHiAwAzgaeqtTmBf4M/Lza7ucBzzrv5wBTRUTcjM9N8RZyqqsdYO7cuQwePJiMjAxuu+22ZonfGGPq4vYVx4NEEkS4UtuNwBuqurvavv2BHQCqGgQOAT2rn1BErheRbBHJzs3NbXSAr2at5OyfPM6JV93P2T95nFezVjb6nBBfIae62gGCwSCzZs1i/vz5bNiwgblz57JixYomidEYYxrCtcQhItOBfaq6vFJbP2Am8HBth9TSVmM9FFV9QlXHqeq49PT0RsX4atZKHnhxIfvziwDYn1/EAy8ubJLkEU8hp2gFnhYtWkRmZiajRo0iOTmZCy+8kDlz5jQ6PmOMaSg3rzgmAOeKyFbg38AUYDUwDNjotHcSkY3O/juBgQAikgB0Aw66GB+zX1+GPxCq0uYPhJj9+jJX+qurkFNd7QA7duygX79+FecYOHBgxQq7xhjTElxbHVdVbwVuBRCRU4Gfqer0yvuISKFzMxzgDeC7wFLgIiBLXV6BsfxKI9b2xqpcyKlTp04VhZzqageo7V9BG771Y4xpB1rTcxyzgZ7OFcjNwC/d7rBXWmpc7U2htkJO0dozMjLYtWtXxfHVr0CMMaa5NUviUNWF1a82nPbOld6XqupMVR2mqieo6ma347p2xngSfd4qbYk+L9fOGO9an7UVcorWPnHiRLZs2cK6desoLS1l7ty5Vo/DGNOiOnQhpwumHANE7nXszy+iV1oq184YX9HeGPEUcgLqbPf5fDzwwANMmzaNUCjE5ZdfzvHHH9/o+IwxpqGskFM70pG/uzEmdlbIyRhjTLOyxGGMMSYuljiMMcbExRKHMcaYuFjiMMYYExdLHMYYY+JiicMYY0xcLHEYY4yJiyUOlzRlIaf+/fszYsQIRo4cyZgxY5olfmOMqUuHXnIE4KVPV/Ho+8vYX1hEr86p3DB5PJd8s/FPX19zzTXcdNNNXH311RVtlQs2JScnM2nSJM4//3xKS0trba+cJBYtWkTfvn0bHZcxxjRWh77ieOnTVfxhwUJyC4tQILewiD8sWMhLn65q9LmbqpCTMca0Nh06cTz6/jLKglULOZUFQzz6fusp5FRu6tSpjB49mvvuu8+V2IwxJlYdeqhqf2EdhZzqaG+shhRyAliyZAmZmZnk5OQwZcoURo8ezbRp01yJ0Rhj6tOhrzh6da6jkFMd7U0h3kJOAJmZmUDkJvn06dNZunSpa/EZY0x9OnTiuGHyeJISqhZySkrwcsPk1lPIqaCggPz8/Ir3WVlZtnS6MaZFdeihqvLZU27MqmqqQk45OTnMmDEDgFAoxEUXXWQVAI0xLcoKObUjHfm7G2NiZ4WcjDHGNCtLHMYYY+JiicMYY0xcLHEYY4yJiyUOY4wxcbHEYYwxJi6WOIwxxsTFEocxxpi4WOJwwaZNmzjxxBMZMmQIw4YN45577qnYNnfuXAYPHkxGRga33XZbve31bTPGmGanqm32dfzxx2t1K1eurNEWzQtrPtNx/3xEB82+V8f98xF9Yc1ncR1fm61bt+rixYtVVTUvL08HDRqky5cv10AgoAMGDNA1a9ZoSUmJjhgxImq7qkbdVl28390Y0zEB2dqI394OvVbVi2s/585P/kdZKFKTY19JEXd+8j8ALh91bIPPO2jQIAYNGgRAWloaw4YNY/v27eTl5ZGZmcmoUaMAuPDCC5kzZ06d7WPHjmXRokV1bjPGmJbQoYeqHvx8SUXSKFcWCvHg50uarI/169ezevVqJk2axI4dO+jXr1/FtoEDB5KTk1NnOxB1mzHGtIQOnThyS2ov2FRXe7wOHTrEBRdcwL333kv37t3RWhaUFJE624Go24wxpiXElDhEpLuIjBaRISLSbpJNekrtBZvqao9HWVkZ06dPZ+bMmVx11VUAZGRksGvXrop9yq8m6mqPdowxxrSUOpOAiHQTkdtE5AtgGfA48DKwTUReEZHJzRWkW2YdO4Ekb7VCTl4vs46d0KjzhsNhLrvsMkaMGMEdd9xR0T5x4kS2bNnCunXrKC0tZe7cuVx44YV1tkc7xhhjWkq0m+NzgOeAU1Q1v/IGETkeuFJEhqjqbDcDdFP5DfAHP19CbkkR6SmpzDp2QqNujAP897//5bXXXmP48OGMHDkSgLvvvpuZM2fywAMPMG3aNEKhEJdffjnHH388QJ3tPp+vzm3GGNMS6izkJCJeVQ3VurGVsEJOVXXk726MiZ2bhZyWi8hJDT2xMcaY9ila4vgB8JCIPCki3RvagYh4ReQzEZnnfH5RRNaLyJci8rSI+Jx2EZG/ishGEVklIvaggjHGtEJ1Jg5V/Rg4EVgBZIvII84P+19F5K9x9HETsLbS5xeBkcBRQApwndN+JjDceV0PPBZHH8YYY5pJfVNrewDfBHKB5dVe9RKRAcDZwFPlbao6v9Jj758AA5xN5wHPOZuWAWki0jeeL2OMMcZ9dc6qEpEfArcAfwau1bruokf3IPBzoEst5/cBVxK5IgHoD+yotMtOp213A/o1xhjjkmhXHKcAJ6nq3xuSNERkOrBPVeu6Ovkb8IGqflh+SC371OhXRK4XkWwRyc7NzY03LGOMMY1U5xWHql4uIoki8j1gNJEf8TXAP1W1LIZzTwDOFZGzgGSgq4i8oKpXiMjtQDqRG/DldgIDK30eAOyiGlV9AngCItNxY4jDGGNME4r25Pg3iCSKU4HtRH7YTwVWO9uiUtVbVXWAqmYClwJZTtK4DjgDuExVw5UOeQO4ypldNR44pKo2TGWMMa1MtKGqh4Efqep3VfWvqvqQqn4X+CHwaCP6/DvQB1gqIp+LyG+d9vnAZmAj8CTw40b00aKaupBT//79GTFiBCNHjmTMmDHN9j2MMaZWdRXqANZF2ba2MUVAmurVFIWcXt/xkc744E495b8/0xkf3Kmv7/goruNr05SFnFRV+/Xrp7t27aq3XyvkZIyJBS4WcvKISJJWu58hIslEX+OqzfjPzqU8/NUb+MNBAA74D/PwV28AcN6Ahj8035SFnIwxprWJNlT1HDBXRDLLG5z3LwPPuxlUc/nHlvcqkkY5fzjIP7a812R9NLaQU7mpU6cyevRo7rvvviaLzRhjGiLarKp7RORG4AMR6URkumwh8BdVfbi5AnTTAf/huNrj1RSFnACWLFlCZmYmOTk5TJkyhdGjRzNt2rQmidEYY+IV9clxVX1EVTOAwUCmqg5qL0kDoGdijecSo7bHo6kKOQFkZmYCkZvk06dPZ+nSpY2OzxhjGipq4hCRI0XkPuDfwL9F5C8iMqJ5QnPf1YNPI9FT9aIr0ZPA1YNPa9R5m7KQU0FBAfn5+RXvs7KybOl0Y0yLirbkyEnAq0QetnuCyFDVccBCEblAI+tJtWnlN8D/seU9DvgP0zOxC1cPPq1RN8ahaQs55eTkMGPGDABCoRAXXXSRVQA0xrSoaIWcFgD3qurCau2TgF+q6pnuhxedFXKqqiN/d2NM7Nws5DS0etIAUNVFwJCGdmiMMaZti5Y4ok0tKmrqQIwxxrQN0R7kG1hHwSYhsty5McaYDiha4rglyrbsKNuMMca0Y9EeAHy2OQMxxhjTNkRbVv0JEal1KVYRSRWRa0TkcvdCM8YY0xpFG6r6G/BbETkK+JJI3fFkYDjQFXgaeNH1CI0xxrQq0YaqPgcuFpHOwDigL1BCZEn19c0UnzHGmFYm6pIjAKpaqKoLVfVfqvq6JY36NaSQ08UXX0yPHj0YPnx4jfNFK/JkjDHNrd7E0d6tODCfh9Zfye9WT+eh9Vey4sD8Rp8zISGB+++/n82bN5Odnc1TTz3FihUrCAaDzJo1i/nz57Nhwwbmzp3LihUrALjmmmt48803a5wr2jHGGNMSOnTiWHFgPu/ufZLCYB4AhcE83t37ZKOTx6BBg5gwYQJQtZDTokWLKgo2JScnVxRsApg2bRq9evWqca5oxxhjTEuIOXGISKqbgbSED/f/i5AGqrSFNMCH+//VZH3EUsgpmoYcY4wxbqo3cYjIt0RkDbDW+XyMiPzN9ciaQfmVRqzt8Yq1kFM0DTnGGGPcFMsVxwPAGcABAFVdCUx0M6jm0jmhe1zt8YinkFM0DTnGGGPcFNNQlaruqNYUciGWZndKr8vwiq9Km1d8nNLrskadN95CTtE05BhjjHFTLIljh4h8C1ARSRSRn+EMW7V1Y3uexel9vl9xhdE5oTun9/k+Y3ue1ajzlhdy+vDDDxk5ciQjR47klVdewefzVRRsGj58OOeff35FwaZzzjmHk08+mS1bttCnTx8efPBBgKjHGGNMS6izkFPFDiK9gIeAbxNZGfdd4CeqetD98KKzQk5VdeTvboyJXWMLOUVbcqTckapaZU0qEZkALGlop8YYY9quWIaqHo6xzRhjTAdQ5xWHiJwEfAtIF5GbK23qCnjdDswYY1pSKBRm4YJVvPufz/B4hDNmjGXiGWPweDr0c9NA9KGqRKCzs0+XSu0FwEVuBmWMMS1JVbn75n/x+cebKC2JPCS8duV2li5cx633XtzC0bW8aKvjLgIWicg/VHVbM8ZkjDEtavVn26okDYDSkgDLFq5jw+ocRozu2NWzY7k5XiwifwZGE6nHAYCqTnEtKmOMaUGfLdtEaWmgRnswGGLlJ5s7fOKIZbDuRWAdMBi4E9gKfOpiTMYY06K6pqWSmFjz72qfz0uXbp1aIKLWJZbE0VNVZwMBVV2kqtcA412OyxhjWsykaUfh8dSyJpwIp5w2uvkDamViSRzl12u7ReRsETkOGOBiTG1eUxdy6t+/PyNGjGDkyJGMGVNrGXhjTBNK65HK7Q9dTpeuKXRKTSIlNZGu3Ttxz6NXkdoluf4TtHeqGvUFTAe6AWOA94HlwLn1Hdccr+OPP16rW7lyZY22aPILntOtOUfrph19dWvO0Zpf8Fxcx9dm69atunjxYlVVzcvL00GDBuny5cs1EAjogAEDdM2aNVpSUqIjRozQ5cuXq6rqggULdPHixTps2LAa5+vXr5/u2rWr3n7j/e7GmOgC/qB+uWKrrv5smwaDoZYOp8kA2dqI3956b46r6jzn7SFgsjvpq2UcOvw8Bw/9BqUMgFB4HwcP/QaAbl2ubPB5Bw0axKBBg4CqhZzy8vIqijIBFUWZxo4dy7Rp01i/3qryGtOaJPi8jD5uUEuH0epEHaoSkckiMldEVjuvOSJyajPF5rr8w3+pSBrllDLyD/+lyfpobCGnclOnTmX06NHcd999TRabMcY0RLQnx88GHgHucl4CjAWeFpEbVbXxxblbWCicG1d7vJqikBPAkiVLyMzMJCcnhylTpjB69GimTZvWJDEaY0y8ol1x3ALMUNVnVHWlqn6uqk8DM4BfNE947vJ60uNqj0dTFXICyMzMBCI3yadPn87SpUsbHZ8xxjRUtMRxhEaq/VWhqquAPrF2ICJeEflMROY5nweLyMci8pWIvCQiiU57kvN5o7M9M76vEr+0Lj9DSKoaL0mkdflZo87blIWcCgoKyM/Pr3iflZVlS6cbY1pUtMRR1MBt1d1E1cJP9wIPqOpwIA+41mm/FshT1WFEytXeG0cfDdKty5X06HY3Xk9vQPB6etOj292NujEOTVvIKScnhxNPPJEjjzySsWPHcsYZZ1gFQGNMi6qzkJOI5AMf1LYJOFlV6y3MLSIDgGeB3wE3A+cAuUSuZoLOCrx3qOoZIvKO836piCQAe4B0rStArJBTdR35uxtjYudmIafzomyLddrRg8DP+Xp13Z5AvqoGnc87gfJFX/oDOwCcpHLI2X9/5ROKyPXA9RC5Z2CMMaZ51bc6boOJyHRgn6ourzSFt7YpRBrDtspxPQE8AZErjsbEaIwxJn6xrI7bUBOAc0XkLCKr6nYlcgWSJiIJzlXHAKB8mtFOYCCw0xmq6ga0eF1zY4wxVblWykpVb1XVAaqaCVwKZGmkdvn7fF0I6rvAf5z3bzifcbZnRbu/EU04HG5w3G1VR/zOxpiWEXPiEJGuItKl/j3r9QvgZhHZSOQexmynfTbQ02m/GfhlQ07u9XrJzc3tUD+k4XCY3NxcvF6r6GuMcV+ds6oqdhAZBzxD5Aa3APnANaq63P3woqttVlVJSQmbN28mFAq1UFQtw+v1MmTIEFJSUlo6FGNMK+fmrKpyTwM/VtUPnQ5PJpJIWuW8z5SUFEaPtvXyjTHGLbEMVR0uTxoAqroYOOxeSMYYY1qzWK44PhGRx4F/EZkeewmwUETGAqjqChfjM8YY08rEkjiOdf55e7X2bxFJJFOaNCJjjDGtWiyFnNpV8SZjjDGNU2/iEJHf1tauqnc1fTjGGGNau1iGqiqvhJtMpAb52jr2NcYY087FMlRVpVapiPyFyFPexhhjOqCGLDnSCRjS1IEYY4xpG2K5x/EFX69S6wXSidQgN8YY0wHFco9jeqX3QWBvpXoaxhhjOph6h6pUdVv5CzjbkoYxxnRs8d7j+KErURhjjGkz4k0ctVXpM8YY04HEmzjOcSUKY4wxbUa9iUNE+ojIbBFZoKo7ReQbInJtcwRnjDGm9YnliuMfwDtAP+fzBmCWWwEZY4xp3WJJHL1U9WUgDODMqupY5fWMMcZUiCVxFIlIT5yHAEVkPHDI1aiMMca0WrE8APh/RNamGioiS4g8OT7T1aiMMca0WrEscrhcRCYBRxKZjrteVQOuR2aMMaZVimVW1SbgOlVdrapfqmpAROY1Q2zGGGNaoVjucQSAySLyjIgkOm39XYzJGGNMKxZL4ihW1UuIFG/6UEQG8fVqucYYYzqYWG6OC4Cq/klElhN5pqOHq1EZY4xptWJJHBU1x1X1fyJyOnC1axEZY4yJmYb2gSQgnub7ez6WoTDk+YkAACAASURBVKp5InKFiPy2UtvbbgVkjDGmfhr4gnDuNDR3CrrvFMIHLkFDu5ql71gSx9+Ak4DLnM+HgUddi8gYY0xUGj6IHrwKQpsBPxCAwEr0wHdQdX9hj1gSx4mqegNQCqCqeUBi9EOMMca4RYtfhRo19cKgh8C/2PX+Y5qOKyJevl5yJB1n3SpjjDEtILQNKKvZriEI7Xa9+1gSx1+B14A+IvI7YDHwe1ejMsaYdu7woWLefvkT5s7+gC3r4/uxl8RxIJ1q2wK+o5omwChiWXLkRWca7tRIVMxQ1bWuR2aMMe3UZx9t5K4fPQtAMBji+YfeZeqMsdx45/mIxFBoNflMKHwUQruI3OMASIbEExDfaNfiLhdrBcBeRB4EfATYLyKDXYzJGGPaLb8/yD03Pk9piZ/SEj/BQIiy0gBZ//mMTxeui+kcIolIz1eg0+XgOQK8GdD5BqR788xbqveKQ0RuB8YRWeTwGcAHvABMcDc0Y4xpf774ZDO1Lb5RWuLnvdeWc8LkUTGdRzzdkK63QtdbmzjC+sVyxXE+cC5QBKCqu4AubgZljDHtlYbrXrEpFGob845ieXLcr6oqIuWzqlJdjskYY9qto04YQriW5JGcksjU88bGdA5VJVD6JmWFT6HhfHzJp5PU+cd4vM3z9HgsVxwvi8jjQJqIfB/4L/Cku2EZY0z7lJTs4xf3X0ZSsg9fYgJIJGmcMHkkJ337GzGdo7TgDxTn/4xQ4DPCoS2UFT3D4dxphMPNU5xVVOtf6FZETgNOJzKr6h1Vfc/twGIxbtw4zc7ObukwjDEmbgf2FrDorc8pLChh3MQjGXXcoJhmVIVD+ynYeyJfz6Yql0xyl5+S3OXH9Z5DRJar6riGRV7PUJXz4N87qvptIK5kISLJwAdAktPPHFW9XUSmAn8mcrVTCFytqhtFJAl4DjgeOABcoqpb4/w+xhjTJvTs05ULrpkY93GhwBcgSaDVE0cpwbJFEEPiaKyoQ1UaWfSkWES6NeDcZcAUVT0GOBaYJiLjgceAy1X1WOCfwK+d/a8F8lR1GPAAcG8D+jTGmHbN4+0DVF9uBMCDeAc2Swyx3BwvBb4QkfdwZlYBqOpPoh2kkTGwQuejz3mp8+rqtHcDypdzPA+4w3k/B3hERERjGUszxpgOwpMwCo93MOHgBqomkESSO1/TLDHEkjjecl5xc4a6lgPDgEdV9WMRuQ6YLyIlQAEw3tm9P7ADQFWDInII6Ansr3bO64HrATIyMhoSljHGtFkiQueeL1B08IeEAisjtThIJCXtXry+2G6uN1adiUNEMlR1u6o+29CTO0Ndx4pIGvCaiIwBfgqc5SSRW4D7getwKg1WP0Ut53wCeAIiN8cbGpsxxrRVHm86XdLnEg7tRsMFeBKGIhLLdUAT9R9l2+vlb0RkbmM6UdV8YCFwJnCMqn7sbHoJ+Jbzficw0Okvgcgw1sHG9GuMMe2Zx9sXr+/IZk0aED1xVL4CGBLviUUk3bnSQERSgG8Da4FuIjLC2e00pw3gDeC7zvuLgCy7v2GMMa1PtDSldbyPVV/gWec+hwd4WVXnOQ8RzhWRMJAHlN/NmQ08LyIbiVxpXNqAPo0xxrgsWuI4RkQKiFx5pDjvcT6rqnat+1BQ1VXAcbW0v0akvkf19lJgZqyBG2OMaRl1Jg5V9TZnIMYYY9qGWOtxGGOMMYAlDmOMMXGyxGGMMSYuzTv51xhjTKPsKs7jf3tXEwqHmdRnFIM7pzd7DNGeHD9MlGm49c2qMsYY07Tmbv+E+9a+RRhFVXli4//47pCJ/GD41GaNI9qsqi4AInIXsAd4nshU3Mux0rHGGNNkVEOUBrfi9XQh0du71n1ySwv4y9q38Ie/XtgwqGGe3fwBU44YzfAuRzRXuDENVZ2hqidW+vyYiHwM/MmlmIwxpsPIL/4vWw7+nLCWohokNelYhvV6FJ+36hDUon1r8dSypF8gHOK93V80a+KI5eZ4SEQuFxGviHhE5HIg5HZgxhjT3hX717PpwI0EwwcJazGKn8KyFazf911iXXFJAKl1jVj3xHLF8R3gIeelwBKnzRhj2rRAIMSr/1rGO/M+JxQKM+WMMVx8xQRSOiU2aT+hcJjn133G82s/pywU5KzBR3LD0eM5WPgMYQ1U2ztIWXArJYE1dEocXdE6qfco7l87v8a5EzxeTut7VJPGW596E4dTvvU890Mxxpjmo6r85v/+yeqVOygri9w3eOWFpSxbvIFHnv4+3oSme1rh/y18k6wdmygJRfp5ZvVy3tn2FY+dvJ3aBnAEL/7QXjrxdeJIT+7KLd+Yzp/XzCOsiqJ4xcM1Q09lWJc+TRZrLOpNHM5Kto8BfVR1jIgcDZyrqve4Hp0xxrhk7Zc7WfPFzoqkAeD3B9m1I49lSzYwYdLIJunnq/wD/G/HJkpDlfoJh9hXXMj2ohH0SlxBZKm+r4W1jNTEMTXOdf7Ab3JSr+Fk7VlNUCPTcQel9mqSOOMRS0p9ErgVCEDF4oW2cq0xpk1btzqHUDBco72kxM/qVTuarJ+VubvxSM17EMXBAFm7vkGCJw3BV9HukRR6d7kKXx2zq45ISeM7gydw1ZBTWiRpQGz3ODqp6idS9YvXVindGGPajF69u5Lg8xIIVB0qSkpOoM8R3ZqsnyNSuyC1JI5Ej5e+qf0ZfcRb7C54jPyS9/B6utKnyzX07DSjyfp3QyyJY7+IDMV5GFBELgJ2uxqVMca4bPzJI0hK8lFa4qfyBCav18uUM5ruZvO3+mbQIymF0mCAUKWOEjweLj3yaHzeLmR0/zUZ3X/dZH26LZahqhuAx4GRIpIDzAJ+6GpUxhjjssTEBB54/GqGDO+DL9FLYlIC/Qf24M+PXkWXrilN1o9HhJfPuoxj0vuS6PGS7E2gf+euPHf6TPqmts1nqaW+ucIi4lXVkIikAh5VPdw8odVv3Lhxmp2d3dJhGGPauAO5hwmFw6T37lrrsFJTyS0pojQYZEBnd/upj4gsV9VxDT0+lqGqLSLyNvASkNXQjowxprXqmd48f/mnp6Q2Sz9ui2Wo6kjgv0SGrLaIyCMicrK7YRljjGmt6k0cqlqiqi+r6gVEaoh3BRa5HpkxxphWKaZHI0Vkkoj8DVgBJAMXuxqVMcaYViuWJ8e3AJ8DLwO3qGqR61EZY4xptaImDhHxAs+o6l3NFI8xxphWLupQlaqGgMnNFIsxxpg2IJbpuB+JyCNEpuNWDFOp6grXojLGGNNqxZI4vuX8s/JwlQJTmj4cY4wxrV0s9ThsqMoYY0yFeqfjikgfEZktIgucz98QkWvdD80YY0w0YVXW5+9jy+EDMZeabQqxDFX9A3gG+JXzeQOR+x2zXYrJGGNMPZbu28pPl71GcTBAWJV+nbry2ISZDO3qfo2OWB4A7KWqLwNhAFUNUlutQ2OMMc1iT3EB13/4EvtLiygO+ikNBdhy+ADfef95/CH3f55jSRxFItKTr+txjAcOuRqVMcaYOs3duoqQVq1eqEBpKMCiPRtd7z+WoaqbgTeAoSKyBEgHLnI1KmOMMXXaU1yAP1zzyiKkyv7SQtf7j2VW1QoRmURklVwB1qtqwPXIjDHG1OqkPpn8Z/sXFAer/xQrx/fKcL3/WGZVzQRSVHU1MAN4SUTGuh6ZMcaYWp3W/0gGd+5Jsvfrv/1TvD5O7z+SEd3SXe8/lqGq36jqK04NjjOAvwCPASe6Gpkxxpha+Txe/j3luzz31ae8sf1LkrwJfGfoWC7IPKZZ+o8lcZQPpJ0NPKaq/xGRO9wLyRhjTH1SEnz8YNS3+MGob9W/cxOLZVZVjog8TqQGx3wRSYrxOGOMMe1QLAngYuAdYJqq5gM9gFvqO0hEkkXkExFZKSKrReROp11E5HciskFE1orITyq1/1VENorIKruPYowxrVMss6qKRWQrcKaITAOWqOq7MZy7DJiiqoUi4gMWO8uWjAIGAiNVNSwivZ39zwSGO68TsfsoxhjTKsUyq+q3wLNAT6AX8IyI/Lq+4zSifEKxz3kp8CPgLlUtfxJ9n7PPecBzznHLgDQR6RvvFzLGGOOuWIaqLgO+qaq3q+rtwHjg8lhOLiJeEfkc2Ae8p6ofA0OBS0QkW0QWiMhwZ/f+wI5Kh+902owxxrQisSSOrUBypc9JwKZYTq6qIVU9FhgAnCAiY5zjS1V1HPAk8LSzu9R2iuoNInK9k3Syc3NzYwnDGGNaBVVlX+lW9pRsIqxtd8m/Ou9xiMjDRH64y4DVIvKe8/k0YHE8nahqvogsBKYRuZKY62x6jcjKuzjtAysdNgDYVcu5ngCeABg3blzzrSNsjDGNsK90K69sv4eiYB4iHhIkkRkDbmFw52NbOrS4RbviyAaWE/lxvw14H1hIZHn1BfWdWETSRSTNeZ8CfBtYB7zO19UDJxFZph0i62Fd5cyuGg8cUtXd8X4hY4xpbYJhPy9svZX8wB4CWoY/XEJx6BCvbL+bw4EDLR1e3Oq84lDVZyEyrRYYRuRqY5OqlsZ47r7AsyLiJZKgXlbVeSKyGHhRRH4KFALXOfvPB84CNgLFwPca8H2MMSaq/cUfsinvrxQHd9DZN4xhPWbRPXmcq31uLPyUkAZrtIcJsyr/f0xIv9jV/ptatKGqBOD3wDXANiI//gNE5BngV/UtdKiqq4DjamnPJ/IUevV2BW6IK3pjjInDnsJ3WL3/l4Sdv3/zy5azYs/3Oa7P3+mR4t7s/8Jgfq33NEIaoDBw0LV+3RJtqOrPRB72G6yqx6vqcURmRKURWa/KGGPaDFXlq4P3ViSNcmEtZcPBP7nad0an0dQ2/8fnSSaznd3jmA58X1UPlzeoagGR5zDOcjswY4xpSkqA0tCeWrcVBWKaKNpgvZMzGdn1JHySVNGWIEn0TspkeJdvutq3G6I9Oa5aS/VzVQ2JiM1mMsa0KYKPBE9nguHDNbYlet1fivzc/jfzZeeFrDi4gJAGGZM2mbHdp+ERb8U+YQ2ztWgHHvGQ0ak/HmmdywJGSxxrROQqVX2ucqOIXEFkdpQxxrQZIsKgrtey5dDjhLWkot0jKQxJ+1Ez9O/hqLQpHJU2pdbtqw+t58ENT+IPR24fd0pI4ZYjf8SQzoNcjy1e0RLHDcCrInINkWm5CnwTSAHOb4bYjDGmSQ1O+z5KkG2HnkYJ4pEkhqTdSP8uF7RoXIcCBdy77lHKwv6KtlJ/GXeveZDHjv8jyd6kKEc3v2jTcXOAE0VkClB+Z2eBqv6vuYIzxpimJOJhaPcbGJx2PcFwAQmebngklrJETWPl59t46smFbNu6n7590/jetZMYf9IwFu//lHDNOwOENcynBz/nlPTWtd5rLKvjZgFZzRCLMcY0C4/4SPT2bNY+P/tsK7/65cuUlUWe59i4cS933fEqt/xiOvlDDxGo5QmHoIYoCNS8J9PSWuedF2OMaWcefyyrImnIgADeaYcJnprHo3PfYnTXI0n21ByO8oqHUV1HNHeo9Wq+azRjjOnAtm3dD4D39MN4Tiiu+PUtDBazYv8mhnbO5KvCLfid+xxJnkTGdj+KIZ0zWirkOlniMMaYZtCzV2f2hHLpM2Ev3ZJLyfGnUao+SIT/7P2Ah469ma8KN7Iodyle8TCl98mt7t5GOUscxhjTDL537TGE0n5MRvoBguohQcIsyB/Fu/nfIKzK8vx1XDhgCqcdMbGlQ62XJQ5jjGkGo456gKLSA3g9IRKJrFs1LW0te/zdWFeaia8ZZ3c1lt0cN8YYlwVDeyjzZ+P1VF3oMMkTYmq39QCc3OuYlgitQSxxGGOMy8LhfKjjiqJrQhmzhl9Kj6RuzRxVw7WdayNjjGmjfAlDEbw1amGHSSCj64X06+luPZCmZlccxhjjMhEfPdN+R6QYavny6on4PN3pnXZzS4bWIHbFYYwxzaBL6kX4EgaRf/jvBEM5pCRPIq3z9/F6e7V0aHGzxGGMMc0kOembHJHU9upvVGdDVcYYY+JiVxzGGNMG7S8t5L7V75G1ex0JHi8zMo7lxpGnkpKQ6HrfljiMMaaNKQn6uXjRE+SWFhLSMAAvbv6YVXk7ee7k7yFSs755U7KhKmOMaWPm53xJgb+kImkA+MMh1uTvZlXeTtf7t8RhjDFtzKqDORSHatbvUJR1h/a63r8lDmOMaWOGdOlFstdXo90rHgamdne9f0scxhjTxszIOBafx0vlOxkJ4qF3chfGpw92vX9LHMYY08Z0S0zhxVOu4eju/fGKhwTxcHKfYTx3yvfwiPs/6zaryhhjXLar4DB/XLiIRZu3kJiQwMyjxnDThJNISoj9J7i0xM/zf5nPe698QjAY4sSpY3j01xeT0jMFr3hI9Dbfz7moVl92q+0YN26cZmdnt3QYxhhTp4LSUk6b/Q/ySkoIO7+3SQlevjlgAP+YeWFM51BVfj7zYdav3E7AqVvu9XpI69WZJxf+ipTUmvXKoxGR5ara4JUV7YrDmHZm+948XvvgCw4WFDPhqMFMGTuMhARvS4fVYc35cjVFfn9F0gAoC4bI3pnDun25jOydXu851n++jY1f7qxIGgDBUJh9ncL88O8vM+zI/px/3DcY0/8IV75DdZY4jGlHspZ/xW9mLyAUChMMhcla8RX/fG85j//8YpJ89p97S1i5ezelwWCNdq8I63JjSxxb1u6m8uiQAvtP6UlxRid25O8n+9MDvPrZan486US+P/GEpgy/VnZz3Jh2wh8Icucz71DmDxIMRR4MKykLsDFnP28s/rKFo+u4RvTqRZK35hWfAoO6p1VpKyguZdmGbWzYtb9Ke7/MXng8X/9clx6RTHFGJ9TnAYGwKqWBII8sXMaeQ4dd+R6VWeIwpp1Ys7X2B79K/UHe/nhdM0djyl1y9FH4qiUOn8fD4O7dObZv34q2x9/9mKl3PMHN/5jHlX/9Fxf9+XlyCwoBOGr8UNL7peF1hhyLB6WgCTWXFfGKh8Ubt7n4bSIscRjTTiQnJlDXZJdOSe4vfGdq1ys1lX9fdglHH3EEHhF8Hg+nDR/Gc5dcVLGm1MLVm3g66xPKgiEKS/2U+INs2nOAm55+EwCPx8OfXvl/nPjtb+BN8OAJKlAzcYhAcjMMSdqgpzHtxJEZvemamkxxWdWlKFKSfFx46tEtFFX7FwqHeWpFNs98voJDZWWM7duXX59yKqPSe1fsM7J3Oq9e+R3KgsFI8qh2BfLCB59R4q96HySkyle797PzwCEG9OxGtx6d+c0T1xIMhPhqz36+88zLNe6dqMLkI4e492UddsVhTDshIjx40/l075JCanIiKUk+En1eZpxyFJOOHdrS4bVbdy7K4q+fLGVfcRFloSBLd+5g5px/szU/r8a+SQkJNZIGQH5hSa3nTvB4KCgurdrm8zJqYB9+Pm0iSQleOiX6SE1KpFOij4cvO4fUZri6tCsOY9qgYCDI4tezWbbgc7r37sqZV08iY2R/hvXvxYI/X8+yNdvILyxl7Ij+9OvVrUVjLSwr47VVa1m8dhNlG3IZuKOUSaeNZeLMk0hMqrneUltysKSYl9d8iT8UqtJeFgzy+PJP+cPU02M6z6ljhrA1Nw9/sOp5EBjWt2etx1x2wjFMGz2CJZu2kej1cvLwTDolNs+/T0scxrQx/rIAP5/2e7as3klpURneBA9vPZXFzY9dx6kzx5OQ4OXko90frojFvsOFXDD7n+QXleDXMHhDSIay7I4XmHP/mzy4+B6SO8X38FprsjU/nySvt0biCKmyau+emM9z5cTjeePTteQVFlMWDCESuTq57YIpJEZ5urx7agrTjx7Z4PgbyoaqjGlj3nvhQzZ/uYPSojIAQsEwZSV+HrhhNmUl/haOrqr73l/CweLiSNIA8HnRpAR2Tctg5/pdzPv7uy0bYCMN7NqtRtIA8IgwomftVwq16ZaazJxbruAHp49n7JD+nHnskcz+8UzOGTeqKcNtMq4lDhFJFpFPRGSliKwWkTurbX9YRAorfU4SkZdEZKOIfCwimW7FZkxbtvCVZZQV10wQHo+HdZ9uaoGI6pa1YTPBcM2ZXoH0FErCId5/aUkLRNV00lNTOX3oMJKrrROV5PXyg+PjexCva0oy1337BP5x48X88cqzOGpQ8zwF3hBuXnGUAVNU9RjgWGCaiIwHEJFxQFq1/a8F8lR1GPAAcK+LsRnTZqV0Tq61XTVMUqfWNe02qa6lThQIKymptX+XtuRP357GpWOOJiUhAY8Iw3r04OlzL2Bkr/qfCG+rXLvHoZEJ5eVXFD7npSLiBf4MfAc4v9Ih5wF3OO/nAI+IiGhbXoXRGBdMv24qKxetpbS4rEp7ardURox1vxZDPC4ZexRPfJRNWeVpo8EwKRvz6ZTo49wfn9FywTWRpIQEfjtpMr+eeCqBUCiuFW/bKlfvcYiIV0Q+B/YB76nqx8CNwBuqurva7v2BHQCqGgQOATUGCUXkehHJFpHs3NxcN8M3plX65hlHc+4Pv40vyUdK52RSuiSTlt6Ve179vyrLUrQGP5hwAidlDiTJ68XjDyH+EEkHSzni7W1M+95kTrlwfEuH2GQ8Ih0iaUAzLasuImnAa8DtwO+BU1U1KCKFqtrZ2Wc1cIaq7nQ+bwJOUNUDdZ23JZdV94cK2FzwHw6UrqZb0lCGdj2flIReLRKL6Zhycw7yxYfr6NI9leOmjCahFS9iuGHfflbv2kvxpv30KRGOOmUUfQa136Gc1q6xy6o3Wz0OEbndefsjoPyJlgxgs6oOE5F3gDtUdamIJAB7gPRoQ1UtlTiKAnt4b8dVBLWYkJbhIRGP+Jg64EnSkoY3ezzGGBOPxiYON2dVpTtXGohICvBtYLmqHqGqmaqaCRQ7N8MB3gC+67y/CMhqrfc3Pt//IGXhQ4Q0MsYcxk9Qi/h03z0tHJkxxrjPzWvbvsCzzs1wD/Cyqs6Lsv9s4HkR2QgcBC51MbZG2VP8ERCu0Z5Xto5QuAyvp+0+0GSMMfVxc1bVKuC4evbpXOl9KTDTrXiakkeSQGuuLSN4ieRJY4xpv1rXFIw2YkjX8/BK1asKDz76dz4Vj7TeG5TGGNMULHE0wJie15OechxeSSJBOuGVFLolDWNc+q0tHZoxxrjO/jxuAK8kMqnfw+SXbeSQfyOdfQPpkfSNiqIsxhjTnlniaIS0pGGkJQ2rf0djjGlHLHE0wvLcrSzes4FETzJjegzkhN4DSfK23X+lm/ccYPu+fIb168mAXtWXEnPXIX8OB/1bSfMNpHtSBiWhIrYWrSfZk8Kg1BF4mnjSQVkgyPKtOQCMG9w/6tLVsQgGQ6xak0MgGGLMyH5s3Lmfw8VlHD28H12rrcekqny1cS/7DxRy5PAj6Nmzc43zbd+wm12b9zFoZD/6ZrbtB+VUlQ1b95F7oJAjh/QmvUeXlg4pJkXBQjYVbiQ1IZXBqUPxiI3sl2u7v3ItqCTgZ+bCBzjEfsofNPFvS6DkcA/+NvEiTunbutYLqk9xqZ9Zj7/Bqi278Xk9+IMhJh41hN9/78xaq5U1pZAGeG/X3Wwr+hiv+AhrkBRvX9YVleIRH6AkeVK4dsiv6JsyqEn6/HD9Vv7vn28hRNbaE+D+y6czYUTDzv/luhx+ec+rBINhQiiHfWF8iV58CV4CwRA3XHwKl54xFoCDBwu55baX2bU7H69H8AdCnDv9OG74wRREhJLCUu666u+s+XQTCT4vQX+Ib542hl88fi2+xLb3n2teQTE//d1ctu/Kw+sRAsEQ0yeP4f+undqqh3bf2f0Wb+x6Da8koCidEzoza8Qt9EluvSvWNidLoQ1w07LnyWc/4gGP80r0BUnslM8PFs4hr6y4pUOMyx9ffp+Vm3dRFghSWOrHHwzx4ZdbeOrtT1zv+9P9z7Kt6BNC6scfLiKoZRwKbKWrN4+ycAll4VIKgnk8ufluwlqz7kG8DhQWM+uFNykq81NY5q/450+ef4O8otrLd0ZTWhbgljvnUHC4lKISPwUJIRTFHwhRVOLHHwjx2CuL+Xx95Ormjt/9h23b91NaGqCo2E8gEOKtBSv5b9YaAP7+q5dZ/clG/KUBig+X4i8L8On/vuRf989v9HdvCbc/9Babt++ntCxQ8e9j/qLVvJn1RUuHVqd1BWt5c/frBDRAqfP/wYP+Azy04S+00meSm50ljgZYW7KB6mvJeTzQKcUPKPO3rWuRuBoiGArz9vL1NUpWlgWCzPlgpev9r85/s+IJ/HIegVRvGfD1f6SBsJ/NhWsa3d/bqzZQ13/773yxIe7zffTpporzqZfI5Uu1v6TL/EHm/m8l+w8cZt363YRCVQMoLQ0w57VsQqEw78/5hEBZsMp2f0mA+f/4IO7YWlp+QTEr1+YQDFV9WLa0LMgrCz5roajq9/6+/+IPV613oiiHgwVsK97aMkG1MpY4GkLq/qsjEA5zOFBW5/bWJhQOEwrVfAoeoNgfcL3/oJbW2l7bIEZpuPFXcoWlZQRqqdgWCIU4XBL//27FxX7CYeffX/nYVzUKHCoqobjYj9db+39yRUVlhENhgtVrTjv+f3vnHmVFcefxz3fmzgwDyENwUUGF8BAJCAF8RHyAEp8cYaMRRTfHd+LqutGTmE3WNequWV2zPjYmGnWNya6rxhUV37oRj4/gAxUFxQAGElEPghEUBJTht39UXaa503eYnpk7cy/5fc65Z25XV1d96ze3+tdV3f2rDSkLN5U7n234gqqq9OmotZ+Vbx9Zt2ltanoVVaxvqKzZhFLhjqMVdLdeqVetX2yqJldVzYE7V849jrqaHMP6N43qWyWx3567l7z+/vVjSHMTn1tuq/QG28Sgbm1fRvOrQ/egNmVxoZrq6lbd4xi79+5sjivcaROpHq9LbY5J44fSf9fe1NU1vU+Ry1Ux4YAh1NTmGDxytyb7JTH6oD0za+tsdu7bgx26NQ2/k6uu4sBxgztBUcsY23s8NVVNX32SsAAADmFJREFUF8RqoIFB3cpXd0fijqMVXDb6JGyzyF9omsHmzbDu0x5M2WMvRvaprBtoF8+YTNe6GnLxarg2V033+lou/PohJa/7wH7nUVvVjSpqAKgih1HFJw29t+SpUR2H7zydbrkeba5v1IB+fG3kEOpra7ak1dfWcMSoYYzo3y9zebvu3IvjpoyjS10NAqo32Fajji61OfbYdUeOnjCC6uoqvn/h0dTV5bZcidfV5ujduxszTgjrUpx/zcl06VZHriY4t5q6HF17dOHsy49vfaM7iaoq8cNzjqRLQXt79ajntDJeh2NC34PpV9ePWgXnIURtVS0n7DaDLtWVv2Jhe9BhYdVLQWeux7F49Qqumj+Ld9a9R0NDjv51u3H6ngcwecCQsn5apBjvf7SGu56ex+L3VzFq4C6ccPBo+vbs1iF1f7bpz8z/+H5WbFhI37ovMaLnsSxZt4jXV8+ha3V39u/zNQZ1b/toI4+Z8dRb7/DAK2+BYNrYLzNpxJda/X8zM16et4yHn5zPho1fsOfwXfjjR6tZs3YDk8YP5cgD9qIu8UTU0mUrue+BV/lgxWrGjx3EMUeNpnviynzFux/x4K2z+cNb7zF83CCmnHYIO/br2eZ2dxbLln/EPY+9xnsrVjNu5O5MPWxvehRZ/rZc+Hzz58xZ9TzzVr9Cj5oeTNzpMAZ1335GGxWzHkcp6EzH4TiOU6mU7XocjuM4zvaJOw7HcRwnE+44HMdxnEy443Acx3Ey4Y7DcRzHyYQ7DsdxHCcT7jgcx3GcTLjjcBzHcTLhjsNxHMfJhDsOx3EcJxPuOBzHcZxMuONwHMdxMuGOw3Ecx8mEOw7HcRwnE+44HMdxnEy443Acx3Ey4Y7DcRzHyYQ7DsdxHCcT7jgcx3GcTLjjcBzHcTLhjsNxHMfJhDsOx3EcJxPuOBzHcZxMuONwHMdxMlEyxyGpi6SXJL0u6U1Jl8X0OyT9XtICSbdJqonpkvQfkpZIekPS2FJpcxzHcVpPKUccG4FDzWw0MAY4UtL+wB3AcGAUUA+cGfMfBQyNn7OBG0uozXEcx2klJXMcFlgbN2vix8zskbjPgJeAATHPVODXcdcLQC9Ju5RKn+M4jtM6SnqPQ1K1pHnAh8CTZvZiYl8N8DfAYzGpP/Bu4vDlMc1xHMcpI3KlLNzMGoAxknoB90kaaWYL4u6fA8+Y2bNxW2lFFCZIOpswlQWwUdKCwjwVRF9gVWeLaAOuv/OoZO3g+jubPdtycEkdRx4zWy3paeBIYIGkHwE7Ad9KZFsO7JbYHgC8n1LWzcDNAJLmmtn4UukuNa6/c6lk/ZWsHVx/ZyNpbluOL+VTVTvFkQaS6oHJwNuSzgSOAE4ys82JQ2YB34xPV+0PrDGzD0qlz3Ecx2kdpRxx7AL8SlI1wUH9xswekrQJ+CMwRxLATDO7HHgEOBpYAnwGnFZCbY7jOE4rKZnjMLM3gK+kpKfWGZ+yOjdjNTe3Qlo54fo7l0rWX8nawfV3Nm3Sr3C+dhzHcZyW4SFHHMdxnEyUreOo9JAlxfQn9v9U0trEdp2ku6P+FyUN7GjNSZqxvyRdIWmRpIWSzk+kl739JR0m6VVJ8yQ9J2lITC8r+0dN1ZJek/RQ3B4UtS2OWmtjetlph1T9FdF38xTqT6SXdd/Nk2L/duu7Zes4qPyQJcX0I2k80Ksg/xnAx2Y2BLgWuKojxaZQTP+phMemh5vZXsBdMX+l2P9G4GQzGwP8D3BxzF9u9gf4e2BhYvsq4FozGwp8TNAM5akdmuqvlL6bp1B/pfTdPIX6T6Wd+m7ZOo5KD1lSTL/CU2ZXAxcVHDIV+FX8/r/AYZLSXorsEIrpB84BLs8/Sm1mH8Y8FWH/+OkR03vS+K5QWdlf0gDgGODWuC3g0KgNgtZp8XtZaYem+gEqpe9Cuv5K6buQrp927Ltl6zig8kOWFNF/HjAr5R2VLfrNbBOwBujTkXoLKaJ/MDBd0lxJj0oaGrNXiv3PBB6RtJzw+7kyZi83+19HOEHl33XqA6yO2mBr+5abdmiqfwuV0HdJ118xfZd0/e3Wd8vacZhZQ5xSGADsK2lkYnerQpZ0JCn6Dwa+Afw0JXsl6B8J1AEb4luztwC3xeyVov8C4GgzGwD8ErgmZi8b/ZKmAB+a2SvJ5JSs1oJ9HU4R/UnKuu+m6Ze0KxXSd5uxf7v13Q4JOdJW2jNkSWeQ0D8JGAIsiSPZrpKWxLnRvP7lknKEaZQ/d5LkrSiw/3Lg3rjrPsLJFyrD/kcBoxMj17tpvOotJ/tPAI6VdDTQhTC1dh1hCiEXr2qT9i0n7ZCiX9J/m9kpFdJ30+z/JuG+WSX03VT7055918zK8kP4cfWK3+uBZ4EphKmG3wH1BfmPAR4leM/9gZfKUX9BnrWJ7+cCN8XvJxLetC87/YSpndNj+kTg5UqyPyEw3bCYfgZwbznaP9GOicBD8fs9wInx+03A35az9hT9FdF3i+kvSC/bvtuM/dut75bziKPSQ5ak6m8m/38C/yVpCeFq5cQO0Ngcxez/HHCHpAuAtTQ+GVMR9pd0FnCvpM2EJ5NOj/nLzf5pfB+4S9K/AK8RNENlaIfg7Cqh72alUux/Je3Ud/3NccdxHCcTZX1z3HEcxyk/3HE4juM4mXDH4TiO42TCHYfjOI6TCXccjuM4TibccWzHSGpQiAKb//xDJ2i4VNJ3U9IHSlrQ0XpaS9S7vsCetZ2tqy1IOlXSSoUIqoslPS7pgBYcN03SiFbUN03SJfH7pZLei3ZcLGlmskxJtZKuk/RO3P+ApAEKPCfpqETeEyQ9Fo95Jr6E55QQN/D2zXoLITec9uGd5uyZeKu7krjbzM4DkDQJmClpkpktbOaYacBDwFsZ67oIODaxfa2Z/STWPR14StIoM1sJ/BjYgfCyZoOk04CZwH7At4F7JM0GqoErgCPN7HNJvwWmEyLxOiXCRxx/gUhaJukyhXUp5ksaHtMPSVxNvyZph5j+PUkvK8Tqz69rMVDS25JuVVhf4Q5JkyU9H68Q901UOVrSUzH9rBQ91ZKuTtTxrcI8Md/tkm6UNFvSH6Le2xTWFrg9ke9wSXNi++6R1D2mXxLrWCDpZsW30CQ9LekqhfU7Fkk6KIMtL41lPQH8urm2FLHjtxM2XxpPhkXbUFD3RCXWipB0g6RT4/dliTa9pLjuSHOY2WzCkqJnxzLOinpfl3SvpK5xRHIscHXUPDh+HpP0iqRn87+nAq3DgI1mtqpI3XcDTwAzJHUlvIR2gZk1xP2/pDFU/gLgQcILkT8iRHZ9JxZ1P3DyttrqtA13HNs39dp6amV6Yt8qMxtLiL2fn0r6LnBuvKo+CFgv6XBCnP59CetajFMI1ggh7tb1wN6EdRZmAAfGcn6YqGtvQliDrwKXKASMS3IGsMbM9gH2Ac6SNKhIm3oTwotfQDh5XAt8GRglaYykvoQ1NibH9s0FLozH3mBm+5jZSEIYkimJcnNmti/wHcLJKI3BCVv+LJE+DphqZjOKtaWYHc3spmjvfQgxg67ZRhuy8Els0w2EWFct4VXC/xLCm937WFjTZCFwhpn9DpgFfM/MxsQT9s3A35nZOML//ucp5U6IZbek7iHAn8zsk4L9cwn/a4DLCL+3o4B/S+RZQLClU0J8qmr7prmpqpnx7yvA1+P35wknrjsIJ43l8YR3OCHEBUB3wgnwT8BSM5sPIOlN4LdmZpLmAwMTdT1gZusJjmg24eQ5L7H/cGBvScfH7Z6xjqUpuh9M1LGioP6BhABtI4Dn44CiFpgTj50k6SKgK7AjIXDdgyn2SGpPUmyqalZsX3NtKWbHZ+L29cBTZvagQnTTYm3Iwp2Jv9e28JhkpNSRCuFNekW9jzfJHEZCBxCmjvLJdSnl7gKsbGHdIj0665Z0M1sn6W5CzKiN+QxxWutzSTuY2afbqM9pJe44/nLJd7YG4u/AzK6U9DAhbs0LkiYTOuu/mtkvkgcrLI+5MZG0ObG9ma1/W4UngcJtEa5YtzoxSbqCMFIhccJO1lFYfy6250kzO6mgrC6EK+HxZvaupEsJkUPzNLFHBta1oC1HkGLHuO9UYA/Ceg/5MtLasB+QP/4SQlyk5KxBsj2wtZ1bGlvoKzSuGnc7MM3MXo8aJ6bkryKsE7Kte2nrCU50W3XPJcRM2iPl5D+WRkcP4X/eZL0PYvjwbdTltAGfqnK2IGmwmc03s6sIHXg44Srz9MR9gv6S/ipj0VMV1gDvQ4zKWbD/ceAcNa5BPUxSNzP7xzgdkuUG/wvABDWuJd41zq/nT6qrYluOL1ZAG0ltC0XsKCk/vXOKxZXZirXBzF7M28PMZhECBo5QWPO6J3BYgZbpib/bHLFIOoRwf+OWmLQD8EFsS/K+wadxH3E6aamkb8QyJGl0SvELCVNQxeo+jjAiu9PM1hFW1LtGIUglkr5JGCk+tY029AFWmtkX22iu0wZ8xLF9U6+wAl6ex8ysuUdyv6PwZE0D4YmZR81so6S9aIxouhY4JeZpKS8BDwO7A/9sZu/HEUueWwnTQ68qVLKSxmVRM2FmK+PV8Z2S8lMmF5vZIkm3APOBZTR1Xu1FalvM7IkidjyPMG02O6bPNbMz09oALCpo67uSfgO8ASymcRosT52kFwkXiCeRznRJBxJOykuB4xJPVP0T8CLBQc0nOgvCWtW3SDqf4IBPBm6UdDFhid67gNcL6nkG+HdJssbIqhdIOgXoRrg3cWh8ogrgB8BPgEUKkYzfBv46cWwxJhGivTolxKPjOs52iKRlhGm51KeYOgNJ1xPuUf1fCeuYCfzAzH5fqjocn6pyHKfj+DFhZFMSFF7IvN+dRunxEYfjOI6TCR9xOI7jOJlwx+E4juNkwh2H4ziOkwl3HI7jOE4m3HE4juM4mXDH4TiO42Ti/wG3Tr9rpoqFiAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = plt.figure(num=1, figsize=(12,4))\n", "ax1 = fig.add_subplot(111)\n", "ax1.plot(years_all, X_all, '0.75', linewidth=1) # all ensemble members\n", "ax1.plot(years_all, X_all.mean(axis=1), 'k-', linewidth=2, label='ensemble mean') \n", "ax1.plot(years_all, Y_all, 'ro', label='observation')\n", "ax1.grid(linestyle='--',color='0.5',linewidth=1.0)\n", "ax1.legend(loc='upper left')\n", "ax1.set_ylabel('Freeze-up Date (DOY)')\n", "\n", "fig = plt.figure(num=2, figsize=(6,6))\n", "ax2 = fig.add_subplot(111)\n", "scatter = ax2.scatter(X_all.mean(axis=1), Y_all, c=years_all, cmap='viridis')\n", "ax2.set_title(\"correlation=%s\" %str(np.around(pearsonr(X_all.mean(axis=1), Y_all)[0],3)))\n", "legend = ax2.legend(*scatter.legend_elements(),\n", " loc=\"upper left\", ncol=1)\n", "ax2.set_xlim((320,460))\n", "ax2.set_ylim((320,460))\n", "ax2.set_xlabel('Ensemble-mean Freeze-up Date (DOY)')\n", "ax2.set_ylabel('Observed Freeze-up Date (DOY)')\n", "ax2.add_artist(legend)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The hindcast freeze-up date typically occurs much later than the observed freeze-up date; However, there is a relatively strong correlation between the ensemble mean and the observation (0.61). \n", "\n", "We'll now create some time-like variables required for calibration:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "a=334, b=455\n" ] } ], "source": [ "im=12 # define the initialization month (December)\n", "\n", "# instantiate the sitdates class and provide the event of interest\n", "si_time = sitdates.sitdates(event='fud')\n", "a = si_time.get_min(im) # minimum date possible\n", "b = si_time.get_max(im) # maximum date possible\n", "print(\"a=%d, b=%d\"%(a,b))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the above section of code, we've made use of the `sidates` module by instantiating it with the event variable (in this case 'fud' for *freeze-up date*; for the retreat season we would set event='ifd'). By doing so, and by calling on `si_time.get_min()` and `si_time.get_max()` with the initialization month as an input argument, the minimum and maximum dates possible are returned. These are printed as `a=334` and `b=455`. Reading the docstring for `sidates.sitdates()` you'll see that default dates are set to those given in [1]." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### *Useful side-note*:\n", "If rather different conventions are used for your data, you can change the dates associated with `si_time.get_min()` and `si_time.get_max()` using `si_time.set_min()` and `si_time.set_max()`. For example, for `im=12`, the convention in [1] is to consider the freeze-up date for the current freeze season. If rather your hindcast and observed data were created so that the freeze-up date is for the following season, you could change the minimum and maximum dates corresponding to `im=12` accordingly with:\n", "\n", "```python\n", "si_time.set_min(im, 538) # corresponding to the start of the freeze season (October 1) for the next year\n", "si_time.set_max(im, 699) # November 30 of the following year (if e.g. the forecast was run for 12 months)\n", "```\n", "Those dates are then set and the minimum and maximum dates can be retrieved at any time within the working script using:\n", "\n", "```python\n", "a = si_time.get_min(im)\n", "b = si_time.get_max(im)\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We'll now define the year of the forecast to be calibrate and split the data in ``X_all`` and ``Y_all`` into the training and testing data. For this example, we'll emulate an operational scenario and only use data from years preceding the year being calibrated to train the calibration model." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "t = 2016 # the forecast initialization year\n", "tau = years_all[years_all" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Observed date: 01/17\n", "Observed category: Later than normal\n" ] } ], "source": [ "#########################################################\n", "######################## PLOTTING #######################\n", "#########################################################\n", "dcnorm = dcnorm_gen(a=a, b=b) # instantiate instance of dcnorm_gen (the DCNORM distribution module) with given a and b \n", "x = np.linspace(a,b,1000) # discretize the support\n", "rv = dcnorm(mu_t, sigma_t) # freeze DCNORM distribution object with calibrated distribution parameters\n", "rv_c = dcnorm(clim_params[0], clim_params[1]) # # freeze DCNORM distribution object with the climatology distribution parameters\n", "\n", "event='fud'\n", "\n", "# A way to convert day-of-year dates to regular mm/dd for x-tick labelling in the plot\n", "xticks = []\n", "xticklabels = []\n", "count = a \n", "while count=clim_terc[1]:\n", " plt.setp(p, 'facecolor', 'b',alpha=0.5)\n", " else:\n", " plt.setp(p, 'facecolor', 'g',alpha=0.5)\n", " \n", " idx+=1\n", " \n", "# insert axes for bars showing probabilities\n", "ax2 = inset_axes(ax,\n", " height=\"45%\", # set height\n", " width=\"45%\", # and width\n", " bbox_to_anchor=(0.0,0.0,0.95,0.95),\n", " axes_kwargs={'frame_on':False},\n", " bbox_transform=ax.transAxes)\n", "center=0.5\n", "offset=1.0\n", "\n", "# compute forecast probabilities for the raw ensemble\n", "p_en_raw = dcnorm.ecdf(clim_terc[0],X_t)[0]\n", "p_nn_raw = dcnorm.ecdf(clim_terc[1],X_t)[0] - p_en_raw\n", "p_ln_raw = 1 - dcnorm.ecdf(clim_terc[1],X_t)[0]\n", "\n", "# creates a bar for each forecast probability in the inset axes\n", "rect1 = plt.bar(center-offset/2,p_en_raw ,0.4, fc='r',alpha=0.5)\n", "rect2 = plt.bar(center,p_nn_raw, 0.4, fc='g',alpha=0.5)\n", "rect3 = plt.bar(center+offset/2,p_ln_raw, 0.4, fc='b',alpha=0.5)\n", "\n", "plt.xlim(-0.5, 1.5)\n", "plt.xticks([center-offset/2,center,center+offset/2])\n", "plt.yticks([])\n", "ax2.set_xticklabels(np.array(['E','N','L']))\n", "\n", "ax2.text(center-offset/2, p_en_raw+0.01,str(int(round(p_en_raw*100,0)))+'%', \n", " fontsize=11, ha='center',fontweight='semibold')\n", "ax2.text(center, p_nn_raw+0.01,str(int(round(p_nn_raw*100,0)))+'%', \n", " fontsize=11, ha='center',fontweight='semibold')\n", "ax2.text(center+offset/2, p_ln_raw+0.01,str(int(round(p_ln_raw*100,0)))+'%', \n", " fontsize=11, ha='center',fontweight='semibold')\n", "\n", "ax.set_xlim((a-10.,b+5))\n", "ax.set_ylim((0.0,rv.pdf(x).max()+0.05))\n", "\n", "ax.set_xticks(xticks)\n", "ax.set_xticklabels(xticklabels,fontsize=11)\n", "ax.set_yticklabels(np.arange(0.0,rv.pdf(x).max()+0.05,0.01),fontsize=11)\n", "ax.set_ylabel('Probability density', fontsize=13)\n", "ax.set_xlabel('Freeze-up date', fontsize=13)\n", "\n", "ax.set_title('Raw '+event.upper()+' forecast from '+\"%02d\" % im+'/01/'+str(int(t))+\n", " ' \\n E=Early, N=Near-normal, L=Late \\n relative to '+str(t-n_clim)+'-'+str(t-1)+\n", " ' climatology', fontsize=14)\n", "\n", "\n", "\n", "######### Now plot the probability distribution for the calibrated for the calibrated forecast #######\n", "ax = fig.add_subplot(122)\n", "# plot the climatological distribution first\n", "ax.plot(x[1:-1], rv_c.pdf(x[1:-1]), color='0.5', linewidth=2)\n", "# now plot end-dates if probability exists for pre/non occurence of the event\n", "if rv_c.pdf(x[0])>1e-4:\n", " ax.plot(x[0], rv_c.pdf(x[0]), 'o', color='0.5', ms=5)\n", "if rv_c.pdf(x[-1])>1e-4:\n", " ax.plot(x[-1], rv_c.pdf(x[-1]), 'o', color='0.5', ms=5)\n", "\n", "# plot main forecast pdf curve (but exclude end-dates)\n", "ax.plot(x[1:-1], rv.pdf(x[1:-1]), 'k',linewidth=2)\n", "# now plot end-dates if probability exists for pre/non occurence of the event\n", "if rv.pdf(x[0])>1e-4:\n", " ax.plot(x[0], rv.pdf(x[0]), 'ko', ms=5)\n", "if rv.pdf(x[-1])>1e-4:\n", " ax.plot(x[-1], rv.pdf(x[-1]), 'ko', ms=5)\n", "\n", "# fill colors under the forecast pdf according to each category\n", "ax.fill_between(x[xclim_terc[0])&(xclim_terc[0])&(xclim_terc[1]], 0.0, rv.pdf(x[x>clim_terc[1]]),color='b',alpha=0.5)\n", "\n", "# plot climatological terciles as dashed lines\n", "ax.vlines(clim_terc[0],0.0,max(rv.pdf(x).max()+0.05,rv_c.pdf(x).max()+0.05),\n", " colors='0.3', linestyles='--', linewidths=1.5)\n", "ax.vlines(clim_terc[1],0.0,max(rv.pdf(x).max()+0.05,rv_c.pdf(x).max()+0.05),\n", " colors='0.3', linestyles='--', linewidths=1.5)\n", "\n", "# insert axis for bars showing probabilities\n", "ax2 = inset_axes(ax,\n", " height=\"45%\", # set height\n", " width=\"45%\", # and width\n", " bbox_to_anchor=(0.0,0.0,0.95,0.95),\n", " axes_kwargs={'frame_on':False},\n", " bbox_transform=ax.transAxes)\n", "\n", "center=0.5\n", "offset=1.0\n", "\n", "# plot bars showing probabilities\n", "p_en_ncgr, p_nn_ncgr, p_ln_ncgr = fcst_probs\n", "\n", "rect1 = plt.bar(center-offset/2, p_en_ncgr, 0.4, fc='r',alpha=0.5)\n", "rect2 = plt.bar(center, p_nn_ncgr, 0.4, fc='g',alpha=0.5)\n", "rect3 = plt.bar(center+offset/2, fcst_probs[2], 0.4, fc='b',alpha=0.5)\n", "\n", "\n", "plt.xlim(-0.5, 1.5)\n", "plt.xticks([center-offset/2,center,center+offset/2])\n", "plt.yticks([])\n", "ax2.set_xticklabels(np.array(['E','N','L']), fontsize=12)\n", "\n", "ax2.text(center-offset/2, p_en_ncgr+0.01,str(int(round(p_en_ncgr*100,0)))+'%', \n", " fontsize=11, ha='center',fontweight='semibold')\n", "ax2.text(center, p_nn_ncgr+0.01,str(int(round(p_nn_ncgr*100,0)))+'%', \n", " fontsize=11, ha='center',fontweight='semibold')\n", "ax2.text(center+offset/2, p_ln_ncgr+0.01,str(int(round(p_ln_ncgr*100,0)))+'%', \n", " fontsize=11, ha='center',fontweight='semibold')\n", "\n", "ax.set_xlim((a-10.,b+5))\n", "ax.set_ylim((0.0,max(rv.pdf(x).max()+0.05,rv_c.pdf(x).max()+0.05)))\n", "\n", "ax.set_xticks(xticks)\n", "ax.set_xticklabels(xticklabels,fontsize=11)\n", "ax.set_yticklabels(np.arange(0.0,rv.pdf(x).max()+0.05,0.01),fontsize=11)\n", "ax.set_ylabel('Probability density', fontsize=13)\n", "ax.set_xlabel('Freeze-up date', fontsize=13)\n", "\n", "ax.set_title('Calibrated '+event.upper()+' forecast from '+\"%02d\" % im+'/01/'+str(int(t))+\n", " ' \\n E=Early, N=Near-normal, L=Late \\n relative to '+str(t-n_clim)+'-'+str(t-1)+\n", " ' climatology', fontsize=14)\n", "\n", "fig.subplots_adjust(top=0.8, bottom=0.1, left=0.1, right=0.99, wspace=0.2)\n", "\n", "plt.show()\n", "\n", "print(\"Observed date:\", si_time.doy_to_date(int(Y_t)))\n", "if Y_t<=clim_terc[0]:\n", " p_en_obs = 1.0\n", " p_nn_obs = 0.0\n", " p_ln_obs = 0.0\n", " print(\"Observed category: Earlier than normal\")\n", "elif (Y_t>clim_terc[0])&(Y_t<=clim_terc[1]):\n", " p_en_obs = 0.0\n", " p_nn_obs = 1.0\n", " p_ln_obs = 0.0 \n", " print(\"Observed category: Near normal\")\n", "else:\n", " p_en_obs = 0.0\n", " p_nn_obs = 0.0\n", " p_ln_obs = 1.0\n", " print(\"Observed category: Later than normal\")\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To compute the CRPS for the raw forecast, we'll use the ``dcnorm.ecdf()`` function, and for the calibrated forecast we'll use the ``dcnorm.cdf()`` function:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CRPS (raw) = 75.88238238238236\n", "CRPS (NCGR, approximate) = 13.487902217819506\n" ] } ], "source": [ "crps_raw = np.trapz((dcnorm.ecdf(x,X_t) - dcnorm.ecdf(x,Y_t))**2.,x) \n", "crps_ncgr = np.trapz((rv.cdf(x) - dcnorm.ecdf(x,Y_t))**2.,x)\n", "print(\"CRPS (raw) =\",crps_raw)\n", "print(\"CRPS (NCGR, approximate) =\",crps_ncgr)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we can also compute the CRPS for the calibrated forecast using the analytic expression for the CRPS with the NCGR library using the following:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CRPS (NCGR, analytic) = 13.442205531634983\n" ] } ], "source": [ "crps_funcs = ncgr.crps_funcs(a,b) # instantiate instance of ``ncgr.crps_funcs``\n", "crps_ncgr_true = crps_funcs.crps_dcnorm_single(Y_t, mu_t, sigma_t) # ``crps_dcnorm_single`` is used to compute CRPS for a single forecast\n", "print(\"CRPS (NCGR, analytic) =\",crps_ncgr_true)" ] } ], "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.6" } }, "nbformat": 4, "nbformat_minor": 4 }