ag6_era5_extreme_quantiles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [],
   "source": [
    "import rasterio as rio\n",
    "import matplotlib.pyplot as plt \n",
    "import numpy as np\n",
    "import numpy.matlib\n",
    "from scipy import interpolate\n",
    "import statsmodels.api as sm\n",
    "import scipy.stats as st\n",
    "import os, sys, pickle, gzip\n",
    "import datetime\n",
    "import geopy.distance\n",
    "import xarray as xr\n",
    "import xesmf as xe\n",
    "import cartopy.crs as ccrs\n",
    "import glob\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The autoreload extension is already loaded. To reload it, use:\n",
      "  %reload_ext autoreload\n"
     ]
    }
   ],
   "source": [
    "run ../util/setupConsole"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dirAgData = '/home/edcoffel/drive/MAX-Filer/Research/Climate-01/Personal-F20/edcoffel-F20/data/projects/ag-land-climate'\n",
    "# dirEra5 = '/home/edcoffel/drive/MAX-Filer/Research/Climate-02/Data-02-edcoffel-F20/ERA5'\n",
    "\n",
    "dirAgData = '/dartfs-hpc/rc/lab/C/CMIG/ecoffel/data/projects/ag-land-climate'\n",
    "dirEra5 = '/dartfs-hpc/rc/lab/C/CMIG/ERA5'\n",
    "\n",
    "years = [1979, 2019]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [],
   "source": [
    "sacksMaizeNc = xr.open_dataset('%s/sacks/Maize.crop.calendar.fill.nc'%dirAgData)\n",
    "sacksStart = sacksMaizeNc['plant'].values\n",
    "sacksStart = np.roll(sacksStart, -int(sacksStart.shape[1]/2), axis=1)\n",
    "sacksStart[sacksStart < 0] = np.nan\n",
    "sacksEnd = sacksMaizeNc['harvest'].values\n",
    "sacksEnd = np.roll(sacksEnd, -int(sacksEnd.shape[1]/2), axis=1)\n",
    "sacksEnd[sacksEnd < 0] = np.nan\n",
    "\n",
    "sacksLat = np.linspace(90, -90, 360)\n",
    "sacksLon = np.linspace(0, 360, 720)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<title>Show/Hide data repr</title>\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<title>Show/Hide attributes</title>\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt, dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><div class='xr-wrap'><div class='xr-header'><div class='xr-obj-type'>xarray.Dataset</div></div><ul class='xr-sections'><li class='xr-section-item'><input id='section-0099c1e8-159c-4d38-b375-7a8d48777808' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-0099c1e8-159c-4d38-b375-7a8d48777808' class='xr-section-summary'  title='Expand/collapse section'>Dimensions:</label><div class='xr-section-inline-details'><ul class='xr-dim-list'><li><span class='xr-has-index'>latitude</span>: 721</li><li><span class='xr-has-index'>longitude</span>: 1440</li><li><span class='xr-has-index'>time</span>: 480</li></ul></div><div class='xr-section-details'></div></li><li class='xr-section-item'><input id='section-699ad4d5-c383-4f11-b434-784c01c991db' class='xr-section-summary-in' type='checkbox'  checked><label for='section-699ad4d5-c383-4f11-b434-784c01c991db' class='xr-section-summary' >Coordinates: <span>(3)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>latitude</span></div><div class='xr-var-dims'>(latitude)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>90.0 89.75 89.5 ... -89.75 -90.0</div><input id='attrs-aa20a44d-70b4-4cdc-b1dc-3e0a392d3957' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-aa20a44d-70b4-4cdc-b1dc-3e0a392d3957' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-fecc792e-36dc-4318-ae5b-0aff627c7cd2' class='xr-var-data-in' type='checkbox'><label for='data-fecc792e-36dc-4318-ae5b-0aff627c7cd2' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_north</dd><dt><span>long_name :</span></dt><dd>latitude</dd></dl></div><pre class='xr-var-data'>array([ 90.  ,  89.75,  89.5 , ..., -89.5 , -89.75, -90.  ], dtype=float32)</pre></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>longitude</span></div><div class='xr-var-dims'>(longitude)</div><div class='xr-var-dtype'>float32</div><div class='xr-var-preview xr-preview'>0.0 0.25 0.5 ... 359.5 359.75</div><input id='attrs-e697bd91-3fe8-4942-861f-87e8b5804226' class='xr-var-attrs-in' type='checkbox' ><label for='attrs-e697bd91-3fe8-4942-861f-87e8b5804226' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-6bd09dca-1672-4fcc-9ba7-e7bccf35a16c' class='xr-var-data-in' type='checkbox'><label for='data-6bd09dca-1672-4fcc-9ba7-e7bccf35a16c' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'><dt><span>units :</span></dt><dd>degrees_east</dd><dt><span>long_name :</span></dt><dd>longitude</dd></dl></div><pre class='xr-var-data'>array([0.0000e+00, 2.5000e-01, 5.0000e-01, ..., 3.5925e+02, 3.5950e+02,\n",
       "       3.5975e+02], dtype=float32)</pre></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>time</span></div><div class='xr-var-dims'>(time)</div><div class='xr-var-dtype'>datetime64[ns]</div><div class='xr-var-preview xr-preview'>1979-01-31 ... 2018-12-31</div><input id='attrs-e440e4b6-df32-461e-8cfa-a14bf87a3d1e' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-e440e4b6-df32-461e-8cfa-a14bf87a3d1e' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-b73b99ca-0f36-4e0d-8a34-02adb5b4a2ef' class='xr-var-data-in' type='checkbox'><label for='data-b73b99ca-0f36-4e0d-8a34-02adb5b4a2ef' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><pre class='xr-var-data'>array([&#x27;1979-01-31T00:00:00.000000000&#x27;, &#x27;1979-02-28T00:00:00.000000000&#x27;,\n",
       "       &#x27;1979-03-31T00:00:00.000000000&#x27;, ..., &#x27;2018-10-31T00:00:00.000000000&#x27;,\n",
       "       &#x27;2018-11-30T00:00:00.000000000&#x27;, &#x27;2018-12-31T00:00:00.000000000&#x27;],\n",
       "      dtype=&#x27;datetime64[ns]&#x27;)</pre></li></ul></div></li><li class='xr-section-item'><input id='section-ada0f106-108b-49ce-9882-51b27f4cf858' class='xr-section-summary-in' type='checkbox'  checked><label for='section-ada0f106-108b-49ce-9882-51b27f4cf858' class='xr-section-summary' >Data variables: <span>(1)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span>tp</span></div><div class='xr-var-dims'>(time, latitude, longitude)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>8.971 8.971 8.971 ... 10.47 10.47</div><input id='attrs-14dd6628-933a-4904-8b9a-cf1a139aa551' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-14dd6628-933a-4904-8b9a-cf1a139aa551' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-04cea408-f4c9-462e-991e-256284e25201' class='xr-var-data-in' type='checkbox'><label for='data-04cea408-f4c9-462e-991e-256284e25201' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><pre class='xr-var-data'>array([[[ 8.97053852,  8.97053852,  8.97053852, ...,  8.97053852,\n",
       "          8.97053852,  8.97053852],\n",
       "        [10.25733271, 10.22031185, 10.22031185, ..., 10.22797134,\n",
       "         10.22797134, 10.22669476],\n",
       "        [11.1088126 , 11.1024297 , 11.1024297 , ..., 11.10370628,\n",
       "         11.1088126 , 11.10625944],\n",
       "        ...,\n",
       "        [ 1.11828543,  1.1208386 ,  1.1208386 , ...,  1.11700885,\n",
       "          1.11700885,  1.11828543],\n",
       "        [ 1.0238184 ,  1.02509498,  1.0238184 , ...,  1.02637156,\n",
       "          1.02637156,  1.02637156],\n",
       "        [ 1.15275314,  1.15275314,  1.15275314, ...,  1.15275314,\n",
       "          1.15275314,  1.15275314]],\n",
       "\n",
       "       [[ 4.3991    ,  4.3991    ,  4.3991    , ...,  4.3991    ,\n",
       "          4.3991    ,  4.3991    ],\n",
       "        [ 4.65696948,  4.65058657,  4.65186315, ...,  4.65952264,\n",
       "          4.65952264,  4.65824606],\n",
       "        [ 4.60335305,  4.60207647,  4.60207647, ...,  4.60079989,\n",
       "          4.60335305,  4.60207647],\n",
       "        ...,\n",
       "        [ 1.96082925,  1.96210583,  1.96210583, ...,  1.96210583,\n",
       "          1.96338242,  1.96082925],\n",
       "        [ 1.89827676,  1.89317043,  1.89317043, ...,  1.89700018,\n",
       "          1.89700018,  1.89700018],\n",
       "        [ 2.28125122,  2.28125122,  2.28125122, ...,  2.28125122,\n",
       "          2.28125122,  2.28125122]],\n",
       "\n",
       "       [[ 5.42802473,  5.42802473,  5.42802473, ...,  5.42802473,\n",
       "          5.42802473,  5.42802473],\n",
       "        [ 6.02291173,  6.00503959,  5.99993326, ...,  6.04078387,\n",
       "          6.0318478 ,  6.02801806],\n",
       "        [ 6.80545622,  6.79524356,  6.78503091, ...,  6.84630683,\n",
       "          6.83226443,  6.82332836],\n",
       "        ...,\n",
       "        [10.2445669 , 10.24584348, 10.24712006, ..., 10.2445669 ,\n",
       "         10.24329032, 10.24329032],\n",
       "        [10.91604879, 10.91987854, 10.91860195, ..., 10.94030384,\n",
       "         10.94030384, 10.93775068],\n",
       "        [11.97816464, 11.97816464, 11.97816464, ..., 11.97816464,\n",
       "         11.97816464, 11.97816464]],\n",
       "\n",
       "       ...,\n",
       "\n",
       "       [[38.90229008, 38.90229008, 38.90229008, ..., 38.90229008,\n",
       "         38.90229008, 38.90229008],\n",
       "        [39.58918417, 39.59655692, 39.60024329, ..., 39.61744636,\n",
       "         39.60270087, 39.59655692],\n",
       "        [40.53043798, 40.53658193, 40.553785  , ..., 40.553785  ,\n",
       "         40.55132742, 40.54886984],\n",
       "        ...,\n",
       "        [13.5093268 , 13.50809801, 13.50809801, ..., 13.50318284,\n",
       "         13.51301317, 13.51424196],\n",
       "        [13.66661203, 13.65555291, 13.65923928, ..., 13.66661203,\n",
       "         13.66784082, 13.66661203],\n",
       "        [13.9123702 , 13.9123702 , 13.9123702 , ..., 13.9123702 ,\n",
       "         13.9123702 , 13.9123702 ]],\n",
       "\n",
       "       [[13.50072526, 13.50072526, 13.50072526, ..., 13.50072526,\n",
       "         13.50072526, 13.50072526],\n",
       "        [12.95022695, 12.96005728, 12.96497244, ..., 12.94039663,\n",
       "         12.93916784, 12.944083  ],\n",
       "        [13.18861238, 13.18984117, 13.20212908, ..., 13.1726381 ,\n",
       "         13.17632447, 13.18492601],\n",
       "        ...,\n",
       "        [ 5.14494736,  5.14617615,  5.14494736, ...,  5.15109132,\n",
       "          5.14863374,  5.14617615],\n",
       "        [ 5.16092164,  5.15846406,  5.15477769, ...,  5.16460802,\n",
       "          5.16460802,  5.16337923],\n",
       "        [ 5.63523492,  5.63523492,  5.63523492, ...,  5.63523492,\n",
       "          5.63523492,  5.63523492]],\n",
       "\n",
       "       [[21.62671927, 21.62671927, 21.62671927, ..., 21.62671927,\n",
       "         21.62671927, 21.62671927],\n",
       "        [22.85673893, 22.86288288, 22.86288288, ..., 22.86534046,\n",
       "         22.86656926, 22.86165409],\n",
       "        [23.32613704, 23.32982341, 23.3310522 , ..., 23.31016276,\n",
       "         23.30770518, 23.31999309],\n",
       "        ...,\n",
       "        [10.9657297 , 10.96818728, 10.96941608, ..., 10.9657297 ,\n",
       "         10.9657297 , 10.96327212],\n",
       "        [10.51353466, 10.51230587, 10.50984829, ..., 10.51353466,\n",
       "         10.51353466, 10.51230587],\n",
       "        [10.47421336, 10.47421336, 10.47421336, ..., 10.47421336,\n",
       "         10.47421336, 10.47421336]]])</pre></li></ul></div></li><li class='xr-section-item'><input id='section-616a13b7-8850-44c9-83b3-ba1cb4c1f7d7' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-616a13b7-8850-44c9-83b3-ba1cb4c1f7d7' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:    (latitude: 721, longitude: 1440, time: 480)\n",
       "Coordinates:\n",
       "  * latitude   (latitude) float32 90.0 89.75 89.5 89.25 ... -89.5 -89.75 -90.0\n",
       "  * longitude  (longitude) float32 0.0 0.25 0.5 0.75 ... 359.25 359.5 359.75\n",
       "  * time       (time) datetime64[ns] 1979-01-31 1979-02-28 ... 2018-12-31\n",
       "Data variables:\n",
       "    tp         (time, latitude, longitude) float64 8.971 8.971 ... 10.47 10.47"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calc era5 quantiles\n",
    "ds_tasmax = xr.open_mfdataset('%s/monthly/tasmax_*.nc'%dirEra5)\n",
    "ds_tasmax.load()\n",
    "ds_tasmax['mx2t'] -= 273.15\n",
    "\n",
    "ds_tp = xr.open_mfdataset('%s/monthly/tp_*.nc'%dirEra5)\n",
    "ds_tp.load()\n",
    "# ds_quantile = ds.chunk({'latitude':500, 'longitude':500, 'time': -1}).quantile(q=[.01, .05, .1, .9, .95, .99], dim='time')\n",
    "# ds_quantile.tw.to_netcdf('era5_tw_mean_quantiles.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [],
   "source": [
    "lat = ds_tasmax.latitude.values\n",
    "lon = ds_tasmax.longitude.values\n",
    "\n",
    "# regrid sacks data\n",
    "regridMesh = xr.Dataset({'lat': (['lat'], lat),\n",
    "                         'lon': (['lon'], lon),})\n",
    "\n",
    "regridder_start = xe.Regridder(xr.DataArray(data=sacksStart, dims=['lat', 'lon'], coords={'lat':sacksLat, 'lon':sacksLon}), regridMesh, 'bilinear')\n",
    "regridder_end = xe.Regridder(xr.DataArray(data=sacksEnd, dims=['lat', 'lon'], coords={'lat':sacksLat, 'lon':sacksLon}), regridMesh, 'bilinear')\n",
    "\n",
    "sacksStart_regrid = regridder_start(sacksStart)\n",
    "sacksEnd_regrid = regridder_end(sacksEnd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_growing_season(month, sacks_start_month, sacks_end_month):\n",
    "    if sacks_start_month < sacks_end_month:\n",
    "        return (month >= sacks_start_month) & (month <= sacks_end_month)\n",
    "    else:\n",
    "        return (month >= sacks_start_month) | (month <= sacks_end_month)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('maize_hw_2005_regrid.dat', 'rb') as f:\n",
    "    maize_ha_regrid = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [],
   "source": [
    "nnLen = len(np.where( (~np.isnan(np.reshape(sacksStart_regrid, [sacksStart_regrid.size, 1]))) & (np.reshape(maize_ha_regrid, [maize_ha_regrid.size, 1]) > 0))[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 % complete\n",
      "1.1 % complete\n",
      "2.2 % complete\n",
      "3.3 % complete\n",
      "4.4 % complete\n",
      "5.5 % complete\n",
      "6.6 % complete\n",
      "7.7 % complete\n",
      "8.8 % complete\n",
      "9.9 % complete\n",
      "11.0 % complete\n",
      "12.1 % complete\n",
      "13.2 % complete\n",
      "14.3 % complete\n",
      "15.4 % complete\n",
      "16.5 % complete\n",
      "17.6 % complete\n",
      "18.7 % complete\n",
      "19.8 % complete\n",
      "20.9 % complete\n",
      "22.0 % complete\n",
      "23.1 % complete\n",
      "24.2 % complete\n",
      "25.4 % complete\n",
      "26.5 % complete\n",
      "27.6 % complete\n",
      "28.7 % complete\n",
      "29.8 % complete\n",
      "30.9 % complete\n",
      "32.0 % complete\n",
      "33.1 % complete\n",
      "34.2 % complete\n",
      "35.3 % complete\n",
      "36.4 % complete\n",
      "37.5 % complete\n",
      "38.6 % complete\n",
      "39.7 % complete\n",
      "40.8 % complete\n",
      "41.9 % complete\n",
      "43.0 % complete\n",
      "44.1 % complete\n",
      "45.2 % complete\n",
      "46.3 % complete\n",
      "47.4 % complete\n",
      "48.5 % complete\n",
      "49.6 % complete\n",
      "50.7 % complete\n",
      "51.8 % complete\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "n = 0\n",
    "\n",
    "growing_season_tasmax_quantiles = np.full([len(lat), len(lon), 4], np.nan)\n",
    "growing_season_tp_quantiles = np.full([len(lat), len(lon), 4], np.nan)\n",
    "\n",
    "for xlat in range(len(lat)):\n",
    "\n",
    "    for ylon in range(len(lon)):\n",
    "\n",
    "        if ~np.isnan(sacksStart_regrid[xlat, ylon]) and ~np.isnan(sacksEnd_regrid[xlat, ylon]) and maize_ha_regrid[xlat, ylon] > 0:\n",
    "\n",
    "            if n % 1000 == 0:\n",
    "                print('%.1f %% complete'%(n/(nnLen)*100))\n",
    "\n",
    "            sacks_start_month = datetime.datetime.strptime('2020%d'%int(sacksStart_regrid[xlat,ylon]), '%Y%j').date().month\n",
    "            sacks_end_month = datetime.datetime.strptime('2020%d'%int(sacksEnd_regrid[xlat,ylon]), '%Y%j').date().month\n",
    "            \n",
    "            # select growing seasons for all years\n",
    "            ds_tasmax_growing = ds_tasmax.mx2t[:,xlat,ylon].sel(time=is_growing_season(ds_tasmax['time.month'], sacks_start_month, sacks_end_month))\n",
    "            ds_tp_growing = ds_tp.tp[:,xlat,ylon].sel(time=is_growing_season(ds_tasmax['time.month'], sacks_start_month, sacks_end_month))\n",
    "\n",
    "            # resample to group growing seasons (1 year frequency, offset starting at 1 month before start of growing season)\n",
    "            ds_tasmax_growing_1y = ds_tasmax_growing.resample(time='1Y', loffset='%dM'%(sacks_start_month-1)).mean()\n",
    "            ds_tp_growing_1y = ds_tp_growing.resample(time='1Y', loffset='%dM'%(sacks_start_month-1)).mean()\n",
    "\n",
    "            # calculate quantiles across all growing seasons\n",
    "            tasmax_q = ds_tasmax_growing_1y.chunk({'time': -1}).quantile(q=[.75, .9, .95, .99], dim='time', skipna=True)\n",
    "            tp_q = ds_tp_growing_1y.chunk({'time': -1}).quantile(q=[.75, .9, .95, .99], dim='time', skipna=True)\n",
    "\n",
    "            growing_season_tasmax_quantiles[xlat, ylon, :] = tasmax_q.values\n",
    "            growing_season_tp_quantiles[xlat, ylon, :] = tp_q.values\n",
    "\n",
    "            n += 1\n",
    "                \n",
    "#                 sys.exit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ecoffel",
   "language": "python",
   "name": "ecoffel"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}