ccarouge/First try.ipynb

## First try.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Correlation coeff with 3D arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray (c: 2, time: 10, z: 100)>\n",
       "array([[[250.31708, 250.31708, ..., 255.19339, 253.62856],\n",
       "        [247.862  , 247.862  , ..., 249.95056, 248.54636],\n",
       "        ...,\n",
       "        [227.74347, 227.74347, ..., 241.1389 , 241.71703],\n",
       "        [238.68315, 238.68315, ..., 247.05644, 245.7788 ]],\n",
       "\n",
       "       [[242.73505, 242.73505, ..., 243.99998, 242.34166],\n",
       "        [230.58101, 230.58101, ..., 241.5458 , 240.07744],\n",
       "        ...,\n",
       "        [209.07805, 209.07805, ..., 226.06282, 224.69025],\n",
       "        [214.15959, 214.15959, ..., 232.0684 , 231.35303]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * time     (time) object 1850-01-16 12:00:00 ... 1850-10-16 12:00:00\n",
       "    height   float64 2.0\n",
       "  * z        (z) MultiIndex\n",
       "  - lat      (z) float64 -90.0 -90.0 -90.0 -90.0 ... -86.21 -86.21 -86.21 -86.21\n",
       "  - lon      (z) float64 0.0 3.75 7.5 11.25 15.0 ... 18.75 22.5 26.25 30.0 33.75\n",
       "Dimensions without coordinates: c\n",
       "Attributes:\n",
       "    standard_name:     air_temperature\n",
       "    units:             K\n",
       "    cell_measures:     area: areacella\n",
       "    associated_files:  baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation..."
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Take 2 files from CMIP5 for testing\n",
    "data = xr.open_dataset('/g/data/al33/replicas/CMIP5/combined/IPSL/IPSL-CM5A-LR/historical/mon/atmos/Amon/r1i1p1/v20110406/tasmax/tasmax_Amon_IPSL-CM5A-LR_historical_r1i1p1_185001-200512.nc')\n",
    "tmax=data.tasmax\n",
    "\n",
    "data2 = xr.open_dataset('/g/data/al33/replicas/CMIP5/combined/IPSL/IPSL-CM5A-LR/historical/mon/atmos/Amon/r1i1p1/v20110406/tasmin/tasmin_Amon_IPSL-CM5A-LR_historical_r1i1p1_185001-200512.nc')\n",
    "tmin = data2.tasmin\n",
    "\n",
    "# We need to create 1 dataarray with all the data in. And we reduce the size to a (10,10,10) cube for testing.\n",
    "# To create a dataarray, we add a dimension ('c') and concat along this dimension.\n",
    "# We also need to transform our 3D arrays into 2D arrays. We use the stack() method of xarray for that. (New MultiIndex dimension z)\n",
    "sample1=tmax[0:10,0:10,0:10]\n",
    "sample1=sample1.stack(z=('lat','lon'))\n",
    "sample1= sample1.expand_dims(dim='c',axis=0)\n",
    "sample2=tmin[0:10,0:10,0:10]\n",
    "sample2=sample2.stack(z=('lat','lon'))\n",
    "sample2= sample2.expand_dims(dim='c',axis=0)\n",
    "sample = xr.concat((sample1, sample2),dim='c')\n",
    "sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can simply use groupby and scipy.stats.pearsonr to calculate what we need\n",
    "# The squeeze() is because the stacking means arrays have degenerated dimensions.\n",
    "def corr2(ar,axis):\n",
    "    ar1=ar[0,...].squeeze()\n",
    "    ar2=ar[1,...].squeeze()\n",
    "    stats = xr.DataArray(np.asarray(scipy.stats.pearsonr(ar1,ar2)),dims='c',coords={'c':['file1','file2']})\n",
    "    return stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "stats = sample.groupby('z').reduce(corr2,dim='time')\n",
    "stats=stats.unstack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7fa5a9affc88>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu4XVV97vHvm4RwC5cgiCHhEjACASFKRKuCiKCgtaitQFoBqYhY8CAPjwq0Hu8eSqFKH6hpjCgcqYgiklo0BVqKelQIEggBNCFQCYmBcBUQQnbe88ecG9ZeWWuvlWTNtdYm7+d55pN5H2MvNuu3xxhz/oZsExERMWhUrysQERH9JYEhIiKGSGCIiIghEhgiImKIBIaIiBgigSEiIoZIYIiuknS/pMPW89qFkg6pupyIjV0CQ4wYtvexfeOG3kfSIZKWdqBKje59vqRFkv4g6R5Jx7c4/y8l/Y+kpyX9UNJ2VdQrYl0kMER01tPAu4FtgBOACyW9sdGJkvYB/gU4DtgReAb45y7VM6KpBIbohWmS7pD0hKTvStps8ICkP5U0X9Ljkv6fpP1qjr3QPSRpc0mXSnpM0t2SPtmgFbBWOZK2BH4M7CTpqXLZqVM/mO3P2L7H9hrbvwJ+CvxJk9P/Cvg32zfZfgr4NPA+SVt1qj4R6yOBIXrhaOAIYDKwH/BBAEmvBS4BPgK8jOKv6TmSNm1wj88AuwG7A4cDH2inHNtPA0cCy2yPK5dl9RdKOqsMTg2Xdn5ISZsDrwMWNjllH+D2wQ3b9wKrgFe1c/+IqiQwRC/8k+1lth8F/g2YVu7/MPAvtn9le8D2pcBzwBsa3ONo4Mu2H7O9FPindSinJdvn2t622dLmbWZSfPHPbXJ8HPBE3b4ngLQYoqcSGKIXfl+z/gzFFyTArsCZdX+Z7ww06urZCXigZvuBBuc0K6dykv4B2Bc42s0zVT4FbF23b2vgD1XWLaKVBIboJw8AX6r763wL299pcO5yYFLN9s7rUE7LlMKSzqkZg1hraXHt5yi6q95u+8lhTl0I7F9z3e7ApsBv2/sxIqqRwBD95OvAKZJer8KWkt7VZDD2SuBsSeMlTQROW4dyVgAvk7RNsxNsf7lmDGKtpdl1ks4G/hI43PYjLepxOfBuSQeVg+KfB35gOy2G6KkEhugbtudRjDNcBDwGLKYcmG7g88BS4D7geuD7FOMR7ZRzD/AdYEnZZdWxp5KALwO7AItqWhjnDB4stw8q67EQOIUiQDxEMbbwNx2sS8R6USbqiZcCSR8FjrX9ll7XJWKkS4shRiRJEyS9SdIoSXsCZwJX97peES8FY3pdgYj1NJbiPYfJwOPAFeSt4YiOSFdSREQMka6kiIgYYqPqStpy/FiPn7h55eU8/EzX3qNCq7oX28dssbprZe246XCP/3fOuFHPd6UcgBWru/d7MZo1XSurm30Oz7t7X1kr735kpe0dNuQe73jrln7k0YGW5916x3NzbR+xIWV10kYVGMZP3JxTr3xT5eXMvuPNlZcxSA9s1vqkDtl+v4e7VtaZr/yPrpRz0GbLu1IOwPkrD+paWVuP/mPXylqDulbWsmfbzUay4b7+usv+Z0Pv8cijA9w8d5eW542esGj7DS2rkzaqwBAR0U0G1nSx9dYpPRljKFMgzy+X+yXNL/cfLulWSQvKfw9tcv1nJT1Yc493dvcniIhozZjnPdBy6Tc9aTHYPmZwXdIFvJhhciXwbtvLJO1LkZVyYpPbfMX2+dXWNCJiw4zEFkNPu5IkiSJ98qEAtm+rObwQ2EzSprbbSnUQEdFPjBkYga8E9Ppx1YOAFbYXNTj258BtwwSF08rZuS6RNL5ZAZJOljRP0rynH13ViTpHRLRtDW659JvKAoOk6yXd2WA5qua0GRTJzOqv3Qf4e4qZvBr5GrAHxcQry4ELmtXD9izb021P33K7sev980RErCsDA7jl0m8q60qyfdhwxyWNAd4HHFC3fxJFzpvjy6kOG917Rc35Xwd+tMEVjoioQD+2CFrp5RjDYcA95bSMAEjaFvh34GzbP292oaQJtgcfQH8vcGelNY2IWA8Gns8Ywzo5lrW7kU4DXgl8uuZR1JcDSJotaXp53nnlI613AG8FzuharSMi2uQ2upE2qq6kVmx/sMG+LwJfbHL+STXrx1VXs4iIDjEM9N/3fkt58zkioiLFm88jTwJDRERlxEAXc0l1SgJDRERFisHnBIaIiCgV7zEkMERERI01aTFERMSgtBgiImIIIwZ6npJu3SUwRERUaCR2JY28UBYRMUIYscqjWy7tkHSEpN9IWizprAbHx0u6usw6fXM5p83gsW0lfV/SPZLulvQnw5WVwBARUZHiBbdRLZdWJI0GLgaOBKYCMyRNrTvtHGC+7f2A44ELa45dCPzE9l7A/sDdw5WXwBARUaGB8iW34ZY2HAgstr3E9irgCuCounOmAjcA2L4H2E3SjpK2Bg4GvlEeW2X78eEK26jGGB55dkv+dfHrKi9nj6+urryMQavGd6+s/zzuqq6V9cWVe3WlnE9cP6Mr5QBM+eivulbW8jPf2rWyVjWdJquCsrbrv/mRh2OLAbf19/f2kubVbM+yPatmeyLwQM32UuD1dfe4nWIqg59JOhDYFZgEDAAPA9+UtD9wK3C67aebVSYthoiICq1BLRdg5eCEYuUyq+42jZoV9en5zgXGS5oPfAy4DVhN0QB4LfA1268BngbWGqOotVG1GCIiuqkYfO7I1+xSYOea7UnAsiFl2U8CJwJIEnBfuWwBLLU92GT9Pi0CQ1oMEREV6dTgM3ALMEXSZEljKeazmVN7Qvnk0eD8xScBN9l+0vbvgQck7Vkeextw13CFpcUQEVGhgQ68x2B7taTTgLnAaOAS2wslnVIenwnsDVwmaYDii/9DNbf4GHB5GTiWULYsmklgiIioSCfffLZ9LXBt3b6ZNeu/AKY0uXY+ML3RsUYSGCIiKrSmvaeS+koCQ0RERYokegkMERFRMuL5NlNe9JOeBAZJ3wUGR8i3BR63PU3SbhSvav+mPPZL26c0uH474LvAbsD9wNG2H6u21hER68am3Rfc+kpPAoPtYwbXJV0APFFz+F7b01rc4izgBtvnlsmkzgI+1fmaRkRsiBdeYBtRehrKypcwjga+s46XHgVcWq5fCrynk/WKiOgEU7QYWi39ptc1OghYYXtRzb7Jkm6T9N+SDmpy3Y62lwOU/768WQGSTpY0T9K8gSef6VzNIyLaMMColku/qawrSdL1wCsaHPpb29eU6zMY2lpYDuxi+xFJBwA/lLRP+ar3eilzjswC2PyVO9XnFomIqIzRiJyop7LAYPuw4Y5LGkORCfCAmmueA54r12+VdC/wKmBe3eUrJE2wvVzSBOChjlY+IqIDDDzfmVxJXdXLNsxhwD22lw7ukLRDOSEFknaneItvSYNr5wAnlOsnANc0OCciosdaz8XQ5nwMXdXLwHAsaw86HwzcIel2igyAp9h+FEDSbEmDr3SfCxwuaRFweLkdEdFXTPHmc6ul3/SsjWP7gw32XQU0nA3G9kk1649QZAiMiOhr/dgiaGXkdX5FRIwQtvqyRdBKAkNEREWKweekxIiIiBe0PedzX9moAsPYMQPssm31KZU+ePnPKy9j0NzHXt21svb6+XFdK2vc5s91pZwtdnyqK+UAPHfdbl0ra6vZa7pW1rI9u1fW2IdH1ldWMficMYaIiKjRj282tzLyahwRMUIMvvncammHpCMk/UbS4jJ5aP3x8ZKulnSHpJsl7Vt3fHSZbuhHrcpKYIiIqNAaRrVcWilf/L0YOBKYCsyQNLXutHOA+bb3A44HLqw7fjrFtAYtJTBERFTEhufXjGq5tOFAYLHtJbZXAVdQZJmuNRW4oSjX9wC7SdoRQNIk4F3A7HYKS2CIiKhI0ZXUkTefJwIP1GwvLffVup0i/xySDgR2BSaVx74KfBJo60mBBIaIiAq1mStp+8HpAcrl5LrbNBqIqM8WfS4wXtJ84GPAbcBqSX8KPGT71nbrnKeSIiIqsg6Pq660PX2Y40uBnWu2JwHLhpRVTE9wIrwwCdp95XIs8GeS3glsBmwt6du2P9CssLQYIiIq07GupFuAKZImSxpL8WU/Z0hJ0rblMYCTgJtsP2n7bNuTbO9WXvefwwUFSIshIqJSnZjz2fZqSacBc4HRwCW2F0o6pTw+E9gbuEzSAHAX8KH1LS+BISKiIsVTSZ3JlWT7WuDaun0za9Z/QTGHzXD3uBG4sVVZCQwRERXJ1J4REbGWTnQldVsCQ0RERZJELyIi1pKJeiIi4gW2WD0CA0NPaizpu5Lml8v95Zt6SPqrmv3zJa2RNK3B9Z+V9GDNee/s/k8REdFap7KrdlNPWgy2jxlcl3QB8ES5/3Lg8nL/q4FrbM9vcpuv2D6/6rpGRKyvjDGsh/K17aOBQxscngF8p7s1iojorJEYGHrd+XUQsML2ogbHjmH4wHBaOSHFJZLGNztJ0smDialWPf7HDa1vRETbOjlRTzdVFhgkXS/pzgZLbQ7xhq0CSa8HnrF9Z5Pbfw3YA5gGLAcuaFYP27NsT7c9fey2m6//DxQRsR7WoJZLv6msK8n2YcMdlzSGInf4AQ0OH8swrQXbK2ru83Wg5VR1ERHdZsPq9ibi6Su9HGM4DLjH9tLanZJGAe8HDm52oaQJtpeXm+8FmrUshnh21RjuXvqK9axu+85Z8d7Kyxg0sLp7v3RrVnUm50s7dn3F77tSzt3XvbIr5QD8fuw2XStrzB5dK4pRz3avrGlv+W3Xylrcofv0Y1dRK70MDM1aBQcDS20vqd0paTYw0/Y84LzyMVYD9wMfqbiuERHrLLmS1pHtDzbZfyPwhgb7T6pZP66yikVEdJATGCIiolY/Di63ksAQEVERO2MMERExhBjIU0kREVFrJI4xjLxQFhExQgzmSurEm8+SjpD0G0mLJZ3V4Ph4SVeXGSFulrRvuX9nSf8l6W5JCyWd3qqsBIaIiKq4GGdotbQiaTRwMXAkMBWYIWlq3WnnAPNt7wccD1xY7l8NnGl7b4onPk9tcO0QCQwRERXqUEqMA4HFtpfYXgVcARxVd85U4AYA2/cAu0na0fZy278u9/8BuBuYOFxhCQwRERVxOfjcagG2H0z2WS4n191qIvBAzfZS1v5yv50izRCSDgR2BSbVniBpN+A1wK+Gq3cGnyMiKtROVxGw0vb0YY43albU3/lc4MJy4rMFwG0U3UjFDaRxwFXAx20/OVxlEhgiIirUoaeSlgI712xPApYNLcdPAifCC3Pd3FcuSNqEIihcbvsHrQpLV1JEREWKwWW1XNpwCzBF0mRJYylyzc2pPUHStuUxgJOAm2w/WQaJbwB32/7HdgpLiyEiokKdePPZ9mpJpwFzgdHAJbYXSjqlPD4T2Bu4TNIAcBfwofLyNwHHAQvKbiaAc2xf26y8BIaIiAq1OcbQxn18LXBt3b6ZNeu/AKY0uO5nNB6jaCqBISKiIkasSUqMiIio1aEGQ1clMEREVMUjM1dSAkNERJVGYJMhgSEiokJpMURExAsMrFmTwNDXRv1xFJst3Lzycp555fOVlzFoyuTlXStr3CarulbWKHWn/X3sX9zYlXIA/vuhV3atrE3O2qZrZT08fauulXXrmu59hh1hYAS2GHryHJWkaZJ+KWl+mTDqwJpjZ5f5xn8j6R1Nrt9O0nWSFpX/ju9e7SMi2teJtNvd1qsHbM8DPmd7GvC/y23KHOHHAvsARwD/XOYhr3cWcIPtKRRpZteatCIioi+4jaXP9CowGNi6XN+GF5NBHQVcYfs52/cBiynykNc7Cri0XL8UeE+FdY2IWE+t8yT14+B0r8YYPg7MlXQ+RXB6Y7l/IvDLmvMa5RwH2NH2cgDbyyW9vFlBZV7zkwHGbJ0ep4josj5sEbRSWWCQdD3wigaH/hZ4G3CG7askHU2R+e8w2ss5vk5szwJmAWw+YecR+J8oIkYsg/NU0otsH9bsmKTLgMEJqb8HzC7XW+YcL62QNKFsLUwAHupAlSMiKjDyAkOvxhiWAW8p1w8FFpXrc4BjJW0qaTJFpsCbG1w/BzihXD8BuKbCukZErL8ROPjcqzGGD1NMQTcGeJZyDKDML34lRS7x1cCptgcAJM0GZtqeRzGF3ZWSPgT8Dnh/D36GiIjW+vCLv5WeBIYyP/gBTY59CfhSg/0n1aw/QjFOERHRv0boC24b1ZvPERHd1o8vsLWSwBARUaUR+FTSyJtaKCJiBJFbL23dRzqiTBW0WNJa2R4kjZd0taQ7JN0sad92r62XwBARUZV2nkhqIzCUqYEuBo4EpgIzyhRCtc4B5tveDzgeuHAdrh0igSEiojIqBp9bLa0dCCy2vcT2KuAKitRAtaZS5I7D9j3AbpJ2bPPaIRIYIiKq1F6LYfsy0/TgcnLdXSYCD9RsN0oXdDvwPoAyY/WuFC8Jt3PtEBl8joio0pq2zlppe/owx9tJF3Quxfth84EFwG0U74Otc6qhtloMkv6+nX0REVFj8D2GDe9KapkuyPaTtk8spzM4HtgBuK+da+u125V0eIN9R7Z5bUTERqtDTyXdAkyRNFnSWIp5a+YMKUfatjwGcBJwk+0n27m23rBdSZI+CvwNsLukO2oObQX8vK0fJyJiY9aBF9xsr5Z0GjAXGA1cUqYQOqU8PhPYG7hM0gBFWqEPDXftcOW1GmP4V+DHwP9h6Cxpf7D96Dr/dBERsV5sXwtcW7dvZs36LygSj7Z17XCGDQy2nwCeAGYAlBPibAaMkzTO9u/aLSgiYmPU7gts/aStp5IkvRv4R2AnirkPdgXuppibecTQuAHGvOGxysuZvOXTlZcx6AMTf9W1snbf5KU37cWiVY3mkqrGx3bt3t9RE7/3eNfK+sTiv+haWU8/tnXrk/qJeUmnxPgi8Abgt7YnU2Q2zRhDREQrI3A+hnYDw/NlqutRkkbZ/i9gWoX1ioh4SehUrqRuavcFt8cljQNuAi6X9BDFixMRETGcPvzib6XdFsNRwB+BM4CfAPcC766qUhERLxkjsCuprRaD7drR1EsrqktExEtKv3YVtdLqBbc/0DieCbDtEfaIQEREl43Ap5JavcewVbcqEhHxUjQSWww9SbstaZqkX0qaX6aYPbDcf7ikWyUtKP89tMn1n5X0YHn9fEnv7O5PEBHRppfqGEMFzgM+Z/vH5Zf6ecAhwErg3baXldPSzaV53vCv2D6/K7WNiFgfL8UxhgoZGByf2IYyBazt22rOWQhsJmlT2891uX4REZ2RwNC2jwNzJZ1P0Z31xgbn/Dlw2zBB4TRJxwPzgDNtN8x1Uc6EdDLAJjtss8EVj4hYF2pvop6+UtkYg6TrJd3ZYDkK+Chwhu2dKd6N+EbdtfsAfw98pMntvwbsQfH29XLggmb1sD3L9nTb08dss0UHfrKIiJe2yloMtg9rdkzSZcDp5eb3gNk1xyYBVwPH2763yb1X1Jz/deBHnahzRETHjcCupJ48lUQxpvCWcv1QYBEUMxAB/w6cbbtpkj5JE2o23wvcWVE9IyLWXxt5kvpxcLpXYwwfppi0egzwLOUYAHAa8Erg05I+Xe57u+2HJM0GZtqeB5wnaRpFLL6f5l1OERG91Ydf/K30JDDY/hlwQIP9X6RI8d3ompNq1o+rrnYRER3UocAg6QjgQorpOWfbPrfu+DbAt4FdKL7bz7f9zfLYGRTzQBtYAJxo+9lmZfWqKyki4iVPFE8ltVpa3kcaDVwMHAlMBWZImlp32qnAXbb3p3gv7AJJYyVNBP4XMN32vhSB5djhyktgiIioSufGGA4EFtteYnsVcAVF1uu60thKkoBxwKO8OD3CGGDzsvt+C8p3x5pJYIiIqFJnUmJMBB6o2V7K2lkhLgL2pvjSXwCcbnuN7QeB84HfUTze/4Tt/xiusASGiIgqtRcYti/zxg0uJ9fdpVGK1vqQ8g5gPrATxTteF0naWtJ4itbF5PLYlpI+MFyVe/VUUk8MDIinn96s8nJWD3Qv3n77wdd3raztN3u69Ukdsue4Fa1P6oCHVr00EwiPH/NM18radVzDpAOVOHLCwq6VdU6H7tNmV9FK29OHOb4U2LlmexJrdwedCJxr28BiSfcBewG7AvfZfhhA0g8osk18u1lhaTFERFSpM11JtwBTJE2WNJZi8HhO3Tm/A94GIGlHYE9gSbn/DZK2KMcf3gbcPVxhG1WLISKiq9yZXEm2V0s6jSLj9GjgEtsLJZ1SHp8JfAH4lqQFFF1Pn7K9Elgp6fvArykGo28DZg1XXgJDRESVOvQeg+1rgWvr9s2sWV8GvL3JtZ8BPtNuWQkMEREV6seUF60kMEREVCmBISIiXtCnU3e2ksAQEVERka6kiIiok8AQERFDJTBERMQQCQwREfGCPp2hrZUEhoiIKiUwRERErU6kxOi2BIaIiAqNxK6knmRXlTRN0i8lzS9zjx9Y7t9N0h/L/fMlzWxy/XaSrpO0qPx3fHd/goiINrSTWbUPA0ev0m6fB3zO9jTgf5fbg+61Pa1cTmly/VnADbanADeU2xER/SeBoW0Gti7Xt6HF/KMNHAVcWq5fCrynQ/WKiOiYwTefOzDnc1f1aozh48BcSedTBKc31hybLOk24Eng72z/tMH1O9peDmB7uaSXV17jiIj1oDV9+M3fQmWBQdL1wCsaHPpbihmEzrB9laSjgW8Ah1FMVL2L7UckHQD8UNI+tp/cgHqcDJwMMGb7bdb3NhER665Pu4paqSww2D6s2TFJlwGnl5vfA2aX1zwHPFeu3yrpXuBVwLy6W6yQNKFsLUwAHhqmHrMoZyvadPeJI/A/UUSMZP3YVdRKr8YYlgFvKdcPBRYBSNpB0uhyfXdgCsWcpfXmACeU6ycA11Ra24iI9TUCB597NcbwYeBCSWOAZym7eoCDgc9LWg0MAKfYfhRA0mxgpu15wLnAlZI+RDHR9fvbKVSrRjHqvs07+5M0sGpU9WUMWjJ669Yndci9o7tWFL942R5dKqkP/6/sBKt7ZT3dvV+Mn27+qq6VBT/uyF061WKQdARwIcWcz7Ntn1t3fBvg28AuFN/t59v+ZnlsW4qemX0pfun/2vYvmpXVk8Bg+2fAAQ32XwVc1eSak2rWH6EYp4iI6G8dCAxlT8rFwOHAUuAWSXNs31Vz2qnAXbbfLWkH4DeSLre9iiKg/MT2X0gaC2wxXHl58zkioiruWEqMA4HFtpcASLqC4rH92sBgYCtJAsYBjwKrJW1N0RvzQYAyUKwarrBejTFERLzkrcN7DNuXWSAGl5PrbjUReKBme2m5r9ZFwN4UY7gLgNNtrwF2Bx4GvinpNkmzJW05XL0TGCIiqmS3XmCl7ek1y6y6uzQaNKrvpHoHMB/YCZgGXFS2FsYArwW+Zvs1wNO0yBaRwBARUaEOvfm8FNi5ZnsSa2eMOBH4gQuLgfuAvcprl9r+VXne9ykCRVMJDBERVelcEr1bgCmSJpeDx8dSPLZf63eUD+VI2hHYE1hi+/fAA5L2LM97G0PHJtaSweeIiAp1YvDZ9mpJpwFzKR5XvcT2QkmnlMdnAl8AviVpAUXX06dsryxv8THg8jKoLKFoXTSVwBARUaFOTdRj+1rg2rp9M2vWlwFvb3LtfGB6u2UlMEREVMUMDi6PKAkMEREVGom5khIYIiKqlMAQERGDBl9wG2kSGCIiqmJnop6IiKgz8uJCAkNERJXSlRQRES8ykK6kiIgYYuTFhQSGiIgqpSspIiKGyFNJERHxovazp/aVBIaIiIoUL7iNvMjQk8AgaRowE9gMWA38je2bJf0V8ImaU/cDXltmBqy9/rPAhymmqwM4p8w8ODzDqOc3vP6tdKOMQXKjiZ0qKqubP9fysd0pqIv/z3Yqy2ZbZXXx5xrYtJtlje5eYZ3Sxf/undKrFsN5wOds/1jSO8vtQ2xfDlwOIOnVwDX1QaHGV2yf353qRkSsn7QY2mdg63J9G9aeog5gBvCdrtUoIqLTMsawTj4OzJV0PsX0om9scM4xwFHD3OM0SccD84AzbT/W6CRJJwMnA4zZZvwGVToiYt2MzFxJlc35LOl6SXc2WI4CPgqcYXtn4AzgG3XXvh54xvadTW7/NWAPYBqwHLigWT1sz7I93fb00Vtu2YkfLSKifXbrpc9U1mKwfVizY5IuA04vN78HzK475ViG6UayvaLmXl8HfrT+NY2IqIg799CBpCOACynmfJ5t+9y649sA3wZ2ofhuP9/2N2uOj6boYXnQ9p8OV1ZlLYYWlgFvKdcPBRYNHpA0Cng/cEWziyVNqNl8L9CsZRER0VsdaDGUX+oXA0cCU4EZkqbWnXYqcJft/YFDgAsk1T7edzpwdztV7lVg+DBFpW8Hvkw5BlA6GFhqe0ntBZJmSxqczPo8SQsk3QG8laI7KiKi/7iNpbUDgcW2l9heRfGHc/0YrIGtJAkYBzxK8ToAkiYB72Lt3pmGejL4bPtnwAFNjt0IvKHB/pNq1o+rrHIRER2kNW31JW0vaV7N9izbs2q2JwIP1GwvBV5fd4+LgDkUPTJbAcfYHiz8q8Any/0t5c3niIiqmHZfcFtpe/owxxu9yVrf1ngHMJ+ie34P4DpJP6XohXnI9q2SDmmnMr3qSoqIeMkTRm69tGEpsHPN9iTWfv/rROAHLiwG7gP2At4E/Jmk+ym6oA6V9O3hCktgiIioUmceV70FmCJpcjmgfCxFt1Gt3wFvA5C0I7AnsMT22bYn2d6tvO4/bX9guMLSlRQRUaUOvKdge7Wk04C5FI+rXmJ7oaRTyuMzgS8A35K0gKLr6VO2V65PeQkMERFVaX+MofWtikSh19btm1mzvgx4e4t73Ajc2KqsBIaIiAq1+VRSX0lgiIioTH+mvGglgSEioiomgSEiIuqMvJ6kBIaIiCplop6IiBgqgSEiIl5gw8DI60tKYIiIqFJaDP3Nm5hnd1rd62pExMYkgSEiIl5gYATO+ZzAEBFRGYMzxhAREYNMBp8jIqJOxhgiImKIBIaIiHjRyEyi15MZ3CTtL+kXkhZI+jdJW9ccO1vSYkm/kfSOJtdvJ+k6SYtNl/cLAAAHTklEQVTKf8d3r/YREW0ysGZN66XP9Gpqz9nAWbZfDVwNfAJA0lSKqef2AY4A/lnS6AbXnwXcYHsKcEO5HRHRfzoztWdX9Sow7AncVK5fB/x5uX4UcIXt52zfBywGDmxw/VHApeX6pcB7KqxrRMR6KlNitFraIOmIsidlsaS1/hiWtE3ZA3O7pIWSTiz37yzpvyTdXe4/vVVZvQoMdwJ/Vq6/H9i5XJ8IPFBz3tJyX70dbS8HKP99ebOCJJ0saZ6keQN/eHqDKx4R0TaDvabl0krZc3IxcCQwFZhR9rDUOhW4y/b+wCHABZLGAquBM23vDbwBOLXBtUNUFhgkXS/pzgbLUcBfl5W7FdgKWDV4WYNbbVA7y/Ys29NtTx+91ZYbcquIiHW3xq2X1g4EFtteYnsVcAVFz0ktA1tJEjAOeBRYbXu57V8D2P4DcDeN/+B+QWVPJdk+rMUpbweQ9CrgXeW+pbzYegCYBCxrcO0KSRNsL5c0AXhoQ+sbEVGJ9sYQtpc0r2Z7lu1ZNduNelNeX3ePi4A5FN+ZWwHHuK45Imk34DXAr4arTK+eSnp5+e8o4O+AmeWhOcCxkjaVNBmYAtzc4BZzgBPK9ROAa6qtcUTEerDbfSpp5WDPRrnMqrtTO70p7wDmAzsB04CL6p74HAdcBXzc9pPDVbtXYwwzJP0WuIciun0TwPZC4ErgLuAnwKm2BwAkzZY0vbz+XOBwSYuAw8vtiIj+05mnktrpTTkR+IELi4H7gL0AJG1CERQut/2DVoX15AU32xcCFzY59iXgSw32n1Sz/gjwtsoqGBHREcYDA5240S3AlLIn5UGKx/r/su6c31F8L/5U0o4UT38uKcccvgHcbfsf2yksbz5HRFSlQ2m3ba+WdBowFxgNXGJ7oaRTyuMzgS8A35K0gKLr6VO2V0p6M3AcsEDS/PKW59i+tll5CQwREVXqUNrt8ov82rp9M2vWl1E+1FN3zs9oPEbRVAJDRERFDDgT9URExAuciXoiIqJOhwafu0ruwwROVZH0MPA/63jZ9sDKCqqzvlKf5vqpLpD6tNJP9WlUl11t77AhN5X0k/Leray0fcSGlNVJG1VgWB+S5tme3vrM7kh9muunukDq00o/1aef6tIPevWCW0RE9KkEhoiIGCKBobX6nCW9lvo01091gdSnlX6qTz/VpecyxhAREUOkxRAREUMkMERExBAJDMNoNcdqD+pzv6QFkubXTerRrfIvkfSQpDtr9m0n6TpJi8p/x/ewLp+V9GD5+cyX9M5u1KUsu+G8ur34fIapS08+H0mbSbq5Zi7iz5X7e/W706w+Pfv96TcZY2iinGP1txTzPSylSHs7w/ZdPazT/cB02z15KUjSwcBTwGW29y33nQc8avvcMniOt/2pHtXls8BTts+vuvwG9ZkATLD9a0lbAbcC7wE+SJc/n2HqcjQ9+HzKtM9b2n6qnBfgZ8DpwPvoze9Os/ocQY9+f/pNWgzNtTPH6kbF9k0U88jWOgq4tFy/lOILqFd16Zlh5tXt+uezPnP8Vlwf236q3NykXEzvfnea1SdKCQzNNZpjtWf/c5UM/IekWyWd3OO6DNrR9nIovpCAl/e4PqdJuqPsaupK10S9unl1e/r5NJjjtyefj6TR5VwADwHX2e7pZ9OkPtAHvz/9IIGhuXbmWO22N9l+LXAkcGrZnRIv+hqwB8V8t8uBC7pdgXWZV7cHdenZ52N7wPY0iikpD5S0b7fKXof69Pz3p18kMDTXzhyrXVVOxIHth4CrKbq7em1F2ac92Lf9UK8qYntF+T/8GuDrdPnzaTKvbk8+n0Z16fXnU9bhceBGiv78nv/u1NanHz6ffpHA0NwLc6xKGksxx+qcXlVG0pblQCKStqSYqenO4a/qijnACeX6CcA1varI4JdM6b108fMpBzQbzavb9c+nWV169flI2kHStuX65sBhwD306HenWX16+fvTb/JU0jDKx9W+yotzrH6ph3XZnaKVAMU8Gv/a7fpI+g5wCEUa4RXAZ4AfAlcCu1BMRv5+25UPCjepyyEU3QAG7gc+MtiH3YX6vBn4KbAAGJyZ5RyKvv2ufj7D1GUGPfh8JO1HMbg8muKP0Sttf17Sy+jN706z+vxfevT7028SGCIiYoh0JUVExBAJDBERMUQCQ0REDJHAEBERQyQwRETEEAkMsVGQ9FTrsyICEhgiIqJOAkNsVFT4B0l3qpjb4phy/yGSbpT0fUn3SLq8fIM4YqMzptcViOiy91G83bo/xVvTt0i6qTz2GmAfipxYPwfeRJGrP2KjkhZDbGzeDHynTJa2Avhv4HXlsZttLy2TqM0HdutRHSN6KoEhNjbDdQ89V7M+QFrUsZFKYIiNzU3AMeVELTsABwM397hOEX0lfxHFxuZq4E+A2ymyaH7S9u8l7dXbakX0j2RXjYiIIdKVFBERQyQwRETEEAkMERExRAJDREQMkcAQERFDJDBERMQQCQwRETHE/wdSmlyhmja4/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stats.isel(c=0).plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
	{
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"# Correlation coeff with 3D arrays."
	]
	},
	{
	"cell_type": "code",
	"execution_count": 4,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"<xarray.DataArray (c: 2, time: 10, z: 100)>\n",
	"array([[[250.31708, 250.31708, ..., 255.19339, 253.62856],\n",
	" [247.862 , 247.862 , ..., 249.95056, 248.54636],\n",
	" ...,\n",
	" [227.74347, 227.74347, ..., 241.1389 , 241.71703],\n",
	" [238.68315, 238.68315, ..., 247.05644, 245.7788 ]],\n",
	"\n",
	" [[242.73505, 242.73505, ..., 243.99998, 242.34166],\n",
	" [230.58101, 230.58101, ..., 241.5458 , 240.07744],\n",
	" ...,\n",
	" [209.07805, 209.07805, ..., 226.06282, 224.69025],\n",
	" [214.15959, 214.15959, ..., 232.0684 , 231.35303]]], dtype=float32)\n",
	"Coordinates:\n",
	" * time (time) object 1850-01-16 12:00:00 ... 1850-10-16 12:00:00\n",
	" height float64 2.0\n",
	" * z (z) MultiIndex\n",
	" - lat (z) float64 -90.0 -90.0 -90.0 -90.0 ... -86.21 -86.21 -86.21 -86.21\n",
	" - lon (z) float64 0.0 3.75 7.5 11.25 15.0 ... 18.75 22.5 26.25 30.0 33.75\n",
	"Dimensions without coordinates: c\n",
	"Attributes:\n",
	" standard_name: air_temperature\n",
	" units: K\n",
	" cell_measures: area: areacella\n",
	" associated_files: baseURL: http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation..."
	]
	},
	"execution_count": 4,
	"metadata": {},
	"output_type": "execute_result"
	}
	],
	"source": [
	"# Take 2 files from CMIP5 for testing\n",
	"data = xr.open_dataset('/g/data/al33/replicas/CMIP5/combined/IPSL/IPSL-CM5A-LR/historical/mon/atmos/Amon/r1i1p1/v20110406/tasmax/tasmax_Amon_IPSL-CM5A-LR_historical_r1i1p1_185001-200512.nc')\n",
	"tmax=data.tasmax\n",
	"\n",
	"data2 = xr.open_dataset('/g/data/al33/replicas/CMIP5/combined/IPSL/IPSL-CM5A-LR/historical/mon/atmos/Amon/r1i1p1/v20110406/tasmin/tasmin_Amon_IPSL-CM5A-LR_historical_r1i1p1_185001-200512.nc')\n",
	"tmin = data2.tasmin\n",
	"\n",
	"# We need to create 1 dataarray with all the data in. And we reduce the size to a (10,10,10) cube for testing.\n",
	"# To create a dataarray, we add a dimension ('c') and concat along this dimension.\n",
	"# We also need to transform our 3D arrays into 2D arrays. We use the stack() method of xarray for that. (New MultiIndex dimension z)\n",
	"sample1=tmax[0:10,0:10,0:10]\n",
	"sample1=sample1.stack(z=('lat','lon'))\n",
	"sample1= sample1.expand_dims(dim='c',axis=0)\n",
	"sample2=tmin[0:10,0:10,0:10]\n",
	"sample2=sample2.stack(z=('lat','lon'))\n",
	"sample2= sample2.expand_dims(dim='c',axis=0)\n",
	"sample = xr.concat((sample1, sample2),dim='c')\n",
	"sample"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 16,
	"metadata": {},
	"outputs": [],
	"source": [
	"# Now we can simply use groupby and scipy.stats.pearsonr to calculate what we need\n",
	"# The squeeze() is because the stacking means arrays have degenerated dimensions.\n",
	"def corr2(ar,axis):\n",
	" ar1=ar[0,...].squeeze()\n",
	" ar2=ar[1,...].squeeze()\n",
	" stats = xr.DataArray(np.asarray(scipy.stats.pearsonr(ar1,ar2)),dims='c',coords={'c':['file1','file2']})\n",
	" return stats"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 17,
	"metadata": {},
	"outputs": [],
	"source": [
	"stats = sample.groupby('z').reduce(corr2,dim='time')\n",
	"stats=stats.unstack()"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 18,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"<matplotlib.collections.QuadMesh at 0x7fa5a9affc88>"
	]
	},
	"execution_count": 18,
	"metadata": {},
	"output_type": "execute_result"
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEWCAYAAABi5jCmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xu4XVV97vHvm4RwC5cgiCHhEjACASFKRKuCiKCgtaitQFoBqYhY8CAPjwq0Hu8eSqFKH6hpjCgcqYgiklo0BVqKelQIEggBNCFQCYmBcBUQQnbe88ecG9ZeWWuvlWTNtdYm7+d55pN5H2MvNuu3xxhz/oZsExERMWhUrysQERH9JYEhIiKGSGCIiIghEhgiImKIBIaIiBgigSEiIoZIYIiuknS/pMPW89qFkg6pupyIjV0CQ4wYtvexfeOG3kfSIZKWdqBKje59vqRFkv4g6R5Jx7c4/y8l/Y+kpyX9UNJ2VdQrYl0kMER01tPAu4FtgBOACyW9sdGJkvYB/gU4DtgReAb45y7VM6KpBIbohWmS7pD0hKTvStps8ICkP5U0X9Ljkv6fpP1qjr3QPSRpc0mXSnpM0t2SPtmgFbBWOZK2BH4M7CTpqXLZqVM/mO3P2L7H9hrbvwJ+CvxJk9P/Cvg32zfZfgr4NPA+SVt1qj4R6yOBIXrhaOAIYDKwH/BBAEmvBS4BPgK8jOKv6TmSNm1wj88AuwG7A4cDH2inHNtPA0cCy2yPK5dl9RdKOqsMTg2Xdn5ISZsDrwMWNjllH+D2wQ3b9wKrgFe1c/+IqiQwRC/8k+1lth8F/g2YVu7/MPAvtn9le8D2pcBzwBsa3ONo4Mu2H7O9FPindSinJdvn2t622dLmbWZSfPHPbXJ8HPBE3b4ngLQYoqcSGKIXfl+z/gzFFyTArsCZdX+Z7ww06urZCXigZvuBBuc0K6dykv4B2Bc42s0zVT4FbF23b2vgD1XWLaKVBIboJw8AX6r763wL299pcO5yYFLN9s7rUE7LlMKSzqkZg1hraXHt5yi6q95u+8lhTl0I7F9z3e7ApsBv2/sxIqqRwBD95OvAKZJer8KWkt7VZDD2SuBsSeMlTQROW4dyVgAvk7RNsxNsf7lmDGKtpdl1ks4G/hI43PYjLepxOfBuSQeVg+KfB35gOy2G6KkEhugbtudRjDNcBDwGLKYcmG7g88BS4D7geuD7FOMR7ZRzD/AdYEnZZdWxp5KALwO7AItqWhjnDB4stw8q67EQOIUiQDxEMbbwNx2sS8R6USbqiZcCSR8FjrX9ll7XJWKkS4shRiRJEyS9SdIoSXsCZwJX97peES8FY3pdgYj1NJbiPYfJwOPAFeSt4YiOSFdSREQMka6kiIgYYqPqStpy/FiPn7h55eU8/EzX3qNCq7oX28dssbprZe246XCP/3fOuFHPd6UcgBWru/d7MZo1XSurm30Oz7t7X1kr735kpe0dNuQe73jrln7k0YGW5916x3NzbR+xIWV10kYVGMZP3JxTr3xT5eXMvuPNlZcxSA9s1vqkDtl+v4e7VtaZr/yPrpRz0GbLu1IOwPkrD+paWVuP/mPXylqDulbWsmfbzUay4b7+usv+Z0Pv8cijA9w8d5eW542esGj7DS2rkzaqwBAR0U0G1nSx9dYpPRljKFMgzy+X+yXNL/cfLulWSQvKfw9tcv1nJT1Yc493dvcniIhozZjnPdBy6Tc9aTHYPmZwXdIFvJhhciXwbtvLJO1LkZVyYpPbfMX2+dXWNCJiw4zEFkNPu5IkiSJ98qEAtm+rObwQ2EzSprbbSnUQEdFPjBkYga8E9Ppx1YOAFbYXNTj258BtwwSF08rZuS6RNL5ZAZJOljRP0rynH13ViTpHRLRtDW659JvKAoOk6yXd2WA5qua0GRTJzOqv3Qf4e4qZvBr5GrAHxcQry4ELmtXD9izb021P33K7sev980RErCsDA7jl0m8q60qyfdhwxyWNAd4HHFC3fxJFzpvjy6kOG917Rc35Xwd+tMEVjoioQD+2CFrp5RjDYcA95bSMAEjaFvh34GzbP292oaQJtgcfQH8vcGelNY2IWA8Gns8Ywzo5lrW7kU4DXgl8uuZR1JcDSJotaXp53nnlI613AG8FzuharSMi2uQ2upE2qq6kVmx/sMG+LwJfbHL+STXrx1VXs4iIDjEM9N/3fkt58zkioiLFm88jTwJDRERlxEAXc0l1SgJDRERFisHnBIaIiCgV7zEkMERERI01aTFERMSgtBgiImIIIwZ6npJu3SUwRERUaCR2JY28UBYRMUIYscqjWy7tkHSEpN9IWizprAbHx0u6usw6fXM5p83gsW0lfV/SPZLulvQnw5WVwBARUZHiBbdRLZdWJI0GLgaOBKYCMyRNrTvtHGC+7f2A44ELa45dCPzE9l7A/sDdw5WXwBARUaGB8iW34ZY2HAgstr3E9irgCuCounOmAjcA2L4H2E3SjpK2Bg4GvlEeW2X78eEK26jGGB55dkv+dfHrKi9nj6+urryMQavGd6+s/zzuqq6V9cWVe3WlnE9cP6Mr5QBM+eivulbW8jPf2rWyVjWdJquCsrbrv/mRh2OLAbf19/f2kubVbM+yPatmeyLwQM32UuD1dfe4nWIqg59JOhDYFZgEDAAPA9+UtD9wK3C67aebVSYthoiICq1BLRdg5eCEYuUyq+42jZoV9en5zgXGS5oPfAy4DVhN0QB4LfA1268BngbWGqOotVG1GCIiuqkYfO7I1+xSYOea7UnAsiFl2U8CJwJIEnBfuWwBLLU92GT9Pi0CQ1oMEREV6dTgM3ALMEXSZEljKeazmVN7Qvnk0eD8xScBN9l+0vbvgQck7Vkeextw13CFpcUQEVGhgQ68x2B7taTTgLnAaOAS2wslnVIenwnsDVwmaYDii/9DNbf4GHB5GTiWULYsmklgiIioSCfffLZ9LXBt3b6ZNeu/AKY0uXY+ML3RsUYSGCIiKrSmvaeS+koCQ0RERYokegkMERFRMuL5NlNe9JOeBAZJ3wUGR8i3BR63PU3SbhSvav+mPPZL26c0uH474LvAbsD9wNG2H6u21hER68am3Rfc+kpPAoPtYwbXJV0APFFz+F7b01rc4izgBtvnlsmkzgI+1fmaRkRsiBdeYBtRehrKypcwjga+s46XHgVcWq5fCrynk/WKiOgEU7QYWi39ptc1OghYYXtRzb7Jkm6T9N+SDmpy3Y62lwOU/768WQGSTpY0T9K8gSef6VzNIyLaMMColku/qawrSdL1wCsaHPpb29eU6zMY2lpYDuxi+xFJBwA/lLRP+ar3eilzjswC2PyVO9XnFomIqIzRiJyop7LAYPuw4Y5LGkORCfCAmmueA54r12+VdC/wKmBe3eUrJE2wvVzSBOChjlY+IqIDDDzfmVxJXdXLNsxhwD22lw7ukLRDOSEFknaneItvSYNr5wAnlOsnANc0OCciosdaz8XQ5nwMXdXLwHAsaw86HwzcIel2igyAp9h+FEDSbEmDr3SfCxwuaRFweLkdEdFXTPHmc6ul3/SsjWP7gw32XQU0nA3G9kk1649QZAiMiOhr/dgiaGXkdX5FRIwQtvqyRdBKAkNEREWKweekxIiIiBe0PedzX9moAsPYMQPssm31KZU+ePnPKy9j0NzHXt21svb6+XFdK2vc5s91pZwtdnyqK+UAPHfdbl0ra6vZa7pW1rI9u1fW2IdH1ldWMficMYaIiKjRj282tzLyahwRMUIMvvncammHpCMk/UbS4jJ5aP3x8ZKulnSHpJsl7Vt3fHSZbuhHrcpKYIiIqNAaRrVcWilf/L0YOBKYCsyQNLXutHOA+bb3A44HLqw7fjrFtAYtJTBERFTEhufXjGq5tOFAYLHtJbZXAVdQZJmuNRW4oSjX9wC7SdoRQNIk4F3A7HYKS2CIiKhI0ZXUkTefJwIP1GwvLffVup0i/xySDgR2BSaVx74KfBJo60mBBIaIiAq1mStp+8HpAcrl5LrbNBqIqM8WfS4wXtJ84GPAbcBqSX8KPGT71nbrnKeSIiIqsg6Pq660PX2Y40uBnWu2JwHLhpRVTE9wIrwwCdp95XIs8GeS3glsBmwt6du2P9CssLQYIiIq07GupFuAKZImSxpL8WU/Z0hJ0rblMYCTgJtsP2n7bNuTbO9WXvefwwUFSIshIqJSnZjz2fZqSacBc4HRwCW2F0o6pTw+E9gbuEzSAHAX8KH1LS+BISKiIsVTSZ3JlWT7WuDaun0za9Z/QTGHzXD3uBG4sVVZCQwRERXJ1J4REbGWTnQldVsCQ0RERZJELyIi1pKJeiIi4gW2WD0CA0NPaizpu5Lml8v95Zt6SPqrmv3zJa2RNK3B9Z+V9GDNee/s/k8REdFap7KrdlNPWgy2jxlcl3QB8ES5/3Lg8nL/q4FrbM9vcpuv2D6/6rpGRKyvjDGsh/K17aOBQxscngF8p7s1iojorJEYGHrd+XUQsML2ogbHjmH4wHBaOSHFJZLGNztJ0smDialWPf7HDa1vRETbOjlRTzdVFhgkXS/pzgZLbQ7xhq0CSa8HnrF9Z5Pbfw3YA5gGLAcuaFYP27NsT7c9fey2m6//DxQRsR7WoJZLv6msK8n2YcMdlzSGInf4AQ0OH8swrQXbK2ru83Wg5VR1ERHdZsPq9ibi6Su9HGM4DLjH9tLanZJGAe8HDm52oaQJtpeXm+8FmrUshnh21RjuXvqK9axu+85Z8d7Kyxg0sLp7v3RrVnUm50s7dn3F77tSzt3XvbIr5QD8fuw2XStrzB5dK4pRz3avrGlv+W3Xylrcofv0Y1dRK70MDM1aBQcDS20vqd0paTYw0/Y84LzyMVYD9wMfqbiuERHrLLmS1pHtDzbZfyPwhgb7T6pZP66yikVEdJATGCIiolY/Di63ksAQEVERO2MMERExhBjIU0kREVFrJI4xjLxQFhExQgzmSurEm8+SjpD0G0mLJZ3V4Ph4SVeXGSFulrRvuX9nSf8l6W5JCyWd3qqsBIaIiKq4GGdotbQiaTRwMXAkMBWYIWlq3WnnAPNt7wccD1xY7l8NnGl7b4onPk9tcO0QCQwRERXqUEqMA4HFtpfYXgVcARxVd85U4AYA2/cAu0na0fZy278u9/8BuBuYOFxhCQwRERVxOfjcagG2H0z2WS4n191qIvBAzfZS1v5yv50izRCSDgR2BSbVniBpN+A1wK+Gq3cGnyMiKtROVxGw0vb0YY43albU3/lc4MJy4rMFwG0U3UjFDaRxwFXAx20/OVxlEhgiIirUoaeSlgI712xPApYNLcdPAifCC3Pd3FcuSNqEIihcbvsHrQpLV1JEREWKwWW1XNpwCzBF0mRJYylyzc2pPUHStuUxgJOAm2w/WQaJbwB32/7HdgpLiyEiokKdePPZ9mpJpwFzgdHAJbYXSjqlPD4T2Bu4TNIAcBfwofLyNwHHAQvKbiaAc2xf26y8BIaIiAq1OcbQxn18LXBt3b6ZNeu/AKY0uO5nNB6jaCqBISKiIkasSUqMiIio1aEGQ1clMEREVMUjM1dSAkNERJVGYJMhgSEiokJpMURExAsMrFmTwNDXRv1xFJst3Lzycp555fOVlzFoyuTlXStr3CarulbWKHWn/X3sX9zYlXIA/vuhV3atrE3O2qZrZT08fauulXXrmu59hh1hYAS2GHryHJWkaZJ+KWl+mTDqwJpjZ5f5xn8j6R1Nrt9O0nWSFpX/ju9e7SMi2teJtNvd1qsHbM8DPmd7GvC/y23KHOHHAvsARwD/XOYhr3cWcIPtKRRpZteatCIioi+4jaXP9CowGNi6XN+GF5NBHQVcYfs52/cBiynykNc7Cri0XL8UeE+FdY2IWE+t8yT14+B0r8YYPg7MlXQ+RXB6Y7l/IvDLmvMa5RwH2NH2cgDbyyW9vFlBZV7zkwHGbJ0ep4josj5sEbRSWWCQdD3wigaH/hZ4G3CG7askHU2R+e8w2ss5vk5szwJmAWw+YecR+J8oIkYsg/NU0otsH9bsmKTLgMEJqb8HzC7XW+YcL62QNKFsLUwAHupAlSMiKjDyAkOvxhiWAW8p1w8FFpXrc4BjJW0qaTJFpsCbG1w/BzihXD8BuKbCukZErL8ROPjcqzGGD1NMQTcGeJZyDKDML34lRS7x1cCptgcAJM0GZtqeRzGF3ZWSPgT8Dnh/D36GiIjW+vCLv5WeBIYyP/gBTY59CfhSg/0n1aw/QjFOERHRv0boC24b1ZvPERHd1o8vsLWSwBARUaUR+FTSyJtaKCJiBJFbL23dRzqiTBW0WNJa2R4kjZd0taQ7JN0sad92r62XwBARUZV2nkhqIzCUqYEuBo4EpgIzyhRCtc4B5tveDzgeuHAdrh0igSEiojIqBp9bLa0dCCy2vcT2KuAKitRAtaZS5I7D9j3AbpJ2bPPaIRIYIiKq1F6LYfsy0/TgcnLdXSYCD9RsN0oXdDvwPoAyY/WuFC8Jt3PtEBl8joio0pq2zlppe/owx9tJF3Quxfth84EFwG0U74Otc6qhtloMkv6+nX0REVFj8D2GDe9KapkuyPaTtk8spzM4HtgBuK+da+u125V0eIN9R7Z5bUTERqtDTyXdAkyRNFnSWIp5a+YMKUfatjwGcBJwk+0n27m23rBdSZI+CvwNsLukO2oObQX8vK0fJyJiY9aBF9xsr5Z0GjAXGA1cUqYQOqU8PhPYG7hM0gBFWqEPDXftcOW1GmP4V+DHwP9h6Cxpf7D96Dr/dBERsV5sXwtcW7dvZs36LygSj7Z17XCGDQy2nwCeAGYAlBPibAaMkzTO9u/aLSgiYmPU7gts/aStp5IkvRv4R2AnirkPdgXuppibecTQuAHGvOGxysuZvOXTlZcx6AMTf9W1snbf5KU37cWiVY3mkqrGx3bt3t9RE7/3eNfK+sTiv+haWU8/tnXrk/qJeUmnxPgi8Abgt7YnU2Q2zRhDREQrI3A+hnYDw/NlqutRkkbZ/i9gWoX1ioh4SehUrqRuavcFt8cljQNuAi6X9BDFixMRETGcPvzib6XdFsNRwB+BM4CfAPcC766qUhERLxkjsCuprRaD7drR1EsrqktExEtKv3YVtdLqBbc/0DieCbDtEfaIQEREl43Ap5JavcewVbcqEhHxUjQSWww9SbstaZqkX0qaX6aYPbDcf7ikWyUtKP89tMn1n5X0YHn9fEnv7O5PEBHRppfqGEMFzgM+Z/vH5Zf6ecAhwErg3baXldPSzaV53vCv2D6/K7WNiFgfL8UxhgoZGByf2IYyBazt22rOWQhsJmlT2891uX4REZ2RwNC2jwNzJZ1P0Z31xgbn/Dlw2zBB4TRJxwPzgDNtN8x1Uc6EdDLAJjtss8EVj4hYF2pvop6+UtkYg6TrJd3ZYDkK+Chwhu2dKd6N+EbdtfsAfw98pMntvwbsQfH29XLggmb1sD3L9nTb08dss0UHfrKIiJe2yloMtg9rdkzSZcDp5eb3gNk1xyYBVwPH2763yb1X1Jz/deBHnahzRETHjcCupJ48lUQxpvCWcv1QYBEUMxAB/w6cbbtpkj5JE2o23wvcWVE9IyLWXxt5kvpxcLpXYwwfppi0egzwLOUYAHAa8Erg05I+Xe57u+2HJM0GZtqeB5wnaRpFLL6f5l1OERG91Ydf/K30JDDY/hlwQIP9X6RI8d3ompNq1o+rrnYRER3UocAg6QjgQorpOWfbPrfu+DbAt4FdKL7bz7f9zfLYGRTzQBtYAJxo+9lmZfWqKyki4iVPFE8ltVpa3kcaDVwMHAlMBWZImlp32qnAXbb3p3gv7AJJYyVNBP4XMN32vhSB5djhyktgiIioSufGGA4EFtteYnsVcAVF1uu60thKkoBxwKO8OD3CGGDzsvt+C8p3x5pJYIiIqFJnUmJMBB6o2V7K2lkhLgL2pvjSXwCcbnuN7QeB84HfUTze/4Tt/xiusASGiIgqtRcYti/zxg0uJ9fdpVGK1vqQ8g5gPrATxTteF0naWtJ4itbF5PLYlpI+MFyVe/VUUk8MDIinn96s8nJWD3Qv3n77wdd3raztN3u69Ukdsue4Fa1P6oCHVr00EwiPH/NM18radVzDpAOVOHLCwq6VdU6H7tNmV9FK29OHOb4U2LlmexJrdwedCJxr28BiSfcBewG7AvfZfhhA0g8osk18u1lhaTFERFSpM11JtwBTJE2WNJZi8HhO3Tm/A94GIGlHYE9gSbn/DZK2KMcf3gbcPVxhG1WLISKiq9yZXEm2V0s6jSLj9GjgEtsLJZ1SHp8JfAH4lqQFFF1Pn7K9Elgp6fvArykGo28DZg1XXgJDRESVOvQeg+1rgWvr9s2sWV8GvL3JtZ8BPtNuWQkMEREV6seUF60kMEREVCmBISIiXtCnU3e2ksAQEVERka6kiIiok8AQERFDJTBERMQQCQwREfGCPp2hrZUEhoiIKiUwRERErU6kxOi2BIaIiAqNxK6knmRXlTRN0i8lzS9zjx9Y7t9N0h/L/fMlzWxy/XaSrpO0qPx3fHd/goiINrSTWbUPA0ev0m6fB3zO9jTgf5fbg+61Pa1cTmly/VnADbanADeU2xER/SeBoW0Gti7Xt6HF/KMNHAVcWq5fCrynQ/WKiOiYwTefOzDnc1f1aozh48BcSedTBKc31hybLOk24Eng72z/tMH1O9peDmB7uaSXV17jiIj1oDV9+M3fQmWBQdL1wCsaHPpbihmEzrB9laSjgW8Ah1FMVL2L7UckHQD8UNI+tp/cgHqcDJwMMGb7bdb3NhER665Pu4paqSww2D6s2TFJlwGnl5vfA2aX1zwHPFeu3yrpXuBVwLy6W6yQNKFsLUwAHhqmHrMoZyvadPeJI/A/UUSMZP3YVdRKr8YYlgFvKdcPBRYBSNpB0uhyfXdgCsWcpfXmACeU6ycA11Ra24iI9TUCB597NcbwYeBCSWOAZym7eoCDgc9LWg0MAKfYfhRA0mxgpu15wLnAlZI+RDHR9fvbKVSrRjHqvs07+5M0sGpU9WUMWjJ669Yndci9o7tWFL942R5dKqkP/6/sBKt7ZT3dvV+Mn27+qq6VBT/uyF061WKQdARwIcWcz7Ntn1t3fBvg28AuFN/t59v+ZnlsW4qemX0pfun/2vYvmpXVk8Bg+2fAAQ32XwVc1eSak2rWH6EYp4iI6G8dCAxlT8rFwOHAUuAWSXNs31Vz2qnAXbbfLWkH4DeSLre9iiKg/MT2X0gaC2wxXHl58zkioiruWEqMA4HFtpcASLqC4rH92sBgYCtJAsYBjwKrJW1N0RvzQYAyUKwarrBejTFERLzkrcN7DNuXWSAGl5PrbjUReKBme2m5r9ZFwN4UY7gLgNNtrwF2Bx4GvinpNkmzJW05XL0TGCIiqmS3XmCl7ek1y6y6uzQaNKrvpHoHMB/YCZgGXFS2FsYArwW+Zvs1wNO0yBaRwBARUaEOvfm8FNi5ZnsSa2eMOBH4gQuLgfuAvcprl9r+VXne9ykCRVMJDBERVelcEr1bgCmSJpeDx8dSPLZf63eUD+VI2hHYE1hi+/fAA5L2LM97G0PHJtaSweeIiAp1YvDZ9mpJpwFzKR5XvcT2QkmnlMdnAl8AviVpAUXX06dsryxv8THg8jKoLKFoXTSVwBARUaFOTdRj+1rg2rp9M2vWlwFvb3LtfGB6u2UlMEREVMUMDi6PKAkMEREVGom5khIYIiKqlMAQERGDBl9wG2kSGCIiqmJnop6IiKgz8uJCAkNERJXSlRQRES8ykK6kiIgYYuTFhQSGiIgqpSspIiKGyFNJERHxovazp/aVBIaIiIoUL7iNvMjQk8AgaRowE9gMWA38je2bJf0V8ImaU/cDXltmBqy9/rPAhymmqwM4p8w8ODzDqOc3vP6tdKOMQXKjiZ0qKqubP9fysd0pqIv/z3Yqy2ZbZXXx5xrYtJtlje5eYZ3Sxf/undKrFsN5wOds/1jSO8vtQ2xfDlwOIOnVwDX1QaHGV2yf353qRkSsn7QY2mdg63J9G9aeog5gBvCdrtUoIqLTMsawTj4OzJV0PsX0om9scM4xwFHD3OM0SccD84AzbT/W6CRJJwMnA4zZZvwGVToiYt2MzFxJlc35LOl6SXc2WI4CPgqcYXtn4AzgG3XXvh54xvadTW7/NWAPYBqwHLigWT1sz7I93fb00Vtu2YkfLSKifXbrpc9U1mKwfVizY5IuA04vN78HzK475ViG6UayvaLmXl8HfrT+NY2IqIg799CBpCOACynmfJ5t+9y649sA3wZ2ofhuP9/2N2uOj6boYXnQ9p8OV1ZlLYYWlgFvKdcPBRYNHpA0Cng/cEWziyVNqNl8L9CsZRER0VsdaDGUX+oXA0cCU4EZkqbWnXYqcJft/YFDgAsk1T7edzpwdztV7lVg+DBFpW8Hvkw5BlA6GFhqe0ntBZJmSxqczPo8SQsk3QG8laI7KiKi/7iNpbUDgcW2l9heRfGHc/0YrIGtJAkYBzxK8ToAkiYB72Lt3pmGejL4bPtnwAFNjt0IvKHB/pNq1o+rrHIRER2kNW31JW0vaV7N9izbs2q2JwIP1GwvBV5fd4+LgDkUPTJbAcfYHiz8q8Any/0t5c3niIiqmHZfcFtpe/owxxu9yVrf1ngHMJ+ie34P4DpJP6XohXnI9q2SDmmnMr3qSoqIeMkTRm69tGEpsHPN9iTWfv/rROAHLiwG7gP2At4E/Jmk+ym6oA6V9O3hCktgiIioUmceV70FmCJpcjmgfCxFt1Gt3wFvA5C0I7AnsMT22bYn2d6tvO4/bX9guMLSlRQRUaUOvKdge7Wk04C5FI+rXmJ7oaRTyuMzgS8A35K0gKLr6VO2V65PeQkMERFVaX+MofWtikSh19btm1mzvgx4e4t73Ajc2KqsBIaIiAq1+VRSX0lgiIioTH+mvGglgSEioiomgSEiIuqMvJ6kBIaIiCplop6IiBgqgSEiIl5gw8DI60tKYIiIqFJaDP3Nm5hnd1rd62pExMYkgSEiIl5gYATO+ZzAEBFRGYMzxhAREYNMBp8jIqJOxhgiImKIBIaIiHjRyEyi15MZ3CTtL+kXkhZI+jdJW9ccO1vSYkm/kfSOJtdvJ+k6SYtNl/cLAAAHTklEQVTKf8d3r/YREW0ysGZN66XP9Gpqz9nAWbZfDVwNfAJA0lSKqef2AY4A/lnS6AbXnwXcYHsKcEO5HRHRfzoztWdX9Sow7AncVK5fB/x5uX4UcIXt52zfBywGDmxw/VHApeX6pcB7KqxrRMR6KlNitFraIOmIsidlsaS1/hiWtE3ZA3O7pIWSTiz37yzpvyTdXe4/vVVZvQoMdwJ/Vq6/H9i5XJ8IPFBz3tJyX70dbS8HKP99ebOCJJ0saZ6keQN/eHqDKx4R0TaDvabl0krZc3IxcCQwFZhR9rDUOhW4y/b+wCHABZLGAquBM23vDbwBOLXBtUNUFhgkXS/pzgbLUcBfl5W7FdgKWDV4WYNbbVA7y/Ys29NtTx+91ZYbcquIiHW3xq2X1g4EFtteYnsVcAVFz0ktA1tJEjAOeBRYbXu57V8D2P4DcDeN/+B+QWVPJdk+rMUpbweQ9CrgXeW+pbzYegCYBCxrcO0KSRNsL5c0AXhoQ+sbEVGJ9sYQtpc0r2Z7lu1ZNduNelNeX3ePi4A5FN+ZWwHHuK45Imk34DXAr4arTK+eSnp5+e8o4O+AmeWhOcCxkjaVNBmYAtzc4BZzgBPK9ROAa6qtcUTEerDbfSpp5WDPRrnMqrtTO70p7wDmAzsB04CL6p74HAdcBXzc9pPDVbtXYwwzJP0WuIciun0TwPZC4ErgLuAnwKm2BwAkzZY0vbz+XOBwSYuAw8vtiIj+05mnktrpTTkR+IELi4H7gL0AJG1CERQut/2DVoX15AU32xcCFzY59iXgSw32n1Sz/gjwtsoqGBHREcYDA5240S3AlLIn5UGKx/r/su6c31F8L/5U0o4UT38uKcccvgHcbfsf2yksbz5HRFSlQ2m3ba+WdBowFxgNXGJ7oaRTyuMzgS8A35K0gKLr6VO2V0p6M3AcsEDS/PKW59i+tll5CQwREVXqUNrt8ov82rp9M2vWl1E+1FN3zs9oPEbRVAJDRERFDDgT9URExAuciXoiIqJOhwafu0ruwwROVZH0MPA/63jZ9sDKCqqzvlKf5vqpLpD6tNJP9WlUl11t77AhN5X0k/Leray0fcSGlNVJG1VgWB+S5tme3vrM7kh9muunukDq00o/1aef6tIPevWCW0RE9KkEhoiIGCKBobX6nCW9lvo01091gdSnlX6qTz/VpecyxhAREUOkxRAREUMkMERExBAJDMNoNcdqD+pzv6QFkubXTerRrfIvkfSQpDtr9m0n6TpJi8p/x/ewLp+V9GD5+cyX9M5u1KUsu+G8ur34fIapS08+H0mbSbq5Zi7iz5X7e/W706w+Pfv96TcZY2iinGP1txTzPSylSHs7w/ZdPazT/cB02z15KUjSwcBTwGW29y33nQc8avvcMniOt/2pHtXls8BTts+vuvwG9ZkATLD9a0lbAbcC7wE+SJc/n2HqcjQ9+HzKtM9b2n6qnBfgZ8DpwPvoze9Os/ocQY9+f/pNWgzNtTPH6kbF9k0U88jWOgq4tFy/lOILqFd16Zlh5tXt+uezPnP8Vlwf236q3NykXEzvfnea1SdKCQzNNZpjtWf/c5UM/IekWyWd3OO6DNrR9nIovpCAl/e4PqdJuqPsaupK10S9unl1e/r5NJjjtyefj6TR5VwADwHX2e7pZ9OkPtAHvz/9IIGhuXbmWO22N9l+LXAkcGrZnRIv+hqwB8V8t8uBC7pdgXWZV7cHdenZ52N7wPY0iikpD5S0b7fKXof69Pz3p18kMDTXzhyrXVVOxIHth4CrKbq7em1F2ac92Lf9UK8qYntF+T/8GuDrdPnzaTKvbk8+n0Z16fXnU9bhceBGiv78nv/u1NanHz6ffpHA0NwLc6xKGksxx+qcXlVG0pblQCKStqSYqenO4a/qijnACeX6CcA1varI4JdM6b108fMpBzQbzavb9c+nWV169flI2kHStuX65sBhwD306HenWX16+fvTb/JU0jDKx9W+yotzrH6ph3XZnaKVAMU8Gv/a7fpI+g5wCEUa4RXAZ4AfAlcCu1BMRv5+25UPCjepyyEU3QAG7gc+MtiH3YX6vBn4KbAAGJyZ5RyKvv2ufj7D1GUGPfh8JO1HMbg8muKP0Sttf17Sy+jN706z+vxfevT7028SGCIiYoh0JUVExBAJDBERMUQCQ0REDJHAEBERQyQwRETEEAkMsVGQ9FTrsyICEhgiIqJOAkNsVFT4B0l3qpjb4phy/yGSbpT0fUn3SLq8fIM4YqMzptcViOiy91G83bo/xVvTt0i6qTz2GmAfipxYPwfeRJGrP2KjkhZDbGzeDHynTJa2Avhv4HXlsZttLy2TqM0HdutRHSN6KoEhNjbDdQ89V7M+QFrUsZFKYIiNzU3AMeVELTsABwM397hOEX0lfxHFxuZq4E+A2ymyaH7S9u8l7dXbakX0j2RXjYiIIdKVFBERQyQwRETEEAkMERExRAJDREQMkcAQERFDJDBERMQQCQwRETHE/wdSmlyhmja4/wAAAABJRU5ErkJggg==\n",
	"text/plain": [
	"<Figure size 432x288 with 2 Axes>"
	]
	},
	"metadata": {
	"needs_background": "light"
	},
	"output_type": "display_data"
	}
	],
	"source": [
	"stats.isel(c=0).plot()"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {},
	"outputs": [],
	"source": []
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {},
	"outputs": [],
	"source": []
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {},
	"outputs": [],
	"source": []
	}
	],
	"metadata": {
	"kernelspec": {
	"display_name": "Python 3",
	"language": "python",
	"name": "python3"
	},
	"language_info": {
	"codemirror_mode": {
	"name": "ipython",
	"version": 3
	},
	"file_extension": ".py",
	"mimetype": "text/x-python",
	"name": "python",
	"nbconvert_exporter": "python",
	"pygments_lexer": "ipython3",
	"version": "3.6.7"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 2
	}