OISST_ENSO_MHWs-empty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9f2f752f-bd01-472d-9b82-f7dbbb04ae15",
   "metadata": {},
   "source": [
    "# Working with Satellite Remote Sensing Data: Ocean Edition\n",
    "## Focus on ENSO and Marine Heat Waves\n",
    "### EDS 220 Fall 2021\n",
    "\n",
    "This week's class will focus on some commonly employed diagnostics used in physical oceanography and climate science, using satellite remote sensing products. \n",
    "\n",
    "During the course of this lab exercise, we will cover the following topics/analysis skills:\n",
    "- Use of OpeNDAP and THREDDS to access data remotely\n",
    "- Spatial subsetting and regional averaging of data\n",
    "- Removal of seasonal cycle using \"groupby\"\n",
    "- Use of threshold exceedances to identify events (El Nino/La Nina, marine heat waves)\n",
    "- Performing computations on all gridpoints within a spatial data field\n",
    "- Creating georeferenced maps of spatial fields using cartopy\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdc03366-9316-4c90-abcc-415fcfd1f913",
   "metadata": {},
   "source": [
    "#### **Datasets: ERSSTv5 and OISSTv2**\n",
    "In this lab, we'll be using data from two sources:\n",
    "\n",
    "- MONTHLY data from the Extended Reconstructed Sea Surface Temperature dataset version 5 (ERSSTv5), to identify El Nino and La Nina periods:\n",
    "- DAILY data from NOAA's 1/4° Daily Optimum Interpolation Sea Surface Temperature (OISST) version 2, or OISSTv2, to identify marine heat waves. \n",
    "\n",
    "\n",
    "#### **OpeNDAP and THREDDS**\n",
    "\n",
    "The method of accessing these files will be via the NOAA Open-source Project for a Network Data Access Protocol (OpeNDAP) server. OpeNDAP is a fairly common method of providing remote access to large data files without the need to download the entire file; more information can be found here:\n",
    "\n",
    "https://earthdata.nasa.gov/collaborate/open-data-services-and-software/api/opendap\n",
    "\n",
    "https://earthdata.nasa.gov/collaborate/open-data-services-and-software/api/opendap/opendap-user-guide\n",
    "\n",
    "The OpeNDAP protocol is implemented on many scientific servers using a Web service called the Thematic Real-time Environmental Distributed Data Services (THREDDS), which is responsible for reading in the various metadata fields associated with the appropriate files. More information on THREDDS can be found here:\n",
    "\n",
    "https://www.unidata.ucar.edu/software/tds/current/"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11c9c014-5a06-4bb7-b0d6-6913058d2e49",
   "metadata": {},
   "source": [
    "### PART 1: ENSO\n",
    "\n",
    "First, we will apply the ERSSTv5 product to identify El Nino and La Nina events. The Python packages we will need for this and subsequent parts of the notebook are numpy, matplotlib, cartopy, and xarray:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8cef0b07-ae72-4095-baa7-5a72a868d08a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure environment\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "716464bc-c29b-4ba4-826e-a1cd4977f9cf",
   "metadata": {},
   "source": [
    "The ERSSTv5 data is stored in a single file, containing information for the entire period. The website describing ERSSTv5 data can be found here:\n",
    "\n",
    "https://psl.noaa.gov/thredds/catalog/Datasets/noaa.ersst.v5/catalog.html?dataset=Datasets/noaa.ersst.v5/sst.mnmean.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "69a5a5d2-9bf9-4419-9549-491980bfc92e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# location of ERSSTv5 data\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "22470693-d4c9-479c-b563-e6b2d0cebf91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# open dataset using remote URL\n",
    "\n",
    "# display dataset to see what it looks like\n"
   ]
  },
  {
   "attachments": {
    "022a5d1d-b083-4383-8c15-caea17b25821.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAFICAYAAABePvdzAAAEGWlDQ1BrQ0dDb2xvclNwYWNlR2VuZXJpY1JHQgAAOI2NVV1oHFUUPrtzZyMkzlNsNIV0qD8NJQ2TVjShtLp/3d02bpZJNtoi6GT27s6Yyc44M7v9oU9FUHwx6psUxL+3gCAo9Q/bPrQvlQol2tQgKD60+INQ6Ium65k7M5lpurHeZe58853vnnvuuWfvBei5qliWkRQBFpquLRcy4nOHj4g9K5CEh6AXBqFXUR0rXalMAjZPC3e1W99Dwntf2dXd/p+tt0YdFSBxH2Kz5qgLiI8B8KdVy3YBevqRHz/qWh72Yui3MUDEL3q44WPXw3M+fo1pZuQs4tOIBVVTaoiXEI/MxfhGDPsxsNZfoE1q66ro5aJim3XdoLFw72H+n23BaIXzbcOnz5mfPoTvYVz7KzUl5+FRxEuqkp9G/Ajia219thzg25abkRE/BpDc3pqvphHvRFys2weqvp+krbWKIX7nhDbzLOItiM8358pTwdirqpPFnMF2xLc1WvLyOwTAibpbmvHHcvttU57y5+XqNZrLe3lE/Pq8eUj2fXKfOe3pfOjzhJYtB/yll5SDFcSDiH+hRkH25+L+sdxKEAMZahrlSX8ukqMOWy/jXW2m6M9LDBc31B9LFuv6gVKg/0Szi3KAr1kGq1GMjU/aLbnq6/lRxc4XfJ98hTargX++DbMJBSiYMIe9Ck1YAxFkKEAG3xbYaKmDDgYyFK0UGYpfoWYXG+fAPPI6tJnNwb7ClP7IyF+D+bjOtCpkhz6CFrIa/I6sFtNl8auFXGMTP34sNwI/JhkgEtmDz14ySfaRcTIBInmKPE32kxyyE2Tv+thKbEVePDfW/byMM1Kmm0XdObS7oGD/MypMXFPXrCwOtoYjyyn7BV29/MZfsVzpLDdRtuIZnbpXzvlf+ev8MvYr/Gqk4H/kV/G3csdazLuyTMPsbFhzd1UabQbjFvDRmcWJxR3zcfHkVw9GfpbJmeev9F08WW8uDkaslwX6avlWGU6NRKz0g/SHtCy9J30o/ca9zX3Kfc19zn3BXQKRO8ud477hLnAfc1/G9mrzGlrfexZ5GLdn6ZZrrEohI2wVHhZywjbhUWEy8icMCGNCUdiBlq3r+xafL549HQ5jH+an+1y+LlYBifuxAvRN/lVVVOlwlCkdVm9NOL5BE4wkQ2SMlDZU97hX86EilU/lUmkQUztTE6mx1EEPh7OmdqBtAvv8HdWpbrJS6tJj3n0CWdM6busNzRV3S9KTYhqvNiqWmuroiKgYhshMjmhTh9ptWhsF7970j/SbMrsPE1suR5z7DMC+P/Hs+y7ijrQAlhyAgccjbhjPygfeBTjzhNqy28EdkUh8C+DU9+z2v/oyeH791OncxHOs5y2AtTc7nb/f73TWPkD/qwBnjX8BoJ98VQNcC+8AAABcZVhJZk1NACoAAAAIAAQBBgADAAAAAQACAAABEgADAAAAAQABAAABKAADAAAAAQACAACHaQAEAAAAAQAAAD4AAAAAAAKgAgAEAAAAAQAAAf+gAwAEAAAAAQAAAUgAAAAAxb2QJgAAAgtpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDUuNC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6dGlmZj0iaHR0cDovL25zLmFkb2JlLmNvbS90aWZmLzEuMC8iPgogICAgICAgICA8dGlmZjpPcmllbnRhdGlvbj4xPC90aWZmOk9yaWVudGF0aW9uPgogICAgICAgICA8dGlmZjpQaG90b21ldHJpY0ludGVycHJldGF0aW9uPjI8L3RpZmY6UGhvdG9tZXRyaWNJbnRlcnByZXRhdGlvbj4KICAgICAgICAgPHRpZmY6UmVzb2x1dGlvblVuaXQ+MjwvdGlmZjpSZXNvbHV0aW9uVW5pdD4KICAgICAgICAgPHRpZmY6Q29tcHJlc3Npb24+MTwvdGlmZjpDb21wcmVzc2lvbj4KICAgICAgPC9yZGY6RGVzY3JpcHRpb24+CiAgIDwvcmRmOlJERj4KPC94OnhtcG1ldGE+CqZd9jAAAEAASURBVHgB7Z2LdeQ40qXVe8qFXiMkI7rGCJUR8xshlRH/GFEyYqqNkIzYMWKWN6UrhUIACYAACCIvzsnCKxCPD0wGH6nuP/67lBsVERABERABERCBqyHwf64mUgUqAiIgAiIgAiJwIaDkrwNBBERABERABK6MgJL/lW24whUBERABERABJX8dAyIgAiIgAiJwZQSU/K9swxWuCIiACIiACCj56xgQAREQAREQgSsjoOR/ZRuucEVABERABERAyV/HgAiIgAiIgAhcGYFvLeL9448/WqiVThEQAREQAREQgTcCe/4bfU2SP/za49QRO8sLlhK/sZbrbDs1jj22U21ITgREoD+BkvNBfy+vz6I/5/q+JxLbR78uJuf11ejTdqmuZsm/1KEj1tXcMFwE5OjjBvLi4Yj4ZVMEREAEro2AP+f6/hYPe54/43n86t/52w3c2mzOPz4+shmseQHAA4JCvo9xyOYedNTXot6KrYXNnjoVX0/a9W1p/+oz7aVx1L0ryQGeGc7hDw8PfjjYh71QLggKNxz8Y3G6+v/YpwbMhjF/Ul3D1zUdfpMb4P4UjzoiIALjEFg7N4zj5fV5YvfFtkMk7Dmc5++tNdBj11m90JGy3q4JtffquOrH/nvhYUO2dNiDJbSBGhMBERABERibgD2Psx3zGDkBZUsutj5lnDZSZGMyUyZ/Cya0AZwPzcVA+XHoyFkPWaypYdv7or4IiIAIiEAaAX/u5rk5dj6PjYesUbc914fkRhib9p1/bMO4ObH5lE3huyvo2io8CCgLu3tsb9nbO8/Y9uoZdb3iG3Vn0vw60/7xXJMW2avUmeLLiQuyM8SGc3dsXxGfnee5Plbn8qM8cwn7pfV07/z9xnhQ2Ii9xduI6UuVi63XuAiIwHkJ6Ps/5t6F9iU0luI91qEw6fs1NfKN10lfWfv51P5Uj/25ETb42vBpg7W1hTbt7d0Yr1d9ERABERCBcQjYc3ysXdvbWN4psTNF8icQJt4SECVrQvboS4k+rREBEZiHAM4PNinME9lckeTuEc/x/vwfG69Ji8dUDZ2nf+fPjfMbUQNOSAfs4O85Y5uA8V6+hPzbOzbDe7k1Bopvjc74c9q/8fco5uEoe2fP3UzYMZ85Djl+7DmeY5BjXuCaFjXsWf/32Dh18ieIPQBK18L2WjnzBcBaXJoTAREQgRkI8Byee67GOn54IZCro4Qf/WVdosOuOf0P/gCiB/hP0NzVV2/71he1RUAExiVwxPlpXBryLJVALMHbXLP32Opy5x8LJBXEiHIzxjQiZ/kkAiIgAtdEIJRbkPRt4q/BI5r8rQNo8wOjfi7mCOXgNNsx2dxx60/u2r3yfLdTezP2+lVj/Sjv5WrEEtKh+EJUzjOm/TvPXnlPtXeeyNd+KE+2yjPBX/vHHLDjaK855ed5AbC25iuK8IjXHZbqM1ozrj4ey4oIiIAIiMBoBJhfa+TIlNiid/5cnJNoIcsAuB41rvisHvTxYfF9jNt59qkfcFLX0AZ1WL2lOrxO/39zsjZo16+xMql+lKzJtWvla/pu9abGu7amVMeazli8tdeU+p6y/17mTL6v+UpmW/Gl6PAyW31v0/axNqWfImN1Qd6u8X0rizaKlWffjtXQQb2oWawNjHk77Fs52+Ya6mPfylCHl9nql+goWUOfWZfosEnfx+v7Nu7S9pcf/IWSN51iArc1DHOeTnA+1sf4mh2uY01Zb4fzvWsfH+yHxnr7JXsiIALjEdC5Ybw9GdGj3Dy397j65iEwwVrFtu3lc/oMDmtox67nfGguNGbXqi0CIiACIiACZyDAnMqaPvfMc5uP/ekUEzP7qEscxZq1daE5jIXsW196teGHf8zfy3YPO/ZxVQ97vW0ovt7E69rT/tXl2VOb9u6Vts9l6OMTyn0t9+fLY/8axmwgNtCt4Oy6kB9b86E1tcdCPoTGatuVPhEQgXMS0PnhnPvWymscDyj2pnYrN4Z82XtcJd/5h4ynjCEoflLk12QsrDW5VnN7YbfyS3pFQATGJXD0eWtcMvIMx8ZRpXnyTw3sjIkVPuOjIgIiIAJrBOwFgM4Za6TmnrN7b9tHRD1M8k8N3n6JUtfUkLMbZd9dwZ8jr95qxGZ12Njs+CxtxXfunTz7/vE8wtrvxtnj8/HY/syxIc6U+JgrWOM4YNuy6tFu9s6fzucE5r8Qa2t7QovZio0zdtUiIAIi4AnwPLd2fvNr1J+HAPPG3uOAekrJfPlTv1JFfh0PbAZo5zlnx9D243uD8/q3+tZX+mLHttZrXgREQARSCPD8kiIrmXkJHHkcNH/sj+AYIOsaWwldrRIzdfPio6bfNWKXDhEQgXMTwLml1fnr3GSuw/sR9r558udBzoSamkixLlW29uFC2zH7eLdDmdq2j9aX8t7qaB/32Fd8e+gdv3aG/cN5Ze3ccjzlNh7MsHdrZFLjW9v/Nf2155o99oejPkHGDvhQUJD160NytcZ62qrls/SIgAickwDONyoicCSBoX7w50FsJeStea8v1ucXMefipJbtmE8aFwERmJuAziFz7+9adCU5x+vbe/w0u/PPSaQ+qCP6Z/P3CEayKQIiUI9A76eb9TyXplICTPql62uua/7Ov6azJboAe+uToxe6+G5npI3MiWFNlrGtyZx5TvGdeffS/pb6TBHyAoA+z3x8zhwb9i8lPuw3P9zzo+pmd/57A0JiBaRY2Zq369b0WLm1tk/0NXSu2dOcCIiACIjAfAR8Ljkqwmbv/Pcmx1hyJ7gU/TEdMdgxeWvTtmN6NC4CIiACOQRi554cHZIdnwDzBzxNyWFrEe09Zoa987dB5wKjfAgu56x+KxeahyxlWNv1aouACIiACIjAGgGbW0bII82SPwLdEyDWElaqHspjA7jGjnFj7Bzbfo591SIgAiLQmgDOQ/Zc5c9Lre1Lf1sC2NvR9nToH/wB1h5g9suEraU+P4650BjGQyXlhx2hdWcYmzk28Fd8ZzgK4z7OvH84Pz08PFzOUzgf8ROnca6ZmfcOOxGLD/u4J4+12uVm7/zhcM+AQ8kb9tfAc451K8jSKwIiIAIlBHhe63kuLfFTaz4T4L5htNXe7c1bzR77I+i9zkHHViFkC5h2WUOHleM4Lw62bGheBERABI4gwPOaPX8d4Yds5hPg3uWv7LOi2WP/HoEziW+h8l8c6xvatr+lS/MiIAIi0JuAzlG9iZfbY74p19BnZbPk38N9fCEImjXshsZrfnli73Z6xNzaxsyxgZ3ia30EtdWv/WvLt6X2a9g75CHkmpr5ptWeNHvnzwTcGoJP+gTFTUDftjmvWgREQATORkDnsnF3rPfe7LXX9M6fFwBwslWBDX+B4aH4+Va+SK8IiIAIiMD1EfA55wwEmiZ/AOiReM8I/gwHh3wUAREYiwDOpy1vpsaKdnxvsBdnzT/Nkz8PVNatt7PHRsz87mrm2HDsKb7W38C2+rV/Yb447/U6x4Y92B6dde94g3u2+Jonf4JhvX2ISEIEREAERCCXAJ8KnOFCIDc2ydcn0PQHf9bdlnfk9opXFxmWutoiIAIzEvDnU99HzKGxGVkcFdPRfPfaD975Qyk/BIs+S6zNeVt7PXauVhsJn59aOqVHBERABEYlwLv8Uf0b2S+bv0b2s7VvweTPu2fWgFVysHEd1lJX64CkXwREQASujUDJ+Xl2RjbJo82PjTs0ZudnbgeTPwMmmFDixtxa4Tx1rMmebe5sP+zI4TtzbOCg+HKOhvFktX/pezLaBUHPvUPeYfxso89cxrwUGksn/FmyZ3yfLZf1osmfkHLUEijWECpq/J+qMMcCSBaU70POzsf6KTK0SR12Ddq2Txm/ZqtfosOuQdv2Y35YmdQ13vctHVa+lx8xO9YXH6/vp+pY05mqwzIsWVPqu7Ub0+Flascbs+vt+L73y/ZjDK0M7foxb8f3vbztx+y21rFml/6hZpvyto8xf15F3xbI2zH0vY6UfoqMt2vXoM0P5ew8xkJ9O4Y2+8wz7GM9xmwfeQf9UPw2N/k1kMeYtWFlrK/QY+1SDjXblLd9jqFmCc3bMbRtn+v21F9+8IdgEBQK22s15CiPNgrlX3uv/4bG7LzaIiACIiACeQRwXkXhOTh0ng2N5Vk5VpoxWi+YeDlm+2TBuZz4rS2rJ6YjNk7bLeu9tr+FnLMAOB8a45xqERABERCB/gSYoJgImAQ53t+j+hZjMfkYKUcPSnKW10ldM9ZfHvsjePtB0KG+Hd8CU7IJWzqPnK/9+OXIWLztmWNDrIrP7/i5+tq/8H7hHM3zrG2HpY8ZLd07Xtikeg15ywLr2E/V4eViTO14aXzeVq/+l8f+NQz7zfL9GjakQwREQARE4IMAExwTEmqWM56DbTyM4+g6xDE01sPPvXa/3Pn3cFo2REAEREAE6hKwyZ4XALTg+xwftWZiszGN6uvZ2JJj8+SPTVQRAREQARHoS8AnJd/v602aNeQLJv60FZIqJdA0+XMTz3D1lgPwbO92FNsHgZn3DlEqvo+9PmOrdP+YNLdi5gUA5HuXrdjOki9i7Lbi6817y16z5B8DtOWQ5kVABERABPII8AaLNVf7PsZDY5RXvU1gFn7NfvAHhLNA2j4cJCECIiACxxLgnbP3Infcr+/Vj/nZy36KnTUf1+ZSdOfK7LXX7M4/NxDJi4AIiIAIlBPAzRYSAj62cNyOoR0b93Lqfybg+X6ePU+vWfKf+a7/bO92cg7HmWMDB8WXczSMJ6v9W98TnHeZ1G2S4lhotZULzdcam2Hv1vKa/U8Jg6n9gGEvzqn71Sz5pzogOREQAREQgboE7EUANYcuAEJylFcdJhDi6CUhg8Laz4/Qb/bOf+SgRwAvH0RABESgBwHccfJ8bNu0HRrjXM96FD9SYg75GhqzurbmrWxKe68+3fmnUJaMCIiACJyUABI/EoVKPQIlTHkBVs+LfZqU/Pfx02oREAEROA2BUNLimC4Q8raR3LBq7114nuU60s2TP6DwU8fl47XM8MOVGMWZY0PMii+28+cY1/612SckMpvMWlhZ27szJk8ygu8oofhGzn1N3/nbDbVtQlMtAiIgAiLQlkDo3Bsaoxdrc5SpXR9hs3YMIX2ICwUXVrXLXmbN7/xrByx9IiACIiAC6QR4R89EtJU0KJ9uQZJrBEbl2Sz5bx1ga7A0JwIiIAIiUI8AEhCTEGteDISsUCY0p7F0AiNzbJb8LZ7ZLgRC73ZsvGduzxwb9kXxnfno1P7t3T0mI9Rsx3RivmbRd68mzf26mr7zp3uzJX/GpVoEREAEzkiA52TWazGkyKytT5nrYSPFj1YyLeLbq7PLnX8roNIrAiIgAiKQT4B3/azzNWhFDoHaT1FybMdkv8Umao3vvTqp5Yf0iIAIiIAIfCaQen5Olfusfa4eGKyVERP8mr+681+jE5mb+d3VzLFhOxVf5KA+ybD2r95GIVmlJCzKbSW/Lc/OuneI2178kIevEd9eRlsMa843u/MnhJSDq2ZA0iUCIiACItCGgE2CtSy00FnDN5/D6CfHYcPnN/TX5mv4VUtH0x/8EVYtZ6VHBERgncBZTjzrUWh2VAItzuktdNbgZ79L0GcTPefsmLeZIuPX5PT3cmt2558ThGRFQAT2E+DJJqRp74kipFNjIlCDAO+W1xJpDTs5Oux3iX5xDH0/Rt0cR59tu45yI9RN3/kz+BECrenDWd9dpTCYOTbEP3N8+L49PDxcTjr2u8eTT8r+jy4z8/6d4fjEcVV6PK3t3R69qccs/A59QuvhDz+heY5RBjW+e14/5CjDOa49utad/9E7IPsiIAIiIALVCSDZ2oIk7IuX8fO+H9JhZey81805jLNt17LNedYcr103fedf21npEwER2E+AJ6W1E9B+K9IwI4GWCamm7lRdqXIlewndKKHvGedC85wLrbsofPtnr+/Bx/5Qyg/ssI2a/UvjbY5t1SIgAuMT4EmF3+fxPZaHIxHQcZO2G/ie4ePzJ1bzO4g259FG4brXXrt/vyR/OBIybp2FOzoA2m2KNItAawL8jut73Jr0fPp9LqgVIfT2PB6Z62r5H9PD75qf5zh59owdvnxJ/t6RHECQDQXgf+iBvh3zfThm52P9FBmsZfF2fB9yKTqtTE8dR9klP/Kp7Qf1eju+X9tujb0bzXfPaI0h5vgjJcrZ9RjzjHyfMlwf61u9qTpCa/xYC7tep7dp+7F413SsraFu1GxT3vY55u34vl2Dtu2X6ChZQ7vWtm1bncg/yCFcw3h8366xMmxz3tpJ1VGyhvZY5+rA9xCF+dP76vsX4Z3/RN/5wwluhL0gsGN0lPP0hWvZVy0CIjAeAfs9te3xPJVHIxFofazU0J+iI0WmFnfYWis2h27JWj12nR1PaQd/7W+h2GSfolAyIiACIiACIlBKoEfOsTmu1M+cdTlJmrJbPuZcJIR8/fLYnwpRs42Ftk1FdJJ91SIgAucjEPpuny8KedyLAJPzqMfNqH7l7g85565Llf+S/GHQfqAo1KcBzF1bse9zZot95tiwV4rv44jFdxcnSn6/P2bGbWn/xtgbHjM5iTZn73hs5kZ75PGcE19uXC3ko+/89xjjBvDAwEaicHyPbq0VAREQAREYh0DL8zpzCKJlHmHkdo5jITk7x3ZLn2mDtfXTx0AZ1pS1ciFfQ3LUkVoH3/mnLl6Tsw7T0TV5zYmACIiACJyPABIVz/E2adWIhPqgnzaol3Psp9ZeT+q6EjnaIiPb9/ow5+VCMn6stN8s+VuHuEkM3M6pLQIiIAIicG4CoXM8x2pEVksXc1AtfVuxwY61aRO89yEkC/0cpzxr6t3yITb/5Z1/TFDjHwTO9m7nw/Pt1syxIXrFt30MjCyh/Rt5dz7/Pswnp1H2jsmzNslYfEzeNvHHbEOW/vk218fW5o53ufPPdUryIiACIiAC5ybApMdkdu5o9nlvGdh2jtbSdTEbzX7wB4MhZ2tfvcQC07gIiIAIiMDxBOwTgFBO6OnhTPlnbyxdH/vvdbbnQSJbIiACIiAC+wkg4TPpIwfgo3I8AT32P34P5IEIiIAITE+AFwAIlBcAdmx6AIMF2O3Of6a7/tgPOwbb2yJ3Zo4NQBRf0WExzCLt3zBbke2I3Ts+DUh5EkAZXjBkG35bsHf9ll0b35bsCPPN3vnbKzpAt/0RApcPIiACIiACYxDwiRn5gmM2d4TGUiLw69AP6fW6rIyfO7rvY8j1R4/9c4lJXgREQAREoCoBn2RjiY1ysfmYU3YdZNCHDhbOs8861w7XnaFW8j/DLslHERABEbgiArFkbBGUJGar17atXtuGDC8SUuTt2tHb3d75jw4ix7+zvdtRbB8EZt47RKn4Pvb6jK2Z969mbEjENjG33OvUpF8zvpbxUHezd/40wA1KBch1qkVABERABERgBAIlTxla+73Xpy53/rwAaA1D+kVABERABK6PAB/NX1/k5RE3S/6xO31sEj/lbmulCIiACIiACLwS4A2mLgLSj4hmyd+6wA1hjY0682ad7d2O3Yut9syxIXbFt3UEjD2v/Rt7f9a8a713Pq8w36z5VHOudXw1fYWuZr/2J3jU2BSU0BjnLgL6RwREQAREQAR2ELA5hTnHju1QPdXSpj/4A3APnxcDHAdNbcxUx5SCEQEREIGhCNh8Q8dy8g7zFteOUO/1qdljfyZ+D9j3R4AoH0RABERABOYlgLxjP4g0dEEwL4GvkTVL/l9NfR6xG8FNYP1Zcrze2d7t5BCcOTZwUHw5R8N4stq/8fYk1aOR9i7nJjT1Dnuk+FL2pGnyB+BUcJQ7ywVAClzJiIAIiIAIjEmA+WlM7z57hbzIz+eZ8l6zd/65V1YIgZuRs7Y8dK0UAREQARG4ZgK86VxjkCKztn7vnLXPNmqUPbmy6Z2/DRrO8mPHGQATv59TXwREQAREQARaEBg97zDJI3a0mexZ72HSLflbZ/1FAPtWZk9QWisCIiACIiACKQSQd2ySTVnTWsbmRPpn82MNf7slfzgL5xkAawaBPtutwe7Vf7YfduTEO3Ns4KD4co6G8WS1f+PtSapHI+9djfyzJz7kPvuBP8yR5OvnOV5aR9/5wxCNo81CSHaObcrYtaEx6qMuyqD2uuyc2iIgAiIgAiLQgkAob8FObLyWD6n6vZzv5/oTvPOHUl9qJGXotbrZrqHb+6u+CIiACIiACOwlELpJ3auT65EDj8p/weRvnYk5x8TNILZq6OQHsradq2vLluZFQAREQAREYBYCLXJkMPmXAoODISftuxDMPzw83GAMbVwEoI8P2ix2DcZC/dAY16eswfoz6bC+ns13vy82lthera0ZJf7RfLdcbZt+hsY8Z9+3a9C2fer1a3zfrknVEVrjx7wd3/fytg/ZlL6VQdv2S3SsraFub8f313RgDsWv8X3KXITf/oGMLTXWUIfVbduwF+rbMbRtP7ZmzfdUHZBjHrNr7Ji3Y/2xaziOMVtS+laGOZI+WF3F7SXhBsui8H2c7VDNsXfhpREa4zznUNuPn2dftQiIQFsC9/d/Xb6Lt7e3nww9P/9+/45ygrLLyYhDzer7+/t3+2inFvgGP3/9+pW6RHIi8IUAc5WdCI35eeY1Ox5rb+njOsqxxrhtUy6nTr7zr3rFsXiNAv/5eR15/beFLatfbREQga8EXl5ebuzdxleJfiM/fny/eXp6ejeI9t3d3Xs/1nh5+fvm58+fy9q/b56fn2NiGheBTQLITT4Xhca8IuY0rPXrKcs5yB5VosnfOsVgOGZrtlMDgDwC9+ssJMqk6pScCIhAHQJInLgIiJVfv35fvruxiwSsXVuP5IzPWsF6JG+U5e79PYm/6l5fe3//zzXVmhOBLAKxXGTzFRM58xr7WBtbDycwl1rW9KTq8HLR5O8Fa/ZjgRBeTVstdMVOfC1s9dY5c2xgqfi2j6g1Rrgjx/eUMraPO3P7sZYgj3V3d98vH7Ttnb2VXV4/LAn/9+V3QMvj/hv0WR4f/5fNLzVsrF14fFkw4AC5DujabpdmiQ35izmMOQt9vJe3fbRRKEuAlGE/tabdVPktuUOSP4MnHNYITkUEROAYAkyySMpbd+feQzwxQOJHskZ5vXt/fWyPkz7mUXCCpMyPHz+idm5v/7rBOujBBQaKvxC4DL79AznYQAyMw86rLQJ7CDBHsYaulGRMeV4AoD9MnlscqV5euXxVi3F+OLsmSxnVIiAC7QjYH/EtifPyHUWd8oM/rkXNQh38USC/8/YHeEsiv9hBvVagg+vhT6xQH2zQJ9qPrdG4COQQYK7i8WjXcs6Ooe3Hfd/L5/T36jrkzn9x+lLWroLAjVdNlFctAiLQlsDT078uBnAn/fj4kGzs9vb7u+zd3Z/vbfsEAXfuLI+P/3Npxh79Uw5PB/DeHwWvDELysIHxJem/P1XgetUiUIvAay7/uONHfjpzjuqW/AmJANE/a4LH48hZy8yxYc8U3/qRi8ftTNL80d36inazuADBI3z4g9cFKKH940XK8/N/Lq8I6DcuCPjKoJ2XdTWH4qtr4Thts8WG/MUchjolPsiNUrolfx/wGgReKPg16ouACLQnwDvtGpZwMcFinwL8+vXvyzCSuy84ieIcwIsQzONCYKtAhomfa2x/a73mRaCUQCxn8eKgVG/LdYck/5TEzycDLYMv1Z1yhVeq++h1M8cGtoov7QjjnXaa9LoUHsej4M/wcDeOD38AGNoPPOpHQTJHGzJYg0J5nB/wwQUF//wQ5xV8aA8xrJ1rLgoH+4fxDeZWFXdmjY3HXSw+zMcuDqqALVTyrXBd9rKcL+GosLKD1gIROCkBnMiQcFPuuLdCRHLGXwIwmVMeydne3XMcTwPw9AGJn0kfc5ANyXOdahEQgXQC0f+lb7qKr5I5d+1eln3WX7VrRAREoCaB1wT/n0Xln5/+TO5j/GYZf717/xh7lfV9+BUawzgSOf+re0jsoUf+kGOBHr6C8PJ8hUC/uAZ1zL6VUVsEehJokc9261zusquXBWqWTsjzg4W567OMVRBe7lgqaBlTxcyxgbjiG/O4S/VK+5dKajy5a967Fjltr87D7/yXAN4LrmRQlsP2fUwNERABERABETgzgd136YHg9+rs9s4/4PuXISX9L0g0cGICvJg9cQhyfRACOjcOshGFbmD/9ibrQtPRZUMl/6iXmuhGIPaL1W4OyJAIiMAXAvpefkGyOiBeq3guk4f8qd+2W2NLzHxgPT39Hhv+Tu9mj28nHi0XgWYEZv/unS0vDPXOv9lRJ8XJBN4P4LvH5DUSDBP4+eP1NyyYffg12e9YnvGngL+XP737fvMYDl+jOwn88fY/Q4Ka5ZeqO7Vdx/LHtzDfz2MDhV37sf9efbrzH+jgkCsiIAIiIALzEdibqFsQUfJvQVU6RUAEREAERGBgAkr+BZsz4iOlgjCCS2Z/Lzf7f+t99vheJv9NyszxzRwbTqZreYG/9g+edA8aVPI/CLzMioAIiIAIXAcBPfafZJ/XrvDOHiJ+wDVz4f/0ZdYYZ4/vdvLjc+b4Zo4N55Oz5QX9nf+sWUBxiYAIiIAIHEoAd/yjFj32L9iZs13h5YSod/45tMaT1Tv/8fYkx6OZ34vPHBv2OJYX8L5/xP9C46mS/8hXUTlfcMmKgAiIgAjMTWDE9/yW+FDJfy25r83ZgHq0Y1d4PWy3tqF3/q0Jt9Wvd/5t+bbWPvN78Zljw3FxtrwwTPJPTe6QS5Vt/UWV/m0CT4/fb/Bfunt6vPsk/PL892Xc/lfw3mWfHj/Jtuw8LbZgF/6oiMAsBO6enm7wXwi8M/+VQMT28vL7Mm7/64GUfezwZ5Q/FhuwjQ/asxbkqBEf9VveTZM/E/VWst4CxflR351YoGqHCby8vCz/OdjH8ORBo0j4L08/lxPi3zfPL/8+yAuZFYF2BF4W1T2SekoEuMh4Wr5rLGhjTOUYAk2TP0JiwrYXAr4NGYyllFS5FF2lMmd7vJMTZ8sf/CHRrpX7x9+X/wb+/f1jUOySrFfu0LfmodT+IO756Z9BO2cetPGdOY6Y77P/aKxFfD+XJIuLgFh5vr+//L8DHiN/Rrlcty8XyLHVr08T8ERhrWA5bgBQft3/dfN8//r/KsDY1trLohP8w7yAHDX6XT9wJv2pn024PlHHgrRrYCgmhznKxmRok/OssfaIwk0+wnZrm63f+eMRO5J8qFwevy8nqtvlxIALANu/eX56P3nc3t4uOp7fVeCJgr+wuH/4fXN799e7DBt8J/6CNWtnNC44Wc34TuZ2sruzvzduFR/u/n9FkjvuvvFdeLj96wYXALb/dPOf9+8Jvne4UGCBTlxYfJS/L4n9/vb7x9Bb63apn5ek/2upv85/f5M6d/VzeZWBz9H5KZVi8p0/ArJBsc/EHTOI+a0PdcV0cHzLFuVUj0cAJw4UPGLPfb9+Sex395eLglcdy93C2xMAm/hx0YAPytPP7+8ylwH3z9PyFAI+0S83ra4ITEHg9Vu3fB/wvdu4O/cBI7Hf3fy5JOvXi2hcIDy96bCJHxcNlPnxFLdzu1wUPC4f3P/zcT/WvZ0avPniPhLwVs5pMQ+HbY4sDqDTwuTkn5J4CZQQmNQflv8dJduo0bdjuJP2d9O2D72Qp17UmPcytk8Z1Cx+je9DbmQd1rdWvvvH/v4xMvp2zPfB0M6jfylI3m/fciRm/47dr7H928vJ5x+Lmn986Hj698UO7/hxp495PDF4vwB4ej1mLvaXf+grnii8lvvlhPjnW/u1snYxEuqHxqwSP+/1YN7LpPStTA0d3q9Yf8uuf1Qd6tsxtG0fdkN9O5a6BrpY/Brfh5y1EetbmRo6YnYwzmJtYqzELtbhPp3fu/vluHv+eGCG6dX4se5xkcGHOrAevvCOH4/wfyzzeKrAC4BHXAAsMize91/Y/7enbvwdsJXHOr/G9ylDG6G+nevZ9udm34cvGLMl1LdjaNu+XVva/mNJxv/dWozka8Vs37apJzTGudIaOlGsH1YX560Mx2Jr7PqcdouNyLHfUvbubvkS4/Hg3WMVM/bR/d3tPy535FCMhI6nACj8f91bWf/Yn78DsDJWH3VAH54K4AIDxY6jzycFsI/XD1YfbUCuRrF/yeD9qKE/pAOJusuj/+fHheXvy7HyGHKk0RiSQKtH441czlJr48Mv4ln++3bzw/5WbR/d/1j+0OZuOS5QkMSZeKnTyvrH/vwdgJWx+qgDuvFkgXbsOOZQGBvu/J8XWTwlQMEFxNdXAZep7H8elxW48z+q1M41a3HszbPJd/4hJ5hcQ3NrY1iXu3YLKub5gW2CwRjtsaZv7KNWaU8A7+B5V87E397qZws4Ab2W/7wn/kt/+U0BLgRURGA2AnjczrtyJv6jYsQNP15FINnjdQHKk3sacRls+A8u/HERZD8NzQ2rOjn5M1EyEibMWFLmPOS5ljV15NZM5GvraMP6hTY/dq2VseNb7dqPX7bs9Zxv/YO/mnfX9gd9/A0AWPGVAh9VWn5vbx4udz/2AgQnRdu3a87U7nLXfyCQme/6gbVVfHgsX6vgYoLF/o7g11sSD33vHvFUYHnSd2/+4gdPAFCelx8Vtiw897N+fv69vP54/vThXEnd0veWupOSvwUCZ3x/zUEkYy/P/tq62BzWxoq1FZPx49SHtSp9CNy//SivhrXX3wMsJ5Dlz/ZwAXD5vJ1g7u7/9cUE/5wQj+Hx4Xo8kej1aP6LUxoQgQ4EHpbH67UKE/z98/IEbUns+PB3AI93f34x8+PtggEJ/8fy9A0/GOTf/IfkvyjQQHUCScm/xCqSqk+oTM4l+rDW6wvpoVyKLNeHLkaohzK2nvnO3//gz8Zdq32LH+XxFnynUiRz6MKdO97z810/krl9MkAz9sdrHJupnj2+j9c2M+3aRywt48Mv7fnr/w+LZS38yR+/d3h3z/f3r7/8//5FKey+P+ZfngDwQgGvI2q97/9iVAOrBL6tzu6c5F011CCZ2n6OapvI1/RY/XZNji3vK/VY3bn6rln+7n75T4wuAHwixt/p20f1YORlfT8kgzHq4uN+/BDQ24NcqIRshOQ0JgJnIvC0/HXNzd3v5Yv3/ZPbT/jh4HKXbouX9X3IhsZwAYA7fv4FwY/lx8Jr1/Q/Fj0/Fvv4tT/KlvxFSP80I5D0a/9c6z5B+36uvhJ5b5NJHLpiiZwyofm1uRL/Rl3z/lSj0q/9R42zh19H/Nq/R1wXGwf92r9bfAMY2vNr/wHcP8SFx8Wq/7V/6Hxe0znmBuhsbcv6Dbt77DW986ejcHCvo9SVWtMm5W0/5EtojGtREzLkWDjGvmoREAEREAEROAOBZu/8ETwSpU2Wtl0DTo4+yoYSNuZC4yEfIccP1vETkq09xhhq65U+ERABERCB6yLQLPn7hJqaXNfwM9Gyhk62UdvCPhM1a8hwnZXPafPRuNVJP3L05MrSb8aWuz5FvscP/lL8aCUz+w/iZo+v5Q/iWh1zOXpnjm/m2HL2eBTZZskfASJJIWHlFiZS1n691WkTsJXjeCxRYp76rT6rI9Zm8uc8baFPnZyrXdPXWFx77bX+O/+9/u1dP/vfwc8eX6u/g997XNVaP3N8M8dWa/976unyzj83ICZmm+h8gqUMdbPPNRy3NRKmnbdtK1fatvqYnO1YqV6/Djp5kdFCv7envgiIgAiIwFwEmt75M0kRme9zPFbHEhwTHmp+mGzXdOXaj+lKGbd+bfmWos/LQD8KGfl59UVABERABEQgRqBp8ofR0oTL5MkkFwuA45Rn39Y2+dIfO2ZlU9r+sf/aGvq1x57X3zLh652/p32uvt75n2u/vLczvxefOTa/j2fod3/sz+SLukcJ2eGYTcgcS/EpJ/lTH/TDXo4drg3V1Bea2zOmd/576B2/Vu/8j9+DPR7M/F585tj27PlRa7slfyZaJr+aibAUHn3Beu9fqc61dUzY1u6avObmIWD/gz/zRIX/oMryn2qdKSDFIgJXQqBL8meyswmW7VE4h3xs4dveCwByQ02fW/h5s/wX3FREQATGIPA4hhvyYiICzd/5h1gxgYXmjh5DQmWCjvlZ8ti/RlxM+C2T/vzv/H/X2ArpEAERyCQw+7klE8fh4l3u/Bklk2rL5EVbazWT6JoMfV2T2TMX8iF2seHtpMhRJpc13vkfdXHj42zTf+wWn/9vjLeJR1qvgcAc38k+3z1979K+EV3+xz5prvSRYtJlvWU1VW5LT2ieCRpzNkmv2Vybszao2+q182qLgAiIwIwEeO5jbK3PgdZea1uMCXVqLrBrbLvrnb81fGR7L7RavtsDxfqEcXtAWXt2jR0PtaknZ01Ij8ZEQAREQATmInDIO/+zI+zxCA4JO/RJZWcTf+xCIqSrR2whu73GFF8v0m3saP/acO2hdfa968Gwpo2rTP45d8JMojWhW11IzDnJ2a7datP3nHi3dGpeBERABETg/ASmfOePZFoz4dXWZw+blrqtHbVFQARE4BoI+JupmrkgxM/aa23L2t+bO4a+87dQbdBr7ZI1MX3QtRcwdcT8wsESm4v5pXEREAEREAER2ENgyOQ/UjKscSVHHbG4dAGw5xDWWhEQAREYh0DsPD+Oh6+eDJH8AcsD8/0ccEy2OWvWZOkf69QfrkCevqBmkg/Fxrk1P3rMpcbWw5cWNhRfC6r9dGr/+rGubWn2vavNq7W+Zu/8reNMbKxDcxyDDIpPkBynXOva2ve2bVIP+cG1fp2VDekIjdk1aouACIiACKwT4PmXUmvnYcrsqa092upxLt9ro+nf+XsQIUiAbuUoY8fY3rNBuWvXbGJuDTznc21Cfk1viT6tEQEREAER6EuAOWAtj/T16Ku1ZsmfwdPkVlLz8+ijjArPxlfiI9fbtRxD3J4HxlREQAREQAREoAaB6Dt/Jl8YQZufUD/mCJKZ/YTkqJdJ0PY5Flp35BjfXTE2+mzrFN8hgzW2hMbsfMs2fEFs3qeWNnvr5t71ttvLnuLrRbqNnZn3b+bYQkfDkefykD9+LPjOnyd/JjD0GYitvTL2/XqO+5p6Oe77HJ+9DsVNhoid+9Cag/XDtlvblX4REAERqEXAnjuhs/X509rztlqeR/fqDj72Z4L3m2HHtwxT1sPwOtnf0ke5a6nJDVx6FmtPe9KTvGyJgAi0IGDPaS30n1Vn9LH/WkBITEzuFiza7LOGHv+4B32OQQ7th4eHTyY5z8FQPzRGedSheTuGtu3H1niddk0NHYgdesiPTGiX/dp2ve/YU/jC/UVtbZKPHfM6KEPfY32rIybjddg1NezW0CHfw9+zGns34n7rmGm33zWOGavjqLY/RnwffmHMllDfjqFt+3ZtcXs5wQfLovB9nG3W7xNLI3XMrkEb6/jxc6P3lwTZ3EWyIV/WrQ372Lwfre231u/ja22vt37F15t4XXsz71+v2Ow5q2c7Fl+rc/devcl3/rj7ZLF3qBzLrZevzPsdZu7aa5AnH8Rq2feOHX6wHOkHfVAtAiIgAmsEeO7sXcd8gh8jnjuDP/iLBZE6jkBt0khdd3Y5bvA1xn72vZP/IiACItCKQIucuFdn8p1/Kyg19AICEy/rGnpzdHAjRr3Ky4lFsiIgAiIgAnMTmCL5Y4uYdFljjBcFtS8Iqv/wAs4OUmaODYgV3yAHWqEb2r9CcAMsu+a9s3lpgK24uDBN8meC5x04ogNw1q0uBC4G9I8IiIAIiIAInIjAdO/8eRHAPeAFAPuovYydC8nb+VDbXnBg3vdDazQmAiIgAiJwXQRq5oa9ur7NhD4Vhk/wWIeC8VQdM3FTLCIgAiIgAu0JjJRjpnrs75N66lZiXc1Noa5U+5ITAREQAREQgZ4Epkj+ve/WZ/7hysyx4Yul+HqeXurb0v7VZ9pK448f3y9PUu/u7i4muHcvL39fxvnEFZOUpUwrn6xe2oQ/uQV+wn98EN/T09NNT99z/Q3KL3ep1ctiqLrOLYV7bGItP1t2QvMh26Gx0FqNnZsAjxvt97n3Ud7XJ3B//9f7edX+1++en3+/j9MqZa0c51rU9/f37z7AH1+en5//i3HUvti19vv/69cvLxrt1zhf7NVx+jt/Xn0tIHaVZZcuj/53KdHiqyZg72SuGoSCFwFH4OfPnzcvLy9u9KP769fvy/k3dveMtWvrcfeecgePO3TeqX9Y/9paEvki9/0GtS/QgbJcqFx8Xi5cLn2OXzon+Gf45J+S3EdM3PBJyeAE34CdLvLYQ42iPd8JVMunJRBL7AiYj+ApY/tI1vZjAUEe3zkkanzQXkvCT0//e7mIWO7erZrkNi5AkPSR8Onr7e33t/X/L1nPCILDJn9sIj48uaLtC+f9eOs+Nz3VTsj31LW95XJj6+3fXnst4+OxutfHPetbxrfHr1prFV8tkv303N7eXowhKafcnVvP8MQAiZ/JGsmXyR3HAuZRXhPya0L/8ePHip3/e7mbD93Rvz5deH2CwKcMfgyxwC6eVLC8vHy0OXaKejlhVS9L4EU6sY6fkAKv1/dDa1LG9upZW8851in+SEYEREAEzk7Avsdfkubl3I465Z0/16JmoQ7+LgDnVHzsu3a+j0e9Vbie7/yhl2OhmnatXtqDb7kFNvaUveuHu/NfYFzu9pfAvhTM4S6aH/a/CGYMtH56AB9RaviaEdYQotynIZyREyIgAocReHr618U27qQfHx+S/fh4pH6zPAH4832dfYLApwKYfHz8n4sMnw68L0hovD5h+OvySH9J5pcVqPGIHx/M24InDLTD2s6P3h4u+acAQyKtlUxr6UnxWzIiIAIicI0Ebm+RQF8fyz895f9pXQ9m8A+P8/Ghr6Ex+IJH/0j4uDhY/iLgUuf6eHTuGSb5p96BAxgL2rbP8ZIa9lMLNn7W0iK2HLatubaIr7XPOfoVXw6t8WRn3j/eTdegjosJFvsU4Nevf1+Ga9qiHdZI+vitgU/8JXvHC4AjzpHDJH+CPaLmRUTrDeBGHxHjUTbJ9ij7sisCIjAOAfwwr1bBo3iU+/t/Xu7CmZQxVpKIsY6FrwD8o37M43E/Cl5hIGfwU/ron+dI6rko7/DPtw42oiYQLBMi6hEKfEJZ82fvgTVCnDEfasZGhtznmM2e4zXj6+l3qi3Fl0pqTLmZ94+xIUny1/R7dgGP55GcoYsJGfpwgcHH9qX6sT6kYy3B298klNjtfb489P/qx0QLUAy8BFqLNakJK1WOPubKc92Za+7zUXtM+2B4lA9n3j/5LgIlBF4T/H+WpX9+eif+MX6zjL/evX+Mvcr6PuyHxjCOhIz37ii4CEh95M/XBfThomDjnw8fQoKf4wxJpI7xnLV2vtqbSw5J/nSadSqQI+S2fNya9z7nyvv16ucRIG/ULGtfKMqoFgEREIEjCfCcFTtf8dxW6mP3d/50mIGVOt5rHcDDV+svH18xlhxfqC9nTU9ZxtbTZktb5I0an5rvHFv6Xap7tv3zHBSfJ3KevvYub694zvL5J09LXLrrnb9NlkymCPAsxfoPn30/NY7Sdan6Jffxn9k90/GlfRMBERCBGAGfN3w/ti423v3O3zpyxhMzgO8tiLuGnr1+zLqeXwpx7rPD4K3juQ9rWbleArXPZ4f92v+MiZ8+80TH/vUejopcBERABESgFwFeANTIPYfe+fcCVtsO3hvXgF/brxr6ZnkvF7sbnSW+2F4fER++C72+D0fEF2PdYnzm+GaODcdCr/jwXeMN6J5jsNs7fzjb6wSxB0ivteLRi7TsiIAIiMBcBJj89+RU3fnPdUwoGhEQAREQARHYJNAl+esuN7wPvHoLz2pUBERABERABNoQCCZ/JCV+YJZtJivb51gb98bUWuPdDh7X8N3NSAxrxDbmrr16pfhG3p1t37R/24xGldDejbUzweQPF/27hFCfCSwWEi8S/NqY/DWOk+FIFwDXuA+KWQREQATOQqBGTg0m/5hijDNJMbFbWHYMbSY2f8WHvh3zfei087F+ioz1z9vx/Zgdr8P3a/iBvyAg25gf1k6p71s6bGy9/IjZsb74eH0/VceazlQdlmHJmlLfrd2YDi9TO96YXW/H971fth9jaGVo1495O77v5W0/Zre1jjW79A8125S3fY55X33frkHb9kt0lKyhXWvbtmM6rQx1+Pi2+iU6StbYGLyvvm9l6b+1yXk7FtLBtaV19Nf+TN5QvNbGvL9YsPKYV0knIHbprCQpAiIgAtdKYG+uCN75x5RyHDWKT/rXsgn2iqx2zPbpSm3dKfpKYsPxwGMixcaRMiXxHelvrm3Fl0tsLPmZ92/m2HAUnS2+YPJHIP6E7k/unL/WCwAwUvkgcPRFy4cnaomACIiACGwRiD7231q4No8LA10UrBHanjsbQ/rLejtCSYiACIiACJQS2Huujd75lzqkdXUInO1Omv6yrkNBWkRABERABFoQUPJvQfVKdTLxs8aVqYoIiIAIiMB4BJT8C/ak5w87eifQvbEx8Rdg7bJkb3xdnNxhRPHtgDfA0pn3b+bYcOicLT698x/gC7/mApM/kurZypl9Pxtr+SsCInBdBHB+3ZMXvl0XrvNFy809YyKl7+ejLo9FQAREYFwCzAd7PNRj/z30Oq5VIu0IW6ZEQAREYHICSv4FG3y2dzs5Ic4cGzgovpyjYTxZ7d94e5LqkfYulVQfOb3z78O5ipW973iqOCElIiACIiAChxLgY/89T4R153/oFuYZx0Zz0/NWSloEREAEREAEPggo+X+wOEWLFwC6CDjFdslJERABERiSgJJ/wbYc/e5qz6OerXCPjm3Lv73zim8vwWPXa/+O5b/HuvZuD736a5X86zPtopFPALoYkxEREAEREIGpCOgHfyfeTj76b/kk4MR45LoIiIAITEmgxrlfd/4nPjSU9E+8eXJdBERABA4koORfAH/md1czx4atVnwFB/xAS7R/A21Gpivau0xgEXHc9de48VPyjwA+y7De/Z9lp+SnCIiACIxDQO/8x9mLYk9qXQkWO6CFIiACIiAC3QjonX831DIkAiIgAiIgAscTqHWzp8f+x++lPBABERABERCBrgSU/Atwz/zDlZljw1YrvoIDfqAl2r+BNiPTFe1dJrDG4nrn3xhwD/W1HgP18FU2REAEREAEygnwfM+6VJPu/EvJDbROv/gfaDPkigiIgAg0IrA34Vu3lPwtDbVFQAREQARE4AoIKPkXbPLM765mjg1brfgKDviBlmj/BtqMTFe0d5nAGovrnX9jwD3V13wk1NNv2RIBERABEdgmYM/xtr298quE7vy/MjnlyN4D4ZRBy2kREAERuBICtc/x30LcYITF/5gs1KesahEQAREQAREQgTYEal4ABO/8keDxsSXU92NWfub2zO+uZo4Nx6TiO/c3U/t33v3T3u3bu9r5Npj84aK9+6fLMM5x1pxTLQIiIAIiIAIicA4C0eS/dpWBOV4I2IsAtG2fCPwVH/p2zPexzs7H+iky9IE67Bq0bZ8yfs1Wv0SHXYO27cf8sDKpa7zvWzqsfC8/YnasLz5e30/VsaYzVYdlWLKm1HdrN6bDy9SON2bX2/F975ftxxhaGdr1Y96O73t524/Zba1jzS79Q8025W2fY95X37dr0Lb9Eh0la2jX2rbtmE4rQx0+vq1+iY6SNTYG76vvW1n6b21y3o6FdHBtaR38tT8SOJN7qLbGKLs1ZufVrksgtAd1LUibCIiACIjAkQT8ed73c337FlsAxbbYvm1bGbVFQAREQAREQATGJxB87M/H+qhR1vqUGT/Ueh7axzH1tI6haebYQFjxjXGclXqh/Ssld/w67d3xe2A9CD72twIl7b2PI0psXvMa8b7m3VfsIiAC10DAn+d9P5dB9LF/rqLe8gh8rVzLE4ktDmuMNCcCIiACInCdBE6X/JnstpI75bCtW7Jn3/rZ4zv7/sh/ERABERiNQPCd/9FOInHHPkh0KcmOcqitrhqxjfLuyl7g1IgLOkaJrVY8Xo/i80TO1df+nWu/rLfaO0vj+PZw7/yR0FKSeym61vpL/SpZN1MsJfFrjQiIgAhcCwF/vvf9XA5DPPZHECwtEz9tpNR7wabY2CNjme3Ro7UiIAIiIALXR+Dw5N87yfI1wNpFBhNrb99SD79R/Ur1X3IiIAIiIALHEhjynX9rJLwAWLOzdnEw87urmWPDfiu+taN+/Dnt3/h7FPNQexcjc8x40zt/3kGvhbaWZNfW1Zhbu4Nem6thu1THqH6VxqN1IiACIiAC/Qk0+8EfQjkysaeg9IkUffrt51L09ZAZ1a8escuGCIiACFwjgdB5PzSWw6bZY//RE7+HRJDWb4yNVOjjSD7JFxEQAREQgfMRaJb8z4ACiR4J1SdVm/RtmzHN/O5qtNhC/LkPJfVo8ZXEsLZG8a3RGX9u5v2bOTYcWWeL76qTf+hUwGRjnwCE5DTWhwAv0PpYkxUREAERGIsAclKLfNTsnX8LZ3tsCZM/bDHxjBBLqwOgB9O9NuyecF/26tR6ERABETgDgdi5PzaeGlPTX/unOjGSXCjR74U8Unxn9YX74i8EzhqP/BYBERCBLQItc4+Sv6EfSyxHPwFoeQCY8IdtWv62PazDckwEREAEBiegd/4bG8Q7TiYdiJ/thx0bIX6aHjU2y/+Tw5mdUePLDCMqrviiaE4xMfP+zRwbDq6zxad3/hmnBD4Z4AVBxtJdotd+5094lr+YkIpqERCBGQlsneO25reY6M5/i5CZ7530YZoJz7hxtU3PH2zE52oPBwUuAiKwg4CS/w54vZb6pNfL7qh2bMIHG9sf1Wf5JQIiIAIjEVDyL9iNs73byQlx9Nh4IYSa7Zniy4klJDv6/oV8zhlTfDm0xpLV3qXtR68nmnrnn7Yfn6SwOSWJ55OSxE5PW4kuDSFmudi2dQ7jKL32ytpWWwREQARKCMTOZ15Xqpxfx77+1I8kMmokk1zwKfKQ8UWJyxP52ud+2BmMhcatjNoiIAIicK0ElPx37DyTNZMM6lChHOZiFwGx8ZA+jX0lEGJvuX9doREREAERGI/AVj6p5bHe+ReQxLsrbBATDhM3ap9wOAczfo5jVqbAnapLzvBeLocX94iQzhAffS2pFV8JtXHWzLx/M8eGI+hs8enOf+f33iYXtm2S5xhrzDF5sd7pgpaLgAiIgAiIQBYB/eAvC1cd4dDFQR3N82shO15MhSJOkQmt05gIiIAIjEAg5cYwRWYtltU7fyrnyRSKcNK1fY6tGdHcZwJrieuzpHqWAI9HO+bbKTJ+jfoiIAIicG0Eou/8txI8Ehg/1wbtbO92cvZn5NhCF542Nn/M2jm2R46PPu6pFd8eesevnXn/Zo4NR87Z4gsm/9jdkz35QsafbENjISiAZEH5fmwNxllqrCnVQR9Ql+rAOpZUHSVraAO1t+P7VpZtazNVR2gN9YV0cMzL+D708hgM+f7w8GCX3LT0wxrKtRPy/QgdiOEIu6nxW99S1/h9GVHH6NzXGJ7Z99Rj6OhjBnkU57ItP+w+lbSD7/x9UodinHRRMGfbGGMfbRQr8zqif0WgPgF/nPl+fYvSKAIiIAJtCaSex1LlYt4G7/yRzJnQWUMBjaFGsXOXAf0TJUBmUQFN7CLAY3OXEi0WAREQgSshsPqDP8vAJy/2dQFgKa2315itza1rnXOWPBgdjjOM8XizbcqoFgEREAERSCMQvPPnUp5oUfODObY5T/lrqe27mJyYyQ2Jy3/sXI7O2rKlsdX2A/rIBG3wipWcC4GR4ovFs2dc8e2hd/zamfdv5thw5JwtvuA7/71fgZyT8V5bZ1mfyiRV7ixxl/ppOYTaGGO51otQxq9aBERgHgI837GORbY1H1vH8eTH/lygWgRaErBJnXZCyT00Rnl55dMjAAAEuUlEQVTVIiACInBGAjz/MbGzbhHL6mP/FgavVSeSFTc2xmBrPrZutnGwUnKfbVcVjwiIwEgElPwLdqPlu52jk17L2ApQf1pS4yp45Pg+BVvYUXyF4AZZNvP+zRwbDp9a8fXKAXrn3/FLv5W8tuY7unqIqbX41+YOcVZGRUAERKAyAZznUHgBsHbeW5tLcUt3/imUKslgQ7m5ViXGQuNWZvb22oG8Njc7F8UnAiJwHQR6n+eU/DsfV/4CgBvOK73O7hxuDvGTQciZtbmQvMZEQARE4KwEcL5jaX3uU/In6Yx677sdXgDYjYb5ES4A9saWgfE96feMu2d8OSxqySq+WiSP0TPz/s0cG46WvfHhPNjzXKg/9TvmO951kw8KMWqWFz1bBzrlWl8Be0dpF+NbPvq16ouACIhAKQF77inVkbpOP/hLJSW53QTsgb2WVCm3JrPbmQ0F8OFI+xvuaVoERGAyAjznpJ7/KF+KQXf+peS0LotA6oGaKpdlvEBYib8AmpaIgAhUIdDj/KN3/gVbtffdToHJbktaxMYr2W5BrBhqEd+Kue5Tiq878qoGZ96/mWPDQVAzvh7nzGaP/at+I6RMBERABERABETgE4E9TwiaPfbf49Sn6AbsjPJougWamWMDL8XX4qjpp1P71491bUvau9pE9+nTY/99/LRaBERABERABE5HoMlj/9NRkMMiIAIiIAIicEUEdt354zEOC9r8YIxtyvg+141a01/rH2PBmJ9n38rYtaO1vZ/0n37aPtusKXOGGj6z0P/QGGXOUFv/4a+Pi33OhWqMjVp8fPDTjs0UH2NhfKG+jZ9yo+6d9ZVtxhTrc9zWaI9cGBP3w/btGGJgf7R4ipO/DQhtvOP37/lDfT82GpCYP7F47XiIQUzfkePWZ/jh98/2KcvYzrR/9N2ytv6H4rSyZ2gzRsZ19pgYj2Vvx2aMj3vHmNm3cXNu9DrkM+Oh777PNdxbyo1a018fB/uoKcN6xFiKkz8DRVBsM1DUHLNBc96OjdoO+b/l61niC8UG3/kJxbk2F5IfYSwW5wi+1fYhduyBAedY17ZdU19oz9b0nz2+WLx2nPuG2o6vcTlqbnT/anBhjNwP1jV099RRnPxjTlowXoZzfnz0Pr988NO2vd9njQ9xwPc1/7fmPYtR+4xxbR9H9T3m1yx7E4rP7pNte1nuK2s/P3IfccX85jjrkeO4Ft+4X9iTtWNydB7Vkj+BIODZDlTEw5hYn3nTRz8oe/t35r3E8ej99337neTx25txqT34S59Zl+oacR32KhRXbHzEGPb4xNhZ79F11Fr4HvrOHeVPqt0qf+fPwFlzI9mnM+xznuNnq7nZiMfGctb4GA/2gfHYWNBm38qcbd/gL+MIxXnGeOgz4rIx2T7HWXPNGWvE4PeQcbE+W1w+HvYZB+NizfGz1D4e32dcrEePC376GOCzHTtDLPpTv9GPNPknAiIgAiIgApUJVHvsX9kvqRMBERABERABEWhEQMm/EVipFQEREAEREIFRCSj5j7oz8ksEREAEREAEGhFQ8m8EVmpFQAREQAREYFQCSv6j7oz8EgEREAEREIFGBJT8G4GVWhEQAREQAREYlcD/BxW0CpDSEE4iAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "968d8ab5-08ab-4a71-ad3e-3b408641d4c7",
   "metadata": {},
   "source": [
    "Now we can subset this data to extract only information in the NINO3.4 region: 5S-5N, 190-240E (the middle of the Pacific):\n",
    "![image.png](attachment:022a5d1d-b083-4383-8c15-caea17b25821.png)\n",
    "\n",
    "This regional filter will be constructed using _logical masking_: assigning \"True\" to values within a specified lat/lon range, and \"False\" to values outside that range. This approach to filtering is nice because it works regardless of the shape or regularity of your data's lat/lon spacing!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "93e5fadf-2d2e-4d03-a5ba-634011b816f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create logical masks for lat and lon variables\n",
    "\n",
    "\n",
    "# Apply lat/lon masks to the ERSST field, then calculate averages over the lat and lon dimensions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d346c50f-42b8-41af-9d4c-246c9ca5e3d0",
   "metadata": {},
   "source": [
    "Let's make a plot to see what this time series looks like!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fc311470-ea68-4dc3-99e0-9151109832f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new figure, set size and resolution (dpi)\n",
    "\n",
    "# add data to the plot\n",
    "  \n",
    "# Specify plot parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e92e21a-d211-4f65-934e-4e6c1aab1991",
   "metadata": {},
   "source": [
    "There is clearly a lot of variability in the time series above, with some years warmer than others: but how do we decide what is really warmer than what would be expected? To do this, we need to remove the effects of seasonal changes in temperature. This is often done by computing the NINO3.4 \"SST anomaly\", where the mean value for each calendar month is removed (so for January, the average of all Januaries, for February the average of all Februaries, and so on). \n",
    "\n",
    "In xarray, the \"groupby\" functionality allows us to do this quite easily! We can get some basic information on groupby using the `?` syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2c715996-0342-4b08-89eb-696b63343380",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take a look at the groupby method in n34_ts\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fda48188-54a1-4353-a245-6022178990e1",
   "metadata": {},
   "source": [
    "We need to pass an argument to groupby, to tell it what to use to construct the groups. There are two possible syntaxes for this: the first uses a feature within xarray called the `DateTimeAccessor`, or .dt for short. This is applied to the time variable within the xarray dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0f0914a9-371a-4ec5-841a-b61b1eab67ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display n34_ts DateTimeAccessor\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd5e394d-028d-4116-ba1a-b9a225b7562c",
   "metadata": {},
   "source": [
    "The DateTimeAccessor knows various things about the time that you can query: for instance, the calendar month:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "63ff39fc-d377-42e1-a66a-1e5e8e31cbde",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Query calendar month\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dad7749b-76bf-4f25-878f-b9a3561ce2f8",
   "metadata": {},
   "source": [
    "Or when variables are already included in the dataframe, you can use a simpler syntax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "94028a79-96dd-4afe-9322-2f7f0b52c49c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Query calendar month more simply\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3a02e9-0825-47b2-bbcd-a1bbbcdef569",
   "metadata": {},
   "source": [
    "The above returned what's called a `GroupBy` object, which is essentially a list of groups and which parts of the dataset belong to which group. we can then calculate the monthly mean using this set of groups. Here I will also use the `sel` functionality to select a time period over which to calculate the climatology: 1971-2000 is the \"base period\" used in the NOAA operational definition of El Nino (more on that below).\n",
    "\n",
    "Stacking up all our commands, we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "43d65017-f0b3-49d2-9744-b6c6a01f3e01",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the 1971-2000 base period, group it by month, and compute a time average for each month\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15b2a095-e008-4d3d-9819-bb8e744ce198",
   "metadata": {},
   "source": [
    "In order to complete our calculation, we next have to _apply_ the GroupBy object to the raw NINO3.4 time series.\n",
    "\n",
    "The most efficient way (that I know of) to do this is to use the Python \"function stacking\" capabilities to simultaneously create the groupings and apply them, in a single line of code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "38556379-b744-44e4-8e42-6aeb511b5331",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the base period climatology and remove it from the NINO3.4 SST time series\n",
    "\n",
    "# Display SSTA dataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e280311c-a3fa-4b89-929a-9dade3df2ee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create a new figure, set size and resolution (dpi)\n",
    " \n",
    "# add data to the plot\n",
    "   \n",
    "# Set plot parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9ed0d67-0a73-451f-86c0-5a46e63ebedb",
   "metadata": {},
   "source": [
    "Now that the seasonal cycle has been removed, it's much more obvious when there are large excursions from what is climatologically \"normal\" for that time of year. \n",
    "\n",
    "There are many definitions for El Nino and La Nina events. Here, let's use the one that NOAA uses in its operational ENSO forecasts:\n",
    "https://www.nws.noaa.gov/ost/climate/STIP/ElNinoDef.htm\n",
    "\n",
    "This definition states that El Nino is \"a positive sea surface temperature departure from normal (for the 1971-2000 base period) in the Niño 3.4 region greater than or equal in magnitude to 0.5 degrees C (0.9 degrees Fahrenheit), averaged over three consecutive months.\" with La Nina being the same thing, but using a threshold of -0.5C. \n",
    "\n",
    "Since we want signals that persist over three consecutive months, we'll use the `rolling` method to isolate periods where the running 3-month average is above/below the threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2ea6289c-ff4f-4b7e-a8ad-e6a8179ab147",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate 3-month rolling mean NINO3.4 SSTA\n",
    "\n",
    "\n",
    "# create a new figure, set size and resolution (dpi)\n",
    "    \n",
    "# add data to the plot\n",
    "   \n",
    "# Set plot parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c35890d5-5357-4f31-ab04-1c98aeba3378",
   "metadata": {},
   "source": [
    "Another powerful data analysis technique is the use of thresholds for _compositing_. In other words: select all parts of a dataset which satisfy some criterion (in this case, exceeding the El Nino/La Nina SSTA threshold) and making an average over those periods.\n",
    "\n",
    "The below code creates logical masks for the NINO 3.4 SSTA time series which are either:\n",
    "- True when El Nino, False other times; or\n",
    "- True when La Nina, False other times:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5ba7d48a-9880-4e9d-b555-4ebab3f0051a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask for El Nino\n",
    "\n",
    "# Mask for La Nina\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b110a342-571f-43ab-85dc-084ed1267d14",
   "metadata": {},
   "source": [
    "Using these logical masks, and the knowledge that the full set of ERSST data has the same size time dimension as the NINO3.4 SSTA time series, we can then use this to construct a composite of the SST anomaly _patterns_ present during El Nino and La Nina events. \n",
    "\n",
    "First, we need to compute a new array, containing spatial fields of SST anomalies. We can do this using the same group by approach we applied to the NINO3.4 time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e0200431-3945-44e3-bb12-a304151d4f04",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate base period climatology and remove it from the ERSST field\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07409899-d809-4f10-a403-6be921ea2323",
   "metadata": {},
   "source": [
    "Now we apply the logical masks to this new SST anomaly field, and average over time to create a composite map:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "a47f6ab1-61ee-4660-8b3b-f3d36a8200a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply El Nino mask and extract the xarray DataArray\n",
    "\n",
    "\n",
    "# Apply La Nina mask and extract the xarray DataArray\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebf81285-f779-442a-b9f1-6c0f8cda942e",
   "metadata": {},
   "source": [
    "The final step is to plot this as a nicely formatted map! This can be accomplished using the functionality of cartopy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a621dcce-9dee-4723-8984-7f519619d236",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an object containing the Plate-Carree map projection, centered at 180E\n",
    "\n",
    "# Create a figure object\n",
    "\n",
    "# Create an axis object using the specified projection\n",
    "\n",
    "# Make a filled contour plot of the El Nino SSTA composite using the specified map projection\n",
    "\n",
    "# Set plot parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "8f2d359f-91aa-4286-b4a4-c89821dd8273",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an object containing the Plate-Carree map projection, centered at 180E\n",
    "\n",
    "# Create a figure object\n",
    "\n",
    "# Create an axis object using the specified projection\n",
    "\n",
    "# Make a filled contour plot of the La Nina SSTA composite using the specified map projection\n",
    "\n",
    "# Set plot parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48a46269-c5b5-429b-b79a-3b83b24ac0d3",
   "metadata": {},
   "source": [
    "### **PART II: Marine Heat Waves**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5a20477-ef55-4fdb-8048-b3143b687ea9",
   "metadata": {},
   "source": [
    "The next example of SST-based data analysis is marine heat waves: these are events where SST is above a (typically very high) local threshold. \n",
    "\n",
    "Here we'll use the OISSTv2 data, since we'll need daily resolution. The OISSTv2 data are stored in many files, each containing information for a single year. The base directory where the information is located can be found here:\n",
    "\n",
    "https://psl.noaa.gov/thredds/catalog/Datasets/noaa.oisst.v2.highres/catalog.html\n",
    "\n",
    "We are also only interested in data within the tropical Pacific - that means that we don't have to read in the entire file, only the subset of it within our region of interest.\n",
    "\n",
    "All of these tasks can be accomplished by combining _spatial subsetting_ with a carefully constructed _for loop_.\n",
    "\n",
    "First, let's specify some parameters appropriate for accessing the OISST files. Since we know that there will be one file per year, we can store the \"base URL\" (first part of the filename) and create a sequence of years of interest using the `range()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d3718979-1b82-4f2a-9fa4-ecfed8059cdf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "range(2013, 2019)\n"
     ]
    }
   ],
   "source": [
    "# NOAA THREDDS server: base URL for all files \n",
    "base_url=\"https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/\"\n",
    "\n",
    "# Naming convention for files: \"sst.day.mean.[year].v2.nc\"\n",
    "yrs=range(2013,2019)\n",
    "print(yrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fbde6f3-baf4-42ca-8285-cb7e20b85950",
   "metadata": {},
   "source": [
    "This combination of information now allows us to open any OISST file. To make things easier, let's work with the anomaly files; these are files which have already had the seasonal climatology subtracted out, so we don't have to do that step again.\n",
    "\n",
    "As an example, the following will open the file for 1981:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d32d4f23-cec0-4d60-9ca9-8dceb7230b54",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "410e1565-e79d-4178-a46a-45f8da0e3878",
   "metadata": {},
   "source": [
    "Now let's extend this to read in data over the entire range of years. \n",
    "\n",
    "The below code chunk will take a few minutes to run, since the data files we're working with are quite large. Run time was roughly 2-3 minutes during initial testing of this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "18dd9c05-7701-4b71-84ee-62fa9fb207cb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2013\n",
      "2014\n",
      "2015\n",
      "2016\n",
      "2017\n",
      "2018\n"
     ]
    }
   ],
   "source": [
    "nepac_sst=[]     # Blank variable to store output data \n",
    "\n",
    "for y in yrs:\n",
    "    print(y)\n",
    "    \n",
    "    # Open the dataset into an xarray object\n",
    "    \n",
    "    # Define logical mask to be true in the Northeast Pacific (30-50N, 180-240E)\n",
    "\n",
    "    # Read in data \n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2573381c-c9cc-4586-a0b6-e77a51564c9a",
   "metadata": {},
   "source": [
    "Since we have appended each year's data to the initially blank output variable, the data returned is a list of xarray Datasets. We can convert this into a single Dataset containing time information for all years using `xr.concat`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "174765dc-ef8f-40a6-9ff3-0e2c61e7a769",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Concatenate list of xarray Datasets into a single Dataset, along dimension \"time\"\n",
    "\n",
    "\n",
    "# Calculate a regionally averaged time series for the NE Pacific\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4145ee72-f157-4435-8ea9-f2d4469b52d2",
   "metadata": {},
   "source": [
    "Now let's try a new analysis: looking at the exceedance of the _daily_ SST anomaly over a given threshold. This is a standard approach used in marine heat wave studies, although the thresholds and persistence times used are not yet as \"standardized\" as those used in ENSO research.\n",
    "\n",
    "There is a paper by Hobday et al. (2016) which proposes numerous marine heat wave definitions - let's use this one:\n",
    "- Daily SST anomaly exceeding the 90th percentile (for 5 or more days)\n",
    "\n",
    "If we were doing this for real, the 90th percentile value would be calculated using a long baseline (roughly 30 years or more), but we'll \"pretend\" for now using the 5-year time series above to shorten the data I/O time needed. To make the computation simpler, we'll also ignore the 5-day persistence requirement for now.\n",
    "\n",
    "The method of calculating percentiles we will use is a numpy function called `percentile`, since xarray itself does not contain a native percentile function. In order to apply this percentile calculation to the xarray dataset, we can use an xarray method called `reduce`: this reduces the size of a dataset by applying a function you pass it as an argument.\n",
    "\n",
    "The below code chunk will return the 90th percentile of the NE Pacific time series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e869857e-daf0-4c59-b924-efc73d275bfb",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "6fd9f521-8b4b-4628-b086-a7e7d1bc23bf",
   "metadata": {},
   "source": [
    "Now we can use the same logical masking method to identify all places within the NE Pacific time series which exceed that threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "0d804e09-e5a1-497d-887e-9179bdc46da3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask for MHW\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c85b6d70-8ab6-4050-9ccf-9ddfe342a2b3",
   "metadata": {},
   "source": [
    "As for the ENSO example, this mask can now be used to composite over all marine heat wave events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "8acfb297-e21a-4dc4-ba33-4678296d7b9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply MHW mask and extract the xarray DataArray\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e79f5123-f8cf-4678-88e3-d90f49475078",
   "metadata": {},
   "source": [
    "Finally, we make another cartopy plot to display the results!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "0737c9a0-1cc8-485d-8690-49da38141dd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an object containing the Plate-Carree map projection, centered at 180E\n",
    "\n",
    "# Create a figure object\n",
    "\n",
    "# Create an axis object using the specified projection\n",
    "\n",
    "# Make a filled contour plot of the La Nina SSTA composite using the specified map projection\n",
    "\n",
    "# Set plot parameters\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48d017ef-26de-4213-b0a6-56150e9a6033",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}