GEVattempt1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-20T19:48:06.365861Z",
     "start_time": "2018-09-20T19:48:04.975795Z"
    }
   },
   "outputs": [],
   "source": [
    "import scipy as scipy\n",
    "import numpy as np\n",
    "import scipy.io as sio\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "from scipy.stats import genextreme as gev\n",
    "from netCDF4 import Dataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-20T19:48:06.471913Z",
     "start_time": "2018-09-20T19:48:06.368686Z"
    }
   },
   "outputs": [],
   "source": [
    "CaMaFile = Dataset('fldare1161.nc')\n",
    "# print(CaMaFile.dimensions.keys())\n",
    "\n",
    "lat =  CaMaFile.variables['lat'][:] \n",
    "lon = CaMaFile.variables['lon'][:] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in Matlab Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-20T19:50:32.368662Z",
     "start_time": "2018-09-20T19:49:09.575037Z"
    }
   },
   "outputs": [],
   "source": [
    "#matData   = sio.loadmat('/Users/meganfowler/Documents/MATLAB/FloodProject/CaMaFlood/Physiology/30yr/control_outflw.mat')\n",
    "\n",
    "matData   = sio.loadmat('control_outflw.mat')\n",
    "maxYearly = matData['maxYearly_full']\n",
    "\n",
    "fullData   = sio.loadmat('full_outflw.mat')\n",
    "fullYearly = fullData['maxYearly_full']\n",
    "\n",
    "physData = sio.loadmat('physiology_outflw.mat')\n",
    "physYearly = physData['maxYearly_full']\n",
    "\n",
    "radData = sio.loadmat('radiation_outflw.mat')\n",
    "radYearly = radData['maxYearly_full']\n",
    "\n",
    "# print(maxYearly.shape)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit GEV \n",
    "The function gev.fit will return the three parameters of the GEV: [0] Shape, [1] Location, and [2] Scale\n",
    "\n",
    "gev.ppf will return the inverted CDF value given a q value that's between 0 and 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-19T00:35:37.819517Z",
     "start_time": "2018-09-19T00:35:37.653199Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-0.8042402812803544, 3.3576293192513486e-05, 1.533450031574482e-05)\n",
      "48.41873261615124\n"
     ]
    }
   ],
   "source": [
    "# TEST CASE \n",
    "\n",
    "i=800;j=400;\n",
    "rec = maxYearly[i,j]\n",
    "gev_fit = gev.fit(rec)\n",
    "K = gev_fit[0]\n",
    "M = gev_fit[1]\n",
    "S = gev_fit[2]\n",
    "print(gev_fit)\n",
    "\n",
    "retPeriod = 100.0\n",
    "\n",
    "R100test = gev.ppf((1-(1/retPeriod)),K,loc=M,scale=S)\n",
    "\n",
    "recFull = fullYearly[i,j]\n",
    "gev_fitFull = gev.fit(recFull)\n",
    "\n",
    "CDFfull = gev.cdf(R100test,*gev_fitFull)\n",
    "T = 1/(1-CDFfull)\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-20T19:59:38.650035Z",
     "start_time": "2018-09-20T19:50:42.187569Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:49: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:48: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:50: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py:2175: RuntimeWarning: invalid value encountered in subtract\n",
      "  -pex2+logpex2-logex2)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:412: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  return m3 / np.power(m2, 1.5)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6\n",
      "7\n",
      "8\n",
      "9\n",
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n",
      "80\n",
      "81\n",
      "82\n",
      "83\n",
      "84\n",
      "85\n",
      "86\n",
      "87\n",
      "88\n",
      "89\n",
      "90\n",
      "91\n",
      "92\n",
      "93\n",
      "94\n",
      "95\n",
      "96\n",
      "97\n",
      "98\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-16aaa62d33f4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     39\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     40\u001b[0m             \u001b[0mgevFull\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgev\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrecFull\u001b[0m\u001b[0;34m)\u001b[0m   \u001b[0;31m#Fit GEV to get parameters in each case\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m             \u001b[0mgevPhys\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgev\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrecPhys\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     42\u001b[0m             \u001b[0mgevRad\u001b[0m  \u001b[0;34m=\u001b[0m \u001b[0mgev\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrecRad\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, data, *args, **kwds)\u001b[0m\n\u001b[1;32m   2209\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Unknown arguments: %s.\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2210\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2211\u001b[0;31m         \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moptimizer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mravel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdisp\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2212\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mrestore\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2213\u001b[0m             \u001b[0mvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrestore\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvals\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfmin\u001b[0;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback, initial_simplex)\u001b[0m\n\u001b[1;32m    406\u001b[0m             'initial_simplex': initial_simplex}\n\u001b[1;32m    407\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 408\u001b[0;31m     \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcallback\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mopts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    409\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    410\u001b[0m         \u001b[0mretlist\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'x'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'fun'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'nit'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'nfev'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'status'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[0;34m(func, x0, args, callback, maxiter, maxfev, disp, return_all, initial_simplex, xatol, fatol, **unknown_options)\u001b[0m\n\u001b[1;32m    551\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mfxr\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mfsim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    552\u001b[0m             \u001b[0mxe\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mrho\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mchi\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mxbar\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mrho\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mchi\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0msim\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 553\u001b[0;31m             \u001b[0mfxe\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mxe\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    554\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    555\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mfxe\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mfxr\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m    290\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 292\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    293\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    294\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36m_penalized_nnlf\u001b[0;34m(self, theta, x)\u001b[0m\n\u001b[1;32m   2018\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mloc\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mscale\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2019\u001b[0m         \u001b[0mn_log_scale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscale\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2020\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_nnlf_and_penalty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mn_log_scale\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2021\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2022\u001b[0m     \u001b[0;31m# return starting point for fit (shape arguments + loc + scale)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py\u001b[0m in \u001b[0;36m_nnlf_and_penalty\u001b[0;34m(self, x, args)\u001b[0m\n\u001b[1;32m   2000\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mn_bad\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2001\u001b[0m             \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0margsreduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m~\u001b[0m\u001b[0mcond0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2002\u001b[0;31m         \u001b[0mlogpdf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_logpdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2003\u001b[0m         \u001b[0mfinite_logpdf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlogpdf\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2004\u001b[0m         \u001b[0mn_bad\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m~\u001b[0m\u001b[0mfinite_logpdf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py\u001b[0m in \u001b[0;36m_logpdf\u001b[0;34m(self, x, c)\u001b[0m\n\u001b[1;32m   2165\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2166\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_logpdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2167\u001b[0;31m         \u001b[0mcx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_lazywhere\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m&\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2168\u001b[0m         \u001b[0mlogex2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog1p\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mcx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2169\u001b[0m         \u001b[0mlogpex2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_loglogcdf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/_lib/_util.py\u001b[0m in \u001b[0;36m_lazywhere\u001b[0;34m(cond, arrays, f, fillvalue, f2)\u001b[0m\n\u001b[1;32m     51\u001b[0m     \u001b[0mtemp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextract\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcond\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0marr\u001b[0m \u001b[0;32min\u001b[0m \u001b[0marrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     52\u001b[0m     \u001b[0mtcode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmintypecode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchar\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0marrays\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m     \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_valarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marrays\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfillvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtypecode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mtcode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m     \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcond\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mtemp\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     55\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mf2\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/_lib/_util.py\u001b[0m in \u001b[0;36m_valarray\u001b[0;34m(shape, value, typecode)\u001b[0m\n\u001b[1;32m     16\u001b[0m     \"\"\"\n\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 18\u001b[0;31m     \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mbool\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     19\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mtypecode\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m         \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtypecode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/numpy/core/numeric.py\u001b[0m in \u001b[0;36mones\u001b[0;34m(shape, dtype, order)\u001b[0m\n\u001b[1;32m    187\u001b[0m     \"\"\"\n\u001b[1;32m    188\u001b[0m     \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mempty\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 189\u001b[0;31m     \u001b[0mmultiarray\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopyto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcasting\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'unsafe'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    190\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    191\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "#Specify return period \n",
    "retPeriod = 100.0\n",
    "\n",
    "# Define empty arrays to save data into\n",
    "K    = np.empty([len(lon),len(lat)])\n",
    "M    = np.empty([len(lon),len(lat)])\n",
    "S    = np.empty([len(lon),len(lat)])\n",
    "R100 = np.empty([len(lon),len(lat)])\n",
    "\n",
    "Tfull = np.empty([len(lon),len(lat)])\n",
    "Tphys = np.empty([len(lon),len(lat)])\n",
    "Trad = np.empty([len(lon),len(lat)])\n",
    "\n",
    "i=0\n",
    "while i<len(lon):\n",
    "    j=0\n",
    "    while j<len(lat):\n",
    "        \n",
    "        # Compute GEV fit for CTRL and compute discharge of \"retPeriod\"-year flood\n",
    "        rec = maxYearly[i,j]\n",
    "        \n",
    "        # Check for NaN \n",
    "        validCount = np.count_nonzero(~np.isnan(rec))\n",
    "        \n",
    "        if validCount>5:\n",
    "        \n",
    "            gev_fit = gev.fit(rec)\n",
    "        \n",
    "            K[i,j] = gev_fit[0]\n",
    "            M[i,j] = gev_fit[1]\n",
    "            S[i,j] = gev_fit[2]\n",
    "        \n",
    "            R100[i,j]=gev.ppf((1-(1/retPeriod)),K[i,j],loc=M[i,j],scale=S[i,j])\n",
    "        \n",
    "            # Compute return period of that magnitude flood in 4xCO2 experiments \n",
    "            recFull = fullYearly[i,j]\n",
    "            recPhys = physYearly[i,j]\n",
    "            recRad  = radYearly[i,j]\n",
    "        \n",
    "            gevFull = gev.fit(recFull)   #Fit GEV to get parameters in each case \n",
    "            gevPhys = gev.fit(recPhys)\n",
    "            gevRad  = gev.fit(recRad)       \n",
    "        \n",
    "            cdfFull = gev.cdf(R100[i,j],*gevFull)    #Compute CDF at the point of R100 \n",
    "            cdfPhys = gev.cdf(R100[i,j],*gevPhys)\n",
    "            cdfRad  = gev.cdf(R100[i,j],*gevRad)\n",
    "        \n",
    "            Tfull[i,j] = 1/(1-cdfFull)    #Return Period \n",
    "            Tphys[i,j] = 1/(1-cdfPhys)\n",
    "            Trad[i,j] = 1/(1-cdfRad)\n",
    "            \n",
    "        j = j+1\n",
    "    i = i+1 \n",
    "            \n",
    "    print(i)  \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-20T19:59:42.739660Z",
     "start_time": "2018-09-20T19:59:42.716242Z"
    }
   },
   "outputs": [],
   "source": [
    "def computeGEV(ctrl,full,phys,rad,i):\n",
    "    \n",
    "    if np.remainder(i,100)==0:\n",
    "        print(i)\n",
    "    \n",
    "    # Define empty arrays to save data into\n",
    "    K    = np.empty([len(lat)])\n",
    "    M    = np.empty([len(lat)])\n",
    "    S    = np.empty([len(lat)])\n",
    "    R100 = np.empty([len(lat)])\n",
    "\n",
    "    Tfull = np.empty([len(lat)])\n",
    "    Tphys = np.empty([len(lat)])\n",
    "    Trad = np.empty([len(lat)])\n",
    "\n",
    "    j=0\n",
    "    while j<len(lat):\n",
    "        \n",
    "        # Compute GEV fit for CTRL and compute discharge of \"retPeriod\"-year flood\n",
    "        rec = ctrl[j,:]\n",
    "        \n",
    "        # Check for NaN \n",
    "        validCount = np.count_nonzero(~np.isnan(rec))\n",
    "        \n",
    "        if validCount>5:\n",
    "        \n",
    "            gev_fit = gev.fit(rec)\n",
    "        \n",
    "            K[j] = gev_fit[0]\n",
    "            M[j] = gev_fit[1]\n",
    "            S[j] = gev_fit[2]\n",
    "        \n",
    "            R100[j]=gev.ppf((1-(1/retPeriod)),K[j],loc=M[j],scale=S[j])\n",
    "        \n",
    "            # Compute return period of that magnitude flood in 4xCO2 experiments \n",
    "            recFull = full[j]\n",
    "            recPhys = phys[j]\n",
    "            recRad  = rad[j]\n",
    "        \n",
    "            gevFull = gev.fit(recFull)   #Fit GEV to get parameters in each case \n",
    "            gevPhys = gev.fit(recPhys)\n",
    "            gevRad  = gev.fit(recRad)       \n",
    "        \n",
    "            cdfFull = gev.cdf(R100[j],*gevFull)    #Compute CDF at the point of R100 \n",
    "            cdfPhys = gev.cdf(R100[j],*gevPhys)\n",
    "            cdfRad  = gev.cdf(R100[j],*gevRad)\n",
    "        \n",
    "            Tfull[j] = 1/(1-cdfFull)    #Return Period \n",
    "            Tphys[j] = 1/(1-cdfPhys)\n",
    "            Trad[j] = 1/(1-cdfRad)\n",
    "            \n",
    "        j = j+1\n",
    "    \n",
    "    Tfull_ret = Tfull\n",
    "    Tphys_ret = Tphys\n",
    "    Trad_ret  = Trad\n",
    "    K_ret     = K\n",
    "    M_ret     = M\n",
    "    S_ret     = S\n",
    "    R100_ret  = R100\n",
    "    \n",
    "    return Tfull_ret,Tphys_ret,Trad_ret,K_ret,M_ret,S_ret,R100_ret\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-20T22:20:08.912356Z",
     "start_time": "2018-09-20T20:01:11.787102Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py:2175: RuntimeWarning: invalid value encountered in subtract\n",
      "  -pex2+logpex2-logex2)\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:49: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:49: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:48: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:50: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/scipy/stats/_continuous_distns.py:2175: RuntimeWarning: invalid value encountered in subtract\n",
      "  -pex2+logpex2-logex2)\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:50: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:412: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  return m3 / np.power(m2, 1.5)\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/scipy/stats/_distn_infrastructure.py:412: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  return m3 / np.power(m2, 1.5)\n",
      "/Users/meganfowler/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:48: RuntimeWarning: divide by zero encountered in double_scalars\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "200\n",
      "300\n",
      "400\n",
      "500\n",
      "600\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-9-bc90ea2209ce>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mretPeriod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m100.0\u001b[0m  \u001b[0;31m#Return period\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mTbig\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mParallel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_jobs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdelayed\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcomputeGEV\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmaxYearly\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfullYearly\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mphysYearly\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mradYearly\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlon\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/parallel.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, iterable)\u001b[0m\n\u001b[1;32m    787\u001b[0m                 \u001b[0;31m# consumption.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    788\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_iterating\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 789\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mretrieve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    790\u001b[0m             \u001b[0;31m# Make sure that we get a last message telling us we are done\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    791\u001b[0m             \u001b[0melapsed_time\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_start_time\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/joblib/parallel.py\u001b[0m in \u001b[0;36mretrieve\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    697\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    698\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'supports_timeout'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 699\u001b[0;31m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    700\u001b[0m                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    701\u001b[0m                     \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mextend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/multiprocessing/pool.py\u001b[0m in \u001b[0;36mget\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    636\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    637\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 638\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    639\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mready\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    640\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mTimeoutError\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/multiprocessing/pool.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    633\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    634\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 635\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_event\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    636\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    637\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtimeout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    549\u001b[0m             \u001b[0msignaled\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_flag\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    550\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0msignaled\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 551\u001b[0;31m                 \u001b[0msignaled\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cond\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwait\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtimeout\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    552\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0msignaled\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    553\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/threading.py\u001b[0m in \u001b[0;36mwait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    293\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m    \u001b[0;31m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    294\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mtimeout\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 295\u001b[0;31m                 \u001b[0mwaiter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    296\u001b[0m                 \u001b[0mgotit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    297\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from joblib import Parallel, delayed\n",
    "\n",
    "retPeriod = 100.0  #Return period\n",
    "Tbig = Parallel(n_jobs=2)(delayed(computeGEV)(maxYearly[i,:,:],fullYearly[i,:,:],physYearly[i,:,:],radYearly[i,:,:],i) for i in range(len(lon)))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import pickle \n",
    "\n",
    "# with open('GEVreturnPeriod_100yr.pkl','wb') as f:\n",
    "#     pickle.dump([R100, Tfull, Tphys, Trad],f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To read in data again: \n",
    "\n",
    "# with open('GEVreturnPeriod_100yr.pkl','rb') as f: \n",
    "#    R100, Tfull, Tphys, Trad = pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save return periods in .mat files for plotting there"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "sio.savemat('pyTfull.mat', mdict={'Tfull': Tfull})\n",
    "\n",
    "sio.savemat('pyTrad.mat', mdict={'Trad': Trad})\n",
    "\n",
    "sio.savemat('pyTphys.mat', mdict={'Tphys': Tphys})\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "fullData = sio.loadmat('pyTfull.mat')\n",
    "Tfull    = fullData['Tfull']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = np.full([len(lon),len(lat)],np.nan)\n",
    "#i1 = np.where((Tfull>=1) & (Tfull<=5))\n",
    "\n",
    "colors[np.where((Tfull>=2) & (Tfull<=5))]     = 1\n",
    "colors[np.where((Tfull>5) & (Tfull<=25))]     = 2\n",
    "colors[np.where((Tfull>25) & (Tfull<=50))]    = 3\n",
    "colors[np.where((Tfull>50) & (Tfull<=75))]    = 4\n",
    "colors[np.where((Tfull>75) & (Tfull<=95))]    = 5\n",
    "colors[np.where((Tfull>95) & (Tfull<=105))]   = 6\n",
    "colors[np.where((Tfull>105) & (Tfull<=125))]  = 7\n",
    "colors[np.where((Tfull>125) & (Tfull<=250))]  = 8\n",
    "colors[np.where((Tfull>250) & (Tfull<=500))]  = 9\n",
    "colors[np.where((Tfull>500) & (Tfull<=1000))] = 10\n",
    "colors[np.where(Tfull>1000)]                  = 11\n",
    "\n",
    "#lonUS  = np.squeeze(np.where((lon >= -99) & (lon <= -82)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<cartopy.mpl.feature_artist.FeatureArtist at 0x1a1ff601d0>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAADqCAYAAAC7kx6uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJztnX90VVd237/HNhKZ8mMmRilIQqhDI5HMrHRCtKYJjqZYL4TJxCWxS9YaIAuTBC9lEbsEV9Txsj2UMl7URYUQT9yoOC1Wa9FptXBC3EwJkUxCh8wkhJlJ7AwiQ5YkBCKDnRUDSZA89u0f9+379t3vnHvv+6X37tP+rKUl6b3747z7zv3effbZex/jeR4URVGU9HFPtRugKIqiFIcKuKIoSkpRAVcURUkpKuCKoigpRQVcURQlpaiAK4qipBQVcEVRlJSiAq4oipJSVMAVRVFSyn2l7Nze3u5NTEyUqy2KoijzBs/zTKnHMKWk0htjPE3FVxRFKQxjTFkEXF0oiqIoKUUFXFEUJaWogCuKoqQUFXBFUZSUUlIUiqIoCmf9SF/w94bll7B5yUUM31qLZz72ehVbVb+ogCtKClk/0ofxqSYAQHvrTes2ZzP9FTn38289hDM31ljfG+gYAgD0Xt6KzUsuovfyVgx0DOH5tx7C5iUXg/cGOobQe3krps+3YOy5PRVp53xABVxRUsDY1eaQ8AGVF+7n33oo+JvE1/879//wrbUhYR6+tRZnbqzBhuWXAOQEnQSfvwcAK9ZdK0tb5ysaB55l7Gpz7DadK68nOg7fjg8p46Abj/YZn2rC+KNPBe9nHjyIkTeeTnw8pT5YP9KHgY6hQByJSlnYAEIWc1LoAZN0u97LWyv6GWqZcsWBq4BnKbTDkhXE4Z1Xvk/+QHo9blsgZ7XQtnQDcytsvt4A8wHZJwv5zteP9IX60dlMf2AY0P7tr7wAADjdczRRe7i1TX/z1+KQ/VzeA9z9It+3tSWpX73zwBGsWHetpu4VFfAyQp2d++x4Z5MdiVtCtk7G96Xj2l4H8m8u+V5cJyZ0kqi+IKubKESw+L4ksL2Xt2LD8kuBC8MlulFWNB2jGMscsN8rvH1JrHdJ58rryDx4EFd23BMarRJRI2v+EKL7cK6oGQH/F7//bwCUZgnKzsrpvbwVd04048LLTxZ9fBfc3UFftM2yBsJCbLNEZJuBfJcIHYeLPw0jpauF3qOhs+umIRHnHZGoJYujXrC5xCrhDug8cAQzbbOBn3t8qgmNkw2xE35JXIEukli+pRyXKOfx5UMlyYig0IdFOV0+1H/+4Mf+U20IeLH70wfhQhWF3GbT8b2YaZvF6Z6jwZcIwPokpXNJ8SR2bX8C765uxKIt10NtkgJM50higUjr2HUM+fTnAsGHlHK4ySeCZJv4jTif/YylIr8Luq42wZAiVew1b3/lBaslKdvF+8b4VBPaW2/ibKa/KN+1FD0uWFFWM/87qfhzgyOuLbZ+z88hjxV37KhzFQLvBzIah76LKM5m+mvHAv/8mz+ZeGjX/soLif1tEj4EBPL9wsUM7Wj/pMPJOEuYv0ediYu3TYRtoVV0I/IZexJ5m4Vlu3m4Na4CnhzbA9T1vUu/bZR42SbApU+au/I2ju4GgECYZR+i/mUzSIpxR8jPQveay0VoM0SkT1yODuU15EZOUguai2cU5bbCo/a3fXbXSJ5G1M9+/P/UhoBfmlwBINkECwl4ki9JdqZjF7qdT7Y4K6FSFPrEB9x+RN6ZCflA4A9KKeS2DsX3VTFPhk3AgXjfcFIxsI3sAOSJtJyD4e60uLkRF4XeE3GfK8nxbCLN93W9H9UmwD4alg+MqGPa7pFCJmPj5r4k/IHYe3lr7bhQSMA5vJNJy6NQd4RtOGcbyvJt+fnLIeRxnVBuW8r5uMBLa4YLOI/R5Te0vAY2NHkin6h5GCCZZZjEHWibHHc91OXIKmpivRLQZ950fC+2PzLiFOBC7gGboQIULt62c9tG4/w+GJ9qKsoDwL8f6Z7k38PG0d2BkRn1fQI15AO3CTiQ7xLgw8M4KyLqZrAJnM0P6RL3KIqxqF0TmqUIeRJ/qs0tQ+eVgi5RKzwMj7Muxipz9a3OldcLnutxHZ+OUWmkQJHlOH2+BSvWXcNAxxA2Hd8b/F2MEVPsfXbsQjce6zoXa2G7wnJt29sMQ8D+cLXpUtJ7nru8Niy/VDsuFNq//ZUXgosLhIeKBA37bRZt1JMtqfjanpDT51uw/ZGR4FyuiSggPmwvjiQdN6kg2CZI5NDbdo1sn4V3TBXvMHzSL+4GfOZjr5cU5VEK3JJ0+Z+B0v3f8pwkltyF6ZpATzryLdeoOCoO3fXQS+pj526sQkZASTTgzI01tWOBu/Ynq2bj6O688Kf1I33BEx1IblkkHaZK90EhvstSJzfiSNKBXLHjQNhXyjM1XdEr9L9Go4SRqelxbo+5sH6LOadrMrxSyBwI10R/Na4ZkXTEsnF0dxDFxh86Ngs+rvYLkTSMseYF3IUMgSokyiOJKMt9yHJwbb/p+F6c2nGokI9QdqImgLkwU4gS/9042RA8CF0TYSrgPjw7MelIiK5ZtSxv3hZXpFW5XHdRx3Ud2/UQLFXAi/lMUdcoyXFc78uHVtT2dWOBu+A3QqGdjdwDgyczkT44IHpYyN0LNPFQaGeLG8IlPUYhGZQ03HdN4tqiTWi0A2DeT1wWIsJxUT/lgvejuHkhTpT16/q/FCisMckkoBRSaZVHfdZi/OM2bPNjxbo96DMBhc1l8ElVeZ/XTBx40v2p7kIxYsmhBB4AaJxsSGw922JX6fVKDfWiOkYpFrFrApOz6fheAL5odx44Evw9nykm0oTH5Vc6LLUSJLESk/jQaTRLAl5IOKPrugL5kSHcogeSu1dLnQewIe/Rrp2HcXuVSaQ5ccZZ6gTcRpKbgn+RcdUA46wk29N9LuLFbdAEazGiSkLkGr7SzXU204+unYfxdvd7wXtxWX5pJi6Lkfu9k2KL207qmy7Uhz0XfuO4CC9OMUKa5PwS6XqRPuckD91CQyyTzKPJsF2pE5uO7y3aKKoLAXclTXCKtVQpaShqVrlcvsNiw6JonyRlajk8lj7OJwrkJmuA0iz/WoQqzdGEeZSAywfftv19eHjPqNU65CJyNtMfnAfIZdPxxAxXYkcSN4HtAUzvl0PU5fHoM8T5jF3HihpV0rGj2sDhsdO8Xa7rybezjQLkd8ffK/Q+A9wP7VLEG6gTAefICmrlKMrE3TZEKfG4hVDsOZJ2Ml6WIGpyyfZeqZ2vlrCVXE3SX3jxMtv3tHF0dygsthIk8Qm7LMukETQSHsJbTLJbISG9hSCFHMgXTdsDwhVBwg22YvWja+dhLNpy3XmNinkgEHUh4LKyG3f2rx/pK4ulSP5faT3JGPBa8m9GdTyqoOhyF7ksLPm6XCwirdiqSBYi4LbRCtUfkdsmJcnwXFrfNtcCMVdGB2+PbCchsx3lvVTIQyHqM1Eflf5xAIlcTnwOqFjaX3kB7a03naOUeSvgUrjnYkhPExA03HbhEvRq+cqJbfv7cOHlJ53zBlGTVXGhXqV0xGrBRZVnBia1uqQoc+G21ShJ4iKR34FrbkIKdlScMVGolV1oX+Wuh0LcKvwzcoElCgnpkw82TlRxLYJchcO31mLwZAaLJ7ySylB3HjiCUzsORVrfXTsPA0DB5ymXgN9T6gGS0HngSPADAHdONONspj/4mSsWT3iJzmfrdJUQ7+Fba0ND2qjtqNStqx2847seQi5Lb+xqc/DdpAEusMO31mKmbRZnM/2hUQavFcPpPHAkcK0RvZe3BsN3WbudW/h0fP6d0eskzsO31gb/28Rmw/JLof3jxHugY8jqSuDfo+xDPNWc72Pbl7eDzmdzx7n2k22l85+5sSb4XLb9XJ/Z5UeXgQe9l7cGoY30GpALd5xpm8Xb3e9h7Goz1o/0ofPAEWefcEHRJvLc3OC58PKTFVmrIClzYoHXShgbuR+iUtGrQSGJCq6SsvQ+EfXZbBZ5miY3ZTKYbUEMwpaIwyd16RhEnE+a72MbysdNSNpwZd6Wqz+WEh9eTIw0/S3nsmyv20Yo/FjyGkgXDp2Pch1IdF3fSSGjTdt9Vq7Raqos8LHn9lRdvIHcxbdZ/rabjsOf9vwm4xaYC7K0bdsVcjNJy5IfK4klT2xYfqmqD6xyQX7Srp2HcedEMzYsv2T9TPJGtMUdc0t385KLsdeGW6p8W7JkqV+QJRplvUZZqnFWrw2bRR71v+1c9De35l1WvetabV5yMbJuNxd8Ovb4VFPof7pv+PnGp5qweUmuZvpAxxBWrLsWEm95nGKuIxfrzpXXa9LVWDNRKNVA3tiF+g6ltZTUAk567CQTYYWmB9v8tERaLHCJbSKTyg3w74VKD5zacSiwoAH7ogj8eDY/NwmDy2csI0xc9UpsfUZatElHAq6wOlvftJ3D5neWnzXOyo5rq+38NIcxfb4Fiye8wF04PtUURAJRG3npi97Lfi6FTKzhIyxueLW33iwo94L6VSWEO9WTmLVC5sGDeGnwxZKOUUjKLsfl1ih0vyTb2xZNpvhzm+shbUI+drU5iAF35RbwVHDpQrEh3UxRQivPFXVM1wSh7Xu1RRHx80VFsLhioXm7XZ/DJbRxkTK2z+Q6j9yn97K/9i3n9iqDmbbZUL2fqPMBuclo/gCngnqndhwqKIQ2yfJ2xaICXkZotjkpLhFNIq7FWuNE0n1tlrasicJrtMvzpUnEKduUJiJLcQ1FhZXa3pMim6QOfVwyTFw8uO0ccbjCSjm2UQKv2MePw48L5Ff2I3iEEImraynBgQ4/uQrwxRtAUDaDizKdl0IFeWSZXHyCxJvXMq8Fd26qfOC1jvxC424I8q/RtvQTFdPKfXmF+OVoW1s0govey1utAsHnIlyTfmmk2CgAWxTQ5iUXI8M0bWnbfHvuTuDfMz9PMRFN3E9PIz7qUwMdQ9iw/FLI38z/pzkPOjdvS9zDwraMoa0vtrfe9EdB7N4A/ElFcrXQqOfMjTXYOLo7z7cNICj7sHjCw6kdh7Ds3AI81nUuZFHfOdEcPBhk+Ch9Rg71+d7LW2tCvMuJWuBZuBUeN8S2ZbPJ910+0ahjljNUMUnKM7eyeY1xIL01UzoPHMFM22xs0bRCijHZcLk8kvjR6fxRBZ6i2l/IvIYkSZQNWepxESFRnyXqc9FvHtvNXR/cZWJzM3Krnv52lQkm/zVlLlM+RbVRF0qZkan8hUwM2nD5/+ImeYr1r7qOY4MXu6pHeH0U23UsxoVB70UluXSuvB5k75FQDp7M5Lnn+ASqS3DGrjY721mKG04KIlFO40L2YxJscuHJZDq+jW3dSnk8+UCR80+2/m0LZa5meLMKeJmRN0ySiS5bSnTS7aP8qbZZ/UIEXfo6+THqVbQ5XESBcIbg4MlM4COVkShysreQ+QZeuY6OQdealhsEwglXsmY2HcMV51/uURrNi5C160qP5xEuQPRyYrLvkp8aQDAyksfYtr8v8Hnzhx3dI3JRFmoT3aP8XnWl0PN1TwdPZqruSlEBLzPypuHDtDi4QNgSMmS4GsHD2uJCwOKGr7xTLzu3IBgm0tCxFmNYi4VnjZ7acSgQAC4QHHIL8cksuh5yxSOKZImKCOEUIggUw3/mxhqrpcmtbxdRrhfbauwkWkC+sWDDNsnJjQqbyNO14etmEtzqJriLaPhWuM44wcP/GicbMNM2G9rGdh1syzfSw5zaDkAFnDWibgTcdtO6LJ4on6McOnLRkNvy7Qk+m88Fnp/LJep07GLrjKcJqu3NsT0oSSzkyIOEmwSGwtUWT/j9+cLLT4ZKErtcD4B9LgGwL2XXtfMwXt0XTvun/iT9s1FC7gohtLUzyiVkI+qz8nMmqeFiixyhOHwZwUJw1xP3cdvyHiiqpGvnYatvm3/PdN5aGIWqgJcZebPE+TtpmygLXZYklTcd3zfJxKcr1pzEnoaTaZ2ALASZvMMjMEhUKK446sbmbg7ppwXcJVxltEVUKVTXw1ZOVsrhf9Tyg/JzSwvZRlSxqCRzMzJix1aYC4D12vDPKAtEbdvfF5RtJVxr1fLRRJIJSQox5aPSWkAFvMzYrJ0oASWirBV6n/sGbW4WPgsvt4mCt6+eXCRxcB8392MDYeubFmKYaZvNe6hJAefHBvIFmVuAz3zs9bxSo0ljwAnu4pATcuSPt5W8lXWzkyzSbOu3HOp/QLI1L/nnozbw45OA2yYl+Rq0NkOE2rLs3ILQSIW3FUDiaBK+ClMt3SMaB15m6MvlsbuyjgOfxKG6G7ZYVk7v5a2h/Xjdi4GOIUyfb8kbQh670O1sJ51n0/G9GJ9qCmLQ5xNSvGX1O87Yc3vQONmQ+Njjjz4VCM/ZTH+eNUnnIBHldU7obx67b0Na0+RaoL5CDxXqk3zb0z1H82LC4x728n2+P+/j/MHA47QBt3EiJ8tpVEHtpFEhXRd6HfDdGjwWf6DDT+Qha3nbfn9B7rOZftz7+fuDdgzfWoulV2YiPzP/HMWGiqaB+6rdgFqCrJAzrf5NeubGmrwbDQhbVrbOIWsXU+cfvrU2FLsKwJ4i3OFuIx3z1I5DQXtrwac3l/jRIv7fLh8sF17bfACtF4qMPamJLHaa1OQPCN/6tsdUxxUVi4okoWOsH8l9p7bJQsCdDSqPL40Dl5iRr3/b/j68uq8/+J9b1b3Y6pyop+PaQgBtVn3v5a1Yse5a3rV4dV8/No7uDlwf9PB9d3UjFiF3zd9d3Wj9HIBfIgMAruy4B42TDWWP3qklVMAZ448+hc4DRzA92RL46QiZMADYLSk+rJQz+ccudANdufRgDs2gb5yKD18k6i26pBBkMgiHXCdR2abjU01YhgXoPHAEY8/1Bwt+2NwtPLwwHFvcgu2PjAT/c5eYbbJZjuq4mPL+NtAxFDygNizPtSOJa02+z4s6ueKr6b3p8y3Aqvz3yKLeOLobm9n29PlcRdVs/m4Z3UVRKHyi83TPUaAndwwgN49B14UKXtmY6lkIAGic9P8fPJnBM8+5V4hPM+oDt8BDj8anmrDs3IJQnGqcy4JHFnB4yBSfeOy97BfyeXVff2jbJFbDfBVwgoT1zonmyJtawqMjChnB8LorfLV6Wg5QPjSSDN9tMdau/sOR2ZjyGNLHzOuE2Cx14nTPUav4ElERJIV8RtfnpOPaJuPbX3khsMpbR+9i5I2n0bXzMB7eMxq6JrzCIVD9tQgk5fKBqwVuYfzRp4LJj9M9R7HtXF/Iaoi7KalTb9vfFwwDT+045Psre/wOtmH5JZwBQu4Umo1Pcg4gm5izsoQPWidsWH4JZ7bkL+UlE0s4tnjxpDRONmB6sgXI+P9HPQCSTk7L79u2nS25BghHoJAf3lYbnUd18OsTWL1AKPGGRNxVqMqGLUOSzrd5yUVsOr4Xw8tz57U9hHLtiY6mIvFetOV63vVaPOHhzkRzMKKoV9QCj4HWvKMZcVusuKuuRFytZO5f5BXfeB1kF2p5x2dMRgk4p1Jx82SZA/m+ehl1YnsvKTJRh4coJj2mDJF01SPhxweio1ZsUVyUxATkEqv4g43XRLFZ4M+/9VAoY1VmLw90DGHX9icw8sbTkZ+32mgUyhxBvjdpWdDvzUvCq7fQLPnG0d3OOhYERaGQeNNwL8o63Di6e16Id9QanaVWUhyfasqLLqG1MuV6mfx18qvLbWzwbbho2eBRIUB+lUT631Y9kYs3f5/cDONTTdZYbcIWpne656g/F8DazC1z2iZKvLlLRkL9e/GEF6yuQ+egkgONkw3WPsBLFlBQAJBbZ7dz5fWaF+9yohZ4AsjvCeQsDjlBaUthpu1s2Zay2hw/ni3NmpgP4h2HS8Bd8eAEL6oko39s21JmJn33SdwuUWLJsa0kw+GRJbbjTZ9vifRn0wiO3CCuxBggP2nGlmSU1Ndt+xxxyWg8g5inzMfV7qbs6TTeE5rIM8dwiyrqRrZNYMroBCBcnwPIVWizZaXRMeZbuKCLQi1wWQ3PVTPFtV8lsbk2ZIKMzc3CBVzuJx9UxcRB84eAzCiOgz9weAq7rcAagJAvX77OLe56QgV8juGRKUD0qiguPzmP4d10fG+obrUtW2z9SF/IGq/nDh3F+pG+yJXnk0DiDeTH3hdKUmHn27n24VmLFMliK8AE5EYVtHYkRUZxtj8yEiSCxdVEd2ELfXRVIpQlc20VOinKiiJFbCV2XaTRuk6CCngVkBNntlKbhC2Zh5J4PvrSt/Cl6V8v6rzzyQovx6pBMq2b/i7HMaOOxc9FlrDcT7pyuNi5QiK5a+f2KhMIIV+KzLbI7+rjH+Dd1Y3W94uF+j/5x+VxubuEwnEB5IX8uY5bz6NODSOsIhuWX8oLkSJIyG2RBRtHd2P8uT3Ac4WfL6qWRT1TiuUtrV4u4qUSpOe3Rr8/PdmS6HUg95DvfDm/1C3B/dAk2gACtxufbOT77Dr+BADfQufhqlHzLXFQv6eSvtJ1yLNXGycbcHuVP3HJR6P0ualWOz3ozrTW74Ij5UQt8CKgus7PfOx1PP/WQ6GYXAB5ljcQHnoWMyyk8/Ca0fU6vOSUIt5AeNK5HBY4t4AB3+K1uWQoNJFWB6IJOoJqXFNb6H2Kwjh2oTsURrd+pC/kBuKTfbZJcQmfhOSrtPO/k0DJQIMn/SB42o8vMMyvweIJL3hY0GfkbsPp8y1oHb2L9599JzQvFDV56WfP1lZiTqGoBV5lNi+5GMSkykqGZKFzK5w6a7Edj6yZsyspxG4vxgq05OcD0sLm/lxpNce5QmREyfT5FjTCF8/FEwtC5+S+bhLnzgNHsHjCwzRy4aF0ThIy2p7aNziZySu1wB9iY8/tQebBg34tkJ6scLbN+uJuqaFDVjF/kLW33gxGAO2tN4PrwbMXKfqGW+e7tj+Bxas9jD7+ABavzgnz8K21eQ+x8akmLMs+5O6caA789Y2TDdjQlbO+V6y7him0YAX8h+3wrbU4hu7I5c7SLt7lRC3wIiGLGPAtar5kE2ELZSvHsJCSi2qpvnElKMekpbR6CRKcQv3Z8m++Dz8WnYuscGoLvS6tb3JDAH7SGFm01H4eRUPHdUEJRLZUfKqjDSAQ1Zm22cA/TaJMort4wsPDe0Yx+vgDmOpZGIw8CHLd0DGpzAAdg6A6MwAiI1qoTCyNXOrVjaKTmDUAd6UACNVv5stmcSuukA7JoxF4+VR+TgDO1UjSTikCLmuASCGXAi5dGrb3AORZ2uSz5WIs4SIO5ESZL6d37EJ3UAOb2szhdT3oQRDnZuCZlFRlEHAvMHHnRHNQ+oG7feTDY8W6a7hzojko6UqTo3LZtGXnFmDRluvBcQGE6grJolZArl/T6KUe+zWgAl6z8JViZEx4VBigLKBFZUwJGdUChIvaU5GfehlelhqBIgWc4OLJ/wYQiLHtPW49c3iiD+3Ho0RoGyBnhZKrgCbteFtlEg4JqYxiSSLg5J+2LTYhkaGDFPZHFvvSKzOY6llobQd9PrK+ueXNo2nINUPXAnCPJDIPHqzrjEoV8BpFrq1pW3XEZoVLwYqy1F0LMNfTcFPG3RebVMN9zDa4UJOY0nwF9wdzFwC5G2SGputcUsD5PjIjk0822rIVaZ8kD2pytQG+W8ZV+S8qW1j2qeffeii0SjyJ/dIrM3j/2XeChRcAv6xr6+jd4DcVnwLq3/0Xh05ipgCeyMChxBTCH+4mPy753OVwXS70m2ZBH3/0qSD7tXGyAY0AxtFk9T0DsAodkPVVI5z1yv242x/JWbu0aAGQE/aZtlnMtGXbIKzNUHtFfDbhH6chz41Cbd50fm8olV1+Lp7wRSQdZYUXSO7PW6AEyE2uywnQ1cc/wMgb+f3HH0Hm5n/u/fz9WIoZvLu6Ea92DGHb6r7AtdI6ejf4PfLG034/33K9bKGcigp42Tmb6Uf7K+E6ErYU4rGruQnPFevWFFSfmqx1srrpf9tqMCT0aXStBEW+RMw0LS6QNH6ZRJ9EefHEAiy94ovO4MkMBrN1YalM7Ip114LvIhDOVt/yJStaJuLQw4HcCDwm+9SOQ+g8gGD1dAAY7849LPiiEONTTdg4tTvk1uEUK3586T/KSaBRB3e1JKnmN3gyg807fAG/suOewNe9cXQ3lrHtXhp8kf33tF+8qi3ZSIpGYPd+/v66dqWUirpQKsjY1WZrarGrSD6RVMDltrZVY9LsG6fPwx9u0l/sEnEef89jlGl/ipogPy9Zutz//VjXuchytLbMSoLHR5NVTt8Brf7D/ck8woVcNLaJQVuJ1STwqCkbNPGe9GHP3Xi8+BW5hGQfT7oIMW/vmRtrcOdEc126W9QHniJIiGRdC+kb5zUn+Lblrn+SlkQIvqK4jTgr3LWsGeAuIiUXI+ALbkikK4cEnKI0KOpj2/6+ICSPxMhVLpeEn8qxFrNqkA05b2KjEJF1zcNwnz6//huWX8JrR3oKOj49GIp9aNUyWg88RZzN9ONspj8kxGcz/aHFjsnNQiup0Aomm5dcTFR/ut6wCY50EUnxpvdJ/JKuoETYVpI5teMQBjqGrK4LKd4r1l0L1dEmy//t7vcC3zktuHtqx6FA6E7tOITFE15QtW/41losO7cA448+FfSdcsEXKh6+tTY4X+fK6wVZup0rrwc/dN3Hp5pCE7L8up+5scZafCuKgY6huhTvcqIWeJXhq7YQtgSMclrhabDAyf1E0IhFrsDuWs1GWufFrkxOYsdjmV0+XF46lf4nNwl3i7SO3kXPF74cfIbRxx8IHadSPl/bQ5F84aX0Bx55FbUWLLlYkkbQUJIQXaN68oVrFEqdQNYV+fykmAPld6HUungD+aFtJHa2eGabhc3FWy5YwBN7kpRcHegYAvb5f28c3W1drmzj6G6cltZnRy7VfdGW6yEf+ejjDwBfCJ/HJlBksVdCvDYd34vT2XooSbGluN850RxcH1ukC18LdtP5/GJbLqbPt6AXfj94bXVP4v3iqKfEN7XAaxA54ZSkaJUt8SXNYYRAfmarLM0btXAG307Ciz9Frdtom5TjixTLVWtk0hb9TT5w7ht2RVeQYE/1LAzqngDlE3A5+UhtT1oYzTVzHqOEAAAfKklEQVR6I5/12Ux/kKEMhEeTdG2T9Eu+Fi1dPyAdxkcSdBKzzuHiFWWBSxeMrS454FvxfDI1LQtDxIkBX4hXRo3Ywje5NS5XSqLSB1Rx77UjPaHFrPmEpoyGWX38AwAydM7tutm1/Yngb0pFl2nrlRArPjFciKDGsX6kLxQxYnPXJC3oxid4aTm4JOUD0oQKeJ3DhYtussyDB/H+s+8A8K1r10SfbfUU7k/mZWlrnfZXXnBGmnARH761Fq8d6QlV0KNJxMUTXkiIObxGyK7tTwQCzF8H8n3qJPJyBGCz+HlEC01mUrIL4As4Pxcdv1JiZes3pfaH9SN9eaMKqsdCtViSuC14+d2x5/bktTUt/TYOFfB5gq0eCmFL2bcl8xDFDmdrgSRhcEQpE5hSuG2Q9WwT3iSQC0Va65xKCpXrWhZ7zsyDB/HS4IvW0WIh/mabf51S93kt83qITFEBn0fI0rXcvcJ9roB7maqo1cHTIOI09OeTmFKYuQVNE3N8PVHa1rYAQlR8uGTX9ifQ84UvA0BopaS4faMsdHowj081ZdPYKxtxwaM8iFIeGtQnyzlyoAcDRwVcHEcFPD1QhTZbx45CirsU/rRY4oVY4RKeLQgkjw0vF+TiIYudhw9KS34u3QR0TcshjKWIOF/Ziqx4PqG7Yt21YASahr4ahybyzEPIKouyzvg6icTmJRcjU8IB31XDf2qRzpXXE68LunF0N7bt78Ou7U+gc6Vfk5qs8enzLUESSyWhdSIB/zvgtb43L7mIlwZfzHsQ88nNuYAScShOvRDGrjaj88CRvKzSmbbZgpLPeH/jiWsjbzyNkTeeRuvoXUyfb8HgyUxdiHc5UQs8pciytQT5f12+XDmZCYTrsPAQvFq9WaiWCE1O0gQZReS4rDQZnlmI2wQIu51oPzmaOXahG0CuEBclAEnLnxYW5hOwQLrKrPIREbmleN16OZKwhSDa3IGVnMCtFdSFMs8pxZ0gcQlZPd9IZDUmXczXBoXF8WPwkEMeP06hcHw5Mb5NIXVCagkp4hTKaQtV5QJuK65Fi1zUa5/jqIDPY4oVb1vii6xXzmOn01qGNil81Rob/LrIFe1dyUK2Alu21/iEa7lLJcw1sj9GVdbsPHAkL/xy2/6+0CglbuRXyezUuUJT6ecpxYh3VGih9APTjUWp4GNXD9VN7K1Erv5uIxBntuABFYDi4W2FQuIdV3u71vHF9AlM9SwMRhO2RUwAv++e2pF/jIf3jIYybTsP5Oqs83VJyZh4d3VjRT5LGtFJzBRRzORi7+WtOHNjDY5d6LZOZA50DMX6gGt1UrMY+IRbEvHl1fvof8B/0Nn2l69t298XWiIt1JaV11Mt3kDOCl484YWqHNoqOEYZEiT6vZe3hkZ97a0384qH3V5lCq5sWK+oCyVFkJCWKwQuqtY2EHa58EShNA/3CVntEHDH0MuY8iioup90zdiSi+rhOhI8Szgqttw1guThrOQf51UdSdRlsg9NXFOdGlrJqNZdf+oDn4dEle2kCaBCJ+WKKbNKYp5m14os6gS4F3lwvQeEY+yjFpjYtr8vWP0HqJ+UcIJqodhKAgA5QY1zAfJENYCt2ZmF0utpXVgg/N3wuQUAzgliW9bnXKICPg+xTbq5KuBJ69plbRcaSueilsMOXUTVWSkE+RDka2HKDNBl5xakMtokCTK0ldcBv3PCF+WlV2ZwZcc9eded9x8ScLkiD6+TIgt/JTVcohKW5lLUNZFnntF54Ahm2maDCTTClV0p64pXQryHb60N2jLQMZQ6X/n4o0/l+aVJfIGcdR3luwWQF9lDq+3IYzdONtSteAN+9AglW9HcC4nshZefxIWXn8TIG09j/NGnsHF0d7Cij7xO5LKizFmZKDT23B6MP/pUUVUbT/ccDRWK48y0zTqXz6tV1AJPCdSJZ9pmrfU9AH/4+P6z78x5mrh8EKQtfnzsanMg3DLEjYibL+DwqoPcpZC261IK3DcNFOYyyjx4MCixSxRy3aLcNLIGetfOw8F3RP34bKa/4qtWqQU+z2gdvYvFEx4aJxuw6fhebBzdjTM31gTp2ly8yUqPsxxLgVuqQK7u+PCttdj+yEiqrHFKX48KK4wTb369SbxlpAS9Ph84m+nH6OMPBNdk7GpzIKyZBw8GCzbYGHnjaVx4+Unf3539SQIZOXKtTg4PcXz+rYdC39H0+ZbUuQHVAk8Z0rqgSRu5QMCiLdeDBQqKXQ8yDlmPGwivXZmmm4EiH6gEgZyQLNTdRIs/8FEJX5V+vsAtXBuljkq4pRyX6cmhSf/W0btVCeXUScx5iquULA9Jo0y1lwZfjK1vXYhrgON6KPDX0xRpwYf8SWqCFwofns9Xxq42511bXnyNXqe6Kry0rlzoAYiedKQHhy37mAwbekhXo5+qgM9jbD6+X/7kwwCAL03/On5ixS8BAGbXtAYiTovqch96FDK9nupbyNBF+T+AVJb85HHMQGkTu3IZNypqdXuVmTc+8GKg8ECynklYu3YeDlLtiynvwBcMJ0OHBF66Aufq+1EBn8dIK5wPByUUskXWRpLaHQTfNqkbJq0JKiTg8tosO7cgZDHStebuEYkUb2LD8ks4dqG7LhYkqCS02AQtI0iU4uqwWet0bD5hmnS9zlKFXgV8nuNaNHZ8qgnLzi0IOmTr6F28/+w7edZzEqLEnZ+z2GiDWqOcFR4lvHBY2kYntYqsKBkV471+pC+x9e6KQFk/0le2706LWSkhyH1xuucohrtyS3S9tCNsJcYl+HBLmyZAgfw0c/JT0ntBKvTKMn+wOcJPQrG/V+w8Aa+yNz7VFCqIZYNbnko8MnnndM9RrB+xVzQc6BhCL5It4OES+Vr8XjSMMKVIS3fzkouBhbd5yUVsWH7JKjq2SnG9l7di4+hubBzdjWMXuoPVakiUyafLw7JO9xzF+FRTEIGSRr838fxbD2HD8kt5has4xazec3uVwbJzC3DnRHNQo0MWeCIorG58qikyxG4+07XzMMauNqNr52Hrij/0Hbmu8UDHkDOJJ62ogKcY2/JplJ1GiwgQw7fWBv+T6G7b79evuHOiGcvOLUDjZAPaW28G1vvgyUzwEJg+35K3gg+vEiczP9PE5iUXI91LAx1DiT7ftv19oWveOnoXS6/MAMjFmD/Wdc65/9vd72HZuQVYtCW9bqhK0XngSDAX8eq+fut1pMqa9LC0wZdsK2fbZLboXKECnmJsMcVUf6P38tYg9R7wOy4NObft7wuW8Hp1Xz8WbbmORVuuY8W6a4H1Mn3er8lMDwk5XOVCFSeAtY5tSTn+nqwxs3F0d2g0Qtfi4T2jQbLOtv19mOpZiCs77gnmI3ovb8XgyUywX+eBI0GCy6v7+tHeehO3V5k5z6Stdbp2Hs7rf1H97dSOQ6GkIc7wrbV55WnnmnI+QFTAU44t64w6N3VUEh3i1X39eLv7Pbzd/R42ju7GhuWXcOdEc0g4qCaEzRoc6BjCqR2HQhEqlV4guJKQW4knJHF4HfWBjiG0t94Mqg9uHN2N7Y+MoPfyVhy70B1YiW93vxcUXOKFl7h/VYoSPTzTPBFcbujhBuSyXTeO7s5z6dn6X9fOw7iy456QS4q+y3JazEkjUjIPHsTY1ebIUVihqIDXCc987PXQjU/DSfqRVge5TCjK5OE9owD8If6G5ZcwPtUUuGO4b3jb/r5Q4X7Af2Dwkp9pRFq9vEgXvUfCQUJ/uucoHus654cGZqN/6Fq1t94M3ahnM/3OG53vo+Qg4ZUT6Y2TDZg+34LRxx8Irh2f26ER0u1VJi+Kilxh5S5alVTENx3fW9YwWw0jrDNooQK+RBWQuwmovCcAZ/SJDX4T2bZNs9XI60oniTiRafVUs3rxhIdFW66H/OVRk7u8bCpfnX6+pdvbsC06zcvTUqVCyqoE8ueEll6ZsS6/RqPKaq75qnHgipO4VU9s1qZ8nVfnA3JDTx7zTW4Zvl5hWuEiTiGSSRJ16NqNTzXhsa5zeO1ID26vMoGYR0XmUIYgkCudmtZInlKhbEtb3gHPbgVyk/CUUAUgqDdONYAaJxuweMILvgsAeX/PtM2WtIIPX1y5/ZUXCkrQUgFXnMgVTaIsSr7iyeIJLyjkRMN/fpMQ0nJPY+alhC9G4CpcJV93FbwCEJs6TzHfVKaAyhzMxyzNqIU1eBmCFeuuBaMbXjSNBHumbRbLzi0AkKurQmUk7pxoDqJ8CDomUFyMd+bBg5jqWRiUpyhkFKrlZBUrthjiuEnGxRNeKA2fJunO3FgTTFbyH2Lb/r5QVEWaOZvpDwm0bUJT+sMJ8sNS7P30+RbcXmUiV4mRrhYAkeFv9UrngSN5vn9+bQc6hvDqvv6gdAFFPG1Yfim0HU0U315l8Hb3e8F3eXuVwYbll6yhh7y8bzGx95TaX62CWIAKeN2xaEs4IoXit6OsyUVbruP9Z98Jwt3IIuduAmLb/r7gB6j9xWMLYXyqKbgulBBFSDE/3XM0sszA4gkvJDA2ztxYE0yA0j7zCVplSibeLDu3IM+fPX2+JW/ynB54NCFP/fZ0z9Fg+1M7DmHwZCYvdBbI1Wsn90sxIj723J6qzv+oC6XOoBhTuTJ6lCuFhIaGlw/vGQ3V+U5rje9CoZKltjrq3A9L/m7+vyuCJOp68RKpQH09DJPA3VZ8PoUeZFSGgBjoGArqrMvX+D6PdZ0LLfBN7y+9MoOXBl8MfWcbll8KtqVV7eeiRri6UBQrjZMNeeIN+JOPNDFJ4XByn9urDF7d1x+y2Mk/27my/mt0RAkoCflAx1Ce5U3/U2kBwLfq4q7X4gkPAx1DWLHumrWSZL0zPtUUquFzuucotj8ygre73wvKEHA3E/XZ6fMtmD7fgvGppkD4b68yobBBmnwfvrUWM22zeLv7vZB4P9Z1DtPnW3DsQnfwQChnfPZcocWs6gyX35VPRkprcfGEh7e73wOQb6lvWH4Jg0h3jHchuOpouOACQ9d1fKoJ4wlCAf1wwSf9AmD1MZVQMGdurAlZy5uXXMTmHn9kMw3f/Td4MhP4uAk5uiSLWvbfzUsu4hi6A2GfPt8CtM3i2IVunM4uPE11cADfl16NFXqKRS3wOkT6a8kXSB2b+7Y3Hd8bWC6ne46G3iP3yXyCRymcubEm71rKpBJumdP783EyshjGH30qyOrlkJieygrsinXX8sR70/G9QTYmRaEQtvkKvi9/TT6wGycbUlVMTH3gdQbVPeYhV68d6QmlI7uWQrPFiNN79RAqmBQZRx+3Hqb0j8/HUMBisSXsELb4b3IDcp81gCCmW/rHOfJ73HR8L7Y/MhL67gDMSTy+xoErefCFXOPiv2kyKGqhh97LW3Hv5+9P1ZCyXIxdbQ78qzKphy9iwf8udYHe+UrngSNWAZcp8kAuYoQiSGhyEgCmehaGBJwMGYLcLJQ1S3HgMqtzLr5HFXAlD6p1TBOPFOFAQ1EgnC7Of1MiCYkVJUDU+8Sli7GrzUESSZIFjjcd31u1Fc7TTueBI8E8DLk2ZLSJC+qrJOhUtvd0z9GgD9v2Iei+4A/kuejzuiKPkgf5q2n1EV5gioSaRDsoB9uBQOSBXIdeemUGI/NUvAG/tssAmoF9dvcSz8IcvrUWK9Zdw8hz8/d6lQJlABNBWOC6fCuaIKuZx4tTsk9Q+34VcDtrmDzWdS64P/jCDyT0x7qbcKZ1TcGT2NVGBbxOGegYAjp8cdm1/Qn0fOHLocI/dFNsOr438Cc+1nUuECm1JHOTYbaJXB5yNvr4Azir16tobq8yeZnCt1cZ3LZUuMyzjln0zvqRvlA0CfnHGycbsLnnYv4cT/b+eHVff1BbJb/0VW2jAl5nSL/35iUXMdiTCcKz+Ha06AOtvk7v79r+BEbemNNm1ySyDrisG/Ma/MqBeq1Kw2ZhcwFOWjWQJ1S1t94EWhHsD2T95/sQmgglNi+5iEFkyur77tp5uOKVJdUHXidQVht3kbhiwqkDy9TtpVdm1PIWlHMlcsVN+ysvlFQZkKCqkhwKM5T3A1nd9LAo98Rl187DzoJmOomphJClSZedWxAUVOKREgBC8d3kVklzPW9F4ZCIc/85n+ex4Y86y2u8kIDbJrdVwJUAKscZFcttK31KHXygY0gFXKlL+MiUXClykYjNSy5WrP9nHjyIlwZfDOVfdK68rrVQlByNkw3ovbzVKd60DRBejHimbRbT51uwa/sTc9JORZlrzmb6sen4XoxPNWGmbTZY9HvT8b2BC6WSxstLgy8CQGj0W07UAq8DeAIP4cq4JEtEE04UxY9Br+S9wBer2Di6G6uPf4CRN55WC1zJIYW69/LWoAYyEH7qD3QMqXgr85qxq83IPHiwrCvTS7p2HsbY1WasPv5B8BqJdznRMMI6hfv5Rh9/AIM9fsCsirei+Gn3VBSr3HQeOIJT2cSklwZfDFyUlYjwUgu8DqCqbIBvbUs/+PvPvlONZilKzRIVkVIKZNXzlZ1G3ni6YuG56gOvA/jKJoC9kNW2/X0VTypQlLRAQluuESmtriSLw7kmSLUWiuJErq5+50SzireiMMqdcXlqX3+oHMBcheWqBZ5ypPXNFxyo9zUsFaUWkPXjgXgB1ygUxcqZG2tw5saawBpYP9JX0dl2RZnvdK68HkyG7tr+xJwmxakFnnKkBU7I6m5qiStK7aAWuJKI6fMtKt6KUqfoJGbKGZ9qwvDyXN1quQSVxn0rSv2iAp5y2ltv5tVAGegYwsbR3ZGL6/KVtzVCRVHSiQp4yjmb6Q/KZ9IEJoCgsH0cKt6Kkl50ErNOIBGXCT3q/1aU2kMnMRUrtHAxsX6kL1itXlGU+kIt8DqELys1PtUU6QtXFGXu0RV5FEVRUoq6UBRFUeY5KuCKoigpRQVcURQlpaiAK4qipJRSE3kmjDGrytISRVGU+cNEOQ5SUhSKoiiKUj3UhaIoipJSVMAVRVFSigq4oihKSlEBVxRFSSkq4IqiKClFBVxRFCWlqIAriqKkFBVwRVGUlKICriiKklJUwBVFUVKKCriiKEpKUQFXFEVJKSrgiqIoKUUFXFEUJaWogCuKoqQUFXBFUZSUogKuKIqSUlTAFUVRUooKuKIoSkpRAVcURUkpKuCKoigpRQVcURQlpaiAK4qipBQVcEVRlJSiAq4oipJSVMAVRVFSigq4oihKSlEBVxRFSSkq4IqiKClFBVxRFCWlqIAriqKklPtK2bmpsc2b/eBuudpSEbyFDdVuQiLeX2iq3YREvJ+Oywmvwat2ExLR0PCdajchEUvuq+37nPjIvX9f7SYk4q0/f++053mfLvU4JQn47Ad3sW7Zz5Tahooyu6a12k1IxLurG6vdhETcXpWOB81M22y1m5CI9tab1W5CIjYsv1TtJiRi85KL1W5CIta0TS8rx3HUhaIoipJSVMAVRVFSigq4oihKSlEBVxRFSSkq4IqiKClFBVxRFCWlqIAriqKkFBVwRVGUlKICriiKklJUwBVFUVKKCriiKEpKUQFXFEVJKSrgiqIoKUUFXFEUJaWogCuKoqQU43nFF743xrwJIA2V3pcBeLvajUiAtrN8pKGNgLaz3KSlnQs9z/t4qQcpaUEHAHc9z+sqtRGVxhhzQdtZPtLQzjS0EdB2lps0tbMcx1EXiqIoSkpRAVcURUkppQr4fylLKyqPtrO8pKGdaWgjoO0sN/OqnSVNYiqKoijVQ10oiqIoKSWRgBtjPm2MGTPGfMsY8yuW9xuNMV/Mvv9VY0x7uRuaoI3/1Rjz7Wxoo+399caYd40xX8/+fG6u25htx0JjzB8bY75hjHnLGLPfsk3Vr2e2HfcaY75mjHnd8t4OY8xNdj13VqmNHzbGDBtjLhljvmmM+RHxvjHG/Fr2Wv6ZMWZtFdrYya7T140xt4wxvyy2qZX+udsY82a2b/6y5f2qXE/b/W2M+W5jzBljzF9mf3/Ese/77LqeqkI7D2X7558ZY14zxnzYsW+kzlrxPC/yB8C9AK4A+CiABgDfAPD9YptdAH4j+/dnAXwx7rjl/gHwKQBrAbzpeH89gNfnul2WdhgAi7J/LwDwVQA/XGvXM3vuJwEM2a4bgB0AvlAD1/MVADuzfzcA+LB4/zMAvpS97j8M4KtVbu+9AG4AWCVer3r/BPBxAG8C+BD8EOPfB/C9tXA9bfc3gP8I4Feyf/8KgBcc+96Zw2toa+ePA7gv+/cLtnYm0VnbTxIL/JMAvuV53l95njcL4H8C+CmxzU9lbyQAGAaQMcaYBMcuG57n/SGAv5nLcxaD53Mn+++C7I+ciKj69TTGtAL4SQAvz+V5C8EYswT+DfObAOB53qzneX8rNvspAIPZ6/4VAB82xqyY46ZyMgCueJ43UcU2uPg+AF/xPO/vPc/7DoA/APCw2KYq19Nxf/P75BUAP13pdsRha6fneb+XvZ4A8BUArZZdk+hsHkkEvAXAVfb/VPY16zbZhr4L4P4Ex55rfiTruviSMeZj1WpE1jXxdQDfBnDG87yvik1q4Xr+KoB/C+CDiG3+VXZYOGyMWTlH7eJ8FMBNAP8t6+p52Rjzj8Q2SfrvXPJZACcc71W7f74J4FPGmPuNMR+Cb23L77WWruc/9jxvGgCyv7/Hsd1CY8wFY8xXjDHVFvmfhz+CkRR1XZMIuM3ykxZjkm2qzUX4w9Z/BuBFAL9VrYZ4nve+53mfgP8k/qQxRqbUVvV6GmMeAvBtz/P+NGKz3wHQ7nneD8Afar8SsW2luA/+cPU/e573gwD+Dv5QmlMzfdMY0wBgE4D/bXm76v3T87xvwh/inwHwf+EP478jNquZ61kAbZ6fnbkVwK8aY1ZXoxHGmGfgX89XbW9bXou9rkkEfArhp3ArgOuubYwx9wFYihpzZ3ied4tcF57n/S6ABcaYZVVu098COAvg0+Ktal/PBwBsMsaMwx/K9Rhj/gffwPO8dzzPm8n+ewzAD81h+4gpAFNsBDMMX9DlNnH9d674CQAXPc/7a/lGrfRPz/N+0/O8tZ7nfQp+n/tLsUktXc+/JvdN9ve3bRt5nnc9+/uv4N9vPzhXDSSMMY8CeAjANi/r9BYUdV2TCPifAPheY8w/yVoQnwUgZ3JPAXg0+/dmAKOORlYNY8xy8iMbYz4J/7O/U4V2NNEstDHmuwD8GIBLYrOqXk/P8572PK/V87x2+N/3qOd5P8u3EX7PTQC+OVftIzzPuwHgqjGmM/tSBsBfiM1OAdiejZ74YQDv0rC7CmyBw31SQ/3ze7K/2wA8gvz21tL15PfJowB+W25gjPmIMaYx+/cy+MaJ7CMVxRjzaQBPAdjked7fOzZLorP5JJxZ/QyAy/BnSZ/Jvvbvsw0CgIXwh4XfAvDHAD46V7O+rI0nAEwDeA/+0+wXAPwigF/Mvv84gLfgDwu/AmDdXLcx244fAPA1AH8G3+f4uVq8nqy965GNjhBtPMiu5xsA1lSpfZ8AcCF7PX8LwEfE924A/Hq27/45gK4qtfND8AV5KXutFvvnOfgC9w0AGUs7q3I9Hff3/QBG4I8SRgB8d3bbLgAvZ/9el23nN7K/f6EK7fwWfP/217M/FGHWDOB32b55Ohv3o5mYiqIoKUUzMRVFUVKKCriiKEpKUQFXFEVJKSrgiqIoKUUFXFEUJaWogCtlh1V/e9MY8zuu6mts+w8bY3ZVsD3txph/yLbpL4wxv2GMKajvG2N+N+5ziO3/nTGmr/DWKkpyVMCVSvAPnud9wvNX3f4bAL8Us/2H4VdgLAhjzL0FbH7F88sX/ACA70fCwkfZhJV7PM/7jJdfKEtRqooKuFJp/gisKI8xZq8x5k+yRbCoFvp/ALA6ayEfMn5t7NfZPl8wxuzI/j1ujPmcMeb/AfgZY8xZY8wLxq+xftkY0x3VGM8vDnYewD91tSdrsX/TGPMS/BolK7PnXZZ9/8ns6OJNw2pmG2OeMX49598H0Jl3ckUpM/dVuwFK/ZK1kDPIlns1xvw4gO+FXzrTADhljPkU/AJUH89ayDDGrI859F3P8340u+0vwq+1/EljzGcA7INfnsDVpg9l2/S5iPZMwhfgn/M8b1d2P9r/hwD8HIB/nt3nq8aYP4BvDH0Wfp2N++ALf1QxMEUpGRVwpRJ8l/HL5bbDF7Ez2dd/PPvztez/i+AL6GSBx/+i+P9k9vefZs9pY3W2TR6A3/Y870vGmP6I9kx4fr1ryY8CeM3zvL8DAGPMSQDd8AX8NS9b68JUeOUXRQFUwJXK8A+e533CGLMUwOvwfeC/Bt9iPeh53gDf2OQvGfcdhN17C8X7fyf+p6qI78Pdp8kHHjp1RHvkOfg+LrQuhTKnqA9cqRie570L4F8D6DPGLABwGsDPG2MWAYAxpiVb/e42gMVs1wkA32/8tUGXwnd5VAJXe6L4QwA/bYz5kPEXj3gYfgGoPwTwsDHmu4wxiwH8ywq1WVEC1AJXKorneV8zxnwDwGc9z/vvxpjvA/BHWZ/yHQA/63neFWPMl42/EOyXPM/ba4z5X/ArDP4lci6Ocrft92ztgW/Ju/a5aIw5Dr9KJOBXvfsaABhjvgi/2twEfFFXlIqi1QgVRVFSirpQFEVRUooKuKIoSkpRAVcURUkpKuCKoigpRQVcURQlpaiAK4qipBQVcEVRlJSiAq4oipJS/j+/LkesPIe6bgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cartopy\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111, projection=cartopy.crs.Robinson(central_longitude=0, globe=None))\n",
    "\n",
    "contour_plot = ax.contourf(lon,lat,np.transpose(colors))\n",
    "cbar = plt.colorbar(contour_plot, orientation='horizontal')\n",
    "cbar.set_label('Return Period')\n",
    "\n",
    "ax.coastlines()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt \n",
    "import cartopy\n",
    "\n",
    "#fig = plt.figure()\n",
    "#ax_cartopy = fig.add_subplot(111, projection=cartopy.crs.PlateCarree())\n",
    "#ax_cartopy.coastlines()"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}