DaymetV4_NGAP_Hyrax_OPeNDAP.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b482542e",
   "metadata": {},
   "source": [
    "# NASA Cloud-Based OPeNDAP Access and Subsetting - Daymet Data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50489f29",
   "metadata": {},
   "source": [
    "*Author:  ORNL DAAC*\n",
    "<br>\n",
    "*Date: March 13, 2023*\n",
    "<br>\n",
    "*Contact for [ORNL DAAC](https://daac.ornl.gov/):  uso@daac.ornl.gov*\n",
    "\n",
    "##### Keywords:  ORNL DAAC, Daymet, THREDDS, OPeNDAP, Earthdata, , `CMR API`, `pydap`, `xarray`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99c311ef",
   "metadata": {},
   "source": [
    "## Tutorial Overview\n",
    "\n",
    "In this tutorial, we're interested in obtaining a **spatial** and **temporal subset** of gridded multidimentional (ie. netCDF) data using OPeNDAP.  We'll demonstrate how to programatically discover, access, and subset using Daymet daily data hosted from NASA's Earthdata Cloud as an example dataset. \n",
    "\n",
    "- **STEP 1: Setting Search and Subset Parameters:** Demonstrates how to use a shape file's bounding box coordinates to obtain `search` and `subsetting` parameters. \n",
    "\n",
    "- **STEP 2: Finding Data Files of Interest:**  Demonstrates how use spatial and temporal parameters to search NASA's Dataset Collections to discover Granules (files) of interest using `NASA's CMR API`. \n",
    "\n",
    "- **STEP 3: Subsetting Data:** Demonstrates how to subset from the North American Daymet dataset netCDF files using `pydap` to inspect data from the NASA Earthdata Cloud Hyrax `OPeNDAP` server.\n",
    "\n",
    "Steps we'll explore: \n",
    "- Store a shape file's (polygon) bounding box parameters to use for searching and subsetting from the North American extent of Daymet Data\n",
    "- Use the Daymet dataset **DOI** or the Daymet **shortname** and NASA's **Common Metadata Repository (CMR) API** to search for Daymet data within a spatial and temporal region of interest\n",
    "- Build a **Request URL** to create a list of files that match our search criteria\n",
    "- From Daymet data on the NASA Earthdata Cloud-based Hyrax **OPeNDAP Server** we will find, subset, and concatenate data as if it were stored locally\n",
    "- **Download subsets** of data NASA Earthdata Cloud to a local computer\n",
    "- **Clip** the subset to the polygon of interest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a673e36",
   "metadata": {},
   "source": [
    "#### Daymet Data: NASA Earthdata\n",
    " - [NASA Earthdata Daymet Collections](https://search.earthdata.nasa.gov/portal/ornldaac/search?fpj=Daymet!Daymet!Daymet)\n",
    "\n",
    "#### Daymet Data: ORNL DAAC \n",
    "- The [Daymet Website](https://daymet.ornl.gov) provides comprehensive information about and access methods for the Daymet datasets under the [Get Data](https://daymet.ornl.gov/getdata) tab. \n",
    "\n",
    "- Landing Page for [Daymet: Daily Surface Weather Data on a 1-km Grid for North America, Version 4 R1](https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=2129).  \n",
    "- Daymet V4 Daily Data dataset DOI **https://doi.org/10.3334/ORNLDAAC/2129**\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "Thornton, M.M., R. Shrestha, Y. Wei, P.E. Thornton, S-C. Kao, and B.E. Wilson. 2022. Daymet: Daily Surface Weather Data on a 1-km Grid for North America, Version 4 R1. ORNL DAAC, Oak Ridge, Tennessee, USA. https://doi.org/10.3334/ORNLDAAC/2129\n",
    "</div>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cbc9a790",
   "metadata": {},
   "source": [
    "### Daymet Overview\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "SKIP this section if you're familiar with Daymet:  The **[Daymet](https://daymet.ornl.gov/)** dataset provides gridded estimates of daily weather data interpolated and extrapoloated from ground-based surface weather observations.  Annual and monthly climatologies are also available.  The spatial resolution of the data is 1km x 1km and the daily data are distributed in a multidimensional, netCDF file format. Daymet data is available for continenal North America (na), with Hawaii (hi), and Puerto Rico (pr) as separate data files.\n",
    "\n",
    "The Daymet daily North American files [Daymet: Daily Surface Weather Data on a 1-km Grid for North America, Version 4 R1](https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=2129) are large and often a researcher needs a spatial or temporal subset of the data based on their region of interest. \n",
    "\n",
    "\n",
    "![title](images/Daymet_V4_Annual_Climatology_Fig1.png)\n",
    "\n",
    "DAYMET Temporal Availability:\n",
    "\n",
    "| Spatial Area    | Years                          \n",
    "|:---------       |:---------------------------------------:|\n",
    "| North America   | 1980 - present\n",
    "| Hawaii          | 1980 - present\n",
    "| Puerto Rico     | 1950 - present\n",
    "\n",
    "DAYMET VARIABLES:\n",
    "\n",
    "| Variable | Description (units)                           \n",
    "|:---------|:---------------------------------------:|\n",
    "| tmax     | Daily maximum 2-meter air temperature (°C)\n",
    "| tmin     | Daily minimum 2-meter air temperature (°C)\n",
    "| prcp     | Daily total precipitation (mm/day)\n",
    "| srad     | Incident shortwave radiation flux density (W/m2)\n",
    "| vp       | Water vapor pressure (Pa)\n",
    "| swe      | Snow water equivalent (kg/m2)\n",
    "| dayl     | Duration of the daylight period (seconds/day)\n",
    "\n",
    "Daymet V4 Daily Data files (or granules) are distributed in netCDF4 format, and each file has one year's worth of data. Data files are organized by variables (dayl, prcp, tmin, tmax, srad, swe, vp) and regions (each for na, pr, hi).\n",
    "\n",
    "The standard naming convention for a file is:  \n",
    "\n",
    "- daymet_v4_daily_`area`_`variable`_`year`.nc\n",
    "- daymet_v4_daily_na_tmax_2022.nc\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ca1257a",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28640b33",
   "metadata": {},
   "source": [
    "## Tutorial Requirements\n",
    "* Python 3.7 or higher\n",
    "* geopandas and dependencies\n",
    "* pydap\n",
    "* xarray\n",
    "* rasterio\n",
    "* shapely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dd5e478",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Running from pyenv, Conda, activate hyrax_daymet\n",
    "from IPython.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:80% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "24350611",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import requests\n",
    "import pyproj\n",
    "import datetime as dt \n",
    "from rasterio.features import Affine\n",
    "from pydap.client import open_url\n",
    "from pydap.cas.urs import setup_session\n",
    "import xarray as xr\n",
    "import rioxarray\n",
    "import shapely\n",
    "import rasterio\n",
    "from shapely.geometry import mapping\n",
    "import json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e6d496fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This Notebook was produced with the following versions\n",
      "geopanda version   :  0.10.2\n",
      "pyproj version     :  3.3.0\n",
      "rioxarray version  :  0.10.3\n",
      "xarray version     :  2022.3.0\n",
      "shapely            :  1.8.0\n",
      "rasterio           :  1.2.10\n"
     ]
    }
   ],
   "source": [
    "print(\"This Notebook was produced with the following versions\")\n",
    "print(\"geopanda version   : \", gpd.__version__)\n",
    "print(\"pyproj version     : \", pyproj.__version__)\n",
    "print(\"rioxarray version  : \", rioxarray.__version__)\n",
    "print(\"xarray version     : \", xr.__version__)\n",
    "print(\"shapely            : \", shapely.__version__)\n",
    "print(\"rasterio           : \", rasterio.__version__)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "dea6a573",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python 3.9.12\n"
     ]
    }
   ],
   "source": [
    "!python --version"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "460917a8",
   "metadata": {},
   "source": [
    "***\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95a1706a",
   "metadata": {},
   "source": [
    "## STEP 1.  Setting Search and Subset Parameters: Spatial Area and Time Range of Interest"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0fe0b9e",
   "metadata": {},
   "source": [
    "## 1.1 Spatial: Search and Subset Parameters\n",
    "You can always hard code subset parameters and if you do, you can skip to STEP 2 of this tutorial. \n",
    "\n",
    "Here we'll demonstrate how to use a polygon boundary file of California to get bounding box coordinates and how to use that for metadata discovery and spatial subsetting. Any polygon area-of-interest can be used.  \n",
    "\n",
    "We'll use Python's [GeoPandas](https://geopandas.org/index.html) to explore the polygon file\n",
    "\n",
    "* **`Geopandas`** is an open-sourced library that enables the use and manipulation of geospatial data in Python. \n",
    "* **`Geopandas`** extends the common datatype used in pandas to allow for many unique geometric operations\n",
    "* We'll use the Geopandas **`GeoDataFrame.to_crs`** to reproject the CA polygon's bounding box coordinates. \n",
    "\n",
    "We'll store bounding box coordinates in 2 reference systems\n",
    "1. **Geographic** WGS84 - for searching NASA Earthdata metadata\n",
    "2. **Lambert Conformal Conic** - Daymet's CRS - for subsetting the Daymet files using OPeNDAP protocols"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8748c00",
   "metadata": {},
   "source": [
    "### Read and store the bounding box coordinates of the shape file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1cd02e0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOUAAAD4CAYAAAATgSFTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAee0lEQVR4nO3deXxU9b3/8ddnZrIQsgFJ2EISdhAMIGFXiyAWUXHD5boAWkvVtg9re2trtb9be7WLba/eXr0qAu610qpoXWq5KK3sJOwIyL6EJWwJCSHrfH9/ZMCwmclkZr5nZj7PxyOPzJw5Z84nw3lzzpzzPd+vGGNQSjmHy3YBSqnTaSiVchgNpVIOo6FUymE0lEo5jCecK8vIyDB5eXnhXKVS1hQVFR0yxmQ2d7mwhjIvL4/CwsJwrlIpa0RkZyDL6eGrUg6joVTKYTSUSjmMhlIph9FQKuUwfodSRNwislJEPvA9/52IbBSRNSLyroikh6xKpWJIc/aUDwAbGj2fC/Q3xuQDXwIPB7MwpWKVX6EUkWzgKmDGyWnGmH8YY+p8T5cA2cEvT6nY42/jgaeBh4CU87x+N/DWuV4QkWnANICcnBy/C/tsUwnvrijG4xLcLsHjbvjtFgGgpt5QV++lzmuorfdSV2+o83qpbfSbk7eKymm/kFPP5fTnZ0x3uQSXgFsEEd9jl+ASOe9rIvJVzS7B7W6ouWGa69TfIXxV3slbWo1vSuNbXE/e73pyWv/sNC7rneX356giT5OhFJGrgRJjTJGIjD7H648AdcAb51reGDMdmA5QUFDg9x3V2w8e5/3Ve/2dPWbcOTxXQxnl/NlTjgImisgEIBFIFZHXjTF3iMgU4GpgrAlyFwZe7RHhnE7uzVX0avI7pTHmYWNMtjEmD7gV+NQXyPHAT4CJxpjKYBd2rKqu6ZmUikItuU75DA3fMeeKyCoReT5INWGM4f1VxcF6u6iiO8ro16y7RIwx84H5vsc9QlAPACLCf1zTj0fnrKO49ESoVhORRI9fo55jW/Rc1ieLJyfl2y5DqbBzbCgBBue2Id7t6BKVCjpHb/GJcW4G5qTbLsNR9Og1+jk6lAAju7ezXYKjiJ7qiXqOD+WIbhrKxnRPGf0cH8qBOem4dEM8RT+K6Of4UAJ4tXGPiiGOD2VpZa3tEhxFD1+jn+NDebC82nYJjqKNB6Kf40M5f1OJ7RIcRSMZ/Rwfyg/W7LNdglJh5ehQbimpYOP+cttlOIvuKqOeo0O5eOsh2yUoFXaODmXXjGTbJSgVdo4O5YXZabZLUCrsHB3KtFZxdExLtF2GUmHl6FBu2HeMfWVVtstQKqwcHcrXlgQ0vJ9SEc2xoTxWVcucldpPj4o9jg3luyuKqaypt12GUmHnyFAaY/TQVcUsR4Zy8bbDbCmpsF2GUlY4MpSv615SxTDHhfLAsSo+WX/AdhlKWeO4UP5p6S7qtasBFcMcFcraei9vLttluwylrHJUKOd+cYAS7WlAxThHhTLB46hylLLCUSkY3TuL7DatbJehlFWOCqXbJQzXzpdVjHNUKAEqdLBYFeMcF8pf3XAhl/bKtF2GUtY4LpRtW8fz8tQh/GhcLx2uQMUkv0MpIm4RWSkiH/ietxWRuSKy2fe7TdCKcgnfH9uT1781jIzk+GC9rVIRoTl7ygeADY2e/xSYZ4zpCczzPQ+qkT0yeHFyge4xVUzxK5Qikg1cBcxoNPla4BXf41eA64Jamc+gnDbcN7p7KN5aKUfyd0/5NPAQ4G00rb0xZh+A73fWuRYUkWkiUigihQcPHgyoyAfG9uIiHdFZxYgmQykiVwMlxpiiQFZgjJlujCkwxhRkZgZ2VjXe4+J/bx+s3y9VTPBnTzkKmCgiO4A/A2NE5HXggIh0BPD9DulIPB3SEnnmtotw6xdMFeWaDKUx5mFjTLYxJg+4FfjUGHMH8D4wxTfbFOC9kFXpM7xbO350Ra9Qr0Ypq1pynfI3wDgR2QyM8z0PuduG5ujeUkU1T3NmNsbMB+b7Hh8Gxga/pK+XnhTPsK5tWbT1cLhX7Qx6/3fUc1yLHn+M6XPOE71KRYWIDGX/zjrwj4peERnKXu1TbJegVMhEZCjTW8URr70UqCgVkVv25pIKauq8Tc+oVASKyFAW7jxiuwRr9ORr9IvIUBbtOGq7BGuM0VhGu4gM5driMtslWKOZjH4RGcoBXdJtl2CNZjL6RWQor87vaLsEa7SBYfSLyFCO6pFBelKc7TKUComIDGWc28X4fh1sl2GFHr5Gv4gMJcDV+Z1sl2CFV8/0RL2IDeXwbm1p1zr2eiIQ/VYZ9SI2lB63i9RWsfe9UjSTUS9iQ1lceoLth47bLiPsPlyzj/V7Y/c6bSyI2FDuLT1huwQr9h+r4qbnF/PpRh2CPlpFbCiH5LXl51dfYLsMKypr6rnnlULeWLrTdikqBCI2lADfurgrP7i8p+0yrPAaeOTddfz27xvxevWMbDSJ6FACPDC2J/dc3NV2GdY8N38rP5y9Sm9liyIRH0oR4ZGr+nLNgNi8bgkwZ9VepsxaRtmJWtulqCCI+FBCQzB/NymfwblBG/gr4izedphJzy1iz9FK26WoFoqKUAIkxrmZfudgctom2S7Fms0lFVz37CJW7y61XYpqgagJJUC75ARmTR1CamKzurONKocqqrll+mL+vm6/7VJUgKIqlAA9spKZMWUIaTHY2uekqlov971RxIv/2qY9FUSgqAslwNCubXn/e6PoHcNdURoDT3y0gf/33nrq6vXMbCSJylAC5LZrzTv3j2RsjPem/tqSnXzntSIqa+psl6L8FLWhBGid4OHx6/vH/IBA8zaWcOv0JRyqqLZdivJDVIcSoGNaKy7vG9t7S4A1e8qY9Nwidh3WSyZOF/WhBJh2aXcu7RXYKNLRZMfhSm54biHrYrg3wEgQE6EcnNuG/739IuLdMfHnfq1DFTXc8sJiPt980HYp6jxiZitNTvAwskc7v+Ztn5pAmyjumOt4TT13vbScOSuLbZeiziFmQgkwsvv5Q3nyjv6heW2pqfOSn51+1jznOmE0rGvbYJUXVnVeww/eWsWL/9pmuxR1hiabvohIIvAvIME3/1+NMf8hIgOB54FEoA643xizLIS1ttjyrxnuIMHj4lsXd2XKiDxcLuHzzQdZv7eMOq9hwoUduSa/E92zWnPbi0vZUlIBwBUXtOcXE/sx/V/bOHK8hvdX7w3XnxI0T3y0gf3HqnhkQl9cMX6W2in8aY9WDYwxxlSISBywQEQ+Bn4JPGaM+VhEJgBPAqNDV2rLlFbWMH9TyXlfr6r1cs2ATmSlJgJw/aBsLuycRm671hyvrmPDvnKWbT/CIxP6sv9YFS6Bawd25tONJby8aAfdMlqH608JupkLtnOwvJrf3ZRPgsdtu5yY12QoTUM7rQrf0zjfj/H9pPqmpwGO3k3M/eIAl/TM5PpBnZmzsph5G08PaOt496kTQaWVNcxauIO1e0rZuL+cfWVVZ71fm6Q43ly2m2JftyTbIry/oPdX7+Xw8Wqev2MwKYnR+306Eog/bSNFxA0UAT2AZ40xPxGRvsAnNPSk7wJGGmPO6p9CRKYB0wBycnIG79xprwsLYwwigtdreHD2Kt5bdfr/I5MGZ/PE9f0Z8etPOXK8xlKVdl3QMZWX7x5CVkqi7VIinogUGWMKmrucXyd6jDH1xpiBQDYwVET6A/cBDxpjugAPAjPPs+x0Y0yBMaYgM9PutULxnc1xuYTfTRrA2D5ZpzVcnzoyj3i3i7Yx2J/sSV/sO8aNzy1i28GKpmdWIdGss6/GmFJgPjAemAK843vpL8DQYBYWavEeFzOmFDDvR99gjK99bNHOo4gIcx+8lPapCZYrtGf3kRNMen4xa/doIwMbmgyliGSKSLrvcSvgcmAjDd8hv+GbbQywOUQ1hoyIkJGcwMwpBTx1ywDys9MoO1GLiMR8nzdHjtfwby8uYem2w7ZLiTn+7Ck7Ap+JyBpgOTDXGPMB8G3gDyKyGvgVvu+NkUhEmHBhR3781zWM+69/cqKmntG9s2K+N/KK6jomz1qmfcyGWZOhNMasMcYMMsbkG2P6G2N+6Zu+wBgz2BgzwBgzzBhTFPpyQ8frhS0lFZSUVzNv4wGeumUgcx+8lHhPTLWvOEt1nZdprxbx3ipt/RMusb3FNVLn/epwtXtmMgA9slL0DhO+av3z2uIdtkuJCRpKn6rar0JZXvXVDcHXDuxsoxzHMQZ+/t56nvl0s3YxEmIaSp+q2vpTjx+ds5ZN+8sBGN07M6b7+znT7//xJb/6aIMGM4Q0lD6Nt7EvD1Qw8ZkFvL5kJ/FuFxMu7GivMAd68fPt/PTttdTrcAkhoaH0yWmXxNO3DDz1vLrOy6Nz1nH/Gyu4fVgON16UTWKcflwnvVW4m+/9aQXVdfVNz6yaRbeyRq4b1Jle7ZNPm/bxuv3c9fJy1hWXEefSj6uxj9ft555XCjlerZ1yBZNuZWcY3//sQ9WD5dVsOlBOuW58Z/l88yFum7GUozHaVjgUNJRnGNQl3XYJEWf17lJuemEx+8picyDfYNNQniHSb8GyZUtJBZOeW8xWbcjeYhrKM+w+ol0wBqq49AQ3a0P2FtNQnuH+0d3pmKb3EgbqsK8h+6Kth2yXErE0lGfISk1kxpQCkuK1W4xAVVTXMXXWch35K0AaynPo1ymN/7p5oO0yIlpNvZf73yhi9vLdtkuJOBrK8xjfvwPTLu1mu4yI5jXw0NtreOGfW22XElE0lF/jx9/szTd0uIMW+/XHG/n1x9pe1l+xO+SxH+LcLl6+awhzVhXz1vLdtEtOoH1KIgu2HOTLA3rqvzle+Oc2So/X8qsbLoz5UdCaoqFsgohw/aBsrh+UfWra+r1lXPvMQuq0QXazvFW4m4qaOp66eWDM3zz+dfSTCUC/TmncP7q77TIi0odr9vGd1wpPu1VOnU5DGaDvj+0Z0z3etcRnmw4yZdYyKrQt8TlpKANUV284VKGNsAO1dPsRbp+xlNJK/QzPpKEM0Jo9pXqTbwut3l3KLS8soaT87GEhYpmGMkArd5faLiEqbDpQzs3PL2bPUW1zfJKGsplmL9/N31bvZcbn222XEjV2HK7k5ucX61AJPnpJpBmMMTzx0QbKTtTaLiXq7C2r4uYXFvPq3cO4oFNq0wtEMd1TNoOIkNcuyXYZUetQRQ23Tl/Mil3nH9w3Fmgom2mA9kwQUseq6rhjxlIWbYndW780lM00IDvddglRr7KmnqkvL+f/vojNMUw0lM2ke8rwqKnzcu/rRby/2tEDhIeEhrKZumW0JiVBz4+FQ53X8MCfV/Lmsl22SwkrDWUzuVzCT67so3c6hIkx8PA7a5nx+TbbpYSNhjIAdwzP5dW7h+oYI2H0+IcbeGrulzFxT6aGMkCjemTw/vdG0S2zte1SYsZ/z9vM4x9G/83S/gyvnigiy0RktYisF5HHGr32fRHZ5Jv+ZGhLdZ7cdq159raL8OihbNjMXLCdh9+J7sGF/DljUQ2MMcZUiEgcsEBEPgZaAdcC+caYahGJydFV+3ZM5b7R3fmfT7fYLiVm/Hn5biqq63jqloHEuaPvYM+f4dWNMeZko8Q4348B7gN+Y4yp9s1XErIqHe57Y3rQXQ9jw+qDNfu497WiqLxZ2q//ZkTELSKrgBJgrjFmKdALuERElorIP0VkSAjrdLQEj5snJ+UjehQbVvM2lnDXS8uprImum6X9CqUxpt4YMxDIBoaKSH8aDn3bAMOBHwOzRc7eLEVkmogUikjhwYMHg1e5wwzObcvk4bm2y4g5i7cdZupLy6NqOL5mHZAbY0qB+cB4YA/wju/wdhngBTLOscx0Y0yBMaYgMzO6u2v88fg+dNIhD8Ju2fYj3BVFwfTn7GumiKT7HrcCLgc2AnOAMb7pvYB4IHZbEQPJCR5mTh2iwbRg2Y4jTH0pOvr98WdP2RH4TETWAMtp+E75ATAL6CYi64A/A1NMtF9A8kPfjqnM+e4o+nRIsV1KzFm+4yhTo6BDLglnjgoKCkxhYWHY1mfTe6uKeeDPq2yXEZMKctvw8t1DSbbcRllEiowxBc1dLvou8jjE2L7tSdAOh60o3HmUKbOWUV4VmT1E6FYTIskJHr59STdaxemQejYURXAwNZQh9O/f7E3Rzy/nu5dpb+o2rNhVGpHB1FCGWFK8hx+N682oHu1slxKTIjGYGsowcLmEP9w0kPQkvdXLhhW7Spk8axnHIiSYGsow6ZCWyG9u0KZ4tqzcVcrkmZERTA1lGI3v34GlPxvLr66/UMNpwardkRFMDWWYZaUkctuwHC7KaWO7lJi0ancpd85c5ugOtTWUlkRLO81ItHp3KZNnLnVsMDWUFhyuqGbj/nLbZcS01XvKHBtMDaUFbpcQ59Yvlbat3lPGnTOXUlbprGBqKC1IT4pn9ndGcEtBF5LitcWPTWv2lHGHw4KpobRkUE4bfjspn2WPXM5jE/vZLiemrS12VjA1lJYlJ3iYPCJX78G0bG1xGbfPXOKI4d41lA5RUx/zt6Jat674GLfPWGo9mBpKByguPcGhimrbZShg/d5j/OzdtVZr0FA6QNHO2B4k1Wl2HKq0un4NpQMcOW7/e4z6StvW8VbXr6F0gKkj8xiSp83unKJw5xEWWhxJWkPpACJCqUNOxyuoqvXyk7fXWBtISEPpAMYY9pdV2S5DNbLn6An+UrjHSjA1lA4gIozr1952GeoMD729htmFu8O+Xg2lQ5Qc00siTvSnZRrKmFRb79XLIg61endp2BsTaCgdoLKmniFd29ouQ51HuMce1VA6QFqrOF69eyjvfXeU9Wtk6mx/X7c/rOuz26+7Os2ALul8f0wPHvvbF7RrHU+n9FZ0Sk+kU3orEuPcHCyvZu2eMjYd0Bukw8kYw7aDFXTLTA7L+nQsEYepq/dS5zUknqdndWMM//jiAM/N30pFdR3JCR5SEj0kJ3jYtL+cbYeOh7ni2CACv70xn5sGZ3OOYVjPs0xgY4loKKOI12tYuPUQT3y4QbsbCZFLembw1C0DyUhOaHJeHeBH4XIJl/TM5J37R3LNgE62y4lKn28+xPf+tCKk69BQRqGkeA9/vHUgj0zoi0u7Agq6HYcqqaqtD9n7ayijlIjw7Uu78fStg2yXElXSk+J4c9rw837nDwYNZZS7+sKOdMtsbbuMqDCwSzqzpg6ha0ZoP88mQykiiSKyTERWi8h6EXnsjNf/XUSMiGSErkwVKJdL+N2kAQzO1VvDAuVxCXeP6srs74wIS8/2/lynrAbGGGMqRCQOWCAiHxtjlohIF2AcsCukVarTGGPYebiSTumtiPdjtOjBuW34670jWLLtCM9+toUFFu8VjEQv3TWES3pmhm19TYbSNFwzqfA9jfP9nLyO8hTwEPBeSKpTAByqqGZdcRl7jp7g9mE5HK+p55E5a1lXfIxv9mvPVfmdKK+q5Z0VxTx8ZR96tk856z1EhBHd2zGieztW7jrKnJXFbNhXzuaScpLiPWSlJpCVkkC8x83fVu+18Fc6U7zbxYAu6WFdp18tekTEDRQBPYBnjTFLRWQiUGyMWf11F1NFZBowDSAnJ6flFceY1btLufbZhaee13sNU0bmcWX/jizccpjZhXt4f/VeLu6RybgL2pPdJqnJ9xyU04ZB5zkMM8ZQ7/Xy0drwNi1zouQED9MnDyY1MbzjivoVSmNMPTBQRNKBd0UkH3gEuMKPZacD06Gh8UDgpcamAV3SuXN4Lq8t2QnA0u2HuXN4LncMz+XWIV2oqfficbn8Ooz1h4jwy2v7s2TbkZjvO2hMnyxGdg//qZJmtX01xpSKyHzgWqArcHIvmQ2sEJGhxhj9LzbIfjGxH/nZaQzskn7aoanH7cLjDv4J9IzkBB6b2I/vv7ky6O8dSY5a6v+1yVCKSCZQ6wtkK+By4LfGmKxG8+wACowxegYhBNwu4aaCLmFd55X9O5CS6KG8KnaH7Mtp2/RXgVDwZ0/ZEXjF973SBcw2xnwQ2rKUbUcra0nwuCkndkNp6z8kf86+rgG+tlmIMSYvWAUpZzhYXk1lTcNG+bMJfejTIZXlO47w8qIdMbP33H/MTmdmej+lOqcLOqVS+OjlfLJ+P+P7daRVvJtLe2VyzYBO3PXScopLT9guMeRsDVOozezUeSXFe7h+UDatGm2cvdqn8O79I+nfOdViZeExpk9W0zOFgIZSNVtWaiJvTRvBZb3D18rFhrRW4b0+eZKGUgWkdYKHFycXcPuw6G0Qkt2mlZX1aihVwDxuF49f15+fXtnHdilBd8/FXRmca6eHQQ2lOqW23stz87eyfm+Z38uICPd+ozvP3DYoaK2KbMtu04oHx/Wytv7o+BRVUMS5XVyd35EP1+xjwn9/zuzl/vcOfnV+J964ZxjpSXa+hwXLVfkdeWxiP1on2LswoaFUp+nSNolpl3bji33HeOjtNczbcMDvZYfkteWNe4aRkhi5V9pG98pkbF+747poKNVZ0pPi+dmEPtwwqDMLtxxu1rL9OqUxa+oQEuMic9PKz063XYJ2MalC47ONJXz71ULqvJF1Y9BV+R159raLgvJe2sWkcpTL+mTx6FV9bZfRbBv2HbNdgoZShc7tw3PJa2fnTotAVdWErutIf2koVcjEuV388IretstolhsHZ9suQUOpQmt8vw7WGnYH4ni1/T1l5J67VkFXVlnLgfIqep2j461AxXtcDO3alvmbDgbtPUPp5UXb8RpD5/RW9GifzGW9w98oXfeU6pSkBDezFmyntt4b1PcdZaGfm0B5Dby8aAdPfLSBe14pZOWu8I+wraFUp8S5XXyjVyY/fXstwbxUNrJHu6C9VzjVew2HK8LfT4+GUp1mfP8O7Cs7wZOfbArae3bLSMYToSMN2ejRT0OpTiMiPHvbRXy8dh8vLdwelPecu+FAxDUigIaOmG205dUTPeosbVrHM/s7I/wesbgpXSzdl9hSj13bjyv6dQj7ejWU6pyyUhOD9l4926cQ73ER73ZRUe3sTrcG57ahuq6enlkpXD+os5UaNJQq5JITPPz+pgHEu138aPYqjltoNZOS4CEzJYGM5ATf73gyUxLOmJZA+9RE3Ja//2ooVVhMHNCJj9bu46LcNny+OTh9difFu78KVXICGSnxZCYnnha6k4EL5SCvwaahVGHTtnU8q3eXNjlfYpyLDqmJtD/1k3DqcYe0RDJ9QbN5I3IoRedfpRypILcNHdNa0S3TfVrQ2qcm+kKYQPu0RFISPEE7yRSJNJQqbDxuF588eKntMhxPr1Mq5TAaSqUcRkOplMNoKJVyGA2lUg6joVTKYTSUSjmMhlIph9FQKuUwYe0hXUQOAjsDXDwDCE5L5sBpDV9xQh1OryHXGNPskXXDGsqWEJHCQLqA1xqit45orUEPX5VyGA2lUg4TSaGcbrsAtIbGnFBHVNYQMd8plYoVkbSnVComaCiVchjHhVJEbhKR9SLiFZGCRtPHiUiRiKz1/R5zjmXfF5F1NuoQkSQR+VBENvqW+024a/C9Ntg3fYuI/FFa2K/G19TQTkQ+E5EKEXnmjGX+zVfDGhH5u4i0aDCRAGuIF5HpIvKl79/kxpbUEGgdjebxf9s0xjjqB+gL9AbmAwWNpg8COvke9weKz1juBuBPwDobdQBJwGW+x/HA58CV4f4sgGXACECAj0NYQ2vgYuBe4JlG0z1ACZDhe/4k8Itw1uB77THgcd9j18l6wl1HINum4/roMcZsAM7qOMkYs7LR0/VAoogkGGOqRSQZ+CEwDZhtqY5K4DPfPDUisgJo0Qikza0BaAukGmMW+5Z7FbiOhnAGu4bjwAIR6XHGIuL7aS0ih4FUYEug6w+wBoC7gT6++bwEoeVPIHUEsm067vDVTzcCK40x1b7n/wn8Aai0XAcAIpIOXAPMC3MNnYE9jV7b45sWNsaYWuA+YC2wF7gAmBnOGnyfP8B/isgKEfmLiLQPZw2NNHvbtLKnFJH/A841SMMjxpj3mli2H/Bb4Arf84FAD2PMgyKSZ6uORtM9wJvAH40x28Jcw7m+PzZ5zaslNZzjveJoCOUgYBvwP8DDwOPhqoGG7TobWGiM+aGI/BD4PXBnUwsG+bMYSADbppVQGmMuD2Q5EckG3gUmG2O2+iaPAAaLyA4a/p4sEZlvjBkd5jpOmg5sNsY87c97BbmGPZx+yJxNw94qJDWcx0Dfe2711Tkb+GmYazhMw57pXd/zvwDf8mfBINcR0LYZMYevvkOSD4GHjTELT043xjxnjOlkjMmj4cv2l/4EMth1+F57HEgDfhCq9X9dDcaYfUC5iAz3nXWdDDR3L9NSxcAFInLy7ohxwIZwFmAazq78DRjtmzQW+CKcNfjqCGzbbOkZqWD/ANfT8D9+NXAA+MQ3/VHgOLCq0U/WGcvmEbyzr82qg4a9kqFhAzw5/Z5wfxZAAbAO2Ao8g6/VVrBr8L22AzgCVPjmucA3/V7f57CGhnC0s1BDLvAvXw3zgJxQbRNfV0cg26Y2s1PKYSLm8FWpWKGhVMphNJRKOYyGUimH0VAq5TAaSqUcRkOplMP8fwvJpTz8lpm+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ca_poly = gpd.read_file('bnds/California.shp')\n",
    "ax = ca_poly.plot() # plot the boundary file of California"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ce55a039",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Geographic 2D CRS: EPSG:4269>\n",
       "Name: NAD83\n",
       "Axis Info [ellipsoidal]:\n",
       "- Lat[north]: Geodetic latitude (degree)\n",
       "- Lon[east]: Geodetic longitude (degree)\n",
       "Area of Use:\n",
       "- name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands.\n",
       "- bounds: (167.65, 14.92, -47.74, 86.46)\n",
       "Datum: North American Datum 1983\n",
       "- Ellipsoid: GRS 1980\n",
       "- Prime Meridian: Greenwich"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ca_poly.crs  # Lists the coordinate reference system (crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "61924efa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>minx</th>\n",
       "      <th>miny</th>\n",
       "      <th>maxx</th>\n",
       "      <th>maxy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-124.392638</td>\n",
       "      <td>32.535781</td>\n",
       "      <td>-114.12523</td>\n",
       "      <td>42.002191</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         minx       miny       maxx       maxy\n",
       "0 -124.392638  32.535781 -114.12523  42.002191"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xy_nad83 = ca_poly.bounds #store the bounds of the shape file\n",
    "xy_nad83"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16907818",
   "metadata": {},
   "source": [
    "As we see above, the boundary file is in the `EPSG:4269` coordinate system.  This differs from `EPSG:4326` only by the datum. We'll see below that differences will be very minor but the process is the same for any file not in EPSG:4326.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fa78b8f",
   "metadata": {},
   "source": [
    "### Geographic Bounding Box - for searching NASA Earthdata CMR Metadata\n",
    "* we'll determine and store, as a python list, the geographic bounding box of the California shape file\n",
    "* data can be reprojected using the **`GeoDataFrame.to_crs()`** command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "52ae71b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-124.39263575019409, 32.53578253612507, -114.1252281747415, 42.002192463378705]\n"
     ]
    }
   ],
   "source": [
    "ca_poly_4269 = ca_poly.to_crs(epsg=4326)\n",
    "xy_wgs84 = ca_poly_4269.bounds # bound of polygon in WGS84\n",
    "xy_wgs84 = ca_poly_4269.bounds.values.tolist()[0] # We'll save this as a Python list to serve as a subsetting parameter\n",
    "print(xy_wgs84)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37c2e371",
   "metadata": {},
   "source": [
    "### Lambert Conformal Conic bounding box - Daymet's CRS - for subsetting the file using OPeNDAP protocols\n",
    "\n",
    "* Again, using **`GeoDataFrame.to_crs()`**, we'll project the polygon geopanda object (`ca_poly`) to LCC\n",
    "* Daymet is natively in a `lambert_conformal_conic (LCC)` projection.  \n",
    "- We'll present two examples to obtain the LCC paramters used to project the California bounding box to LCC."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de484afa",
   "metadata": {},
   "source": [
    "### Example 1:  Daymet's LCC Projection\n",
    "Use a Daymet file and **`rioxarray`** to obtain the LCC projection parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "32b932b5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CONVERSION[\"unnamed\",METHOD[\"Lambert Conic Conformal (2SP)\",ID[\"EPSG\",9802]],PARAMETER[\"Latitude of 1st standard parallel\",25,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8823]],PARAMETER[\"Latitude of 2nd standard parallel\",60,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8824]],PARAMETER[\"Latitude of false origin\",42.5,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8821]],PARAMETER[\"Longitude of false origin\",-100,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8822]],PARAMETER[\"Easting at false origin\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8826]],PARAMETER[\"Northing at false origin\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8827]]]\n"
     ]
    }
   ],
   "source": [
    "# Obtain the Daymet LCC projection parameters from a Daymet file\n",
    "daymet_x = xr.open_dataset('daymet_v4_tmax_annavg_pr_1950.nc',  decode_coords=\"all\")\n",
    "daymet_x.rio.crs\n",
    "daymet_x_proj = daymet_x.rio.crs\n",
    "\n",
    "# Project the California polygon to LCC\n",
    "ca_poly_lcc_1 = ca_poly.to_crs(daymet_x_proj) # to_crs re-projects from 4269 to LCC\n",
    "print(ca_poly_lcc_1.crs.coordinate_operation.to_wkt(pretty=False))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ac030ee",
   "metadata": {},
   "source": [
    "### Example 2:  Daymet's LCC Projection\n",
    "The Daymet LCC projection parameters are found presented as a `Proj.4 String` in the Daymet [User Guide](https://daac.ornl.gov/DAYMET/guides/Daymet_Daily_V4R1.html).\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "PROJ.4: +proj=lcc +ellps=WGS84 +a=6378137 +b=6356752.314245 +lat_1=25 +lat_2=60 +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +units=m +no_defs\n",
    "</dev>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "ebf31e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the Daymet projection using the Proj.4 String\n",
    "daymet_proj = \"+proj=lcc +ellps=WGS84 +a=6378137 +b=6356752.314245 +lat_1=25 +lat_2=60 +lon_0=-100 +lat_0=42.5 +x_0=0 +y_0=0 +units=m +no_defs\"\n",
    "# Project the California polygon to LCC \n",
    "ca_poly_lcc_2 = ca_poly.to_crs(daymet_proj) # to_crs re-projects from 4269 to LCC\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75583f5f",
   "metadata": {},
   "source": [
    "### Save the CA Bounding Box in LCC Projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bc87818a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Projected bounding box of the California polygon file in LCC projection\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>minx</th>\n",
       "      <th>miny</th>\n",
       "      <th>maxx</th>\n",
       "      <th>maxy</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-1949774.11</td>\n",
       "      <td>-924691.39</td>\n",
       "      <td>-1249285.39</td>\n",
       "      <td>222385.26</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         minx       miny        maxx       maxy\n",
       "0 -1949774.11 -924691.39 -1249285.39  222385.26"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lccbounds = ca_poly_lcc_1.bounds # Bounds in LCC projection\n",
    "print('')\n",
    "print('Projected bounding box of the California polygon file in LCC projection')\n",
    "lccbounds.round(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "223fd49e",
   "metadata": {},
   "source": [
    "### Plot the California Shape Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "311e0720",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEYCAYAAABhi+CNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABDeUlEQVR4nO3dd5xV1bXA8d+aDghDlTKUoSqCQiIiYgNERAVRYyHqg0QT8vyYxDxjijGJz5L3ouZFY4xJjCb2ggKKShGQIggKCIgIhCoMTXofhmHW+2Pvwcsw5c7Mvffcsr6fz/nMnVP2WXeUu+45Z++1RVUxxhhjIi0t6ACMMcYkJ0swxhhjosISjDHGmKiwBGOMMSYqLMEYY4yJCkswxhhjosISjDHGmKiwBGOMiTkRuUlEFojIARHZIiITReSCMvt8R0RURG6oZtvPiUiRiOz3y+ci8r8ikhvZd1F7PtaHqthHRaRTBdtaisiz/m+4X0RWiMj9IlLPbxcR+bH/GxwUkQIReUNEzozG+ynLEowxJqZE5C7gceB/gOZAW+ApYFiZXUcCu/zP6npEVesDzYDvAn2AOaUfvMlARBoDc4E6wHn+/V4KNAQ6+t3+BNwJ/BhoDHQB3gKujEmQqmqLLbbYEpMFyAUOANdXsV87oAT4FlAMNK/GOZ4DHiqzrj6wBfih/70j8AGwE9gBvAw09Nt+Bowpc/yfgcf96xnAQ8BH/r28AzTxbewD5gP5IceeDkzBJcuVwA1+/SjgKFBU2k4F70eBTuWsfwhYCqRVcFxn4BjQO6j/3nYFY4yJpfOAHGBcFfuNABao6hhgOXBz6QYRuUBE9lTnpKq6H/chf2FpM8D/Aq2ArkAb4L/9tpeAwSLS0J8vA7gReDGkyeHAfwB5uGQ1F/gX7iphOXCfP7aeP+8rwKnAt4GnRKSbqj6NS0qPqOopqjq0Ou8JGAiMVdWSCrZfAhSo6ifVbDdiLMEYY2KpCbBDVYur2G8E7kMZ//P4bTJVna2qDWtw7s24BICqrlbVKap6RFW3A38ELvbbtgCzgOv9cYN9zAtD2vqXqq5R1b3ARGCNqk717+sN4Bt+vyHAelX9l6oWq+qnwBjguhrEX1YT3FVZTbdHnSUYY0ws7QSa+quCconI+UB74DW/6hXgTBHpWctz5+FuUyEip4rIayKySUT24a5amobs+zxwi399CydevQBsC3l9uJzfT/Gv2wHnisie0gV3Ndailu8F3N+yZS22R50lGGNMLM0FCoGrK9lnJO4W1mIR2Qp87NePqOlJReQU3C2lD/2q/8U92zhLVRvgkoiEHPIWcJaIdMddhbxcw1NvBGaqasOQ5RRVvd1vr005+6nANSJS0ef4NKC1iPSqxTlqxRKMMSZm/C2l3wJ/EZGrRaSuiGSKyOUi8oiI5AA34B6A9wxZfgTcXNmVT3lEJFtEzsYljN245yTgHvofAPaISB7uwX5onIXAm7irp09UdUMN3i7Au0AXEfkP/z4zReQcEenqt28DOoTRTpaI5IQs6bjbeg2A50WknX+/eSLyRxE5S1VX4XrnvSoi/USktI3hIvLLGr6farEEY4yJKVX9I3AX8GtgO+5b/g9xSeBq3C2mF1R1a+kCPAuk4x6+XygiB6o4zc9FZD/ultgLwEKgr6oe9NvvB74J7AXeA8aW08bzwJmcfHssbL5zwSBcp4DNwFbgYSDb7/IscIa/ffZWJU0tw/1dSpfvquouoC+uJ9rH/v1O8+9ptT/ux8CTwF+APcAa4Bpcz7eoE9+dzRhjTAgRaQusAFqo6r6g40lEdgVjjDFl+OcadwGvWXKpuWrdzzTGmGTnx65sA77EdVE2NWS3yIwxxkSF3SIjZoX3DoQsS/y2fN9m6fr1ob07RGSYiCwWkX0iskNEpolIfsj2zr4v/3a/zyoR+bOItK7F3+I+H9PAkHU/88Xy9ovIOhH5WTnH3em3HRSR5SLSxa8XEblXRDb4GF8TkQY1jc8YkzhS/haZL7z3S+A/gcm4ukCDcYX3ZofsGlp4b3Q1T/OIqv66ku0NVbVYRM4DponIYlwvkBeAa3E1k07B9UYp8XF3wo0PeA74hqoWiMipwE3ABXw9SC1sItIRN8K47OhfwY1B+AxXFuN9Edmoqq/5474H3IYroLcc1+1ytz92BK6kxvl+3cu4uk41KWBoAtK0aVPNz88POgwThxYuXLhDVZuVty2lb5GJK9+9Cdfl741K9msHrMOVjngdyFPVbRXtX+bY53D1gE5KMP5qZB2QWVo6Q0Tm+3OsB36tqj0raPclILcG9Ysqi3Ui7sP/KeB7qjq1gv2ewP2/8yP/MPRL4DuqOq2cfd8EPlbVR/3vfXEJs7GqHopU7Ca6evXqpQsWLAg6DBOHRGShqpY7mDPVb5EFUnivPP5W0vlAN2AR8Clwuog8JiL9xY1EDjUQV9OosjbbSkiJinKWm0L2vR4oUtUJVcWJKxi4zK9q7ZfuIrLR3ya7X74eXSycOEJacGMAOld2HmNM4kv1BBOrwnt3l/lgf77M9h2422/PAL9U1Wmquhboh6ufNBrY4Z/nlCaaprhBWwCIyA992wdE5B8+tg1lSlSUXV7xx56Cm5vjJ1W8D3AVZ9P4ekR06fOeQbhBaf1xFWNv8+snAt/zz5tygV/49XXDOJcxJoGleoKJVeG9P5T5YC/7/KGpqjZS1a6q+kTpSlWdp6o3+PubFwIXAfeGxN4yZN8nfaJ7HMisRmzgRjW/qKrrKttJRH6IS7ZXquoRv/qw//mIqu5R1fXA34Er/Pp/Aq/i5tBYBkz36wuqGaMxJsGkeoIJpPBeTajqfFw5i+5+1TRcB4AK+VtkBypZSm/1XQL8WES2+vfYBhgtIr8IaetWXGeIS1Q1NDmsxHWMKPdhnqqWqOp9qpqvqq1xSWaTX4wxSSylE0ysC+9Vh3+2833fMwwROR24Cpjnd/lv4EJxhe3y/D5NcZMnlb6/Db5ya0VLaYXYS3CJq/T9bQZ+gKtfhE9E/wNc6m/dEXKOQ7hOCT8Xkfq+i/T3cUX+EJHGItLRP2M6A1eg74FKJkkyxiSJlE4wENPCe6FXDjvCCG0PLqEs9e1PwnVGeMTH/W/cPOOtgSXiCt3NwSWH34T7/n1bO8u8v2PAblUtfV8P4Z5XzQ95D38LaeKHuMq0m3FXha/gbo2Be1Y0ATiIex7zTz+TnzEmyaV0N2VjTHism7KpiHVTNsYYE3MpP5LfGBMbq1evpqSkhA4dOpCRYR89qcCuYIwJgIi0EZHp4uq2LRORO/36xiIyRVxduSki0ijkmHtEZLWIrBSRy0LWny0iS/22J/xg2NLZHF/36z+WE+vYjfTnWCUiMSnb8+GHH/Lqq6/y6KOP8sYbb7B06VIKCwtjcWoTEPsaYUwwioGfquqnIlIfWCgiU4DvANNU9ffiCp/+EviF74E3HFfpoRUwVUS6qOox4K+4no7zcB0qBuM6VNyG66zRSUSG42ZSvFFEGgP3Ab1w3csXish4Vd1NlBw8eJCNGzcCUFRUxBdffMEXX3xBWloarVu3pkOHDnTo0IG8vDzS0ux7b7KIaYKxgnkmnlVWtC/SVHULvqioqu4XkeW4qg3DcBUcwE3ZOwNX/WAYbvKrI8A6EVkN9BaR9UADVZ0LICIv4Ho/TvTH/Ldv603gSX91cxkwxU+5i09sg3EDYqNi1qxZlNehqKSkhA0bNrBhwwZmzJhBdnY2+fn5xxNOkyZN8BdkJgHFNMHk5+djPVFMvBKRLwM6bz7wDdwg3uY++aCqW0rHQeGSz7yQwwr8uqOcWBWhdH3pMRt9W8UishfX3fz4+nKOCY1rFO7KiLZt29b4/ZWUlFBQEF7hhiNHjrBy5UpWrlwJQIMGDejYsSOdOnWiQ4cO5OTk1DgOE3t2i8yYAPk6cGOAn6jqvkq+rZe3QStZX9Njvl7hxis9Da6bckWBVSUtLY3rrruOJ554ouqdy9i3bx+LFi1i0aJFiAitW7emU6dOdOrUiZYtW9rVTZyzBGNMQEQkE5dcXlbVsX71NhFp6a9eWgJf+fUFuBI+pVrjBrYW8HXB0dD1occU+KoTubiiqgV8fRuu9JgZEXpb5crMrG55vJOpKhs3bmTjxo1Mnz6dunXr0qlTJzp27EjHjh2pV69eBCI1kWQJxpgA+GchzwLLfTWJUuNx9e9+73++HbL+FRH5I+4hf2fgE1U9Jm6m0T64W2wjcHP6hLY1FzeR3AeqqiIyGfifkB5qg4B7ovRWAXfrK9IOHTrEZ599xmeffQZAq1atjl/dWGeB+GAJxphgnI+b6XOpuBlMAX6FSyyjReQ2YANukjtUdZmIjAa+wPVAu8P3IAO4HTezaR3cw/2Jfv2zwIu+Q8AuXC80VHWXiDwIzPf7PVD6wD9aopFgytq8eTObN29m1qxZ5OTk0KFDh+NXN7m5uVE/vzlZ2AlGRNKBBcAmVR0iIo8CQ3GVdNfgZoXcE5UojUkyqjqb8p+FgCs+Wt4xvwN+V876BXxdZTt0fSE+QZWz7Z98XS8u6oqKimJ1KgAKCwuPd4UGaNas2fHOAu3atbOBnjFSnb/ynbjZHBv436cA9/jeKQ/jLrF/UdHBxpjUdfDgwUDPv337drZv3868efPIyMggPz//+NVN06ZNrbNAlISVYHwJ9itx357uAlDV90N2mYe7x2uMMSdZvnx50CEcV1xczOrVq1m9ejUAubm5dOjQgU6dOtG+fXvq1KkTcITJI9wrmMeBnwP1K9h+K25OkJPUpC/9/ffD9OmQkQHp6e5n6RXt0aOVL6VjuUq/kIT7My3t6yU9/cSfFa0rjat0CY01I8O1XRpPuD//8Y+vYzImGRw+fJgVK1YEHUaF9u7de0JX6Ly8vONXN61atSI9PT3oEBNWlQlGRIYAX6nqQhHpV872e3EPHV8uuw1q1pd++XKYOTOcPZPPP/4RdATGRNbs2bM5duxY1TvGAVWloKCAgoICZs6cSVZWFm3atKFdu3bk5+dbwqmmcK5gzgeuEpErgByggYi8pKq3+CJ5Q3DT6EZsYplUnqLGrl5MMtm4cSMfffRR0GHUWFFREWvWrGHNmjWAG89TNuFYh4GKVfmXUdV78H3k/RXM3T65DMY91L/YT5sbEaqwaFGkWjPGBKWoqIhx48YFHUZEHT16lLVr17J2rZs5PCMj43jCOe2002jRokXAEcaX2qTeJ4FsYIrvgTFPVf+ztgHt2weNGlW9nzEmvq1du5bdu6NWoDkuFBcXs27dOtatW0daWpolmDKqlWBUdQa+pISqdopCPOTmwkcfwZAhMGlSNM4Qv+z2mEkmBw4cCDqEmLJK8SeLy1oK6ekwdGjQUcSeJRiTTFIpwWRkZNCqVaugw4g7cZlgAPr3DzoCY0xN7du3L6Wm5mjTpo31LitH3CaY00+HVLudaVcwJhkUFxczevTowEfvx5LdHitf3CYYkdS7irEEY5LBpEmT2LRpU9BhxFS7du2CDiEuxW2CATjzzKAjMMZUx6effsrChQuDDiOmMjIyyMs7aUJQQ5wnmL17g44gtuwKxiSyTZs2MWHChKDDiLnWrVvbYMsKxHWC2bkz6AhiyxKMSVQHDx5k9OjRCVMSJpLs9ljF4jrBpFAnFMASjElMJSUljBkzhn379gUdSiDsAX/F4jbBrFoFixcHHYUxpirTpk1j3bp1QYcRiPT0dHv+Uom4TTBvvBF0BLFnVzAm0Wzfvj2hi1nWVl5eHpmZmUGHEbfiNsFMmxZ0BMaYqqT6w227PVa5uE0wp50WdATGmKqkeoKxB/yVi9sEc+65QUdgjKlKPE2FHGtpaWm0adMm6DDiWtwmmPoVTc5sjIkbnTt3JisrK+gwAmHPX6oWtwnmmWeCjsAYU5VGjRpx+eWXBx1GINq2bRt0CHEvLhPMmjWpNxeMMYmqR48edO3aNegwYs4e8FctLhPM3//upk42xsQ/EWHIkCEpdatMROz5SxjiLsEcPgzPPht0FMaY6qhbty6dO3cOOoyYadmyJdnZ2UGHEffiLsGMHg27dgUdhTHJT0QGi8hKEVktIr+sbXtdunSJRFgJwZ6/hCfuEsxTTwUdgTHJT0TSgb8AlwNnAN8WkTNq0+bpp5+eMr2q7PlLeOIqwSxYAJ98EnQUxqSE3sBqVV2rqkXAa8Cw2jSYlZVFt27dIhJcvLMrmPDEVYKxqxdjYiYP2Bjye4Ffd5yIjBKRBSKyYPv27WE12rNnz4gFGK+aN29OnTp1gg4jIcRNgtm1C159NegojEkZ5ZVWPaHvpqo+raq9VLVXs2bNwmq0bdu21K1bNxLxxS0rDxO+uEkws2ZBYWHQURiTMgqA0H62rYHNtW1URJL+OYwlmPDFTYIZNAgaNgw6CmNSxnygs4i0F5EsYDgwPhINS5LPO2EJJnxxk2Dq1oWhQ4OOwpjUoKrFwA+BycByYLSqLotE20VFRZFoJi41bdqUevXqBR1GwoirWtt2i8yY2FHVCcCESLZ59OhRDh06FMkm44pdvVRP3FzBADz+OFx4YdBRGGNqaleSj5K2BFM9YScYEUkXkUUi8q7/vbGITBGRVf5no9oG06oVfPAB/OIXtW3JGBOEOXPmBB1CVFmCqZ7qXMHcibtXW+qXwDRV7QxM87/XWkYG/P73MH485OREokVjTCxs2rSJpUuXBh1G1DRq1IgGDRoEHUZCCSvBiEhr4EogdJaWYcDz/vXzwNWRDGzoUJdojDHxT1WZMmVK0GFElV29VF+4VzCPAz8HSkLWNVfVLQD+56mRDQ1+9CMYMCDSrRpjIq24uJgvv/wy6DCiyhJM9VWZYERkCPCVqi6syQlqUm7ieHBp8MorkJdX9b7GmOBkZMRVh9SosARTfeFcwZwPXCUi63EF8QaIyEvANhFpCeB/flXewTUpNxGqeXN44w1I8sHBxiQ0EUnqJNOgQQMa2kjwaqsywajqParaWlXzcaN9P1DVW3Cjfkf63UYCb0cryPPOgz/+MVqtG2MiIZkTTH5+ftJXKIiG2oyD+T1wqYisAi71v0fND34AubnRPIMxpjaSuQaZleevmWolGFWdoapD/OudqnqJqnb2P6M6wiozE4YMieYZgqda9T7GxKtkv4Ix1RdXI/mrcvnlQUdgjCnPjh072Lt3b9BhREW9evVo3Lhx0GEkpIRKMD16BB2BMaY8kyZNoqSkpOodE1C7du3s+UsNJVSC6dLFjfQ3xsSX9PT0oEOIGuueXHMJlWAOHoTi4qCjMMaUdfHFFwcdQtTY85eaS6gEM3du0BEYY8rTqlUrunTpEnQYEVenTh1qMn7POAmVYJK8UKv1IjMJLRmvYtq2bWvPX2ohoRLMJ58EHUF0WYIxiaxVq1ZJdzvJnr/UTkIlmDPPDDqC6LIEYxLdeeedF3QIEWUJpnYSKsHceGPQEUSXXYmbRNe5c2eaNGkSdBgRkZWVRYsWLYIOI6ElVILp3RvsC4Ux8UtE6NOnT9BhRESbNm1IS0uoj8i4k1B/PRG44Yago4geu0VmkkGPHj2SovKw3R6rvYRKMJDct8mSdCC0STGZmZlcffXVQYdRa5Zgai/hEsw3vwnJOi22XY2bZNGuXTv69u0bdBg1lp6eTqtWrYIOI+El3EfamjWwb1/QUURHcTF89FHQURgTGeeff37QIdRYmzZtkro6dKwkXIJZtCjoCKJrwAB49dWgozCm9urUqZOwH9LWeywyEi7BfOtbcPPNQUcRPUeOwE03we9/bw/9TWITkYQtc7948eKknX4glhIuwaSlwb/+BcOGBR1JdN1zD9xxhxX3TEYi8qiIrBCRz0RknIg0DNl2j4isFpGVInJZyPqzRWSp3/aE+PolIpItIq/79R+LSH7IMSNFZJVfRoasb+/3XeWPzYrWex01alS0mo6qwsJC3nrrLdS+5dVKwiUYcLNbvvaau52UzP76V7j2WldF2iSVKUB3VT0L+DdwD4CInAEMB7oBg4GnRKS0Dv5fgVFAZ78M9utvA3araifgMeBh31Zj4D7gXKA3cJ+INPLHPAw8pqqdgd2+jahI5HEk69evZ968eUGHkdAS9r9+Tg6MHQvdugUdSXS98w707w/btgUdiYkUVX1fVUuvTecBrf3rYcBrqnpEVdcBq4HeItISaKCqc9V9pX4BuDrkmOf96zeBS/zVzWXAFFXdpaq7cUltsN82wO+LP7a0rahI5GKR06ZN46uvvgo6jISVsAkGIDcX3n0XTj016Eiia/586NMHli8POhITBbcCE/3rPGBjyLYCvy7Pvy67/oRjfNLaCzSppK0mwJ6QBBfa1glEZJSILBCRBdu3b6/Rm5s3b15C32Y6duwY48aN49ixY0GHkpASOsEA5OfD+PHuiiaZrV8P550HH3wQdCQmHCIyVUQ+L2cZFrLPvUAx8HLpqnKa0krW1+SYyto6caXq06raS1V71WROlK1btzJt2rRyt5Ve1aSnp1OnTp24Hvm/detWZs6cGXQYCSnhEwzAuefCSy8l/3TKe/fCZZe5Tg4mvqnqQFXtXs7yNrgH8MAQ4Gb9+it+AdAmpJnWwGa/vnU56084RkQygFxgVyVt7QAa+n3LthVR27Zt41vf+hY9evQ4aZuqctFFFzFq1CjatWvHKaecEo0QImb27Nls3Lix6h3NCZIiwYDrvjxlCjRtGnQk0VVcDLfeCr/5jXVjTlQiMhj4BXCVqh4K2TQeGO57hrXHPcz/RFW3APtFpI9/hjICeDvkmNIeYtcBH/iENRkYJCKN/MP9QcBkv2263xd/bGlbEdWjRw+6du3KsGHDTiq7kpeXx0UXXcSpp57KtddeS/v27QFXwTg3N5e6detGI6QaU1XGjh1LYWFh0KEklKT6zt+vHyxYAEOHwtKlQUcTXQ89BOvWwbPPQnZ20NGYanoSyAam+FtF81T1P1V1mYiMBr7A3Tq7Q1VLb/7fDjwH1ME9syl9bvMs8KKIrMZduQwHUNVdIvIgMN/v94Cq7vKvfwG8JiIPAYt8G1EjIvTr14/nn3/++LrDhw9z7Ngx0tLSOHToEL179+biiy8mPd11mtu2bRt/+9vfju/fuHFj9uzZQ0mABfv27NnDhAkTuPbaawOLIdFILB/A9erVSxcsWBD183z4IVx0UdRPExf69YO33nIdHkztiMhCVe0VdBzxKBL/dhcsWMDEiRMpKSmhf//+nH/++fz5z38+PqAxLS2NevXqUa9ePU4//XRycnJYtWoVzZs357zzzuO9995jy5YtgQ+AvOaaazjrrLMCjSGeVPbvJmlukYW64ALo3j3oKGJjxgyXTDdH5S66MZHTq1cvrrjiCsDV+jp69OjxKxaAkpIS9u/fz9atW5kxYwbvv/8+devWpWvXrkydOpWCggKy4+By/b333mP37t1Bh5EQkjLBiMDTT7sH4qngs89cD7MVK4KOxJjKnX322VxzzTU0adKEnJycSq8ESkpKWLp0Kc8++yxLlizhwIEDcTEmpaioiHHjxgV6uy5RJGWCAfeBe+edQUcROxs2wPnnw9y5QUdiTOU2b97MCy+8wNGjR1mRoN+KNm7cyIcffhh0GHEvaRMMuBHw9epVvd+gQeHNxRLv3aB37YJLLnGDT42JRyUlJSxZsoSdO3fy1Vdf8c1vfjNhR/rPnDmTgoKCqndMYVV+rIpIjoh8IiJLRGSZiNzv1/cUkXkistiP9u0d/XCrR8TVLStPkybQt6+7vZSWBgMHuvX16rlqzT/96Yn7P/CA2/+CC6Ibc20dPgxXX+16lxkTbzZs2EBhYSFNmzalVatWnHPOOfTqlZj9KlSVcePGUVRUFHQocSuc7+RHgAGqekBEMoHZIjIReAC4X1UnisgVwCNAv+iFWn0TJ8KePeVvO3jQjYrPzoYJE2D2bLjtNhgyxCWcAwfgtNNcj7ROneBXv3JdnxOhgvexY/C978GWLXDvvS7RGhMPSm+J5eXlHb9yqV+/fpAh1cquXbt4//33GTJkSNChxKUqE4wfmHXA/5rpF/VL6eTFuURpNHBtvPSSew5TUOAKY4b2yD799K9vi33+OcyZA0uWwP33w8qV7kO6QQOXXAoLYcQIeOONYN5HTf3mN7BpEzz5JIR01jEmMP/+978BWLlyJbt376ZRo0ZxN6iyuhYuXEiXLl3o0qVL0KHEnbCeKviS4QuBTsBfVPVjEfkJMFlE/oC71VbuBNwiMgpXZpy2bdtGIuawPf00NGrkvsE/+ST86Edfb1u8GMaMcc9VbrwRyusQsm8ffPqpWxLV3/7mKjG//DLUqRN0NCbVdevWjdmzZ1NYWMjrr7/OrbfemvAJBmD8+PHcfvvt1AvnoW8KCeshv6oeU9WeuLpFvUWkO25k8X+pahvgv6hgNHBtC+bVRuPGX98e+uEP3dWJr0hBTg5cd527SknQSffCNm6c68hgXfdN0Dp16nT89bZt2xg/fnxSJJiDBw/yzjvvJHTl6GioVi8yVd0DzMBNdjQSGOs3vYGb1Ciu/fa38MUX8JOfuNten37qPnh/85ugI4u+2bPdgMwtW4KOxKSyrKwTJ89ctmwZq1evpk+fPgk9ORm4236LFi0KOoy4Ek4vsmalU7qKSB1gILAC98zlYr/bAGBVlGKMqJwceOwxl2h27XLPWw4fDjqq2Pj8c9cLbu3aoCMxqaq8YpGzZ89m3rx5STFwcdKkSezatavqHVNEOF8ZWgLTReQzXOG8Kar6LvB94P9EZAnwP/jnLIli+3a4/HJ3NTNyJNx1V2oUjVy71iWZZcuCjsSkorlJPhL46NGjNso/RDi9yD4DvlHO+tnA2dEIKhZKb5Vu3QrNm8P//R8cOuQeiie7LVvc7bIJE9xcOsbEQklJSUoMTCwoKODDDz/k4osvrnrnJJfYNz1rofRKvU6drzsC3HJLcPHEWumo/womHDQm4tLS0jjnnHOCDiMmZs6cyaZNm4IOI3Apn2A+/xwmT3av+/Z1UzCnioMH4YorXLl/Y2Khd+/eZMR7zaUIsFH+TsommKNH3c/9+2HwYPjZz9y6m28ONq5YKypys4GGzAVlTNTUq1ev3CmUk9HOnTuZOnVq0GEEKmUTzNVXu9H5pf7wB/fw+7LLXMmYVBovVVIC3/kO/OlPQUdiUkGi1h6rifnz57NmzZqgwwhMyiaYjAy4774T182f765mJk6EI0eCiStIP/mJ+5vYWDETTS1atKBz585BhxEzb7/9NodTZSxEGSmbYAA6dIBvfvPEdYcOudkhi4uDiSloDzzgqh5YL0sTTddddx3tS8tqJLn9+/czceLEoMMIRPI/batCixZBRxB/nnrK9TJ7/nkoM/DamIjIyspixIgRFBUVsWPHDiZMmJDUva6WLl3KaaedRrdu3YIOJaZS+goG3Hww5mSvvQZXXeV6mhkTLVlZWbRq1YqePXsGHUrUvffee+zfvz/oMGIq5RNMinRoqZHJk91EbFb5wkRb165dE74WWVUOHz6ccgUxk/u/aBheeQW6dw86ivg1b54b9Z/Edy9MHKhbty4NGzYMOoyoW7VqVUoVxEz5BNOgAbzzDjRtGnQk8WvZMteFe1VClDM1iWjFihUpUyRy8uTJ7E6RuTNSPsGAG73/yis2tXBl1q93SSaFvnyZGFBVSkpKmD59etChxExRURFvvfVWShTEtATjXXqpG2xpUwtX7KuvoF8/mDUr6EhMMlixYgXjx49n3LhxbN++PehwYmrDhg3MmTMn6DCizhJMiLvuct/Q77wTHnoInnnGbp2VtW+fq3bwzjtBR5L4RORuEVERaRqy7h4RWS0iK0XkspD1Z4vIUr/tCRF3vS0i2SLyul//sYjkhxwzUkRW+WVkyPr2ft9V/thAOqMvWrSIxYsX8/nnnwdx+sDNmDEjqbtmgyWYk5x5Jjz+ONx7rysZ8+STQUcUfwoL4dprXVdmUzMi0ga4FNgQsu4MYDjQDTdr7FMiUnpN/VfcnEud/TLYr78N2K2qnYDHgId9W42B+4BzcbPN3icijfwxDwOPqWpnYLdvI+aOpGK5jBAlJSWMHTs2qQtiWoKpwg03wHXXBR1F/Ckuhptugn/8I+hIEtZjwM+B0D6rw4DXVPWIqq4DVgO9RaQl0EBV56rr4/oCcHXIMaWlSt8ELvFXN5fhJgfcpaq7gSnAYL9tgN8Xf2xpWzGVzB+s4dq1a1dSj/K3BFMFERg0KOgo4pMqjBoFf/xj0JEkFhG5CtikqkvKbMoDNob8XuDX5fnXZdefcIyqFgN7gSaVtNUE2OP3LdtW2ThHicgCEVkQjWckqX4FU2rx4sV88cUXQYcRFSlfKiYc8+YFHUF8++lP3bQHv/2t9cQrJSJTgfIKEd0L/Aoo72tLeX89rWR9TY6prK0TV6o+DTwN0KtXr4iPDrQE87V33nmHvLw8cnNzgw4louwKJgyWYKr23/8Nd99tlZhLqepAVe1edgHWAu2BJSKyHmgNfCoiLXBXE21CmmkNbPbrW5ezntBjRCQDyAV2VdLWDqCh37dsWzHVsWPHIE4blwoLC5Oy67IlmEq89ho89xwk6dVrxP3xj+6W2bFjQUcSv1R1qaqeqqr5qpqPSwTfVNWtwHhguO8Z1h73MP8TVd0C7BeRPv4Zygjgbd/keKC0h9h1wAf+Oc1kYJCINPIP9wcBk/226X5f/LGlbcXUFVdcQaNGjareMUWsX7+euXPnBh1GRFmCqcBXX8G3vw3f/W7QkSSWZ55xs4KWzhhqwqeqy4DRwBfAJOAOVS1N17cDz+Ae/K8BSp8MPws0EZHVwF3AL31bu4AHgfl+ecCvA/gFcJc/polvI+ays7P51re+lfQ1yKrjgw8+YPPmQC4oo0JiWXitV69eumDBgpidrzZKSqBxY9i7N+hIEtOQITB6NNSpE3Qk4RORhaqaOtMtVkM0/+2OHz8+pepzVaVJkyaMGjWKrASZK6Oyfzf21aECaWlwzjlBR5G43n0XrrzSPfw3pjKpNkdKVXbu3MnkyZODDiMiLMFUwhJM7Uyf7krwpEgNQ1ND7du3p04iXerGwKeffsry5cuDDqPWLMFUonfvoCNIfB9/DP37w7ZtQUdi4lVaWhpdu3YNOoy4884777Bv376gw6gVSzCVOO889xzG1M5nn8GFF8KGDVXva1KT3SY72eHDhxk3blxCd122BFOJ5s1h/nyw//drb9Uqm1PGVCw/P5+6desGHUbcWb9+PR999FHQYdSYJZgqdOgAc+fC0KFBR5L4Nm50VzJLlwYdiYk3aWlp9OzZM+gw4tL06dMTtupylQlGRHJE5BMRWSIiy0Tk/pBtP/JlxZeJyCPRDTU49evDG2/AGWcEHUni27YNLr4YPvkk6EhMvOnXrx9NbX6Mk5SUlDBmzJiELK0TzhXMEWCAqvYAeuIqsvYRkf64Sq5nqWo34A/RCzN42dnw7LNWaysSdu+GSy6BGTOCjsTEk8zMTAZZZdly7d69OyGrLleZYNQ54H/N9IviRhb/XlWP+P2+ilqUcaJPH/jxj4OOIjkcOACXXw4TJgQdiYknyVbsMZKWLFnC0gS7vxzWMxgRSReRxcBXuDkmPga6ABf6mfFmiki5o0aiXfI71h56CNq1CzqK5FBYCMOGwZgxQUdi4oU96K/ce++9x549e4IOI2xhJRhVPaaqPXGVV3uLSHdcqf9GQB/gZ8Do0mlcyxz7tKr2UtVezZo1i1zkATnlFHj66aCjSB7FxXDjjfDmm1Xva5JfvXr1uPLKKxOmTEqsHTlyhLFjxyZM1+Vq9SJT1T3ADNx0rQXAWH8L7ROgBEiJJ3SDBsEPfxh0FMnj2DEYPtx1pDCpTUTo1asXw4cPDzqUuLVx40ZmzZoVdBhhCacXWTMRaehf1wEGAiuAt3BTryIiXYAs3FwTKeGJJ+DBB4OOInkcO+aqV7/+etCRmHiQn59PkyZNgg4jbs2aNYsNCTByOZwrmJbAdBH5DFf2e4qqvgv8E+ggIp8DrwEjNZalmQMmAr/+tasabCLj2DFX6t+SjBERzj777KDDiFuqytixYyksLAw6lEqF04vsM1X9hqqe5Wfle8CvL1LVW/y6b6rqB9EPN/7ceGPQESSXY8fgppvcZG8mtfXo0YPMzMygw4hbe/fu5d133yWev9fbSP5auuoq6Ns36CiSS0mJu5J59dWgIzFBqlu3LrfffjsdOnQIOpS4tWzZMpYsWRJ0GBWyBFNLDRrAnDmwZg3k5QUdTfIoKYFbboFXXgk6EhOkRo0a0de+wVVqwoQJ7IrTOTEswURIhw7w4os20j+SSkrgP/4DXn456EhMkAoKCoIOIa4dPXqUMWPGcOzYsap3jjFLMBHUvz/87GdBR5FcSkpgxAh46aWgIzFBOHjwYEJXE46VzZs3M3369KDDOIklmAh78EG4805o0yboSJJHSQmMHGlJJhUdOXKEoqKioMNICHPmzGHdunVBh3ECSzARlpUFjz/unsk0bBh0NMmj9ErmxReDjsTEUuPGjRkxYgQNGjQIOpSEMG7cOA4dOhR0GMdZgomSzEy7iok0VXcl88ILQUdiYql9+/bcfvvtdO/ePehQ4t7+/ft555134qbrsiWYKNm50ybWigZV+M534Pnng47ExFJOTg4DBw4MOoyEsGLFChYtWhR0GIAlmKg5fDjoCJKXKnz3u/Dcc0FHYmIpGaqxx8qkSZPiouuyJZgoad0aZs50YzlycoKOJvmowq23wr/+FXQkJhZUlZkzZwYdRsI4evRoXFRdtgQTRRdd5B5Kb94M118fdDTJRxVuuw3++c+gIzHR9u9//9vGw1TTpk2b+PDDDwONwRJMDDRqZAkmWlThe99LvCQjIj8SkZUiskxEHglZf4+IrPbbLgtZf7aILPXbniide0lEskXkdb/+YxHJDzlmpIis8svIkPXt/b6r/LFxP/lKvBd1jFczZ85k06ZNgZ3fEkyMrF8fdATJK9GSjIj0B4YBZ6lqN+APfv0ZwHCgG27OpadEJN0f9ldgFNDZL4P9+tuA3araCXgMeNi31Ri4DzgX6A3cJyKN/DEPA4+pamdgt28jruVZHaYaUVXGjBnDkSNHAjm/JZgYmTcv6AiSmyr86U9BRxG224Hfq+oRAFX9yq8fBrymqkdUdR2wGjeDbEuggarO9VNivABcHXJMaZ+6N4FL/NXNZbipNXap6m5gCjDYbxvg98UfW9pW3LLR/DW3e/duJk2aFMi5LcHEyNy5QUeQ/NIS5//mLsCF/jbVTBE5x6/PAzaG7Ffg1+X512XXn3CMqhYDe4EmlbTVBNjj9y3b1glEZJSILBCRBUH34Dp69Gig5090ixcvZtmyZTE/b+L8k0xwv/td0BEkv4KC+LkVKSJTReTzcpZhQAbQCOgD/AwY7a8syiuVqpWspwbHVNbWiStVn1bVXqraq1mzZuXtEhPHjh1j8+bNgZ0/Wbz77rvs3bs3pue0BBMjp54adATJb8cO1y08HqjqQD8ZX9nlbdxVw1h1PgFKgKZ+fWj9h9bAZr++dTnrCT1GRDKAXGBXJW3tABr6fcu2FZc+/fTTuBjTkegKCwsZN25cTLsuW4KJkY0bq97H1N5HH0FAt5ur4y3ccxBEpAuQhfvgHw8M9z3D2uMe5n+iqluA/SLSx1/pjADe9m2NB0p7iF0HfOCf00wGBolII/9wfxAw2W+b7vfFH1vaVlzKysoiLYHuf8azL7/8kjlz5sTsfPZfLUYGD656H1N7qnD55fDxx0FHUql/Ah1E5HPgNWCkv5pZBowGvgAmAXeoaukkH7cDz+Ae/K8BJvr1zwJNRGQ1cBfwSwBV3QU8CMz3ywN+HcAvgLv8MU18G3GrR48ejBw5ErHJliJixowZMeu6nFH1LiYSVq0KOoLU8txzcO65QUdRPlUtAsq9maeqvwNOemKnqguAk6o9qmohUO4oK1X9Jy6ZlV2/Ftd1OWG0bduWvLw8G2wZASUlJYwZM4Yf/OAHZGdnR/VclmBiZNasoCNILVOmwJEjEOV/PyaGrrrqKjZv3szu3buPL9u3b7dBmDVQ2nV52LBhUT2PJZgYueMO94H3l79AHE3XkLTWrHF1yv7zP4OOxERKs2bNKNub7dixY6xdu5ZZs2bZ1U01LV68mE6dOtGtW7eoncOewcRIixbwyCPw1FMnb2vWDHJzYx9TstuzJ+gITLSlp6fTuXNnRo4cSc+ePYMOJ+FEu+uyXcHE2C23uAGBeXnQtq2rulxabXnfPnjmGXj4YfjKj+3OzIR69ezDsib+93/d1eIDDwQdiYm2jIwMrrrqKvr27cvWrVvZunUr27Zt48svv6S4uLjqBlJUadflESNGRKWnnsRy5rNevXrpggULYna+RHX4MOzfDw0afP0M4eOP4dFHYezYYGNLRP37u8TdoUPl+4nIQlXtFZuoEkui/tvdvn07Y8aMYdu2bUGHEtcGDBjAhRdeWKNjK/t3Y7fI4lCdOm5gZk4OiLilTx94803429/cVY0J3/TpMHAgHDgQdCQm1po1a8att95KixYtgg4lrkWr67IlmAQiAj/4gZvIrGXLoKNJLLt3u3l5TOrJyspi+PDh1KtXL+hQ4lZp1+VIV122BJOAzjvPFc+0K5nwZGS4W4tdugQdiQlKbm4uN9xwg1UEqEQ0qi5X+dcWkRwR+URElvjJke4vs/1uEVERaRrRyEyl2rWDG28MOor4d9ZZ7oqvf/+gIzGRUFxczMGDBykuLqa6z4/btm3LufE6+jZOLF68mJUrV0asvXDS+RFggKr2AHri5pToAyAibYBLgQ0Ri8iE7R//gL/+FfLzg44kPv3qVzB7NvTtG3QkJlKKi4t5/PHH+d3vfsejjz5KUVERhYWF/Otf/2LBggVVJp02bdpUut3AgQg+rKwywfgaSaVnzPRL6X/Fx4CfU0G5bxM+VSgudj8ffRR++1sonb5hzx749NOTj8nJcQMJ//1vVxrl9NNjGXF8O/NMN0VC/fpBR2IiKScnh+7dXcWcw4cP89lnn5GTk0NRURHvvfcezz77LCtXrmThwoXs2LHjpOM7duxI165dYx12QolkEg7rhqSIpIvIYuAr3Cx5H4vIVcAmVV0SsWhS0MSJMGAANGrkpv0VcYniwQehe3fo2hV69HBXKxV9OcvMhJEjYfly2LoV3n8f3n0XFi50D7ZXr4ZTTont+wpSejrcfXfQUZhoueKKKxg6dCjf+c53jiebm2++mWHDhnHqqaeycuVKtm/fzp5yBo9lZWVx/fXXc+mll1rxzApsjGDp92qNgxGRhsA44E7gH8AgVd0rIuuBXqp60lcGERmFm0uctm3bnv3ll19GIOzksWsXnHEGbNvmPhjnz4eePV0trQMHoGlTd4sno5ZDYp9+2vVASwW//rVL0NVl42AqlqjjYCqzfv163nzzTQ4ePBh0KHGlYcOG3HnnnWHvH7FxMKq6B5iBmwe8PbDEJ5fWwKciclJn83iZFS9eNW4MU6fCvfe6+lnf+Ia7ihk0CK69Fi66qPbJBeD733djQVLB1q1BR2ASQX5+PnfccQeDBg2iUaNGQYcTN/bu3cuWLVsi0laVH10i0gw4qqp7RKQOMBB4WFVPDdlnPRVcwZiqde8ODz0U3XOIuIKbU6dG9zzxwJ5FmXDVqVOH8847jz59+nD48OHjVZrXrFnD4sWLgw4vEKrKSy+9xE9/+tNad+sO57txS+B5EUnHXfGMVtV3a3VWE4jQahlpaRDDmVNjqpxnu8ZUSkSoW7cudevWJS8vj86dO/PFF19QVFQUdGiBOHToEDt37jypenV1hdOL7DNV/YaqnuXnFD+pdKCq5tvVS/xr0sQV2Rw82PVMe/99uOCCoKOKvNJCocbUVHZ2NmeeeWbQYQRq3759tW7DhrWmkOuugy+/hJdfdt13L73UdSYYPjzoyCLLqrabSOjVK7X7e3z++ee1bsMSTIpJT3cdC0rl5LiE88tfBhdTpLVrF3QEJhk0adIk6BACdfTo0Vq3YQnGkJbm5k75+99dAkp0lmBMJHzyySdBhxCovXv3ljtYtToswZjjRo2Cd95xE5wlqssuc+OKjKmtBg0akJmZSfPmzenbty9ZWVlBhxRTBQUFjB49mv3799e4DZtwLMkVFblR/Zdf7uaZCceiRXDllRChrvAx06ABLF3qZgqtCRtoWbFU/be7d+9eVJWGDRtSWFjIq6++yoYNqVV6sUGDBnz3u9+lYcOG5W63CcdSWFYWNG8O3/2uG8m/bl3Vx3zjGzBvHnTrFv34IqVXL/jDH2qeXIwpT25u7vEP1pycHK6//npOSaW6S7jeZFNrOIDOEkwKOP98uOUWVy7m2mshnDmF2rZ1Ze4T5XZTnz6uWoEx0XTKKaek5LwyNe3wkFp/pRR28cUuwVx6KYQ7M2qTJm6sTCJMB3DWWUFHYBLdrl27wuo51aZNGy677LIYRBS83Nxchg8fTv8aTqgUgSpXJhHUr1+zb/h5eW6szAUXnFgJIN7ccw/cemty9IIzwcjMzGTixIl0796dDh06VLrvOeecQ0FBAUuXLo1RdLEjInTq1Imzzz6b0047rVZt2RWMqVKnTvC3vwUdReV27nQdGhKBiPQUkXkislhEFohI75Bt94jIahFZKSKXhaw/W0SW+m1PiK81LyLZIvK6X/+xiOSHHDNSRFb5ZWTI+vZ+31X+2NTqHlWB+vXrc/nll/Piiy8yYcKESvcVEYYOHUrz5s1jFF3s9OrVi5tuuqnWyQUswZgwDRsG8TzbbJs2bl6cBPEIcL+q9gR+639HRM4AhgPdgMHAU74GIMBfcdNedPbLYL/+NmC3qnbCTQD4sG+rMXAfcC7QG7hPREpLBj8MPKaqnYHdvg2DSxyZmZnMnz+frVWU5c7MzOSGG24gJycnRtFFX8eOHekbwSlgLcGYsIjAz38edBQV698/MtMaxIgCDfzrXGCzfz0MeE1Vj6jqOmA10FtEWgINVHWuunEFLwBXhxzzvH/9JnCJv7q5DDc54C5V3Q1MwU13LsAAvy/+2NK2Ul5GRgZ33303V199dViFLhs3bsw111wTg8iir0WLFtx8880VdkeuicT5J2lqpLDQLZH4f2bAgPitwvzKK9CvH5x2GnTs6Lpmx7GfAJNF5A+4L3mlXxnzgHkh+xX4dUf967LrS4/ZCKCqxSKyF2gSur7MMU2APapaXE5bJygzWWB132PCysrKokePHmHv36VLF84//3zmzJkTxaiir2HDhhGf5dOuYJJcYSE8/nhk2mrY0I03iUfFxe4h//nnuyUChWBrRUSmisjn5SzDgNuB/1LVNsB/Ac+WHlZOU1rJ+pocU1lbJ660yQLD1q9fPxqHFvlLQJ06dYp4m5ZgklzDhm6E+9//Hpn2LrkkMu1E05o1cOhQsDGo6kA/vUXZ5W1gJDDW7/oG7hkJuKuJNiHNtMbdPivwr8uuP+EYEcnA3XLbVUlbO4CGft+ybZkaysjIYOjQoUGHUSvZ2dkRb9MSTAq44w6XYN56q/ZtVdF7M25s3x50BJXaDFzsXw8AVvnX44HhvmdYe9zD/E9UdQuwX0T6+GcoI4C3Q44p7SF2HfCBf04zGRgkIo38w/1BwGS/bbrfF39saVspae/evZRE4L5vfn4+PRN4rohoVCiwZzApIDvbJZcrroAePaB9+5q3NW1axMKKmnPOgaZNg46iUt8H/uSvIgrxzzlUdZmIjAa+AIqBO1T1mD/mduA5oA4w0S/gbq+9KCKrcVcuw31bu0TkQWC+3+8BVd3lX/8CeE1EHgIW8fUtupS0detW9u7dG5HnTBdccEFCTrXcoUOHqDxnswSTItq2hQjMHxT3pWOys+G556Bly6AjqZiqzgbOrmDb74DflbN+AdC9nPWFwPUVtPVP4J/lrF/L17flUl6HDh0YPXo0N910U60fcjdp0oSmTZvWusx9rLVt2zYq5W8swZhq6d7dVQR4+eVgn3PUqwctWrjeYi1auB5kTZq4uWDiPQma+JKZmUl2djZLly7lrAjUHOrYsWNCJZi0tLSIDKosjyUYUy0DB7pqy+3bw69+Fdm2c3JOTBqhr0N/Nm8OKVbQ1kRZv379ePHFF+nSpUutB04mUrXlzMxMbrnlFlq0aBGV9i3BmGqpX999uM+eHf4xjRpBq1aurlmrVm5p2fLk5FG/vhvQaUysNW3alCuuuIL9+/fXOsEkwsj+jIwMGjZsSP/+/aM6xskSjKk2ETegcfr0r5NGaPIIfd2qVfgTnRkTpEjdJqqTAP/D16tXjzvuuCPq57EEY2rk0Ufh//7PrjiMKSsRplbuFqPZBG0cjKmR9HRLLsaUZ+bMmUGHUKXly5fjhkRFlyUYY4yJoH1B1ykKQ3p6esTrjpXHEowxxkSQqpIe5zPftWvXLibnsQRjjDERdN1110V0TpVoGDhwYEzOYwnGGGMiqF27duTm5gYdRqWWLFkSk/NYLzJjjImwFStWBB1CpSZNmsSCBQto06YNQ4cOjdrzGLuCMcaYCGvVqhXNmjWL60GXO3bsYNGiRWzYsCFq57ArGGOMibD+/fvTv39/AI4ePcr+/fsrXfbt20dxcXEVrUZHo0aNotZ2TBPMwoULd4jIl7E8ZxWa4iZgiicWU3iiEVNsutaYlJKZmUnjxo0rnfFSVTly5Aj79u07IfEcOHCAAwcOnPB7pBNRNAeGxjTBqGpczbsqIgtUNa4mAbaYwhOPMRlTUyJCTk4OOTk5nHrqqRXup6oUFhaelHTK+3n06NFKz5mbm8sNN9wQ1dt4dovMGGMShIhQp04d6tSpQ7NmFX9fV1WKiopOSDj79+8HXKHLL7/8kr59+9KqVauoxmsJxhhjkoyIkJ2dTXZ2Nk3Lmd61d+/YzDeX6r3Ing46gHJYTOGJx5iMMSFSOsGoatx9SFlM4YnHmIwxJ0rpBGOMMSZ6LMEYY4yJiqRLMCJyvYgsE5ESEekVsv5SEVkoIkv9zwF+fV0ReU9EVvjjfl9Bu/kiclhEFvvlb9GKyW87269fLSJPSAW1HETkHr/PShG5LNyYqoiriYhMF5EDIvJkyPr6Ie9/sYjsEJHHy2k3Gn+rcmPy22b49196vnL7edbmb2WMqb5k7EX2OXAt8Pcy63cAQ1V1s4h0ByYDeX7bH1R1uohkAdNE5HJVnVhO22tUtWeMYvorMAqYB0wABgMnxCQiZwDDgW5AK2CqiHRR1WO1jKsQ+A3Q3S8AqOp+oGfI+RcCYytoO9J/q3JjCnGzqi6oqNEI/K2MMdWUdFcwqrpcVVeWs36Rqm72vy4DckQkW1UPqep0v08R8CnQOsiYRKQl0EBV56qbdu4F4Opymh4GvKaqR1R1HbAaCLv/YSVxHVTV2bgP9XKJSGfgVODDcM8X7ZiqUKu/lTGm+pIuwYTpW8AiVT0SulJEGgJDgWkVHNdeRBaJyEwRuTCKMeUBBSHbCvj6yiZUHrAxjP2i4dvA61rxvKvR/FuV51/+9thvKridGOTfypiUlJC3yERkKtCinE33qurbVRzbDXgYGFRmfQbwKvCEqq4t59AtQFtV3SkiZwNviUg3Vd0XhZjK+4As74O8yv1qE1cVhgP/UcG2qP2tKnCzqm4SkfrAGB/XC2X2CfdvaoyJkIRMMKpao+nYRKQ1MA4Yoaprymx+Glilqo9XcM4jwBH/eqGIrAG6AAuiEFMBJ96maw1sLnus369NZfvVNK7KiEgPIENVF5a3PVp/q4qo6ib/c7+IvIK79VU2wVT5tzIVi5NCtfFUdDVeYomHOCosEpuQCaYm/O2v94B7VHVOmW0PAbnA9yo5vhmwS1WPiUgHoDNQ3pVOrWNS1S0isl9E+gAfAyOAP5fTxHjgFRH5I+7BdWfgk9rEFKZv4672yhWNv1Ul58oAGqrqDhHJBIYAU8vZNai/VVKIh0K1EkcFTuMllniJo0KqmlQLcA3u2+oRYBsw2a//NXAQWByynIr7JqvA8pD13/PHXAU84F9/C/cgfgmuI8DQaMXkt/XC9ahaAzwJSNmY/O/3+n1WApdH4m/lt60HdgEH/D5nhGxbC5xepq2o/q0qigmoBywEPvPn/BOQHum/lS3BL8CCoGOIt1jiJY6KltIPLWOMiWvx9G09XmKJlzgqkqq9yIwxiSee6s/FSyzxEke57ArGGGNMVNgVjDHGmKiwBGOMMSYqLMEYY+JSRYVPy9nvv/x+n4vIqyIS0UnmqxFHQxF5U1zh3OUicl4k46hOLH7fdF9N491IxxEuSzDGmHhVWvh0VkU7iEge8GOgl6p2B9JxVSZiGof3J2CSqp4O9MANfYi0cGMBuDNKMYTNEowxJi5pBYVPy5EB1PGDbusS4QoN4cQhIg2Ai4Bn/TFFqronknGEG4uPpzVwJfBMpGOoDkswxpiEpa5M0B+ADbgaeHtV9f0AQukAbMcVXV0kIs+ISL0A4ij1OPBzoCTAGCzBGGOCIyJT/bOTssuwMI9vhJuKoT2uBFA9Ebkl1nHgrqK+CfxVVb+Bq9Dxy+rGEYlYRGQI8JVWUCswllKmFpkxJv5o7QufDgTWqep2ABEZC/QFXopxHAVAgap+7H9/kxommAjEcj5wlYhcAeQADUTkJVWtduKtLbuCMcYksg1AH3FTnwtwCQE82FbVrcBGETnNr7oE+CLWcfhY7lHV1qqaj+vw8EEQyQUswRhj4pSIXCMiBcB5wHsiMtmvbyUiEwD8FcObuKKqS3GfaREtnxJOHN6PgJdF5DPc1OL/E8k4qhlLXLBSMcYYY6LCrmCMMcZEhSUYY4wxUWEJxhhjTFRYgjHGGBMVlmCMMSYFicg/ReQrEfk8zP1vEJEvfLHNV8I6xnqRGWNM6hGRi4ADwAu+UGhl+3YGRgMDVHW3iJyqql9VdQ67gjHGmBSkqrOAXaHrRKSjiEwSkYUi8qGInO43fR/4i6ru9sdWmVzAEowxxpivPQ38SFXPBu4GnvLruwBdRGSOiMwTkcHhNGa1yIwxxiAip+DquL3hqu4AkO1/ZgCdgX5Aa+BDEele1ZQElmCMMcaAu6O1R1V7lrOtAJinqkeBdSKyEpdw5lfVoDHGmBSnqvtwyeN6AHF6+M1vAf39+qa4W2Zrq2rTEowxxqQgEXkVmAucJiIFInIbcDNwm4gsAZbh5toBmAzsFJEvgOnAz1R1Z5XnsG7KxhhjosGuYIwxxkSFJRhjjDFRYQnGGGNMVFiCMcYYExWWYIwxxkSFJRhjjDFRYQnGGGNMVPw/bVzeR7MCPiUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the CA polygon in EPSG:4269 and the Lambert Conformal Conic projection\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(6, 4))\n",
    "ca_poly.plot(ax=ax1, facecolor='blue');\n",
    "ax1.set_title(\"CA: EPSG=4269\");\n",
    "ca_poly_lcc_2.plot(ax=ax2, facecolor='grey');\n",
    "ax2.set_title(\"CA: Daymet LCC\");\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71d3f3c1",
   "metadata": {},
   "source": [
    "## 1.2 Temporal: Search and Subset Parameters\n",
    "### Set and format a time range of interest to discover and subset data\n",
    "\n",
    "For this example, we'll access 1 year of Daymet daily data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5c8ad184",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2010-01-01 00:00:00\n",
      "2011-12-31 00:00:00\n",
      "2010-01-01T00:00:00Z,2011-12-31T00:00:00Z\n"
     ]
    }
   ],
   "source": [
    "import datetime as dt \n",
    "\n",
    "start_date = dt.datetime(2010, 1, 1) # specify your own start date\n",
    "end_date = dt.datetime(2011, 12, 31)  # specify your end start date\n",
    "print(start_date,end_date,sep='\\n')\n",
    "\n",
    "dt_format = '%Y-%m-%dT%H:%M:%SZ' # format requirement for datetime search\n",
    "temporal_str = start_date.strftime(dt_format) + ',' + end_date.strftime(dt_format)\n",
    "print(temporal_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "342414b7",
   "metadata": {},
   "source": [
    "----\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6f69c02",
   "metadata": {},
   "source": [
    "## STEP 2. Finding Data Files of Interest : NASA's CRM API\n",
    "\n",
    "If you know the OPeNDAP URL's for your data of interest, you can skip this step.\n",
    "\n",
    "In this step, we'll leverage NASA's Common Metadata Repository (CMR) API to search NASA Earthdata for Daymet data that are within a spatial and temporal region of interest.  \n",
    "\n",
    "If you search [NASA's Earthdata Search](https://search.earthdata.nasa.gov/search) for [Daymet Daily data](https://search.earthdata.nasa.gov/search?q=daymet%20daily), you'll find about 1,000 files associated with the most recent version of that dataset.  \n",
    "\n",
    "In STEP 1, we defined a spatial and temporal area-of-interest. We'll see how we can programatically use those as a search parameters to discover which specific Daymet daily data files are within NASA Earthdata holdings that match those parameters.  \n",
    "\n",
    "**At the end of this step, we will have a list of Daymet files that are within our search parameters (spatial and temporal)**\n",
    "\n",
    "* **This methodology can be applied to any NASA Earthdata Dataset**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bda5405",
   "metadata": {},
   "source": [
    "A little backgroud ....\n",
    "\n",
    "- Searching [NASA's Earthdata Holdings](https://search.earthdata.nasa.gov/search) can lead to quite a few files for a user to consider. \n",
    "![title](images/NASAEarthdataCollectionSearch_narrow_LI.jpg.PNG)\n",
    "<br>\n",
    "\n",
    "Some helpful links:\n",
    "- NASA has a metadata system that catalogs all data and metadata records:  [Common Metadata Repository (CMR)](https://earthdata.nasa.gov/eosdis/science-system-description/eosdis-components/cmr) \n",
    "\n",
    "- Using the [CMR API](https://cmr.earthdata.nasa.gov/search), we can search metadata records to discovered data granules (files).\n",
    "\n",
    "\n",
    "**`Request URL's`**\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "https://cmr.earthdata.nasa.gov/search/collections\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "https://cmr.earthdata.nasa.gov/search/granules\n",
    "</div>\n",
    "\n",
    "- The [Daymet V4R1 Daily](https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=2129) data DOI, **`doi.org/10.3334/ORNLDAAC/2129`**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a382f912",
   "metadata": {},
   "source": [
    "## Finding Data Files in NASA's EarthData Cloud: or What's a ConceptID\n",
    "\n",
    "NASA Earthdata holds thousands of datasets and finding the files that fit within a spatial area of interest and time frame can be quite daunting.  \n",
    "* We'll demonstrate how to use NASA's [**CMR API**](https://cmr.earthdata.nasa.gov/search) and either the Daymet dataset `DOI` or the `shortname` to obtain the NASA EarthData's unique ID (called the **`Concept ID`**) for **Daymet Daily** data.  With this information, we can search all of NASA's Earthdata holdings for matching data files (granules) based on our spatial area and time period of interest.\n",
    "\n",
    "Search NASA Earthdata Holdings based on...\n",
    "* *`DOI`* - The Daymet daily data Digital Object Identifier\n",
    "* *`ConceptID`* - remember, this is NASA Earthdata's unique ID for this dataset\n",
    "* Area of Interest\n",
    "* Time Range of Interest\n",
    "\n",
    "## Step 2.1 Obtaining the Daymet V4 Daily data *Concept ID*\n",
    "Using Python's **`requests`** library, we will build and provide the URL `doisearch` to search **NASA's Common Metadata Repository (CMR)** to acquire the **`concept_id`** unique to the Daymet V4 Daily data in NASA's Cloud holdings.  We'll use the **concept_id** in a next step when we search NASA Cloud Holdings for our spatial and temporal area of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f17f82a6",
   "metadata": {},
   "source": [
    "*  We're going to build a **CMR URL** for the Daymet V4R1 Daily Data\n",
    "* **`requests`** is a simple HTTP libray for Python.  It allows you to easily send HTTP requests\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "https://cmr.earthdata.nasa.gov/search/collections\n",
    "</div>\n",
    "\n",
    "* Basically, we're going to programmatically build the CMR URL for the Daymet Daily data and search for and store the ConceptID\n",
    "`https://cmr.earthdata.nasa.gov/search/collections.json?doi=10.3334/ORNLDAAC/2129`\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "cdc053d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://cmr.earthdata.nasa.gov/search/collections.json?doi=10.3334/ORNLDAAC/2129\n",
      "C2532426483-ORNL_CLOUD\n"
     ]
    }
   ],
   "source": [
    "# CMR API base url\n",
    "cmrurl='https://cmr.earthdata.nasa.gov/search/' # define the base url of NASA's CMR API as the variable `cmrurl`\n",
    "doi = '10.3334/ORNLDAAC/2129'                   # Daymet V4R1 DOI - define the DOI as the variable `daymet_doi`\n",
    "doisearch = cmrurl + 'collections.json?doi=' + doi\n",
    "print(doisearch)\n",
    "\n",
    "import requests\n",
    "concept_id = requests.get(doisearch).json()['feed']['entry'][0]['id']\n",
    "print(concept_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2086317c",
   "metadata": {},
   "source": [
    "**`C2532426483-ORNL_CLOUD`** is the unique NASA-given concept ID for the Daymet daily dataset. We can use this to search relevant granules and their names. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bfd98d0",
   "metadata": {},
   "source": [
    "- Another way to get to the Daymet ConceptID from CMR's Collection's URL `https://cmr.earthdata.nasa.gov/search/collections.json` and knowing the shortname `Daymet_Daily_V4R1_2129`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4824d4c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Response [200]>\n",
      "['C2532426483-ORNL_CLOUD']\n",
      "C2532426483-ORNL_CLOUD\n"
     ]
    }
   ],
   "source": [
    "params = {\n",
    "    'short_name': 'Daymet_Daily_V4R1_2129',\n",
    "} # parameter dictionary with known CMR short_name\n",
    "\n",
    "cmr_collections_url = 'https://cmr.earthdata.nasa.gov/search/collections.json'\n",
    "cmr_response = requests.get(cmr_collections_url, params=params)\n",
    "print(cmr_response)\n",
    "cmr_results = json.loads(cmr_response.content) # Get json response from CMR collection metadata\n",
    "#print(cmr_results)\n",
    "collectionlist = [el['id'] for el in cmr_results['feed']['entry']]\n",
    "print(collectionlist)\n",
    "concept_id = collectionlist[0]\n",
    "print(concept_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe9d6cab",
   "metadata": {},
   "source": [
    "## Step 2.2 Searching for Daymet Files \n",
    "Daymet daily data files (or granules) are in netCDF4 format, and each file has one year's worth of data. Data files are organized by variables (each for dayl, prcp, tmin, tmax, srad, swe, vp) and regions (each for us, pr, hi).  Daymet filenames can be used to identify the files from continental North America (`*_na_*.nc`). The files from Puerto Rico and Hawaii are named as (`*_pr_*.nc`) and (`*_hi_*.nc`) respectively.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "https://cmr.earthdata.nasa.gov/search/granules\n",
    "</div>\n",
    "\n",
    "We will first search all the granules for our area (i.e., California) and time period of interest (2010, 2011).  For this tutorial, we set up the search parameters and download maximum temperature data (`tmax`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "50d50640",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Variables we have defined to pass to the CMR API\n",
      "-----------------------------------------------\n",
      "time range :      2010-01-01T00:00:00Z,2011-12-31T00:00:00Z\n",
      "Daymet variable : tmax\n",
      "bounding box :    [-124.39263575019409, 32.53578253612507, -114.1252281747415, 42.002192463378705]\n",
      "concept ID :      C2532426483-ORNL_CLOUD\n",
      "cmr url :         https://cmr.earthdata.nasa.gov/search/\n"
     ]
    }
   ],
   "source": [
    "#var = 'prcp' # select your own variable\n",
    "var = 'tmax' # select your own variable\n",
    "\n",
    "print('Variables we have defined to pass to the CMR API')\n",
    "print('-----------------------------------------------')\n",
    "print('time range :     ', temporal_str)\n",
    "print('Daymet variable :', var)\n",
    "print('bounding box :   ', xy_wgs84)\n",
    "print('concept ID :     ', concept_id)\n",
    "print('cmr url :        ', cmrurl)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "180f470b",
   "metadata": {},
   "source": [
    "### Here, we're building the CMR URL for the files (granules) that match our parameters of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "44f2206f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "https://cmr.earthdata.nasa.gov/search/granules.json?collection_concept_id=C2532426483-ORNL_CLOUD&page_size=1000&temporal=2010-01-01T00:00:00Z,2011-12-31T00:00:00Z&bounding_box[]=-124.39263575019409,32.53578253612507,-114.1252281747415,42.002192463378705\n"
     ]
    }
   ],
   "source": [
    "granulesearch = cmrurl + 'granules.json?collection_concept_id=' + concept_id + \\\n",
    "                '&page_size=1000' + '&temporal=' + temporal_str + \\\n",
    "                '&bounding_box[]=' + ','.join(map(str, xy_wgs84))\n",
    "print(granulesearch)                "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a02f914",
   "metadata": {},
   "source": [
    "### Searching NASA's Earthdata holdings and getting a list of files that match our parameters\n",
    "Now, we'll get a list of files using Python's **`requests`** library. We can provide the URL created above, `granulesearch`, created above to search **NASA's Earthdata Metadata** holdings for Daymet data that fit our spatial, temporal, and variable request."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "cb30dcca",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2010.nc\n",
      "Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2011.nc\n"
     ]
    }
   ],
   "source": [
    "response = requests.get(granulesearch)\n",
    "granules = response.json()['feed']['entry']  \n",
    "granule_names = []\n",
    "for g in granules:\n",
    "    granule_name = g['title']\n",
    "    if var in granule_name:\n",
    "        granule_names.append(granule_name)\n",
    "        print(granule_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6599ba19",
   "metadata": {},
   "source": [
    "The array `granule_names` contains the names of all the granules (files) that match our search parameters.\n",
    "\n",
    "**Note:** A typical Daymet V4 daily file has the name **daymet_v4_daily_`area`_`var`_`year`.nc**.  The beginning of the filename (**\"Daymet_Daily_V4R1.\"**) is specific to the cloud-based naming convention and because we are searching CMR, the cloud-based files will be discovered.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f08124e",
   "metadata": {},
   "source": [
    "## STEP 3. Subsetting data from NASA Earthdata Cloud\n",
    "\n",
    "Using the granule list we programmatically obtained above, we'll subset the North American files to our area of interest using the bounding box from the California polygon."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f31d08e",
   "metadata": {},
   "source": [
    "Daymet data is available through NASA's EarthData Cloud *`Hyrax Data Server`*.\n",
    "\n",
    "For example ...\n",
    "\n",
    "[Daymet V4R1 Daily Data in the NASA Earthdata Cloud Hyrax Data Server](https://opendap.earthdata.nasa.gov/collections/C2532426483-ORNL_CLOUD/granules/Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2010.nc.dmr.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95ea4c73",
   "metadata": {},
   "source": [
    "### **`OPeNDAP`** allows users to use data files that are stored on remote computers with their favorite analysis and visualization client software.\n",
    "\n",
    "NASA Earthdata has great informaion on OPeNDAP: [What is OPeNDAP](https://earthdata.nasa.gov/collaborate/open-data-services-and-software/api/opendap)\n",
    "\n",
    "###  We'll use Pythons **`pydap`** library\n",
    "\n",
    "* **[pydap](https://www.pydap.org/en/latest/client.html)** also **[pydap](https://github.com/pydap/pydap)** can be used to inspect data from any scientific dataset available on the internet from OPeNDAP servers. This allows inspection and manipulation of scientific data as if it were stored locally.  Data can be downloaded on-the-fly as needed.\n",
    "\n",
    "###  We'll use Python's **`Xarray`** Library to \n",
    "* open the data files we are searching for, \n",
    "* subset the North American file based on the Lambert Conformal Conic **`lccbounds`** bounding box we defined earlier, \n",
    "* concatenate the time range of files into one file, and \n",
    "* save the concatened file to a netCDF file.\n",
    "\n",
    "Pydap\n",
    "* **open_url**\n",
    "\n",
    "Xarray\n",
    "* **[open_dataset](http://xarray.pydata.org/en/stable/generated/xarray.open_dataset.html)** \n",
    "* **slice**\n",
    "* **concat**\n",
    "* **to_netcdf**\n",
    "\n",
    "We will subset and download Daymet V4 Daily Data based on\n",
    "- **`granual_names`** array we saved above \n",
    "- our polygonal region of interest (LCC bounding box)\n",
    "- our time range of interest\n",
    "- Daymet variable\n",
    "\n",
    "###  A NASA Earthdata Login is required\n",
    "Get an account here https://urs.earthdata.nasa.gov/\n",
    "\n",
    "- We'll use the Central Authentication Service (CAS) as a single sign-on protocol for the web.\n",
    "- NASA Earthdata Authenticaion is done through `username` and `login`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "57c5f29e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "granule_name ----> Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2010.nc\n",
      "https://opendap.earthdata.nasa.gov/collections/C2532426483-ORNL_CLOUD/granules/Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2010.nc\n",
      "hyrax_ds -----> <DatasetType with children 'y', 'lon', 'lambert_conformal_conic', 'time', 'yearday', 'tmax', 'lat', 'time_bnds', 'x'>\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:397: RuntimeWarning: overflow encountered in int_scalars\n",
      "  count = response_dtype.itemsize * n\n",
      "C:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:397: RuntimeWarning: overflow encountered in int_scalars\n",
      "  count = response_dtype.itemsize * n\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:                  (y: 8075, x: 7814, time: 365, nv: 2)\n",
      "Coordinates:\n",
      "  * y                        (y) float32 4.984e+06 4.983e+06 ... -3.09e+06\n",
      "    lon                      (y, x) float32 ...\n",
      "    lambert_conformal_conic  int16 ...\n",
      "  * time                     (time) datetime64[ns] 2010-01-01T12:00:00 ... 20...\n",
      "    lat                      (y, x) float32 ...\n",
      "    time_bnds                (time, nv) datetime64[ns] ...\n",
      "  * x                        (x) float32 -4.56e+06 -4.559e+06 ... 3.253e+06\n",
      "Dimensions without coordinates: nv\n",
      "Data variables:\n",
      "    yearday                  (time) int16 ...\n",
      "    tmax                     (time, y, x) float32 ...\n",
      "Attributes:\n",
      "    build_dmrpp_metadata.build_dmrpp:    3.20.10-462\n",
      "    build_dmrpp_metadata.bes:            3.20.10-462\n",
      "    build_dmrpp_metadata.libdap:         libdap-3.20.9-84\n",
      "    build_dmrpp_metadata.configuration:  \\n# TheBESKeys::get_as_config()\\nAll...\n",
      "    build_dmrpp_metadata.invocation:     build_dmrpp -c /tmp/bes_conf_NSyi -f...\n",
      "    start_year:                          2010\n",
      "    source:                              Daymet Software Version 4.0\n",
      "    Version_software:                    Daymet Software Version 4.0\n",
      "    Version_data:                        Daymet Data Version 4.0\n",
      "    Conventions:                         CF-1.6\n",
      "    citation:                            Please see http://daymet.ornl.gov/ f...\n",
      "    references:                          Please see http://daymet.ornl.gov/ f...\n",
      "tmax_h ====\n",
      "<xarray.DataArray 'tmax' (time: 365, y: 47, x: 83)>\n",
      "[1423865 values with dtype=float32]\n",
      "Coordinates:\n",
      "  * y                        (y) float32 -5.66e+05 -5.67e+05 ... -6.12e+05\n",
      "    lon                      (y, x) float32 ...\n",
      "    lambert_conformal_conic  int16 ...\n",
      "  * time                     (time) datetime64[ns] 2010-01-01T12:00:00 ... 20...\n",
      "    lat                      (y, x) float32 ...\n",
      "  * x                        (x) float32 1.388e+06 1.389e+06 ... 1.47e+06\n",
      "Attributes:\n",
      "    long_name:     daily maximum temperature\n",
      "    units:         degrees C\n",
      "    cell_methods:  area: mean time: maximum\n",
      "granule_name ----> Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2011.nc\n",
      "https://opendap.earthdata.nasa.gov/collections/C2532426483-ORNL_CLOUD/granules/Daymet_Daily_V4R1.daymet_v4_daily_na_tmax_2011.nc\n",
      "hyrax_ds -----> <DatasetType with children 'y', 'lon', 'lambert_conformal_conic', 'time', 'yearday', 'tmax', 'lat', 'time_bnds', 'x'>\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:397: RuntimeWarning: overflow encountered in int_scalars\n",
      "  count = response_dtype.itemsize * n\n",
      "C:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:397: RuntimeWarning: overflow encountered in int_scalars\n",
      "  count = response_dtype.itemsize * n\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:                  (y: 8075, x: 7814, time: 365, nv: 2)\n",
      "Coordinates:\n",
      "  * y                        (y) float32 4.984e+06 4.983e+06 ... -3.09e+06\n",
      "    lon                      (y, x) float32 ...\n",
      "    lambert_conformal_conic  int16 ...\n",
      "  * time                     (time) datetime64[ns] 2011-01-01T12:00:00 ... 20...\n",
      "    lat                      (y, x) float32 ...\n",
      "    time_bnds                (time, nv) datetime64[ns] ...\n",
      "  * x                        (x) float32 -4.56e+06 -4.559e+06 ... 3.253e+06\n",
      "Dimensions without coordinates: nv\n",
      "Data variables:\n",
      "    yearday                  (time) int16 ...\n",
      "    tmax                     (time, y, x) float32 ...\n",
      "Attributes:\n",
      "    build_dmrpp_metadata.build_dmrpp:    3.20.10-462\n",
      "    build_dmrpp_metadata.bes:            3.20.10-462\n",
      "    build_dmrpp_metadata.libdap:         libdap-3.20.9-84\n",
      "    build_dmrpp_metadata.configuration:  \\n# TheBESKeys::get_as_config()\\nAll...\n",
      "    build_dmrpp_metadata.invocation:     build_dmrpp -c /tmp/bes_conf_BPjs -f...\n",
      "    start_year:                          2011\n",
      "    source:                              Daymet Software Version 4.0\n",
      "    Version_software:                    Daymet Software Version 4.0\n",
      "    Version_data:                        Daymet Data Version 4.0\n",
      "    Conventions:                         CF-1.6\n",
      "    citation:                            Please see http://daymet.ornl.gov/ f...\n",
      "    references:                          Please see http://daymet.ornl.gov/ f...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:397: RuntimeWarning: overflow encountered in int_scalars\n",
      "  count = response_dtype.itemsize * n\n",
      "C:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:397: RuntimeWarning: overflow encountered in int_scalars\n",
      "  count = response_dtype.itemsize * n\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "buffer size must be a multiple of element size",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Input \u001b[1;32mIn [19]\u001b[0m, in \u001b[0;36m<cell line: 18>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     37\u001b[0m         \u001b[38;5;28mprint\u001b[39m(tmax_h)\n\u001b[0;32m     38\u001b[0m         \u001b[38;5;66;03m#tmax_h.to_netcdf('test.nc')\u001b[39;00m\n\u001b[0;32m     39\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m     40\u001b[0m         \u001b[38;5;66;03m#print('elsey')\u001b[39;00m\n\u001b[1;32m---> 41\u001b[0m         tmax_h \u001b[38;5;241m=\u001b[39m \u001b[43mxr\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mconcat\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[43mtmax_h\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtmax_sel\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdim\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mtime\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m     43\u001b[0m     cnt \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[0;32m     45\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mconcat tmax_h ----->\u001b[39m\u001b[38;5;124m'\u001b[39m, tmax_h)    \n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\concat.py:238\u001b[0m, in \u001b[0;36mconcat\u001b[1;34m(objs, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs)\u001b[0m\n\u001b[0;32m    233\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m    234\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTypeError\u001b[39;00m(\n\u001b[0;32m    235\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcan only concatenate xarray Dataset and DataArray \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m    236\u001b[0m         \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mobjects, got \u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mtype\u001b[39m(first_obj)\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m    237\u001b[0m     )\n\u001b[1;32m--> 238\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mf\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m    239\u001b[0m \u001b[43m    \u001b[49m\u001b[43mobjs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdata_vars\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcoords\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcompat\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpositions\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mfill_value\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mjoin\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcombine_attrs\u001b[49m\n\u001b[0;32m    240\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\concat.py:578\u001b[0m, in \u001b[0;36m_dataarray_concat\u001b[1;34m(arrays, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs)\u001b[0m\n\u001b[0;32m    575\u001b[0m             arr \u001b[38;5;241m=\u001b[39m arr\u001b[38;5;241m.\u001b[39mrename(name)\n\u001b[0;32m    576\u001b[0m     datasets\u001b[38;5;241m.\u001b[39mappend(arr\u001b[38;5;241m.\u001b[39m_to_temp_dataset())\n\u001b[1;32m--> 578\u001b[0m ds \u001b[38;5;241m=\u001b[39m \u001b[43m_dataset_concat\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m    579\u001b[0m \u001b[43m    \u001b[49m\u001b[43mdatasets\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    580\u001b[0m \u001b[43m    \u001b[49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    581\u001b[0m \u001b[43m    \u001b[49m\u001b[43mdata_vars\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    582\u001b[0m \u001b[43m    \u001b[49m\u001b[43mcoords\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    583\u001b[0m \u001b[43m    \u001b[49m\u001b[43mcompat\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    584\u001b[0m \u001b[43m    \u001b[49m\u001b[43mpositions\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    585\u001b[0m \u001b[43m    \u001b[49m\u001b[43mfill_value\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mfill_value\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    586\u001b[0m \u001b[43m    \u001b[49m\u001b[43mjoin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mjoin\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    587\u001b[0m \u001b[43m    \u001b[49m\u001b[43mcombine_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcombine_attrs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    588\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    590\u001b[0m merged_attrs \u001b[38;5;241m=\u001b[39m merge_attrs([da\u001b[38;5;241m.\u001b[39mattrs \u001b[38;5;28;01mfor\u001b[39;00m da \u001b[38;5;129;01min\u001b[39;00m arrays], combine_attrs)\n\u001b[0;32m    592\u001b[0m result \u001b[38;5;241m=\u001b[39m arrays[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39m_from_temp_dataset(ds, name)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\concat.py:517\u001b[0m, in \u001b[0;36m_dataset_concat\u001b[1;34m(datasets, dim, data_vars, coords, compat, positions, fill_value, join, combine_attrs)\u001b[0m\n\u001b[0;32m    515\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mKeyError\u001b[39;00m:\n\u001b[0;32m    516\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mk\u001b[38;5;132;01m!r}\u001b[39;00m\u001b[38;5;124m is not present in all datasets.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m--> 517\u001b[0m combined \u001b[38;5;241m=\u001b[39m \u001b[43mconcat_vars\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mvars\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpositions\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcombine_attrs\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcombine_attrs\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    518\u001b[0m \u001b[38;5;28;01massert\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(combined, Variable)\n\u001b[0;32m    519\u001b[0m result_vars[k] \u001b[38;5;241m=\u001b[39m combined\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\variable.py:2999\u001b[0m, in \u001b[0;36mconcat\u001b[1;34m(variables, dim, positions, shortcut, combine_attrs)\u001b[0m\n\u001b[0;32m   2997\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m IndexVariable\u001b[38;5;241m.\u001b[39mconcat(variables, dim, positions, shortcut, combine_attrs)\n\u001b[0;32m   2998\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m-> 2999\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mVariable\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mconcat\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvariables\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdim\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpositions\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mshortcut\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcombine_attrs\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\variable.py:1897\u001b[0m, in \u001b[0;36mVariable.concat\u001b[1;34m(cls, variables, dim, positions, shortcut, combine_attrs)\u001b[0m\n\u001b[0;32m   1894\u001b[0m variables \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(variables)\n\u001b[0;32m   1895\u001b[0m first_var \u001b[38;5;241m=\u001b[39m variables[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m-> 1897\u001b[0m arrays \u001b[38;5;241m=\u001b[39m [v\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;28;01mfor\u001b[39;00m v \u001b[38;5;129;01min\u001b[39;00m variables]\n\u001b[0;32m   1899\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m dim \u001b[38;5;129;01min\u001b[39;00m first_var\u001b[38;5;241m.\u001b[39mdims:\n\u001b[0;32m   1900\u001b[0m     axis \u001b[38;5;241m=\u001b[39m first_var\u001b[38;5;241m.\u001b[39mget_axis_num(dim)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\variable.py:1897\u001b[0m, in \u001b[0;36m<listcomp>\u001b[1;34m(.0)\u001b[0m\n\u001b[0;32m   1894\u001b[0m variables \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mlist\u001b[39m(variables)\n\u001b[0;32m   1895\u001b[0m first_var \u001b[38;5;241m=\u001b[39m variables[\u001b[38;5;241m0\u001b[39m]\n\u001b[1;32m-> 1897\u001b[0m arrays \u001b[38;5;241m=\u001b[39m [\u001b[43mv\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdata\u001b[49m \u001b[38;5;28;01mfor\u001b[39;00m v \u001b[38;5;129;01min\u001b[39;00m variables]\n\u001b[0;32m   1899\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m dim \u001b[38;5;129;01min\u001b[39;00m first_var\u001b[38;5;241m.\u001b[39mdims:\n\u001b[0;32m   1900\u001b[0m     axis \u001b[38;5;241m=\u001b[39m first_var\u001b[38;5;241m.\u001b[39mget_axis_num(dim)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\variable.py:339\u001b[0m, in \u001b[0;36mVariable.data\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    337\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_data\n\u001b[0;32m    338\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 339\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvalues\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\variable.py:512\u001b[0m, in \u001b[0;36mVariable.values\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    509\u001b[0m \u001b[38;5;129m@property\u001b[39m\n\u001b[0;32m    510\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mvalues\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m    511\u001b[0m     \u001b[38;5;124;03m\"\"\"The variable's data as a numpy.ndarray\"\"\"\u001b[39;00m\n\u001b[1;32m--> 512\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43m_as_array_or_item\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_data\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\variable.py:252\u001b[0m, in \u001b[0;36m_as_array_or_item\u001b[1;34m(data)\u001b[0m\n\u001b[0;32m    238\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_as_array_or_item\u001b[39m(data):\n\u001b[0;32m    239\u001b[0m     \u001b[38;5;124;03m\"\"\"Return the given values as a numpy array, or as an individual item if\u001b[39;00m\n\u001b[0;32m    240\u001b[0m \u001b[38;5;124;03m    it's a 0d datetime64 or timedelta64 array.\u001b[39;00m\n\u001b[0;32m    241\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    250\u001b[0m \u001b[38;5;124;03m    TODO: remove this (replace with np.asarray) once these issues are fixed\u001b[39;00m\n\u001b[0;32m    251\u001b[0m \u001b[38;5;124;03m    \"\"\"\u001b[39;00m\n\u001b[1;32m--> 252\u001b[0m     data \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    253\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m data\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m    254\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m data\u001b[38;5;241m.\u001b[39mdtype\u001b[38;5;241m.\u001b[39mkind \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mM\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\indexing.py:552\u001b[0m, in \u001b[0;36mMemoryCachedArray.__array__\u001b[1;34m(self, dtype)\u001b[0m\n\u001b[0;32m    551\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__array__\u001b[39m(\u001b[38;5;28mself\u001b[39m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m--> 552\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_ensure_cached\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    553\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39masarray(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray, dtype\u001b[38;5;241m=\u001b[39mdtype)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\indexing.py:549\u001b[0m, in \u001b[0;36mMemoryCachedArray._ensure_cached\u001b[1;34m(self)\u001b[0m\n\u001b[0;32m    547\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_ensure_cached\u001b[39m(\u001b[38;5;28mself\u001b[39m):\n\u001b[0;32m    548\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray, NumpyIndexingAdapter):\n\u001b[1;32m--> 549\u001b[0m         \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray \u001b[38;5;241m=\u001b[39m NumpyIndexingAdapter(\u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masarray\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m)\u001b[49m)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\indexing.py:522\u001b[0m, in \u001b[0;36mCopyOnWriteArray.__array__\u001b[1;34m(self, dtype)\u001b[0m\n\u001b[0;32m    521\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__array__\u001b[39m(\u001b[38;5;28mself\u001b[39m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m--> 522\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masarray\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\indexing.py:423\u001b[0m, in \u001b[0;36mLazilyIndexedArray.__array__\u001b[1;34m(self, dtype)\u001b[0m\n\u001b[0;32m    421\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__array__\u001b[39m(\u001b[38;5;28mself\u001b[39m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m    422\u001b[0m     array \u001b[38;5;241m=\u001b[39m as_indexable(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray)\n\u001b[1;32m--> 423\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39masarray(\u001b[43marray\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\conventions.py:62\u001b[0m, in \u001b[0;36mNativeEndiannessArray.__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m     61\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[1;32m---> 62\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masarray\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\coding\\variables.py:70\u001b[0m, in \u001b[0;36m_ElementwiseFunctionArray.__array__\u001b[1;34m(self, dtype)\u001b[0m\n\u001b[0;32m     69\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__array__\u001b[39m(\u001b[38;5;28mself\u001b[39m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m---> 70\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfunc\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43marray\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\coding\\variables.py:137\u001b[0m, in \u001b[0;36m_apply_mask\u001b[1;34m(data, encoded_fill_values, decoded_fill_value, dtype)\u001b[0m\n\u001b[0;32m    133\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_apply_mask\u001b[39m(\n\u001b[0;32m    134\u001b[0m     data: np\u001b[38;5;241m.\u001b[39mndarray, encoded_fill_values: \u001b[38;5;28mlist\u001b[39m, decoded_fill_value: Any, dtype: Any\n\u001b[0;32m    135\u001b[0m ) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m np\u001b[38;5;241m.\u001b[39mndarray:\n\u001b[0;32m    136\u001b[0m     \u001b[38;5;124;03m\"\"\"Mask all matching values in a NumPy arrays.\"\"\"\u001b[39;00m\n\u001b[1;32m--> 137\u001b[0m     data \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43masarray\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdtype\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    138\u001b[0m     condition \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mFalse\u001b[39;00m\n\u001b[0;32m    139\u001b[0m     \u001b[38;5;28;01mfor\u001b[39;00m fv \u001b[38;5;129;01min\u001b[39;00m encoded_fill_values:\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\indexing.py:423\u001b[0m, in \u001b[0;36mLazilyIndexedArray.__array__\u001b[1;34m(self, dtype)\u001b[0m\n\u001b[0;32m    421\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__array__\u001b[39m(\u001b[38;5;28mself\u001b[39m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m    422\u001b[0m     array \u001b[38;5;241m=\u001b[39m as_indexable(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray)\n\u001b[1;32m--> 423\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39masarray(\u001b[43marray\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m, dtype\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m)\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\backends\\pydap_.py:37\u001b[0m, in \u001b[0;36mPydapArrayWrapper.__getitem__\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m     36\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[1;32m---> 37\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mindexing\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mexplicit_indexing_adapter\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m     38\u001b[0m \u001b[43m        \u001b[49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mindexing\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mIndexingSupport\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mBASIC\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_getitem\u001b[49m\n\u001b[0;32m     39\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\core\\indexing.py:712\u001b[0m, in \u001b[0;36mexplicit_indexing_adapter\u001b[1;34m(key, shape, indexing_support, raw_indexing_method)\u001b[0m\n\u001b[0;32m    690\u001b[0m \u001b[38;5;124;03m\"\"\"Support explicit indexing by delegating to a raw indexing method.\u001b[39;00m\n\u001b[0;32m    691\u001b[0m \n\u001b[0;32m    692\u001b[0m \u001b[38;5;124;03mOuter and/or vectorized indexers are supported by indexing a second time\u001b[39;00m\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    709\u001b[0m \u001b[38;5;124;03mIndexing result, in the form of a duck numpy-array.\u001b[39;00m\n\u001b[0;32m    710\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m    711\u001b[0m raw_key, numpy_indices \u001b[38;5;241m=\u001b[39m decompose_indexer(key, shape, indexing_support)\n\u001b[1;32m--> 712\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mraw_indexing_method\u001b[49m\u001b[43m(\u001b[49m\u001b[43mraw_key\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtuple\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    713\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m numpy_indices\u001b[38;5;241m.\u001b[39mtuple:\n\u001b[0;32m    714\u001b[0m     \u001b[38;5;66;03m# index the loaded np.ndarray\u001b[39;00m\n\u001b[0;32m    715\u001b[0m     result \u001b[38;5;241m=\u001b[39m NumpyIndexingAdapter(np\u001b[38;5;241m.\u001b[39masarray(result))[numpy_indices]\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\backends\\pydap_.py:45\u001b[0m, in \u001b[0;36mPydapArrayWrapper._getitem\u001b[1;34m(self, key)\u001b[0m\n\u001b[0;32m     41\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m_getitem\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[0;32m     42\u001b[0m     \u001b[38;5;66;03m# pull the data from the array attribute if possible, to avoid\u001b[39;00m\n\u001b[0;32m     43\u001b[0m     \u001b[38;5;66;03m# downloading coordinate data twice\u001b[39;00m\n\u001b[0;32m     44\u001b[0m     array \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mgetattr\u001b[39m(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124marray\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39marray)\n\u001b[1;32m---> 45\u001b[0m     result \u001b[38;5;241m=\u001b[39m \u001b[43mrobust_getitem\u001b[49m\u001b[43m(\u001b[49m\u001b[43marray\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkey\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcatch\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;167;43;01mValueError\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[0;32m     46\u001b[0m     \u001b[38;5;66;03m# in some cases, pydap doesn't squeeze axes automatically like numpy\u001b[39;00m\n\u001b[0;32m     47\u001b[0m     axis \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mtuple\u001b[39m(n \u001b[38;5;28;01mfor\u001b[39;00m n, k \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28menumerate\u001b[39m(key) \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(k, integer_types))\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\xarray\\backends\\common.py:64\u001b[0m, in \u001b[0;36mrobust_getitem\u001b[1;34m(array, key, catch, max_retries, initial_delay)\u001b[0m\n\u001b[0;32m     62\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(max_retries \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m):\n\u001b[0;32m     63\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m---> 64\u001b[0m         \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43marray\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\n\u001b[0;32m     65\u001b[0m     \u001b[38;5;28;01mexcept\u001b[39;00m catch:\n\u001b[0;32m     66\u001b[0m         \u001b[38;5;28;01mif\u001b[39;00m n \u001b[38;5;241m==\u001b[39m max_retries:\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\model.py:323\u001b[0m, in \u001b[0;36mBaseType.__getitem__\u001b[1;34m(self, index)\u001b[0m\n\u001b[0;32m    321\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, index):\n\u001b[0;32m    322\u001b[0m     out \u001b[38;5;241m=\u001b[39m copy\u001b[38;5;241m.\u001b[39mcopy(\u001b[38;5;28mself\u001b[39m)\n\u001b[1;32m--> 323\u001b[0m     out\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_data_index\u001b[49m\u001b[43m(\u001b[49m\u001b[43mindex\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    324\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m out\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\model.py:353\u001b[0m, in \u001b[0;36mBaseType._get_data_index\u001b[1;34m(self, index)\u001b[0m\n\u001b[0;32m    351\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m np\u001b[38;5;241m.\u001b[39mvectorize(decode_np_strings)(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_data[index])\n\u001b[0;32m    352\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 353\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_data\u001b[49m\u001b[43m[\u001b[49m\u001b[43mindex\u001b[49m\u001b[43m]\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:175\u001b[0m, in \u001b[0;36mBaseProxy.__getitem__\u001b[1;34m(self, index)\u001b[0m\n\u001b[0;32m    173\u001b[0m \u001b[38;5;66;03m# Parse received dataset:\u001b[39;00m\n\u001b[0;32m    174\u001b[0m dataset \u001b[38;5;241m=\u001b[39m build_dataset(dds)\n\u001b[1;32m--> 175\u001b[0m dataset\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m=\u001b[39m \u001b[43munpack_data\u001b[49m\u001b[43m(\u001b[49m\u001b[43mBytesReader\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdata\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdataset\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    176\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m dataset[\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mid]\u001b[38;5;241m.\u001b[39mdata\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:452\u001b[0m, in \u001b[0;36munpack_data\u001b[1;34m(xdr_stream, dataset)\u001b[0m\n\u001b[0;32m    450\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21munpack_data\u001b[39m(xdr_stream, dataset):\n\u001b[0;32m    451\u001b[0m     \u001b[38;5;124;03m\"\"\"Unpack a string of encoded data, returning data as lists.\"\"\"\u001b[39;00m\n\u001b[1;32m--> 452\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43munpack_children\u001b[49m\u001b[43m(\u001b[49m\u001b[43mxdr_stream\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdataset\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:387\u001b[0m, in \u001b[0;36munpack_children\u001b[1;34m(stream, template)\u001b[0m\n\u001b[0;32m    383\u001b[0m         out\u001b[38;5;241m.\u001b[39mappend(\u001b[38;5;28mtuple\u001b[39m(unpack_children(stream, col)))\n\u001b[0;32m    385\u001b[0m     \u001b[38;5;66;03m# unpack arrays\u001b[39;00m\n\u001b[0;32m    386\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m--> 387\u001b[0m         out\u001b[38;5;241m.\u001b[39mextend(\u001b[43mconvert_stream_to_list\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstream\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcol\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mdtype\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcol\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[43m,\u001b[49m\n\u001b[0;32m    388\u001b[0m \u001b[43m                                          \u001b[49m\u001b[43mcol\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mid\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[0;32m    389\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m out\n",
      "File \u001b[1;32mC:\\Programs\\Anaconda\\envs\\hyrax_daymet\\lib\\site-packages\\pydap\\handlers\\dap.py:413\u001b[0m, in \u001b[0;36mconvert_stream_to_list\u001b[1;34m(stream, parser_dtype, shape, id)\u001b[0m\n\u001b[0;32m    410\u001b[0m stream\u001b[38;5;241m.\u001b[39mread(\u001b[38;5;241m4\u001b[39m)  \u001b[38;5;66;03m# read additional length\u001b[39;00m\n\u001b[0;32m    411\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m    412\u001b[0m     out\u001b[38;5;241m.\u001b[39mappend(\n\u001b[1;32m--> 413\u001b[0m         \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfrombuffer\u001b[49m\u001b[43m(\u001b[49m\n\u001b[0;32m    414\u001b[0m \u001b[43m            \u001b[49m\u001b[43mstream\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mread\u001b[49m\u001b[43m(\u001b[49m\u001b[43mcount\u001b[49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mresponse_dtype\u001b[49m\u001b[43m)\u001b[49m\n\u001b[0;32m    415\u001b[0m         \u001b[38;5;241m.\u001b[39mastype(parser_dtype)\u001b[38;5;241m.\u001b[39mreshape(shape))\n\u001b[0;32m    416\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m \u001b[38;5;28;01mas\u001b[39;00m e:\n\u001b[0;32m    417\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mstr\u001b[39m(e) \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mtotal size of new array must be unchanged\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[0;32m    418\u001b[0m         \u001b[38;5;66;03m# server-side failure.\u001b[39;00m\n\u001b[0;32m    419\u001b[0m         \u001b[38;5;66;03m# it is expected that the user should be mindful of this:\u001b[39;00m\n",
      "\u001b[1;31mValueError\u001b[0m: buffer size must be a multiple of element size"
     ]
    }
   ],
   "source": [
    "# THIS IS ACCESSING DAYMET DAILY DATA in NGAP HYDRAX via OPeNDAP.  \n",
    "# for loop based on response granules: granule_names  \n",
    "\n",
    "from pydap.client import open_url\n",
    "from pydap.cas.urs import setup_session\n",
    "import xarray as xr\n",
    "import time\n",
    "\n",
    "username = \"YourEarthdataLogin_UserName\" # Earthdata URS Login (only needed for NGAP Hyrax) Get account here https://urs.earthdata.nasa.gov/\n",
    "password= \"YourEarthdataLogin_pw\"\n",
    "\n",
    "hyrax_url = 'https://opendap.earthdata.nasa.gov/collections/C2532426483-ORNL_CLOUD/granules/'\n",
    "\n",
    "session = setup_session(username, password, check_url=hyrax_url)  \n",
    "\n",
    "before = time.time()\n",
    "cnt = 0 \n",
    "for g_name in granule_names:\n",
    "    print('granule_name ---->', g_name)\n",
    "    granule_dap = hyrax_url + g_name \n",
    "    print(granule_dap)\n",
    "    \n",
    "    hyrax_ds = open_url(granule_dap, session=session)\n",
    "    print('hyrax_ds ----->', hyrax_ds)\n",
    "    ds = xr.open_dataset(xr.backends.PydapDataStore(hyrax_ds), decode_coords=\"all\")\n",
    "    print(ds)\n",
    "    \n",
    "    #tmax_sel=ds['tmax'].sel(x=slice(lccbounds.minx[0],lccbounds.maxx[0]), y=slice(lccbounds.maxy[0],lccbounds.miny[0]))\n",
    "    tmax_sel=ds['tmax'].sel(x=slice(1386957,1470587), y=slice(-565105,-612937)) #GSMNP boundary\n",
    "    #print(tmax_sel)\n",
    "    if cnt==0:\n",
    "        tmax_h = tmax_sel\n",
    "        print('tmax_h ====')\n",
    "        print(tmax_h)\n",
    "        #tmax_h.to_netcdf('test.nc')\n",
    "    else:\n",
    "        #print('elsey')\n",
    "        tmax_h = xr.concat([tmax_h, tmax_sel], dim=\"time\")\n",
    "    \n",
    "    cnt += 1\n",
    "\n",
    "print('concat tmax_h ----->', tmax_h)    \n",
    "# save to netcdf\n",
    "#tmax_h.to_netcdf(var + '_cali_hyrax_subset.nc')\n",
    "print(\"Processing Time: \", time.time() - before, 'seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5750a98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# SKIP THIS - TESTING WITH ONE SLICE\n",
    "# THIS IS ACCESSING DAYMET DAILY DATA in NGAP HYDRAX via OPeNDAP.  \n",
    "# for loop based on response granules eg.) Daymet_Daily_V4.daymet_v4_daily_na_tmax_2010.nc, Daymet_Daily_V4.daymet_v4_daily_na_tmax_2011.nc \n",
    "\n",
    "from pydap.client import open_url\n",
    "from pydap.cas.urs import setup_session\n",
    "import xarray as xr\n",
    "import time\n",
    "\n",
    "username = \"mthornton\" # Earthdata URS Login (only needed for NGAP Hyrax) Get account here https://urs.earthdata.nasa.gov/\n",
    "password= \"Earthdata4Michele\"\n",
    "\n",
    "hyrax_url = 'https://opendap.earthdata.nasa.gov/collections/C2532426483-ORNL_CLOUD/granules/'\n",
    "\n",
    "session = setup_session(username, password, check_url=hyrax_url)  \n",
    "\n",
    "before = time.time()\n",
    "cnt = 0 \n",
    "for g_name in granule_names:\n",
    "    print('granule_name ---->', g_name)\n",
    "    #granule_dap = thredds_url + g_name.replace('Daymet_Daily_V4R1.','') #The name suffix is needed for cloud file name, not here for TDS\n",
    "    granule_dap = hyrax_url + g_name \n",
    "    print(granule_dap)\n",
    "    \n",
    "    hyrax_ds = open_url(granule_dap, session=session)\n",
    "    ds = xr.open_dataset(xr.backends.PydapDataStore(hyrax_ds), decode_coords=\"all\")\n",
    "    \n",
    "    subset=ds['tmax'].sel(x=slice(lccbounds.minx[0],lccbounds.maxx[0]), y=slice(lccbounds.maxy[0],lccbounds.miny[0]))\n",
    "    print(subset)\n",
    "    #dummy=ds['tmax'].sel(x=slice(1386957,1470587), y=slice(-565105,-612937))\n",
    "    subset.to_netcdf(g_name.replace('Daymet_Daily_V4R1.','CAsubset_'))\n",
    "    #subset.to_netcdf(CA_test)\n",
    "    \n",
    "    #if cnt==0:\n",
    "    #    tmax_h = dummy\n",
    "    #else:\n",
    "    #    tmax_h = xr.concat([tmax_h, dummy], dim=\"time\")\n",
    "    \n",
    "    #cnt += 1\n",
    "\n",
    "# save to netcdf\n",
    "#tmax_h.to_netcdf(var + '_cali_hyrax_subset.nc')\n",
    "#granule_name ----> Daymet_Daily_V4.daymet_v4_daily_na_tmax_2010.nc\n",
    "#subset.to_netcdf(g_name.replace('Daymet_Daily_V4R1.','CAsubset_'))\n",
    "print(\"Processing Time: \", time.time() - before, 'seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32f18864",
   "metadata": {},
   "outputs": [],
   "source": [
    "# THIS IS ACCESSING DAYMET DAILY DATA in NGAP HYDRAX via OPeNDAP.  \n",
    "# for loop based on response granules eg.) Daymet_Daily_V4.daymet_v4_daily_na_tmax_2010.nc, Daymet_Daily_V4.daymet_v4_daily_na_tmax_2011.nc \n",
    "\n",
    "from pydap.client import open_url\n",
    "from pydap.cas.urs import setup_session\n",
    "import xarray as xr\n",
    "import time\n",
    "\n",
    "username = \"mthornton\" # Earthdata URS Login (only needed for NGAP Hyrax) Get account here https://urs.earthdata.nasa.gov/\n",
    "password= \"Earthdata4Michele\"\n",
    "\n",
    "hyrax_url = 'https://opendap.earthdata.nasa.gov/collections/C2532426483-ORNL_CLOUD/granules/'\n",
    "\n",
    "session = setup_session(username, password, check_url=hyrax_url)  \n",
    "\n",
    "before = time.time()\n",
    "cnt = 0 \n",
    "for g_name in granule_names:\n",
    "    print('granule_name ---->', g_name)\n",
    "    #granule_dap = thredds_url + g_name.replace('Daymet_Daily_V4R1.','') #The name suffix is needed for cloud file name\n",
    "    granule_dap = hyrax_url + g_name \n",
    "    print(granule_dap)\n",
    "    \n",
    "    hyrax_ds = open_url(granule_dap, session=session)\n",
    "    ds = xr.open_dataset(xr.backends.PydapDataStore(hyrax_ds), decode_coords=\"all\")\n",
    "    print(ds)\n",
    "    \n",
    "    dummy=ds['tmax'].sel(x=slice(lccbounds.minx[0],lccbounds.maxx[0]), y=slice(lccbounds.maxy[0],lccbounds.miny[0]))\n",
    "    print('dummy')\n",
    "    print(dummy)\n",
    "    #dummy=ds['prcp'].sel(x=slice(1386957,1470587), y=slice(-565105,-612937))\n",
    "    if cnt==0:\n",
    "        tmax_h = dummy\n",
    "        print('tmax_h ====')\n",
    "        print(tmax_h)\n",
    "        tmax_h.to_netcdf('test.nc')\n",
    "    else:\n",
    "        tmax_h = xr.concat([tmax_h, dummy], dim=\"time\")\n",
    "    \n",
    "    cnt += 1\n",
    "\n",
    "print(tmax_h)    \n",
    "# save to netcdf\n",
    "tmax_h.to_netcdf(var + '_cali_hyrax_subset.nc')\n",
    "print(\"Processing Time: \", time.time() - before, 'seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "403eb22e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax_h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07ceadbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "tmax = tmax_h.isel(time=200) # isel = xarray index selection (python index start from 0, so day 7) \n",
    "ax = plt.subplot()\n",
    "tmax.plot(ax=ax, robust=True)\n",
    "# Note xarray uses dimension and attribute information image - time:2010-01-07, projection is recognized"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf62c55c",
   "metadata": {},
   "source": [
    "##### Let's take a look at the Daymet Temperature Subset for 2011"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b03f5085",
   "metadata": {},
   "outputs": [],
   "source": [
    "#subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec758d47",
   "metadata": {},
   "outputs": [],
   "source": [
    "###SKIP THIS\n",
    "import matplotlib.pyplot as plt\n",
    "#tmax_ca = subset.isel(time=6) # isel = xarray index selection (python index start from 0, so day 7) \n",
    "tmax_ca = tmax_h.isel(time=6) # isel = xarray index selection (python index start from 0, so day 7) \n",
    "ax = plt.subplot()\n",
    "tmax_ca.plot(ax=ax, robust=True)\n",
    "# Note xarray uses dimension and attribute information image - time:2010-01-07, projection is recognized"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb437875",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "tmax_ca = tmax_h.isel(time=6)\n",
    "fig, ax = plt.subplots(figsize = (20,10))\n",
    "tmax_ca.plot(ax=ax, robust=True, cbar_kwargs={'label': 'degrees C'})\n",
    "ca_poly_lcc_1.plot(ax = ax, color = 'none', edgecolor = 'red')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e94e398",
   "metadata": {},
   "source": [
    "### Step 4. Clipping the subset area by polygon\n",
    "In this step, we demonstrate how to clip the subset to our polygon of interest.  We'll use Python's *rasterio* library.\n",
    "\n",
    "* **`rasterio`** provides many capabilities found in GDAL\n",
    "\n",
    "First, we'll map the subset precipitation data matrix to real world coordinates using the **'translation'** and **'scale'** methods of **'Affine'**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34ec5a09",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.features import Affine\n",
    "res = 1000 # daymet spatial resolution is 1km x 1km\n",
    "#transform = Affine.translation(min(subset.x).values - res / 2, min(subset.y).values - res / 2) * Affine.scale(res, res)\n",
    "transform = Affine.translation(min(tmax_h.x).values - res / 2, min(tmax_h.y).values - res / 2) * Affine.scale(res, res)\n",
    "transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e33d30e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from rasterio.features import geometry_mask\n",
    "#ca_mask = geometry_mask(ca_poly_lcc_1.geometry, out_shape=(len(subset.y), len(subset.x)), transform=transform, invert=True)\n",
    "ca_mask = geometry_mask(ca_poly_lcc_1.geometry, out_shape=(len(tmax_h.y), len(tmax_h.x)), transform=transform, invert=True)\n",
    "polymask = xr.DataArray(ca_mask , dims=(\"y\", \"x\"))\n",
    "polymask = polymask.sel(y=slice(None, None, -1)) \n",
    "tmax_masked = tmax_h.where(polymask == True)\n",
    "# write out to a netCDF file\n",
    "tmax_masked.to_netcdf(var + '_subset_poly.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b4df5a",
   "metadata": {},
   "source": [
    "Now let's plot the subset for the day"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e0cf751",
   "metadata": {},
   "outputs": [],
   "source": [
    "##SKIP THIS\n",
    "import matplotlib.pyplot as plt\n",
    "tmax_1 = tmax_masked.isel(time=6)\n",
    "ax = plt.subplot()\n",
    "tmax_1.plot(ax=ax, cmap=\"turbo\", cbar_kwargs={'label': 'degrees C'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9725d5ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "tmax_1 = tmax_masked.isel(time=6)\n",
    "fig, ax = plt.subplots(figsize = (20,10))\n",
    "tmax_1.plot(ax=ax, robust=True, cmap=\"turbo\", cbar_kwargs={'label': 'degrees C'})\n",
    "ca_poly_lcc_1.plot(ax = ax, color = 'none', edgecolor = 'red')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b346cee5",
   "metadata": {},
   "source": [
    "#### Saving to a GeoTIFF\n",
    "The subset data can also be clipped by a polygon of interest and converted to GeoTiff using the rioxarray and shapely libraries.\n",
    "\n",
    "* **`rioxarray`** is an extention of xarray.  It can be used to clip, merge, and reproject rasters\n",
    "* **`shapely`** is a Python package for manipulation and analysis of planar geometery objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "080269cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rioxarray\n",
    "from shapely.geometry import mapping\n",
    "print(var)\n",
    "\n",
    "#prcp_sub = rioxarray.open_rasterio(var + '_subset_poly.nc', decode_coords=\"all\")\n",
    "#tmax_sub = rioxarray.open_rasterio(var + '_tdssubset.nc', decode_coords=\"all\", decode_times=False)\n",
    "\n",
    "#tmax_sub = rioxarray.open_rasterio(var + '_subset_poly.nc', decode_coords=\"all\", decode_times=False)\n",
    "tmax_sub = rioxarray.open_rasterio(var + '_cali_hyrax_subset.nc', masked=True)\n",
    "tmax_clip = tmax_sub.tmax.rio.clip(ca_poly_lcc_1.geometry.apply(mapping), ca_poly_lcc_1.crs, invert=False, from_disk=True)\n",
    "\n",
    "# saving one time step to geotiff\n",
    "tmax_clip_day6 = tmax_clip.isel(time=6)\n",
    "tmax_clip_day6.rio.to_raster(\"tmax_clip_day6.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20d641ed",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69f56f18",
   "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}