ccarouge/Correlation coef.ipynb

## Correlation coef.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import scipy.stats as stats\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "DMI = xr.open_dataset(\"/g/data/w97/sl7808/DMI/9Dec/all20_1/output/test2.nc\")\n",
    "DMI_p = DMI.sel(time=DMI.time.dt.month.isin([8,9,10]))\n",
    "x=DMI_p.tos.groupby('time.year').mean(dim=\"time\")\n",
    "y=x.where(x>0.5)\n",
    "z=y.rename({\"year\":\"time\"})\n",
    "del z[\"time\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "year_HWX = xr.open_dataset(\"/g/data/w97/sl7808/corr/hwx/2_hwx_exp1.nc\",decode_times=False)\n",
    "year_HWX[\"time\"] = y[\"year\"]\n",
    "HWA = year_HWX.HWA_EHF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask with the land-sea mask\n",
    "land_sea = xr.open_dataset(\"/g/data/w97/sl7808/ppt/landmask_AUS_2.nc\")\n",
    "topo=land_sea.topo\n",
    "# Rename dimensions to be the same as HWA\n",
    "topo = topo.rename({\"lat_v\":\"lat\",\"lon_u\":\"lon\"})\n",
    "HWA_mask = HWA.where(topo == 1,drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to select only the valid data at a point and calculate the correlation coefficient.\n",
    "# Return NaN if no valid data at the point.\n",
    "def corr_miss(data, tos):\n",
    "    '''data and tos should both be 1D arrays'''\n",
    "    # select the non-NAN:\n",
    "    mask=~np.isnan(data) & ~np.isnan(tos)\n",
    "    if (mask.any(\"time\") == False):\n",
    "        return np.nan\n",
    "    else:\n",
    "        res = stats.linregress(data[mask],tos[mask])\n",
    "        return res.rvalue**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate correlations between tos and HWA_mask:\n",
    "rvals=np.apply_along_axis(corr_miss,0,HWA_mask,z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "rvals=xr.DataArray(rvals,dims=[\"lat\",\"lon\"],coords={\"lon\":HWA_mask.lon,\"lat\":HWA_mask.lat})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.QuadMesh at 0x7f0e62ad21f0>"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEJCAYAAACT/UyFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqOUlEQVR4nO3de5wcVZ338c83M4RLCCAkyCUJCWxAgiKXEVEUZeUSUMkKusKqK16eLCqyuqsixgddlF0RFa+IUXFdBREXo4hASPbRsIoIiUJIuAYIEAKScAsIJMz09/mjakgz6ZmuM9M9ffu9edUrXVW/OnVqEn5dc+rUObJNCCGEzjGm0RUIIYQwuiLxhxBCh4nEH0IIHSYSfwghdJhI/CGE0GEi8YcQQoeJxB9CCA0maaak2yWtkPTJCvvfIWlpvlwr6eVl+1ZKulnSjZIWFzpfu/TjnzBhgqdOndroaoTAnUvvLxzr556rY03SaMstko+ZPmOXOtSk/pYsWbLW9sSRlHHUYeP8yKN9xc63dP182zMr7ZPUBdwBHAGsAm4ATrR9S1nMq4FbbT8m6Wjgs7Zfme9bCfTYXlu07t1FA0eTpM8C/wdYk2/6lO0rhjpm6tSpLF5c6MsuhLo6etKphWN7Vz+YVrgSf0l3qXBo1557p5UNXLX4zORjmoGke0daxiOP9nH9/CmFYrt2vnPCELsPAlbYvjuv28XALOD5xG/72rL464BJyRUu05SJP3eu7S81uhIhhFCJgRKFv1gnDGiGmWt7bv55V6D818RVwCuHKOt9wJUDqnK1JAPfKSt3UM2c+EMIoWkZ85yLNfUAa233DLJPFYuvFCgdRpb4X1O2+RDbqyXtCCyQdJvta4aqTDM/3D0lf5BxgaQXVQqQNFvSYkmL16xZUykkhBDqplTwvypWAZPL1icBqwcGSdoX+B4wy/Yj/dttr87/fBiYR9Z0NKSGJX5JCyUtq7DMAr4N7AHsBzwIfLlSGbbn2u6x3TNx4oie04QQQhJj+lxsqeIGYLqkaZLGAicAl5UHSJoC/Bx4l+07yraPkzS+/zNwJLCs2gkb1tRj+/AicZK+C1xe5+qEEEKyUuUWmSS2eyWdAswHuoALbC+XdHK+/3zgDGAH4DxJAL1509GLgXn5tm7gIttXVTtnU7bxS9rZdn93h7dQ4BsshDD6jtriHUnx85+9sE41GX0G+mqQ+AHyXotXDNh2ftnn9wPvr3Dc3cDLB26vpikTP/BFSfuR/WxXAv/U0NqEkODKVV+vW9kzt3lPUnzp2fWFY/uW3ppU9pixY5Pi21Et7vgboSkTv+13NboOIYQwFAPPtegLsE2Z+EMIodkZ16ypZ7RF4g8hhOEw9LVm3o/EH0IIw5G9uduaIvGHEMKwiL6KL902v7ZJ/HcsuZsjNzuhUOyY8eOTy7/q0e8mHxNCaF8GStHUE0KoN2/YkBbfW3zYZ3V1JZVdSqxLuzGwoalHvRlcJP4QQhimkqOpJ4QQOkb25m4k/hBC6BhG9EVTTwghdJZo6gkhhA5ixAanPRBvFpH4QwijZuaE2YVjr1pbdQbBhspe4Iqmnoba88DduXrxxY2uRgh1lTqs8RFj3lY41n19LCj9LLVKhaUk/VYRD3dDCKGD2KLPcccfQggdpRR3/CGE0Dmyfvxxxx9CCB3DiOfcmim0NWsdQghNoC/68YcQQudo5Td3G1ZrSW+TtFxSSVLPgH2nS1oh6XZJRzWqjiGEMJSSxxRamk0j7/iXAccB3ynfKGkGcAKwD7ALsFDSnrb7Rr+KIYRQWTzcHQbbtwJIm7SRzQIutr0euEfSCuAg4A+jW8MQQhicUbTx19CuwHVl66vybZuQNBuYDTBlypT61yyEEHI20aunEkkLgZ0q7Jpj+5eDHVZhW8UJzmzPBeYC9PT0tOgkaCHUz9PHH9zoKjyv2cfeSad4gasS24cP47BVwOSy9UnA6trUKIQQasPQskM2NGOtLwNOkLS5pGnAdOD6BtcphBA20ceYQkuzaVgDlaS3AN8AJgK/lnSj7aNsL5d0CXAL0At8KHr0hBCajVFMxJLK9jxg3iD7zgLOGt0ahRBCcSYe7oYQQodRjMcfQgidxNCUb+UWEYk/hBCGKe74Qwihg9iKO/4QGunovU8vHNt7+4qksrsnVXxxfEhX3vfV5GNC66lVP35JM4GvAV3A92x/YcD+dwCn5atPAR+wfVORYyuJxB9CG3tmQmvekbaCbCKWrhGXI6kL+BZwBNkLrDdIusz2LWVh9wCvs/2YpKPJRix4ZcFjNxH/KkIIYRiyh7sqtFRxELDC9t22NwAXkw1WufFc9rW2H8tXryMb0aDQsZVE4g8hhGGq0Zu7uwL3l60POjBl7n3AlcM8FoimnhBCGJbEN3cnSFpctj43H2QSEgamlHQYWeJ/Teqx5SLxhxDCMJWKN5qstd0zyL5CA1NK2hf4HnC07UdSjh0omnpCCGEY7Gyy9SJLFTcA0yVNkzSWbAbCy8oDJE0Bfg68y/YdKcdWMuQdv6TjqhUAPGv7igJxIYTQNozoLY28V4/tXkmnAPPJumRekA9WeXK+/3zgDGAH4Lx81sJe2z2DHVvtnNWaer4L/JLK7Uj9DgUi8YeGSumb3z11t6Sy/cijqdVpGjv+bm2jq9DWavXmbn7zfMWAbeeXfX4/8P6ix1ZTLfFfafu9QwVI+nHKCUMIoR30d+dsRUMmftvvrFZAkZgQQmg/HTBkg6RXA1PLj7H9X3WoUwghtIS2nnNX0o+APYAbgf7ZsAxE4g8hdCQbnqvBw91GKHrH3wPMsF31xYAQQugEnTD14jJgJ+DBOtYltJjDX/fvSfG9Wxa/O3pu67Q7qa2SosNgjuh6e1L8gr6f1qkmraEtm3ok/YqsSWc8cIuk64H1/fttHzvcE0t6G/BZYG/gINuL8+1TgVuB2/PQ62yfPNzzhM7Qte22hWN7V96bVPaYLbdMrU7T8D2rCseWnnm2jjVpP23bqwf4Uh3PvQw4DvhOhX132d6vjucOIYQRa8tePbYXAUg62/Zp5fsknQ0sGu6Jbd+alzPcIkIIoXGKDbnclIp+XR1RYdvRtazIANMk/VnSIkmvHSxI0mxJiyUtXrNmTR2rE0IIL2Sg12MKLc2mWhv/B4APAntIWlq2azzw+2qFS1pI9lB4oDm2fznIYQ8CU2w/IulA4BeS9rG9bmBgPqzpXICenp7ocRRCGDXt3MZ/EdmA//8BfLJs+5O2qw5gYvvw1ArZXk/+ANn2Ekl3AXsCi4c8MIQQRllbJn7bT0h6EniZ7bSuEMMkaSLwqO0+SbsD04G7R+Pcne5VJ345KX5cnerRDl75rq8Ujv3jj/6ljjUJ9dLW/fhtlyTdJGmK7ftqdWJJbwG+AUwEfi3pRttHkY32eaakXrK3hE8u8ttFGH1r903r5vjiH9xYOHbzl0xLKrtv3ZOFY7unTKoeVKb3vuJdIgHWveNVSfH1NGbyLoVjSwkjnIZMW/bjL7MzsDzvx//X/o0j6cdvex4wr8L2S4FLh1tuCCGMCrdpU0+Zf6trLUIIocUY6C01X4+dIgolftuLJL0YeEW+6XrbD9evWiGE0NxauY2/0NeVpL8HrgfeBvw98EdJb61nxUIIodnZKrQ0m6JNPXOAV/Tf5ec9bxYC/12vioUQQrNr94e7YwY07TxC8bd+QwhDOHrv05Pir7z1P+pUk5DCHfBw9ypJ84Gf5OtvJyZYbzvbLbwzKX7rPdO6RbLX1MKhuiutCyUuFQ9d91Ra2Ym2ueiPhWO799y9jjVJm4Qexb1cGtHX5g93Py7peOAQQMDcvDtmCCF0rGZsvy+i8Jy70b8+hBA2auWxeor26jlO0p2SnpC0TtKTkjYZNC2EEDqGs3b+IkuzKXrH/0Xgzf1j6IcQQmj/Xj1/iaQfQggbmfZv418s6afAL3jhnLs/r0elQu0cOfYfCsd2bbtNHWsSauWIMW9rdBWed9hRZyfF/2b+adWDWoboK7V34t8GeBo4smybgUj87WTzsUnh3avTBk3148UfC5WeSuty2b1H8dE8S6sfSio7ebL1hOlEe+9IHHE8odtqqq6t0wbaXn/wXnWqSeto6zt+2+8Zar+k023HWyUhhI6RPbhtzcRfq7cPmud3zxBCGCWlfML1akuzKdyPv4rmu7IQQqizZuyqWUStEn+LXn4IIQyPEaV2HrKhgLjjDyF0nFa9461V4v9ZjcoJVbzujV9Mit+8TvUIYTiO6Hp74dgFfT+tY01qoIUf7hZK/JK+CHweeAa4Cng58BHbPwaw/e+pJ5Z0DvBmYANwF/Ae24/n+04H3kc22fqptuenlh8yzxyzf+HYrW9Zm1R2710r0yqT0BWx+8U7JhXdd2/x0TzV1ZVUdmnDhqT4lOscs/kWSUVrs7R7tb6X7VE4tnTDLUlldy/8U1J8W6rRLb+kmcDXgC7ge7a/MGD/S4AfAAcAc2x/qWzfSuBJsnzZa7un2vmKNlAdaXsd8CZgFbAn8PGCxw5mAfBS2/sCdwCnA0iaAZwA7APMBM6TlPZ/agghjIJazMCV57dvAUcDM4AT8zxY7lHgVOBLVHaY7f2KJH0onvg3y/88BviJ7bQ3dyqwfbXt3nz1OqB/cPdZwMW219u+B1gBHDTS84UQQq3VaJC2g4AVtu+2vQG4mCwPlp3HD9u+AXiuFvUumvh/Jek2oAf4n3zqxWdrUYHce4Er88+7AveX7VuVb9uEpNmSFktavGbNmhpWJ4QQhmaDS2MKLVUUznmDVQW4WtISSbOLHFD0zd1PSjobWGe7T9LTDPhGqkTSQmCnCrvm2P5lHjMH6AUu7D+sUhUGqddcYC5AT09Pqz5gDyG0qIR+/BMkLS5bn5vnL0jIeYM4xPZqSTsCCyTdZvuaoQ4o+nB3K+BDwBRgNrALsBdw+VDH2T68SrnvJntu8Ab7+R/hKmByWdgkYHWReoYQwqgqnp7XDtH+PqKcZ3t1/ufDkuaRNR2NPPGTPU1eAry6rKI/o0riH0r+FPs04HW2ny7bdRlwkaSvkH3BTAeuH+55GqHnfV8pHDvhv5elFf6alyTWJoRQH9Uf3BZ0AzBd0jTgAbLOLYWG1ZU0Dhhj+8n885HAmdWOK5r497D9dkknAth+RkoYgrCyb5J1M1+QF3Wd7ZNtL5d0CXALWRPQh2z3jfBcbWOrOx9JOyChK2LfTi9KKrrbuyXFe4vNqgf1e+KvSWVrTPF/jqX900aVHLM0bRL60tNPVw/qlzra5rS0Ce67Hyj+76U0NuHvB/DTac8Zx4xNG/21JdSggdl2r6RTgPlk3TkvyPPgyfn+8yXtBCwmGym5JOkjZD2AJgDz8hzaDVxk+6pq5yya+DdI2pL8MiXtQdm4/MNh+2+G2HcWcNZIyg8hhLqq4Qtctq8Arhiw7fyyzw+xsedjuXVk71UlKZr4P0P24tZkSRcChwAnpZ4shBDaSju/uWt7gaQ/AQeTPYH+Z9tpr3mGEEK7adG+hIX68eft+UcDB9q+HNhKUrxUFULobC64NJmiL3CdB7wKODFff5LsFeMQQuhMJmvqKbI0maJt/K+0fYCkPwPYfkxSGz6iDyE0i5nbvjcp/qonLqhTTQbX7hOxPJcPJNTfq2ciUL9Zn1vchCWPF45VYhe3DbtsmxQ/9q6/FI4dc+vKpLL7nkkbtaNr4g6FY3sfTJsQPakey+5Kii+tH1EHtiG5L62nsm9Nq3vKzzypGyqg7rTunymjnHaNH59UdsOUmu9uvoiiTT1fB+YBO0o6C/gdkDwUcwghtBO52NJsqt7xSxoD3AN8AngDWa+ev7N9a53rFkIIzatJH9wWUTXx2y5J+rLtVwG3jUKdQgihBTTng9siijb1XC3p+BoM0xBCCO2jRbtzFn24+y/AOKBX0rNkzT22vU3dahZCCM2uCZN6EUXf3G2RR+z1MXOfOWkHbBYzRYbQ9kzL9uopOh7/ARU2PwHcWzZ9Ysg9MaN4l8vtHnkiqezN1qSNWpnS0Vi7vDip6L6d037h61p+X+HY1JEcU7oK9j31VFLZ3VPSRsTs23n7wrH+49Kksrv2KT55OkDp7vurBz1f9p5JZfctvyMpPkVqN+dGacYeO0UUbeo5j2x295vz9ZcBNwE7SDrZ9tX1qFwIITS1Fk38RR/urgT2t32g7QOB/YBlwOHAF+tTtRBCCPVQ9I7/JbaX96/YvkXS/rbvjo4+IYRO1e5NPbdL+jZwcb7+duAOSZsDadPwhBBCu2jRfvxFE/9JwAeBj5B15fwd8DGypH9YPSoWQghNzbTsiGVFu3M+I+k84HLbtw/YndZFIoQQ6uDoSaeO+jnbuqlH0rHAOcBYYJqk/YAzbR873BNLOgd4M7ABuAt4j+3HJU0FbgX6v2Cus33ycM9TC1vOfTwpft3XinfnKz36WFLZTx+Y1rVw908VH+Xyph9NSyp7y7VptzvPTBh0muVNbHv9qqSy+UvChHCpE5z/9ZmkcCdM/N619dZpdRlTtD9GpvRM8bp3P554D6e0uijh59L7SPFJ4gG6d9k5Kb5mWjTxF/2b+wxwEPA4gO0bgakjPPcC4KW29wXuAE4v23eX7f3ypaFJP4QQBtWiQzYUTfy9ttPeNKrC9tVlL39dR+UZ5EMIoSkVHZK5GZuDiib+ZZL+AeiSNF3SN4Bra1iP9wJXlq1Pk/RnSYskvXawgyTNlrRY0uI1a9bUsDohhFBAScWWJlM08X8Y2AdYD/wEWEfWw2dIkhZKWlZhmVUWMwfoBS7MNz0ITLG9P9ngcBdJqjg2gO25tnts90ycOLHgpYQQQm206h1/0V49TwNz8qUw24cPtV/Su4E3AW+ws0FlbK8n+4LB9hJJdwF7AotTzh1CCHXXhEm9iCETv6RfMcSljbBXz0zgNOB1+RdL//aJwKO2+yTtDkwH7h7ueUIIoS6a9G6+iGp3/F/K/zwO2An4cb5+Itn4PSPxTWBzYEE+7EN/t81DgTMl9QJ9wMm2Hx3huUbkpJ1/nxT/8X2Ld1tc+/L9k8qe/D9pE5zf/vV9CsfuvGhlUtmPHrZbUvz2v7u3eHBXWlfBMVN2LR78VNoIp70PPZwU33Vz8cnZtdOOSWX7/rRJ6LtfXLz8vof+klR2PXVPSvj7BEqPNChFtGPit70IQNLnbB9atutXkq4ZyYltV8yOti8FLh1J2SGEMCpaNPEXva2amDe7ACBpGhBPU0MIHa2tH+4CHwV+K6m/rX0qMLsuNQohhFbRhEm9iKK9eq6SNB14Sb7ptrz3TQghdKYmvZsvYsimnvIpF22vt31TvqyvFBNCCB2lRYdsqHbH/wNJrycbinkw3wfSuqaEEEI7aMKkXkS1xL8tsIShE3/bj5Vw2o3HJcVPPqT4BNd8Lu0Z+WaPPF09qEzXDsUnrV5zZFr3zImX35kUX0roRpkyqiRA947Ff44p9QDonlB8tFWA3jXFRwrtfu5FSWU/t9/u1YPK48cXfYwH4+7cNqlsP5wwIirQlzASbWlN2uicpfVp3ZxrQbRpU4/tqbZ3tz1tiOWg0apsCCE0DYNKxZZqJM2UdLukFZI+WWH/SyT9QdJ6SR9LObaStLdkQgghbFSDNn5JXcC3gKOBGcCJkmYMCHsUOJWNL9WmHLuJSPwhhDBctXm4exCwwvbdtjeQzW0+qzzA9sO2b2DTOc6rHltJJP4QQhimGr3AtStQ/mBwVb6tiGEdWyjxK/NOSWfk61MkRdt+CKGzFb/jn9A/d0i+lL8AW6nzTNHHxsM6tugj//PI5pP/W+BM4Emy8XReUfD4EEJoLy724Da31nbPIPtWAZPL1icBqwuWO6xjiyb+V9o+QNKfAWw/Jql4P8EWd9txZyTFv+xfzy0cW0r8venmj/wgKf4VZ3ygcOwO/3lDUtnecouk+NKzCaNWdm+WVPZTBxefKH6Ly65PKjtV9847FY4tPZTWG3qzbcYlxY99qK948JNp3VzpS5y0PoHGbZkU35UYT1pv0cHVpjvnDcD0fAy0B4ATgH+o57FFE/9z+dNjw/Nj5tfvbz2EEFpALfrx2+6VdAowH+gCLrC9XNLJ+f7zJe1ENhnVNkBJ0keAGbbXVTq22jmLJv6vA/OAHSWdBbwV+HTa5YUQQpup0Qtctq8Arhiw7fyyzw+RNeMUOraaooO0XShpCfAGsocJf2f71pQThRBCW2nScXiKqDb1Yvm76g+TTbT+/L5Gz4wVQgiNIoYey6aZVbvjX0L2nSZgCvBY/nk74D6g+BO1EEJoMwm9eppKtbF6ptnenezBwZttT7C9A/Am4OejUcEQQmhabTosc79X5BOhA2D7SkmfG8mJ8+NnkfUOehg4yfbqfN/pwPvIJls/1fb8kZxrtI19snisSmn/Kva+9p1J8dP+kDCC4l5pv8Dd9onxSfF7n160azJ4+22Syt5q4bLCsak3aaWn00ZETenmWNqwIa3sm29LCh8zto69rpX24n/KBOq9qx5IKnvM5mldi2umCZN6EUX/5tZK+rSkqZJ2kzSHkfeEPcf2vrb3Ay4H+t8KnkHWF3UfYCZwXt6VNIQQmkfB4Rqacejmoon/RLLJ1ecBvwB2zLcNm+11Zavj2PjdOQu4OJ/x6x5gBdlARCGE0Fzauakn773zz7U+ef5OwD8CTwCH5Zt3Ba4rCxt00KF8vIvZAFOmTKl19UIIYUit+nC3UOKX9BsqfG/Z/tsqxy0EKr2/Psf2L23PAebkbfqnAJ8hYdAh23OBuQA9PT1N+L0aQmhnzdiMU0TRh7vlM75sARwP9FY7yPbhBcu/CPg1WeIfyYBFIYQwOpq0GaeIok09SwZs+r2kRSM5saTptvsnbT0W6O+ucBlwkaSvALsA04H6jqoVQgjD0c6Jf8AbvGOAA6nchJPiC5L2Iutddy/QPyDRckmXALeQ/VbxIdsJQww23jM7Fo+dcl7ayBcTr02bnL33tuIToo/Zv+qMbS8w/qbNk+JXnVC8u+ikS+5NKru0YeDERLXTvesuSfG9DxT/BVVdqR3W0uJLz1X9xXxjyVunjfzpvsT/LbuLT/ye+nMZs/12aXWpQRtCK0+2XvRvovwN3l7gHrJ+9sNm+/gh9p0FnDWS8kMIoe7aPPHvbfvZ8g2S0m73QgihnTj9BcxmUbQf/7UVtv2hlhUJIYRW06ovcFUbnXMnsj70W0ran41dLbcBtqpz3UIIobk1YVIvolpTz1HASWRdKr9Stv1J4FN1qlMIIbSEZrybL2LIxG/7h8APJR1v+9JRqlMIIbSGdkz8kt5p+8fAVEn/MnC/7a9UOKzj7XbJQ8WDxyW2mD36eFJ41/YvKhxbWnp7UtmTHqs4E9yQrrz7S8UCC4aNhkOPPScpflyp+Hv8vQ89nFR2926TqweVKa1+sHBs35MJw8oOw5iu4qN5diVMWA/Q92DC/3O10qTt90VUa+rp79i7dYV9LXrJIYQwcqJNx+qx/Z3840Lbvy/fJ+mQutUqhBBagVvz/rfo717fKLgthBA6Rrt253wV8Gpg4oA2/m1IfXc8hBDaSRsP0jaWrH2/GyifZ28d8NZ6VSqEEFpBu7bxLwIWSfpP22mjZoUQQptry8Rf5mlJ55DNg/v8rMbVJmLpVL0Tik9Cvvq0LZPKnvLxZ5Li+1bcnRS/oPSzpPhOcM1lH0+KP+yoswvH9r4ibea4lEnlIW10zu4ddkgq2+vXJ8X33nVP8bokTMwOWffPK+/7auF46adJ5Vdk2v7h7oVk4+VPA/4NWAncUKc6hRBCS2jVh7tFE/8Otr8PPGd7ke33AgfXsV4hhND82nmydaB/losHJb2RbBqD9Nc2QwihTXTCRCyfl7Qt8K9k/fe3AT5Sr0qFEELTs1u2jb/onLuX5x+fAA4DkPSROtUphBBaQqv26ik+atKmNhm0LYQQOkmrPtwtPvvxplQ9ZIiDpc8Bs8gmW38YOMn2aklTgVuB/qEir7N98kjONdoWXPvp+hU+q35Fh9r4zfzTCsfu+9Fzk8reYr89k+J7t0r4X/ypxAnr/3hzUnh3yoibY0aUXkaHgTaferGSkV7xObb3tb0fcDlwRtm+u2zvly8tlfRDCB2kHXv1SHqSytUWkPbm0QC215WtjhvkPCGE0LSasRmniCHv+G2Pt71NhWW87ZE0EwEg6SxJ9wPv4IV3/NMk/VnSIkmvHeL42ZIWS1q8Zs2akVYnhBDS9PfsqbZUIWmmpNslrZD0yQr7Jenr+f6lkg4o27dS0s2SbpS0uEi1R9LUU5WkhZKWVVhmAdieY3sy2ZvBp+SHPQhMsb0/2QPkiyRtU6l823Nt99jumThxYj0vJYQQXshZr54iy1AkdQHfAo4GZgAnSpoxIOxoYHq+zAa+PWD/YXnTeE+Rqo/4rn0otg8vGHoR8GvgM7bXA+vz45dIugvYEyj0TRZCCKMhe4GrJm09BwErbN8NIOlism4ct5TFzAL+y7aB6yRtJ2ln28Xn1ixT1zv+oUiaXrZ6LNlYQEiamH8DIml3sm+4tJHGQghhNJQKLjChv1k6X2aXlbIrcH/Z+qp8GwVjDFwtacmAcgdV1zv+Kr4gaS+yH8u9QH/vnUOBMyX1An3AybYfbVAdQ6irped+NCn+6EmnJsWX/qZ4F0pvlnYfuNkeU5Pie+9amRS/oK8GI2jWWcId/9ohmmEq9V0dWPBQMYfkXeF3BBZIus32NUNVpmGJ3/bxg2y/FLh0lKsTQghpatdVcxUwuWx9Etl4aIVibPf/+bCkeWRNR0Mm/oY19YQQQmszKhVbqrgBmC5pmqSxwAnAZQNiLgP+Me/dczDwhO0HJY2TNB5A0jjgSKDqpA2NbOoJIYTWVoOHu7Z7JZ0CzCeby/wC28slnZzvPx+4AjgGWAE8DbwnP/zFwDxJkOXzi2xfVe2ckfhDCGE4XLtB2mxfQZbcy7edX/bZwIcqHHc38PLU80XiDyGE4WrnYZlDCCFU0Jp5PxJ/CK3kylVfT4o/eqcPFo5d9c7p1YPKOW24rm1Xtt/b9TV6gWvUReIPIYThMNAXiT+EEDqGcNzxhxBCx4nEH0IIHSYSfwghdBDTPwBby4nEH0IIwxRt/CGEpnPlQ+fVrewjXv355GMWXPvpOtSkUQyl1rzlj8QfQgjDYaKNP4QQOk5r3vBH4g8hhOGKNv4QQug0kfhDCKGD2NDXmm09kfhDCGG4WvSOv+FTL0r6mCRLmlC27XRJKyTdLumoRtYvhBAGZRdbmkxD7/glTQaOAO4r2zaDbM7JfYBdgIWS9rTd15hahhBCBQaqz6fblBp9x38u8AleOJ3BLOBi2+tt30M2x+RBjahcCCEMzuBSsaXJNCzxSzoWeMD2TQN27QrcX7a+Kt9WqYzZkhZLWrxmzZo61TSEEAYRTT2bkrQQ2KnCrjnAp4AjKx1WYVvFn5ztucBcgJ6enub76YbQxsYsvbPRVWgsE716KrF9eKXtkl4GTANukgQwCfiTpIPI7vAnl4VPAlbXs54hhDAsTXg3X0RDmnps32x7R9tTbU8lS/YH2H4IuAw4QdLmkqYB04HrG1HPEEIYXMFmnib8cmi6fvy2l0u6BLgF6AU+FD16QghNx8TonCOR3/WXr58FnNWY2oQQQkFNeDdfRFMk/hBCaEmR+EMIoYPYuK81W6Ej8YcQwnC16Ju7kfhDCGG4oqknhBA6iGPO3RBChyk9/XSjq9B4cccfQgidJB7uhhBCZ2nhYZkj8YcQwnA14ZDLRUTiDyGEYTDguOMPIYQOYscdfwghdJpWveOXW7Q70kCS1gD3AhOAtQ2uTiPEdXeWTr1uqM2172Z74kgKkHRVXpci1tqeOZLz1VLbJP5+khbb7ml0PUZbXHdn6dTrhs6+9lpp9GTrIYQQRlkk/hBC6DDtmPjnNroCDRLX3Vk69bqhs6+9JtqujT+EEMLQ2vGOP4QQwhAi8YcQQodpucQv6QJJD0taVrbtbZKWSypJ6hkQf7qkFZJul3TU6Ne4Nga57nMk3SZpqaR5krYr29fO1/25/JpvlHS1pF3K9rXtdZft+5gkS5pQtq1tr1vSZyU9kP993yjpmLJ9bXHdo852Sy3AocABwLKybXsDewG/BXrKts8AbgI2B6YBdwFdjb6GGl73kUB3/vls4OwOue5tyj6fCpzfCdedb58MzCd/WbETrhv4LPCxCrFtc92jvbTcHb/ta4BHB2y71fbtFcJnARfbXm/7HmAFcNAoVLPmBrnuq2335qvXAZPyz+1+3evKVseRjZcFbX7duXOBT7DxmqEzrruStrnu0dZyiT/RrsD9Zeur8m3t6L3Alfnntr9uSWdJuh94B3BGvrmtr1vSscADtm8asKutrzt3St68d4GkF+XbOuG666LdE78qbGu7/quS5gC9wIX9myqEtdV1255jezLZNZ+Sb27b65a0FTCHjV9yL9hdYVtbXHfu28AewH7Ag8CX8+3tft110+6JfxVZm2i/ScDqBtWlLiS9G3gT8A7nDZ90wHWXuQg4Pv/czte9B1k79k2SVpJd258k7UR7Xze2/2K7z3YJ+C4bm3Pa+rrrqd0T/2XACZI2lzQNmA5c3+A61YykmcBpwLG2y2e+bvfrnl62eixwW/65ba/b9s22d7Q91fZUsqR3gO2HaOPrBpC0c9nqW4D+Hj9tfd311HLj8Uv6CfB6YIKkVcBnyB4GfQOYCPxa0o22j7K9XNIlwC1kTSEfst2SsyMPct2nk/VoWCAJ4DrbJ3fAdR8jaS+gRNa75WSAdr9u29+vFNvu1w28XtJ+ZM04K4F/gva67tEWQzaEEEKHafemnhBCCANE4g8hhA4TiT+EEDpMJP4QQugwkfhDCKHDROIPIYQOE4m/A0h6qs7lXyFpu3z54DCOf72kyxPjn5B0xSD7/1PSW1Pr0Yryn8Wry9Y/Kuk+Sd9sZL1Cc4vEH0bM9jG2Hwe2A5IT/zD9r+1jqocNn6RWeMHx9cDzid/2uVQezyeE50Xi71CS9pN0XdkkLi/Kt/9W0tmSrpd0h6TX5tu3knRJHv9TSX9UPumNpJX5pCBfAPbIJ8s4Z+CdvKRvSjop/zxT2SQyvwOOK4sZl4/AeIOkP0uaVeBalJd9i6RfAzuW7TtQ0iJJSyTN73/9X9Ir8mv5Q17XZfn2kyT9TNKvgKsHq4+krvy4G/Jy/infvrOka/KfwbL+n98g9T4yP/+f8nNunW8/Iy93maS5yl/LlnRqfo1LJV0saSrZW8sfzc836LlCeIFGTwgQS/0X4KkK25YCr8s/nwl8Nf/8W+DL+edjgIX5548B38k/v5TsFfmefH0lMAGYygsn0Hg9cHnZ+jeBk4AtyIbTnU42wuIl/XHAvwPvzD9vB9wBjBtQ94HlHgcsALqAXYDHgbcCmwHXAhPzuLcDF+SflwGvzj9/ob/eef1WAdsPVR9gNvDpfPvmwGKyQdT+FZiTb+8Cxg/ydzIBuKb/2sjGXDoj/7x9WdyPgDfnn1cDm/fXJf/zswyYpCS/hm82+t9dLM27tMKvsqHGJG1LljgW5Zt+CPysLOTn+Z9LyJI5wGuArwHYXiZp6Qiq8BLgHtt35vX5MVkihWxWsWMlfSxf3wKYAtw6RHmHAj9xNk7Lakn/L9++F9mXVP9YRl3Ag8qmqBxv+9o87iKyEU77LbDdPxnIYPU5Eti37FnCtmRfZDcAF0jaDPiF7RsHqfPBZDNI/T6v21jgD/m+wyR9AtgK2B5YDvyK7Mv6Qkm/AH4xxM8jhCFF4g+VrM//7GPjv5FKY59X08sLmxO3KPs82CBRAo535RnVhlKpPAHLbb/qBRs3TuQxmL9Wq0/e/PJh2/M3Oal0KPBG4EeSzrH9X4PUbYHtEwccuwVwHtlvU/dL+iwbf25vJPuSOxb4v5L2qXIdIVQUbfwdyPYTwGNlbcLvAhYNcQjA74C/B5A0A3hZhZgngfFl6/cCM5QNm7st8IZ8+23ANEl75OvlyW8+8OGydu39C1zSNWTD83blbfiH5dtvByZKelVe1maS9rH9GPCkpIPzuBOGKHuw+swHPpDf2SNpz/x5wG7Aw7a/C3yfbP7YSq4DDpH0N/nxW0nak41Jfm3e5v/WfP8YYLLt35BNvbgdsDWb/sxDqCru+DvDVsqGuO33FeDdwPnKZna6G3hPlTLOA36YN/H8mazZ4YnyANuPSPp9/qD0StsfVzZs7lLgzvw4bD8raTbZENpryb5UXpoX8zngq8DSPNmu5IXNMJXMA/4WuJmsDX5Rfp4NeVPM1/Mvnu687OXA+4DvSvor2XONJzYtdsj6fI+sGexP+fY1wN+RPX/4uKTngKeAf6xUqO01yh50/0TS5vnmT9u+Q9J382tZSdZ0BFkz1Y/z6xBwru3H84fQ/50/dP6w7f+t8rMKIYZlDsVI6gI2y5P2HsD/AHva3tCAurye7IFmtS+EocrY2vZT+edPAjvb/ufa1LCx8i+UHtunVIsNnSnu+ENRWwG/yZs2BHygEUk/twF4qaQrPPy+/G+UdDrZ/wP3kvWEaXmSPkrWxfPSRtclNK+44w+hziT9kazLZ7l32b65EfUJIRJ/CCF0mOjVE0IIHSYSfwghdJhI/CGE0GEi8YcQQof5/wmEFJqsk+V1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "rvals.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:analysis3-21.01] *",
   "language": "python",
   "name": "conda-env-analysis3-21.01-py"
  },
  "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 1,
	"metadata": {},
	"outputs": [],
	"source": [
	"import xarray as xr\n",
	"import scipy.stats as stats\n",
	"import numpy as np"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 85,
	"metadata": {},
	"outputs": [],
	"source": [
	"DMI = xr.open_dataset(\"/g/data/w97/sl7808/DMI/9Dec/all20_1/output/test2.nc\")\n",
	"DMI_p = DMI.sel(time=DMI.time.dt.month.isin([8,9,10]))\n",
	"x=DMI_p.tos.groupby('time.year').mean(dim=\"time\")\n",
	"y=x.where(x>0.5)\n",
	"z=y.rename({\"year\":\"time\"})\n",
	"del z[\"time\"]"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 86,
	"metadata": {},
	"outputs": [],
	"source": [
	"year_HWX = xr.open_dataset(\"/g/data/w97/sl7808/corr/hwx/2_hwx_exp1.nc\",decode_times=False)\n",
	"year_HWX[\"time\"] = y[\"year\"]\n",
	"HWA = year_HWX.HWA_EHF"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 87,
	"metadata": {},
	"outputs": [],
	"source": [
	"# Mask with the land-sea mask\n",
	"land_sea = xr.open_dataset(\"/g/data/w97/sl7808/ppt/landmask_AUS_2.nc\")\n",
	"topo=land_sea.topo\n",
	"# Rename dimensions to be the same as HWA\n",
	"topo = topo.rename({\"lat_v\":\"lat\",\"lon_u\":\"lon\"})\n",
	"HWA_mask = HWA.where(topo == 1,drop=True)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 133,
	"metadata": {},
	"outputs": [],
	"source": [
	"# Function to select only the valid data at a point and calculate the correlation coefficient.\n",
	"# Return NaN if no valid data at the point.\n",
	"def corr_miss(data, tos):\n",
	" '''data and tos should both be 1D arrays'''\n",
	" # select the non-NAN:\n",
	" mask=~np.isnan(data) & ~np.isnan(tos)\n",
	" if (mask.any(\"time\") == False):\n",
	" return np.nan\n",
	" else:\n",
	" res = stats.linregress(data[mask],tos[mask])\n",
	" return res.rvalue**2"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 134,
	"metadata": {},
	"outputs": [],
	"source": [
	"# Calculate correlations between tos and HWA_mask:\n",
	"rvals=np.apply_along_axis(corr_miss,0,HWA_mask,z)"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 135,
	"metadata": {},
	"outputs": [],
	"source": [
	"rvals=xr.DataArray(rvals,dims=[\"lat\",\"lon\"],coords={\"lon\":HWA_mask.lon,\"lat\":HWA_mask.lat})"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 136,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"<matplotlib.collections.QuadMesh at 0x7f0e62ad21f0>"
	]
	},
	"execution_count": 136,
	"metadata": {},
	"output_type": "execute_result"
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEJCAYAAACT/UyFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqOUlEQVR4nO3de5wcVZ338c83M4RLCCAkyCUJCWxAgiKXEVEUZeUSUMkKusKqK16eLCqyuqsixgddlF0RFa+IUXFdBREXo4hASPbRsIoIiUJIuAYIEAKScAsIJMz09/mjakgz6ZmuM9M9ffu9edUrXVW/OnVqEn5dc+rUObJNCCGEzjGm0RUIIYQwuiLxhxBCh4nEH0IIHSYSfwghdJhI/CGE0GEi8YcQQoeJxB9CCA0maaak2yWtkPTJCvvfIWlpvlwr6eVl+1ZKulnSjZIWFzpfu/TjnzBhgqdOndroaoTAnUvvLxzr556rY03SaMstko+ZPmOXOtSk/pYsWbLW9sSRlHHUYeP8yKN9xc63dP182zMr7ZPUBdwBHAGsAm4ATrR9S1nMq4FbbT8m6Wjgs7Zfme9bCfTYXlu07t1FA0eTpM8C/wdYk2/6lO0rhjpm6tSpLF5c6MsuhLo6etKphWN7Vz+YVrgSf0l3qXBo1557p5UNXLX4zORjmoGke0daxiOP9nH9/CmFYrt2vnPCELsPAlbYvjuv28XALOD5xG/72rL464BJyRUu05SJP3eu7S81uhIhhFCJgRKFv1gnDGiGmWt7bv55V6D818RVwCuHKOt9wJUDqnK1JAPfKSt3UM2c+EMIoWkZ85yLNfUAa233DLJPFYuvFCgdRpb4X1O2+RDbqyXtCCyQdJvta4aqTDM/3D0lf5BxgaQXVQqQNFvSYkmL16xZUykkhBDqplTwvypWAZPL1icBqwcGSdoX+B4wy/Yj/dttr87/fBiYR9Z0NKSGJX5JCyUtq7DMAr4N7AHsBzwIfLlSGbbn2u6x3TNx4oie04QQQhJj+lxsqeIGYLqkaZLGAicAl5UHSJoC/Bx4l+07yraPkzS+/zNwJLCs2gkb1tRj+/AicZK+C1xe5+qEEEKyUuUWmSS2eyWdAswHuoALbC+XdHK+/3zgDGAH4DxJAL1509GLgXn5tm7gIttXVTtnU7bxS9rZdn93h7dQ4BsshDD6jtriHUnx85+9sE41GX0G+mqQ+AHyXotXDNh2ftnn9wPvr3Dc3cDLB26vpikTP/BFSfuR/WxXAv/U0NqEkODKVV+vW9kzt3lPUnzp2fWFY/uW3ppU9pixY5Pi21Et7vgboSkTv+13NboOIYQwFAPPtegLsE2Z+EMIodkZ16ypZ7RF4g8hhOEw9LVm3o/EH0IIw5G9uduaIvGHEMKwiL6KL902v7ZJ/HcsuZsjNzuhUOyY8eOTy7/q0e8mHxNCaF8GStHUE0KoN2/YkBbfW3zYZ3V1JZVdSqxLuzGwoalHvRlcJP4QQhimkqOpJ4QQOkb25m4k/hBC6BhG9EVTTwghdJZo6gkhhA5ixAanPRBvFpH4QwijZuaE2YVjr1pbdQbBhspe4Iqmnoba88DduXrxxY2uRgh1lTqs8RFj3lY41n19LCj9LLVKhaUk/VYRD3dDCKGD2KLPcccfQggdpRR3/CGE0Dmyfvxxxx9CCB3DiOfcmim0NWsdQghNoC/68YcQQudo5Td3G1ZrSW+TtFxSSVLPgH2nS1oh6XZJRzWqjiGEMJSSxxRamk0j7/iXAccB3ynfKGkGcAKwD7ALsFDSnrb7Rr+KIYRQWTzcHQbbtwJIm7SRzQIutr0euEfSCuAg4A+jW8MQQhicUbTx19CuwHVl66vybZuQNBuYDTBlypT61yyEEHI20aunEkkLgZ0q7Jpj+5eDHVZhW8UJzmzPBeYC9PT0tOgkaCHUz9PHH9zoKjyv2cfeSad4gasS24cP47BVwOSy9UnA6trUKIQQasPQskM2NGOtLwNOkLS5pGnAdOD6BtcphBA20ceYQkuzaVgDlaS3AN8AJgK/lnSj7aNsL5d0CXAL0At8KHr0hBCajVFMxJLK9jxg3iD7zgLOGt0ahRBCcSYe7oYQQodRjMcfQgidxNCUb+UWEYk/hBCGKe74Qwihg9iKO/4QGunovU8vHNt7+4qksrsnVXxxfEhX3vfV5GNC66lVP35JM4GvAV3A92x/YcD+dwCn5atPAR+wfVORYyuJxB9CG3tmQmvekbaCbCKWrhGXI6kL+BZwBNkLrDdIusz2LWVh9wCvs/2YpKPJRix4ZcFjNxH/KkIIYRiyh7sqtFRxELDC9t22NwAXkw1WufFc9rW2H8tXryMb0aDQsZVE4g8hhGGq0Zu7uwL3l60POjBl7n3AlcM8FoimnhBCGJbEN3cnSFpctj43H2QSEgamlHQYWeJ/Teqx5SLxhxDCMJWKN5qstd0zyL5CA1NK2hf4HnC07UdSjh0omnpCCGEY7Gyy9SJLFTcA0yVNkzSWbAbCy8oDJE0Bfg68y/YdKcdWMuQdv6TjqhUAPGv7igJxIYTQNozoLY28V4/tXkmnAPPJumRekA9WeXK+/3zgDGAH4Lx81sJe2z2DHVvtnNWaer4L/JLK7Uj9DgUi8YeGSumb3z11t6Sy/cijqdVpGjv+bm2jq9DWavXmbn7zfMWAbeeXfX4/8P6ix1ZTLfFfafu9QwVI+nHKCUMIoR30d+dsRUMmftvvrFZAkZgQQmg/HTBkg6RXA1PLj7H9X3WoUwghtIS2nnNX0o+APYAbgf7ZsAxE4g8hdCQbnqvBw91GKHrH3wPMsF31xYAQQugEnTD14jJgJ+DBOtYltJjDX/fvSfG9Wxa/O3pu67Q7qa2SosNgjuh6e1L8gr6f1qkmraEtm3ok/YqsSWc8cIuk64H1/fttHzvcE0t6G/BZYG/gINuL8+1TgVuB2/PQ62yfPNzzhM7Qte22hWN7V96bVPaYLbdMrU7T8D2rCseWnnm2jjVpP23bqwf4Uh3PvQw4DvhOhX132d6vjucOIYQRa8tePbYXAUg62/Zp5fsknQ0sGu6Jbd+alzPcIkIIoXGKDbnclIp+XR1RYdvRtazIANMk/VnSIkmvHSxI0mxJiyUtXrNmTR2rE0IIL2Sg12MKLc2mWhv/B4APAntIWlq2azzw+2qFS1pI9lB4oDm2fznIYQ8CU2w/IulA4BeS9rG9bmBgPqzpXICenp7ocRRCGDXt3MZ/EdmA//8BfLJs+5O2qw5gYvvw1ArZXk/+ANn2Ekl3AXsCi4c8MIQQRllbJn7bT0h6EniZ7bSuEMMkaSLwqO0+SbsD04G7R+Pcne5VJ345KX5cnerRDl75rq8Ujv3jj/6ljjUJ9dLW/fhtlyTdJGmK7ftqdWJJbwG+AUwEfi3pRttHkY32eaakXrK3hE8u8ttFGH1r903r5vjiH9xYOHbzl0xLKrtv3ZOFY7unTKoeVKb3vuJdIgHWveNVSfH1NGbyLoVjSwkjnIZMW/bjL7MzsDzvx//X/o0j6cdvex4wr8L2S4FLh1tuCCGMCrdpU0+Zf6trLUIIocUY6C01X4+dIgolftuLJL0YeEW+6XrbD9evWiGE0NxauY2/0NeVpL8HrgfeBvw98EdJb61nxUIIodnZKrQ0m6JNPXOAV/Tf5ec9bxYC/12vioUQQrNr94e7YwY07TxC8bd+QwhDOHrv05Pir7z1P+pUk5DCHfBw9ypJ84Gf5OtvJyZYbzvbLbwzKX7rPdO6RbLX1MKhuiutCyUuFQ9d91Ra2Ym2ueiPhWO799y9jjVJm4Qexb1cGtHX5g93Py7peOAQQMDcvDtmCCF0rGZsvy+i8Jy70b8+hBA2auWxeor26jlO0p2SnpC0TtKTkjYZNC2EEDqGs3b+IkuzKXrH/0Xgzf1j6IcQQmj/Xj1/iaQfQggbmfZv418s6afAL3jhnLs/r0elQu0cOfYfCsd2bbtNHWsSauWIMW9rdBWed9hRZyfF/2b+adWDWoboK7V34t8GeBo4smybgUj87WTzsUnh3avTBk3148UfC5WeSuty2b1H8dE8S6sfSio7ebL1hOlEe+9IHHE8odtqqq6t0wbaXn/wXnWqSeto6zt+2+8Zar+k023HWyUhhI6RPbhtzcRfq7cPmud3zxBCGCWlfML1akuzKdyPv4rmu7IQQqizZuyqWUStEn+LXn4IIQyPEaV2HrKhgLjjDyF0nFa9461V4v9ZjcoJVbzujV9Mit+8TvUIYTiO6Hp74dgFfT+tY01qoIUf7hZK/JK+CHweeAa4Cng58BHbPwaw/e+pJ5Z0DvBmYANwF/Ae24/n+04H3kc22fqptuenlh8yzxyzf+HYrW9Zm1R2710r0yqT0BWx+8U7JhXdd2/x0TzV1ZVUdmnDhqT4lOscs/kWSUVrs7R7tb6X7VE4tnTDLUlldy/8U1J8W6rRLb+kmcDXgC7ge7a/MGD/S4AfAAcAc2x/qWzfSuBJsnzZa7un2vmKNlAdaXsd8CZgFbAn8PGCxw5mAfBS2/sCdwCnA0iaAZwA7APMBM6TlPZ/agghjIJazMCV57dvAUcDM4AT8zxY7lHgVOBLVHaY7f2KJH0onvg3y/88BviJ7bQ3dyqwfbXt3nz1OqB/cPdZwMW219u+B1gBHDTS84UQQq3VaJC2g4AVtu+2vQG4mCwPlp3HD9u+AXiuFvUumvh/Jek2oAf4n3zqxWdrUYHce4Er88+7AveX7VuVb9uEpNmSFktavGbNmhpWJ4QQhmaDS2MKLVUUznmDVQW4WtISSbOLHFD0zd1PSjobWGe7T9LTDPhGqkTSQmCnCrvm2P5lHjMH6AUu7D+sUhUGqddcYC5AT09Pqz5gDyG0qIR+/BMkLS5bn5vnL0jIeYM4xPZqSTsCCyTdZvuaoQ4o+nB3K+BDwBRgNrALsBdw+VDH2T68SrnvJntu8Ab7+R/hKmByWdgkYHWReoYQwqgqnp7XDtH+PqKcZ3t1/ufDkuaRNR2NPPGTPU1eAry6rKI/o0riH0r+FPs04HW2ny7bdRlwkaSvkH3BTAeuH+55GqHnfV8pHDvhv5elFf6alyTWJoRQH9Uf3BZ0AzBd0jTgAbLOLYWG1ZU0Dhhj+8n885HAmdWOK5r497D9dkknAth+RkoYgrCyb5J1M1+QF3Wd7ZNtL5d0CXALWRPQh2z3jfBcbWOrOx9JOyChK2LfTi9KKrrbuyXFe4vNqgf1e+KvSWVrTPF/jqX900aVHLM0bRL60tNPVw/qlzra5rS0Ce67Hyj+76U0NuHvB/DTac8Zx4xNG/21JdSggdl2r6RTgPlk3TkvyPPgyfn+8yXtBCwmGym5JOkjZD2AJgDz8hzaDVxk+6pq5yya+DdI2pL8MiXtQdm4/MNh+2+G2HcWcNZIyg8hhLqq4Qtctq8Arhiw7fyyzw+xsedjuXVk71UlKZr4P0P24tZkSRcChwAnpZ4shBDaSju/uWt7gaQ/AQeTPYH+Z9tpr3mGEEK7adG+hIX68eft+UcDB9q+HNhKUrxUFULobC64NJmiL3CdB7wKODFff5LsFeMQQuhMJmvqKbI0maJt/K+0fYCkPwPYfkxSGz6iDyE0i5nbvjcp/qonLqhTTQbX7hOxPJcPJNTfq2ciUL9Zn1vchCWPF45VYhe3DbtsmxQ/9q6/FI4dc+vKpLL7nkkbtaNr4g6FY3sfTJsQPakey+5Kii+tH1EHtiG5L62nsm9Nq3vKzzypGyqg7rTunymjnHaNH59UdsOUmu9uvoiiTT1fB+YBO0o6C/gdkDwUcwghtBO52NJsqt7xSxoD3AN8AngDWa+ev7N9a53rFkIIzatJH9wWUTXx2y5J+rLtVwG3jUKdQgihBTTng9siijb1XC3p+BoM0xBCCO2jRbtzFn24+y/AOKBX0rNkzT22vU3dahZCCM2uCZN6EUXf3G2RR+z1MXOfOWkHbBYzRYbQ9kzL9uopOh7/ARU2PwHcWzZ9Ysg9MaN4l8vtHnkiqezN1qSNWpnS0Vi7vDip6L6d037h61p+X+HY1JEcU7oK9j31VFLZ3VPSRsTs23n7wrH+49Kksrv2KT55OkDp7vurBz1f9p5JZfctvyMpPkVqN+dGacYeO0UUbeo5j2x295vz9ZcBNwE7SDrZ9tX1qFwIITS1Fk38RR/urgT2t32g7QOB/YBlwOHAF+tTtRBCCPVQ9I7/JbaX96/YvkXS/rbvjo4+IYRO1e5NPbdL+jZwcb7+duAOSZsDadPwhBBCu2jRfvxFE/9JwAeBj5B15fwd8DGypH9YPSoWQghNzbTsiGVFu3M+I+k84HLbtw/YndZFIoQQ6uDoSaeO+jnbuqlH0rHAOcBYYJqk/YAzbR873BNLOgd4M7ABuAt4j+3HJU0FbgX6v2Cus33ycM9TC1vOfTwpft3XinfnKz36WFLZTx+Y1rVw908VH+Xyph9NSyp7y7VptzvPTBh0muVNbHv9qqSy+UvChHCpE5z/9ZmkcCdM/N619dZpdRlTtD9GpvRM8bp3P554D6e0uijh59L7SPFJ4gG6d9k5Kb5mWjTxF/2b+wxwEPA4gO0bgakjPPcC4KW29wXuAE4v23eX7f3ypaFJP4QQBtWiQzYUTfy9ttPeNKrC9tVlL39dR+UZ5EMIoSkVHZK5GZuDiib+ZZL+AeiSNF3SN4Bra1iP9wJXlq1Pk/RnSYskvXawgyTNlrRY0uI1a9bUsDohhFBAScWWJlM08X8Y2AdYD/wEWEfWw2dIkhZKWlZhmVUWMwfoBS7MNz0ITLG9P9ngcBdJqjg2gO25tnts90ycOLHgpYQQQm206h1/0V49TwNz8qUw24cPtV/Su4E3AW+ws0FlbK8n+4LB9hJJdwF7AotTzh1CCHXXhEm9iCETv6RfMcSljbBXz0zgNOB1+RdL//aJwKO2+yTtDkwH7h7ueUIIoS6a9G6+iGp3/F/K/zwO2An4cb5+Itn4PSPxTWBzYEE+7EN/t81DgTMl9QJ9wMm2Hx3huUbkpJ1/nxT/8X2Ld1tc+/L9k8qe/D9pE5zf/vV9CsfuvGhlUtmPHrZbUvz2v7u3eHBXWlfBMVN2LR78VNoIp70PPZwU33Vz8cnZtdOOSWX7/rRJ6LtfXLz8vof+klR2PXVPSvj7BEqPNChFtGPit70IQNLnbB9atutXkq4ZyYltV8yOti8FLh1J2SGEMCpaNPEXva2amDe7ACBpGhBPU0MIHa2tH+4CHwV+K6m/rX0qMLsuNQohhFbRhEm9iKK9eq6SNB14Sb7ptrz3TQghdKYmvZsvYsimnvIpF22vt31TvqyvFBNCCB2lRYdsqHbH/wNJrycbinkw3wfSuqaEEEI7aMKkXkS1xL8tsIShE3/bj5Vw2o3HJcVPPqT4BNd8Lu0Z+WaPPF09qEzXDsUnrV5zZFr3zImX35kUX0roRpkyqiRA947Ff44p9QDonlB8tFWA3jXFRwrtfu5FSWU/t9/u1YPK48cXfYwH4+7cNqlsP5wwIirQlzASbWlN2uicpfVp3ZxrQbRpU4/tqbZ3tz1tiOWg0apsCCE0DYNKxZZqJM2UdLukFZI+WWH/SyT9QdJ6SR9LObaStLdkQgghbFSDNn5JXcC3gKOBGcCJkmYMCHsUOJWNL9WmHLuJSPwhhDBctXm4exCwwvbdtjeQzW0+qzzA9sO2b2DTOc6rHltJJP4QQhimGr3AtStQ/mBwVb6tiGEdWyjxK/NOSWfk61MkRdt+CKGzFb/jn9A/d0i+lL8AW6nzTNHHxsM6tugj//PI5pP/W+BM4Emy8XReUfD4EEJoLy724Da31nbPIPtWAZPL1icBqwuWO6xjiyb+V9o+QNKfAWw/Jql4P8EWd9txZyTFv+xfzy0cW0r8venmj/wgKf4VZ3ygcOwO/3lDUtnecouk+NKzCaNWdm+WVPZTBxefKH6Ly65PKjtV9847FY4tPZTWG3qzbcYlxY99qK948JNp3VzpS5y0PoHGbZkU35UYT1pv0cHVpjvnDcD0fAy0B4ATgH+o57FFE/9z+dNjw/Nj5tfvbz2EEFpALfrx2+6VdAowH+gCLrC9XNLJ+f7zJe1ENhnVNkBJ0keAGbbXVTq22jmLJv6vA/OAHSWdBbwV+HTa5YUQQpup0Qtctq8Arhiw7fyyzw+RNeMUOraaooO0XShpCfAGsocJf2f71pQThRBCW2nScXiKqDb1Yvm76g+TTbT+/L5Gz4wVQgiNIoYey6aZVbvjX0L2nSZgCvBY/nk74D6g+BO1EEJoMwm9eppKtbF6ptnenezBwZttT7C9A/Am4OejUcEQQmhabTosc79X5BOhA2D7SkmfG8mJ8+NnkfUOehg4yfbqfN/pwPvIJls/1fb8kZxrtI19snisSmn/Kva+9p1J8dP+kDCC4l5pv8Dd9onxSfF7n160azJ4+22Syt5q4bLCsak3aaWn00ZETenmWNqwIa3sm29LCh8zto69rpX24n/KBOq9qx5IKnvM5mldi2umCZN6EUX/5tZK+rSkqZJ2kzSHkfeEPcf2vrb3Ay4H+t8KnkHWF3UfYCZwXt6VNIQQmkfB4Rqacejmoon/RLLJ1ecBvwB2zLcNm+11Zavj2PjdOQu4OJ/x6x5gBdlARCGE0Fzauakn773zz7U+ef5OwD8CTwCH5Zt3Ba4rCxt00KF8vIvZAFOmTKl19UIIYUit+nC3UOKX9BsqfG/Z/tsqxy0EKr2/Psf2L23PAebkbfqnAJ8hYdAh23OBuQA9PT1N+L0aQmhnzdiMU0TRh7vlM75sARwP9FY7yPbhBcu/CPg1WeIfyYBFIYQwOpq0GaeIok09SwZs+r2kRSM5saTptvsnbT0W6O+ucBlwkaSvALsA04H6jqoVQgjD0c6Jf8AbvGOAA6nchJPiC5L2Iutddy/QPyDRckmXALeQ/VbxIdsJQww23jM7Fo+dcl7ayBcTr02bnL33tuIToo/Zv+qMbS8w/qbNk+JXnVC8u+ikS+5NKru0YeDERLXTvesuSfG9DxT/BVVdqR3W0uJLz1X9xXxjyVunjfzpvsT/LbuLT/ye+nMZs/12aXWpQRtCK0+2XvRvovwN3l7gHrJ+9sNm+/gh9p0FnDWS8kMIoe7aPPHvbfvZ8g2S0m73QgihnTj9BcxmUbQf/7UVtv2hlhUJIYRW06ovcFUbnXMnsj70W0ran41dLbcBtqpz3UIIobk1YVIvolpTz1HASWRdKr9Stv1J4FN1qlMIIbSEZrybL2LIxG/7h8APJR1v+9JRqlMIIbSGdkz8kt5p+8fAVEn/MnC/7a9UOKzj7XbJQ8WDxyW2mD36eFJ41/YvKhxbWnp7UtmTHqs4E9yQrrz7S8UCC4aNhkOPPScpflyp+Hv8vQ89nFR2926TqweVKa1+sHBs35MJw8oOw5iu4qN5diVMWA/Q92DC/3O10qTt90VUa+rp79i7dYV9LXrJIYQwcqJNx+qx/Z3840Lbvy/fJ+mQutUqhBBagVvz/rfo717fKLgthBA6Rrt253wV8Gpg4oA2/m1IfXc8hBDaSRsP0jaWrH2/GyifZ28d8NZ6VSqEEFpBu7bxLwIWSfpP22mjZoUQQptry8Rf5mlJ55DNg/v8rMbVJmLpVL0Tik9Cvvq0LZPKnvLxZ5Li+1bcnRS/oPSzpPhOcM1lH0+KP+yoswvH9r4ibea4lEnlIW10zu4ddkgq2+vXJ8X33nVP8bokTMwOWffPK+/7auF46adJ5Vdk2v7h7oVk4+VPA/4NWAncUKc6hRBCS2jVh7tFE/8Otr8PPGd7ke33AgfXsV4hhND82nmydaB/losHJb2RbBqD9Nc2QwihTXTCRCyfl7Qt8K9k/fe3AT5Sr0qFEELTs1u2jb/onLuX5x+fAA4DkPSROtUphBBaQqv26ik+atKmNhm0LYQQOkmrPtwtPvvxplQ9ZIiDpc8Bs8gmW38YOMn2aklTgVuB/qEir7N98kjONdoWXPvp+hU+q35Fh9r4zfzTCsfu+9Fzk8reYr89k+J7t0r4X/ypxAnr/3hzUnh3yoibY0aUXkaHgTaferGSkV7xObb3tb0fcDlwRtm+u2zvly8tlfRDCB2kHXv1SHqSytUWkPbm0QC215WtjhvkPCGE0LSasRmniCHv+G2Pt71NhWW87ZE0EwEg6SxJ9wPv4IV3/NMk/VnSIkmvHeL42ZIWS1q8Zs2akVYnhBDS9PfsqbZUIWmmpNslrZD0yQr7Jenr+f6lkg4o27dS0s2SbpS0uEi1R9LUU5WkhZKWVVhmAdieY3sy2ZvBp+SHPQhMsb0/2QPkiyRtU6l823Nt99jumThxYj0vJYQQXshZr54iy1AkdQHfAo4GZgAnSpoxIOxoYHq+zAa+PWD/YXnTeE+Rqo/4rn0otg8vGHoR8GvgM7bXA+vz45dIugvYEyj0TRZCCKMhe4GrJm09BwErbN8NIOlism4ct5TFzAL+y7aB6yRtJ2ln28Xn1ixT1zv+oUiaXrZ6LNlYQEiamH8DIml3sm+4tJHGQghhNJQKLjChv1k6X2aXlbIrcH/Z+qp8GwVjDFwtacmAcgdV1zv+Kr4gaS+yH8u9QH/vnUOBMyX1An3AybYfbVAdQ6irped+NCn+6EmnJsWX/qZ4F0pvlnYfuNkeU5Pie+9amRS/oK8GI2jWWcId/9ohmmEq9V0dWPBQMYfkXeF3BBZIus32NUNVpmGJ3/bxg2y/FLh0lKsTQghpatdVcxUwuWx9Etl4aIVibPf/+bCkeWRNR0Mm/oY19YQQQmszKhVbqrgBmC5pmqSxwAnAZQNiLgP+Me/dczDwhO0HJY2TNB5A0jjgSKDqpA2NbOoJIYTWVoOHu7Z7JZ0CzCeby/wC28slnZzvPx+4AjgGWAE8DbwnP/zFwDxJkOXzi2xfVe2ckfhDCGE4XLtB2mxfQZbcy7edX/bZwIcqHHc38PLU80XiDyGE4WrnYZlDCCFU0Jp5PxJ/CK3kylVfT4o/eqcPFo5d9c7p1YPKOW24rm1Xtt/b9TV6gWvUReIPIYThMNAXiT+EEDqGcNzxhxBCx4nEH0IIHSYSfwghdBDTPwBby4nEH0IIwxRt/CGEpnPlQ+fVrewjXv355GMWXPvpOtSkUQyl1rzlj8QfQgjDYaKNP4QQOk5r3vBH4g8hhOGKNv4QQug0kfhDCKGD2NDXmm09kfhDCGG4WvSOv+FTL0r6mCRLmlC27XRJKyTdLumoRtYvhBAGZRdbmkxD7/glTQaOAO4r2zaDbM7JfYBdgIWS9rTd15hahhBCBQaqz6fblBp9x38u8AleOJ3BLOBi2+tt30M2x+RBjahcCCEMzuBSsaXJNCzxSzoWeMD2TQN27QrcX7a+Kt9WqYzZkhZLWrxmzZo61TSEEAYRTT2bkrQQ2KnCrjnAp4AjKx1WYVvFn5ztucBcgJ6enub76YbQxsYsvbPRVWgsE716KrF9eKXtkl4GTANukgQwCfiTpIPI7vAnl4VPAlbXs54hhDAsTXg3X0RDmnps32x7R9tTbU8lS/YH2H4IuAw4QdLmkqYB04HrG1HPEEIYXMFmnib8cmi6fvy2l0u6BLgF6AU+FD16QghNx8TonCOR3/WXr58FnNWY2oQQQkFNeDdfRFMk/hBCaEmR+EMIoYPYuK81W6Ej8YcQwnC16Ju7kfhDCGG4oqknhBA6iGPO3RBChyk9/XSjq9B4cccfQgidJB7uhhBCZ2nhYZkj8YcQwnA14ZDLRUTiDyGEYTDguOMPIYQOYscdfwghdJpWveOXW7Q70kCS1gD3AhOAtQ2uTiPEdXeWTr1uqM2172Z74kgKkHRVXpci1tqeOZLz1VLbJP5+khbb7ml0PUZbXHdn6dTrhs6+9lpp9GTrIYQQRlkk/hBC6DDtmPjnNroCDRLX3Vk69bqhs6+9JtqujT+EEMLQ2vGOP4QQwhAi8YcQQodpucQv6QJJD0taVrbtbZKWSypJ6hkQf7qkFZJul3TU6Ne4Nga57nMk3SZpqaR5krYr29fO1/25/JpvlHS1pF3K9rXtdZft+5gkS5pQtq1tr1vSZyU9kP993yjpmLJ9bXHdo852Sy3AocABwLKybXsDewG/BXrKts8AbgI2B6YBdwFdjb6GGl73kUB3/vls4OwOue5tyj6fCpzfCdedb58MzCd/WbETrhv4LPCxCrFtc92jvbTcHb/ta4BHB2y71fbtFcJnARfbXm/7HmAFcNAoVLPmBrnuq2335qvXAZPyz+1+3evKVseRjZcFbX7duXOBT7DxmqEzrruStrnu0dZyiT/RrsD9Zeur8m3t6L3Alfnntr9uSWdJuh94B3BGvrmtr1vSscADtm8asKutrzt3St68d4GkF+XbOuG666LdE78qbGu7/quS5gC9wIX9myqEtdV1255jezLZNZ+Sb27b65a0FTCHjV9yL9hdYVtbXHfu28AewH7Ag8CX8+3tft110+6JfxVZm2i/ScDqBtWlLiS9G3gT8A7nDZ90wHWXuQg4Pv/czte9B1k79k2SVpJd258k7UR7Xze2/2K7z3YJ+C4bm3Pa+rrrqd0T/2XACZI2lzQNmA5c3+A61YykmcBpwLG2y2e+bvfrnl62eixwW/65ba/b9s22d7Q91fZUsqR3gO2HaOPrBpC0c9nqW4D+Hj9tfd311HLj8Uv6CfB6YIKkVcBnyB4GfQOYCPxa0o22j7K9XNIlwC1kTSEfst2SsyMPct2nk/VoWCAJ4DrbJ3fAdR8jaS+gRNa75WSAdr9u29+vFNvu1w28XtJ+ZM04K4F/gva67tEWQzaEEEKHafemnhBCCANE4g8hhA4TiT+EEDpMJP4QQugwkfhDCKHDROIPIYQOE4m/A0h6qs7lXyFpu3z54DCOf72kyxPjn5B0xSD7/1PSW1Pr0Yryn8Wry9Y/Kuk+Sd9sZL1Cc4vEH0bM9jG2Hwe2A5IT/zD9r+1jqocNn6RWeMHx9cDzid/2uVQezyeE50Xi71CS9pN0XdkkLi/Kt/9W0tmSrpd0h6TX5tu3knRJHv9TSX9UPumNpJX5pCBfAPbIJ8s4Z+CdvKRvSjop/zxT2SQyvwOOK4sZl4/AeIOkP0uaVeBalJd9i6RfAzuW7TtQ0iJJSyTN73/9X9Ir8mv5Q17XZfn2kyT9TNKvgKsHq4+krvy4G/Jy/infvrOka/KfwbL+n98g9T4yP/+f8nNunW8/Iy93maS5yl/LlnRqfo1LJV0saSrZW8sfzc836LlCeIFGTwgQS/0X4KkK25YCr8s/nwl8Nf/8W+DL+edjgIX5548B38k/v5TsFfmefH0lMAGYygsn0Hg9cHnZ+jeBk4AtyIbTnU42wuIl/XHAvwPvzD9vB9wBjBtQ94HlHgcsALqAXYDHgbcCmwHXAhPzuLcDF+SflwGvzj9/ob/eef1WAdsPVR9gNvDpfPvmwGKyQdT+FZiTb+8Cxg/ydzIBuKb/2sjGXDoj/7x9WdyPgDfnn1cDm/fXJf/zswyYpCS/hm82+t9dLM27tMKvsqHGJG1LljgW5Zt+CPysLOTn+Z9LyJI5wGuArwHYXiZp6Qiq8BLgHtt35vX5MVkihWxWsWMlfSxf3wKYAtw6RHmHAj9xNk7Lakn/L9++F9mXVP9YRl3Ag8qmqBxv+9o87iKyEU77LbDdPxnIYPU5Eti37FnCtmRfZDcAF0jaDPiF7RsHqfPBZDNI/T6v21jgD/m+wyR9AtgK2B5YDvyK7Mv6Qkm/AH4xxM8jhCFF4g+VrM//7GPjv5FKY59X08sLmxO3KPs82CBRAo535RnVhlKpPAHLbb/qBRs3TuQxmL9Wq0/e/PJh2/M3Oal0KPBG4EeSzrH9X4PUbYHtEwccuwVwHtlvU/dL+iwbf25vJPuSOxb4v5L2qXIdIVQUbfwdyPYTwGNlbcLvAhYNcQjA74C/B5A0A3hZhZgngfFl6/cCM5QNm7st8IZ8+23ANEl75OvlyW8+8OGydu39C1zSNWTD83blbfiH5dtvByZKelVe1maS9rH9GPCkpIPzuBOGKHuw+swHPpDf2SNpz/x5wG7Aw7a/C3yfbP7YSq4DDpH0N/nxW0nak41Jfm3e5v/WfP8YYLLt35BNvbgdsDWb/sxDqCru+DvDVsqGuO33FeDdwPnKZna6G3hPlTLOA36YN/H8mazZ4YnyANuPSPp9/qD0StsfVzZs7lLgzvw4bD8raTbZENpryb5UXpoX8zngq8DSPNmu5IXNMJXMA/4WuJmsDX5Rfp4NeVPM1/Mvnu687OXA+4DvSvor2XONJzYtdsj6fI+sGexP+fY1wN+RPX/4uKTngKeAf6xUqO01yh50/0TS5vnmT9u+Q9J382tZSdZ0BFkz1Y/z6xBwru3H84fQ/50/dP6w7f+t8rMKIYZlDsVI6gI2y5P2HsD/AHva3tCAurye7IFmtS+EocrY2vZT+edPAjvb/ufa1LCx8i+UHtunVIsNnSnu+ENRWwG/yZs2BHygEUk/twF4qaQrPPy+/G+UdDrZ/wP3kvWEaXmSPkrWxfPSRtclNK+44w+hziT9kazLZ7l32b65EfUJIRJ/CCF0mOjVE0IIHSYSfwghdJhI/CGE0GEi8YcQQof5/wmEFJqsk+V1AAAAAElFTkSuQmCC\n",
	"text/plain": [
	"<Figure size 432x288 with 2 Axes>"
	]
	},
	"metadata": {
	"needs_background": "light"
	},
	"output_type": "display_data"
	}
	],
	"source": [
	"rvals.plot()"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"metadata": {},
	"outputs": [],
	"source": []
	}
	],
	"metadata": {
	"kernelspec": {
	"display_name": "Python [conda env:analysis3-21.01] *",
	"language": "python",
	"name": "conda-env-analysis3-21.01-py"
	},
	"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.8.6"
	}
	},
	"nbformat": 4,
	"nbformat_minor": 4
	}