hrrr_zarr_one_day_analysis_not_mf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d8da0b2f-7cfa-4b2b-9102-cf855d0a9d55",
   "metadata": {},
   "source": [
    "# HRRR Zarr Example\n",
    "## Download an analysis variable, whole grid, 24 hours\n",
    "\n",
    "This example shows a fast and convenient way to download and merge multiple hours of analysis data from the HRRR Zarr S3 archive. It should be easy to adapt for cloud computing using Dask as well.\n",
    "\n",
    "If you only want data at a single grid point or a small region at many different time points, please refer to other examples on how to download data by chunk. Doing so is much faster, especially as you look at more hours or variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4efb4228-1505-4074-b1fd-c5f764e8aa31",
   "metadata": {},
   "outputs": [],
   "source": [
    "import s3fs\n",
    "import xarray as xr\n",
    "import metpy\n",
    "import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd6258d1-9b45-4911-87ca-b40107f5d41b",
   "metadata": {},
   "source": [
    "Because the HRRR Zarr archive contains metadata (such as the grid information) in a group and the actual data in a subgroup, the xarray open_dataset function doesn't work if all you want to do is open data for a single hour––for that, call xr.open_mfdataset instead. \n",
    "\n",
    "But xr.open_dataset is faster and has the capability to lazy-load, so since we're already going to be combining multiple datasets, let's go ahead and use that function, which will work equally well if we feed it the group url (it'll just load the metadata) as if we pass the subgroup url (it'll load the data variable but not the metadata)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1e7f2ab4-e2b3-456f-9a65-129e875509b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load(url, time=None):\n",
    "    # Open the dataset from S3 - this is lazy (doesn't download the actual data)\n",
    "    ds = xr.open_dataset(s3fs.S3Map(url, s3=fs), engine=\"zarr\")\n",
    "    ds[\"time\"] = time\n",
    "    ds = ds.set_coords(\"time\")\n",
    "    # Later on we use metpy functions that expect the grid variables to be x and y\n",
    "    ds = ds.rename(projection_x_coordinate=\"x\", projection_y_coordinate=\"y\")\n",
    "    return ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bb7c342-7e04-4b1f-a06c-1797e0d16cea",
   "metadata": {},
   "source": [
    "Now we actually create the dataset we want in 3 steps:\n",
    "- Download the grid data\n",
    "- Download the actual data for the variable of interest\n",
    "- Add the grid data to the dataset and enhance with metpy functions that understand the coordinate reference system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e4cb1cce-c7b2-4a28-8d01-81efa85e4138",
   "metadata": {},
   "outputs": [],
   "source": [
    "projection = ccrs.LambertConformal(central_longitude=262.5, \n",
    "                               central_latitude=38.5, \n",
    "                               standard_parallels=(38.5, 38.5),\n",
    "                                globe=ccrs.Globe(semimajor_axis=6371229,\n",
    "                                                 semiminor_axis=6371229))\n",
    "\n",
    "\n",
    "def load_combined_dataset(start_date, num_hours, level, param_short_name):\n",
    "    group_url_format = f\"s3://hrrrzarr/sfc/%Y%m%d/%Y%m%d_%Hz_anl.zarr/{level}/{param_short_name}\"\n",
    "    \n",
    "    # Get data at the group level just once so we have the grid data\n",
    "    group_url = start_date.strftime(group_url_format)\n",
    "    grid = load(group_url)\n",
    "    \n",
    "    # Load and combine the data for each hour\n",
    "    datasets = []\n",
    "    for i in range(num_hours):\n",
    "        time = start_date + datetime.timedelta(hours=i)\n",
    "        subgroup_url = time.strftime(group_url_format) + f\"/{level}\"\n",
    "        datasets.append(load(subgroup_url, time=time))\n",
    "    ds = xr.concat(datasets, dim=\"time\", combine_attrs=\"override\")\n",
    "    \n",
    "    # Add the grid data to the combined dataset\n",
    "    ds[\"x\"] = grid[\"x\"]\n",
    "    ds[\"y\"] = grid[\"y\"]  \n",
    "    \n",
    "    ds = ds.metpy.assign_crs(projection.to_cf())\n",
    "    ds = ds.metpy.assign_latitude_longitude()\n",
    "    \n",
    "    return ds"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c770b017-2adf-4483-bc5c-ba02e05e6a2d",
   "metadata": {},
   "source": [
    "Running this for 24 hours takes about half a minute."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e2ba87c6-d6c8-417c-a4bc-1ade03c0323e",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.7 s, sys: 433 ms, total: 11.1 s\n",
      "Wall time: 33.2 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "ds = load_combined_dataset(datetime.datetime(2021, 4, 1), 24, \"1000mb\", \"TMP\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0d63513-f407-4874-b842-2f903d9e31cc",
   "metadata": {},
   "source": [
    "At this point, only metadata has been downloaded––the actual temperature data isn't downloaded until we access it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d177c2a4-eba0-47a4-8136-626cd8863403",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.78 s, sys: 47.5 ms, total: 1.83 s\n",
      "Wall time: 1.84 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[0.6787, 0.6807, 0.6787, ..., 0.4114, 0.4114, 0.4114],\n",
       "       [0.6675, 0.6675, 0.6714, ..., 0.4114, 0.4114, 0.4114],\n",
       "       [0.6675, 0.6675, 0.6675, ..., 0.4114, 0.4114, 0.4114],\n",
       "       ...,\n",
       "       [0.9087, 0.903 , 0.9043, ..., 0.7217, 0.7603, 0.7515],\n",
       "       [0.9087, 0.8955, 0.903 , ..., 0.7754, 0.792 , 0.7554],\n",
       "       [0.9155, 0.8955, 0.8955, ..., 0.7617, 0.8086, 0.7954]],\n",
       "      dtype=float16)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "std_dev = ds.TMP.std(dim=\"time\")\n",
    "std_dev.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "aaf49458-cf5c-493c-822a-8b6607359ba5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAACfCAYAAABa3wq8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABRIElEQVR4nO2dd3xb1fn/31fLkrzlbSeOndiJ4+xBNgkJWRBGoRQoqy0t5Vu+dJd+6aa00EJpC2X0B6VQyh6FUlYCKYSEhGyynOnEe+8ly9a4vz+MFFm+kq6WJTl6v155xbr33HOPdM/93Oc+5znPEURRJEaMGDFijA6KcDcgRowYMc4lYqIbI0aMGKNITHRjxIgRYxSJiW6MGDFijCIx0Y0RI0aMUSQmujFixIgxiqg87Vy/fr3Y2trqdr/NZqOmpgZ7mezsbHJyclAo3Gu52Wzm0KFDkvsKCwtJSkri2LFjDA4OkpKSwqRJk+R8jxijRFl9k9/H2rTuwxMVJkFWuUhkhiE73E0YhtVq5fTp0wiCQFFREYIgeD/oc+rr62lvb2fatGk+HRdjOPv27dskiuJ6qX2Cpzjd+fPni3v37vV6gu3bt7Ns2TIApk+fzuHDhz2W7+7u5tNPP2X9+vWsX7+ejRs3ui373//+l1WrVnltQwx5TP3Fn4d9NpWYhn3WHteOOMa1jFxc6/JUj72snHMtLy5n66kij/u94en4YCCnDXaeWfD3oJ//i1/8Iq+//jrPPvssN9xwg0/HCoLAq6++ylVXXRX0dp0rCIKwTxTF+ZL7giG6AEePHuWWW25hx44dvPvuu1x00UVej2lra+PgwYOsXr2aBQsWsGvXLh566CFuv/12uru7mTJlCs3NzQwODqJWq2W141zDVURBnrjJKeev2PqLp3N6EjFXAQ2n6Mo5921ZHwLwWJN7Y8Jd+7THtRz7zfe9nuPee++lpaWFN954g4yMDL75zW9y8803o1Qqh5W77777uPPOO7n11lt58MEH+clPfsJLL73EwYMHyczM9HqeGNJ4El2P7gVfKC0t5ZNPPuGOO+6goKBA1jFpaWmsWrWK73znOzz00EMALF26lObmZpqbm9m4cSMNDQ2oVEFrZtQgJaaumEpMUOJbvaYSU1gEVQ6e2rT1VJFbQfNm+Y4lTCUmt30j7ajV6ZOGT/71Z+6//362bNnCb37zG55++mleeeUVxo0b5yjV3t5OcnIy11xzDT/72c84cOBATHBDTFDVTBAEHnjgAZ+PmzVrFgB6vZ4f/OAHbN26ddj+hoYGsrMjy28WCBcVDFkqPfPyJPe3lSoltzsTiGjKOda1jF3wwiVucl/XI93KBc8WrhzsD05X2kqVw4R32RfP3oti+iXUHvuA0lnLmLHiNgRhaNzl+T/dRlVVFe+88w6bNm0iJSWF1NTUgNoXwzMRYUIuW7aMhIQEJk+ezNatW/n973/PLbfcQmJiYlS6Feyi6gl3giuHSLRSg4mrwPviHw0XwWpjoA8EZ+Ed/vBWoi1cQc+fNnNi9/OMm7IKbXwaU2csICmtEOv2E7TWlgGw5PK7idMlA/DJv34UUHtijCQiRLe4uJienh7++Mc/kpyczA9/+MNhLoW+vj7i4+PD2MKRyBHWUBAOwQ2WoPhqLYdCbMeCG0LKynVG6k3J1FhL5ZN/AsBmNXNi93OYeltRqrQUzR0aMGutPUj2xMUOwYXh1rIrMUH2j6ANpAVCQ0MDx44dw2Qy8cADD9De3s5XvvIV9Ho9NTU13HPPPZw4cYLJkyePOLatrY3u7m4KCwuD3q5AhNVuySbuq/NaRgrJGyeMFq6zALobuJISNF+E05PfVs6x3s4XSsEdTUvXm+hKYWqqp/vIPtJXXIRCpcJQZqHu5EdUlb2HUq1j4SV3IYo2QECh8O7e8sa5LsijEr3gL9u3b2flypUsXryYhoYGTp06xdq1a9FoNCQlJVFYWMg999wDQH5+Pt/97ne58cYbUSgUpKWlIQgCs2bN4sCBA363QUpce+bleRRMuXgTX2fhDbUvNxCkREVK6JxFY7TcAuEWXHf48/1DJbqu2F0QXa1nOLL1r0w+7zoyxs8JuF5vnCtiHNGiW15ezvnnn8/06dPZvHmzY/vChQs5cuQI8+fP5+OPPx5xnMFgYP369bzwwgvccsstPPHEE7LOJ9d6dRVdOZaru3qccT3eV9F1ZTRFOBp8q1JEmkshGA+HYAgvDImvxWxCoVQHxcKVi2tflxMGF01EtOgCtLS0MH/+fKqrq5kyZQo/+MEP+OY3v+mYRPHZZ5+xZ88eTpw4QVVVFf39/aSmpjpmwh04cMARAQFQW1vLtd99yfE5mBarP3VKHRtI5IIzo235RpvwRprg2nH3O/ra3kDFd3iY2eghp5+HQohPnDhBWVkZRqORDz/8kOTkZKqrq2lrayMvLw+tVktzczNqtRqDwUBVVRVxcXHMmzcPrVbLwYMHMZlMKJVKDAYDhYWFpKens3LlSgoKCjhy5Agmk4klS5ZEtugCnDx5kilTpqDT6bBYLMydO5ft27ePCOYGMBqNLFmyhIMHDwKgUyVxfu5XUQgKBqx9fFT7N/JL1zG+ZDVwVuT8FV93Aim3Pl9EF+QLbyS5GiKZcLoW/PFzj3Xh9eeNzo6/Qnzs2DHuv/9+XnvtNZYtW4ZCoWDt2rUMDg6Sn59PamoqLS0tGI1GMjIyGBwcpK6ujqlTp2I0Gjlw4AB9fX3MmDGDpKQkrFYr7e3tlJeX09rayn/+8x8GBgaYMGECCQkJ7N69O/JFF4YG1LKysnjmmWe4+eabUavVaLVaXnvtNdauXQucHU0dNPVyaOMD6JQpCAow20xMSl5Is/E0DcYTzFlzB/pE6QBvu1jKEWI5oV1y6rB3NDkdPBS+XTmzvXy52aNFeCPN0pX7u42m8EaT6ErhTYhramo477zzuPXWW/n+979PSkqKx/Jbt25lxYoVAKSnpzNt2jQmTpxIYWEh5513HsuXL0ev1w87prGxEb1eT1JSEhAF7gUpzpw5Q3HRFGyixbEtp+h8Js68zPF50NTL/vfuITd+GjW9B0fUUTT3ahorPqVo7tXEJ8ubXOHJ5yr3GFd8Ed1ABFfuFN9g4a/wypkGGyycxcuX/A5SeIrg8KcOT/hTv7/CG02i63ztPH1fuxCLosjy5cu5+OKL+clPfiLrHB0dHfzkJz/h/fffp7GxkRUrVnD55ZdTWVnpcHdu2LCB2267jWXLlkkmBooa0XUd5Bq0Gqno3k9F91AbNNokJs35IoacUkeZjsbjnNjzPGmqcbQP1CAgYLYNXRiVoMEiDjK+ZDXxKXmk5U6X3RZna9gXpATYWXTtuOvowbRy5b7iBmIN+iq8dsGF0Iiup5wFEBzB9XQeX+sJtD5nokV0IbhjF1LfO+VQP6beVno6aqg5/gHz1/+U7a//2KdziqLIvn37+N73vsf48eP58Y9/zOzZs2lvb+eFF17gkUceIT4+nhtuuIGbb755mAUdkaIrN4pAFEXq+o5S13uUjoEhQTPklFKy6KsIgsCAsZPakx+i1RvobD5Fb0sluQnTSI3LZtBm4lj7RwCkj5vNlAXXezyXO8H0FbkDZf4K71gRXQi+8IZKdMF3N4ynLGvByhcRbf5cO76Iruu0Z6lrqD2uxWrqp+GtF+k9UYZCoUSbkE6cLoXc4hWkZBZhHuijs/kULYpaBlubmJVloKGhgba2NtLT0ykoKKCoqIgNGzawYcMGx3hSZ2cn9913H6+88gqCILBq1SpKS0spLi6mtraWZ555hp07d1JVVcX48eOBCBHdYMzgstos7Gx8mR5zC5n58xg09dDZfBKApVf+ARhyOZzc8zxdLUM3SLq2gFZTJdr4dGZecDvqOPcz24Ilut5wdjW4czt46pRyOqEcpKbbhtqv6yq6MFJ4A7GGR2MCRCBhXcF0+wQzXnc08dfKdTcuEVempnP/Dlq2vAeCQPZFV1HAdERRpK+zlgZFBb3lxxloqkM/oQhd3gTiMrNRxOlQJSSSUa3DPNjL73+0mrKyMl588UUGBwd58sknmTdvnuM8oihy8OBBtm3bxsmTJzl27BgnTpxArVYjiiJvvfUW06cPvU2HRXRDOU22tb+alox+mqv30d/TzKxV3yMh5aw4DvR38dnG+7CKZqanrUGBgpreIwxqRaYu/hpx+hS3dfvi05Vjtco5zjVRiVQ5Z7x1QrkEksQmkIE0KeF1RySJLsh/MHkTxECuWzBjdEebYMeiC590UfvKU6gSEslafyVx6dn011bSc+Iw3Uf2oUpIQp8/CX3hZOILi1GoNY5j3X3/hjM7OHPgDWBIaP1hVEQ3HLkIegfbOKY4SGfTCeJTxpGSWYwheyr65Gza6g5zYvdzJKuz6TI3AUPfc9LsK8meuNhtnb74ct2JrjsRdT3GFbnWrj83rJQ7wdMsM7l1+YMvouuMNwGOlEgFZ1H0ZeDU+WHqus1d/f4ip3/6K8rBik6wWcyY6qqIyx6HeZZEgS1t1Lz4OBmrLiF5xjxaGz6h880tKBL06OeWkLB0Dold44Yd4vyd3BlYoihScehN+ntbmbb0G479vsymC4nohivhiys98/KoOPwWTRW70SdnM2DsxGYdRJeQQU97FalxeRgtnSDCgK2PhNRxzFr5Xcm6pGaggXsRldO5XMVYCuebzVv5YA0EeSPQPApy8Ed8/Un8PZrYBdFXwfLmTpI6R6D44tIKB/311VQ99eCwbcmzFmBYvBJrv5H6N54lffl6+pVN9GzZgzonHcP1G4grPHvvunsAuos4qiswc/qzfzFg7KBk4U2kZI3M92I//r1K9zmvgyK6kSKyrvTMy+PknhewWgbQ6JJpry/DZrOiVGlQa+JRKNXo+1W09ldhEQcBmLv2x+gSMkbUJSW6Uh3R2RcrB2+d253oeirrDW8DR5ESZ+ur8EaD6LpeQ7lvT+EQPek0kKOPu35qM5oYrGvGXNtEx8sbsfX2f75HwLD4AnpqylDnZpJy1Wo0uWfj8qUeTAn72jm24yn6uupJjcuj29pKWu50UnNKUShUNFbspLutksKZl2LILvU4/gMjhdtZhP0W3eS4LHFJznUeTxxueublUXvyI7paTlOy8CYUShWDph4sg0Ysg0ZsNgvG7iYqD7/lOGZ8yRryS9cOq8fdq0agndEfCzlYwusJb0lq3GX8CiQTmLv43GD5eMMlup7cCb5mmQu3+IULuX3a2t1L7bfvJ3XhcrqO7iH12vXEL509IlZWSnSt7+3i1L6XnbYIxGXnoUpIxHjmJAqlhqkLb+LYzn9gyJnO5POuldUmqWu8serB0C/XE05yJi6lt6OW8s9eY8p51xGnSx6WE7Sr5fSw8uOmrHRbVzCjFdy9YjpbF96CvaVuQn8iF+SGiLmLZAhWpILr9seaVjmE1F9fbzjw5reVM0VcygKO9Ff+cGJubKP5oedBraDr2F6yf/FN1FlpI8q5c78oL1pIyUULSTtqxWaz0NNezdHd/wCbFW18Ov09TRz55HEAMgvO86ltFtsAxzu2kqUvIkPnOc2s+7XSowB7Z1WqNBTN/RLt9WVYLYPDyrQ3HqPu5EfDttlsozNq6+mGCfbN5Iufz5sA2/fLHXyTItj+2khBe1zrcCHY/zmTuK/O5xwfnsqHK452NJFjOPT8dxcNd/0Vc00j6qw0sn95q6TgykWhUJGcPpEZS26lqPgSlqdey+wLf4A6LpHzNtxFSsYk2XXV9x5nc81fqe0to9/S7bW8R0u3eMZ4aJV97rCiUmtJySymseJTtJefvXm7RbPXY91NZrALo9yph3JxFlypsC/ngTtv/mO5Qf+BuAVCRSBuhtHG0wCZnFzJ3oTY04QaX8cQxhrdH3xK578/AoWAdmohmXd8FYXG92W8pK5dQspQ/pXaggHKd71KbvH5aLz4cp0RRRuH2jYCsHr8bagUGt6r/DOC8KDbY7xaup5G6MKNa0fOmbSMysNvU/+fF7EODIlQ0rQ5qFOGnoiqpFQm//h3qNTuhbOtVOn4Fyie3AvO++zWk6/t8PVhINfnaS/nbN36ItjexNNfizZc/lxfBbdnXp5Pgutczl3Zc8HidaX/SDnV3/ot7c+8ha23H/38aWT95OtBE1w7JksvR7b+FVNvKy3Hd5C4r46EvbXUndrKiT0veKlZYKphJdPT1qBSaLyUHWJM+HTtJBryAeg+tIeMlRejjBu6WSbd/jNHGV8GOgIlGDGPdtxZsr4Oqsl1E4R7UcjHmlaF1fL15F93vpZyxgB8TYDvLfXnWEW02bB292E6eoaOF9/D2jH0qq4tnUjypSvQzSgOyXm1qgRW5N1Mh6meE52fIIoi5V07qazehUI5XEilrmFS4tkg4vcq/4zNZvN4Plmi+17lnyM2ZMwZpUqDOi4R80APKn3CiP2ehE9uR3edgisnDtd1nxwBDnWWMKl1ztytfeYLcq1cKavVmwV8W9aHkmWCaeXKeWNwjdt2d91dr3OwloAaa4iiSMujL2HceXjY9uTLV5K4fgmqxMAWpfV238NQTm5NvJ7yrp0cTzxFVfUuAGzWQco3PkFqXC4G7TjQZCCKIoO2fsw2E6JoQ6dKRqVQO7wCf/6zZ+9A1Fu6rmJpHugBoO/MCdQpBuIyzqZ09DRTzJ4JTE4guqvwjgVCvbaZP/kV5Fq7wUoe4wty3E9SmeXcCa+cB/5Y8+2auzvp3LeDtt8OLdMVv3Q2qdeuR5mcgKAIzhi/NwPH+XooBRWlhlUcLN/E3HV3UrnledoHamgynkIhKKno3suAtQ8AtUKLTbRiFc0ka7JYnPNlAHbt2sV9993n8ZyyRTdarN2ktEK62yqofflJx7bJ//f7YXOu7UiF7LgObPliccq9KeSUCXWyFDuB5F9wxZ1AehJV1+P89ff6I7yj/eC0P/TPRdeBM6Joo/blv9NXfgwA7fQiMr9zHQp98K+HJ0NLijTteFJ6Mtm/6feoBA0CCpSCmrb+agZtRke5KannU9b2IRMSZ1OYNI/3Kv/MJ598whVXXMFTTz3FZZdd5vYcPlm60SC8E6ZdRNmef5B/423YBgYYaKxFULn/mpIxkkfVklavHU83ayjiLH0RW3+S4PgjuFLREL74Yd2Jq6cE54GGlDn/NuF6U/FVBFyJlBlk/iKKImce/R3mzjaS5ywifelqrIv03g+Ugb/X1NnatYoWNAodC7K+hNHcwZH2zSRq0snUTUQpqOkYqKPBeIIjbR+QGpdHSeoK2gdqMRqNPPzww9xzzz0eBRfGgHthBAJoDOloM3MB0I8fHqjs3Omlkot7Q2ok291UYfv/wcqUH+zyctwI7kLN3B0rZb2688V6IhR5d13fYNyF24VakAMVXoiuSRRpR600F1upffEJ+qvPADDumq+TUDwNACuBvcEF63pZbIPU9x0byt/dd5RF2deyPPerVHTvp3OwgZb+Cvjc7aHRJpEWl8+OhufpMbdy220J1NXVsXPnTq/n8Vl0I93aVSg10NXvvaCPBGMNqmDcGFJC4SlZtidC6Qf1RSSDNSPNVxdDuBb2DEX4VyRawGlHrbTVHWb7rn8O255+wcUOwQ31NXD3Ww8YOxAUajTaBISdhznSvZ/GvpOO/CxFyYs/99lamJA4C+OMHMYrVCjVWrpbK6k5/gGd8f28+vJzTJ8+nXvvvReDwcDTTz/ttU1jwtJ1HkyzDPahiRsZueDuOHBv7Xqa6ukcPhSqGEpflyhx3SdlwXnLxetqwQYrb64cEXbnnvBFwH0RXqnfINRWrrdY30D9vXLE17WMvxn05LRDFM+GTwkKFXE5eRgWLg+objverpW9DTarGYu5H3VcIgDGrgYOfPhnVBo9CanjMbZUMyFxNlMNF3C47X0Ayrs+hS6nyhogM38+xfOvwZAzFUPOVBL31bFu3ToAHn30UX73u9/x1ltv4Q2/RDeSrV1jTzOKCTlBr9ddJ/RllV9vN4Q3QfRVEKQGBT3VEeyoBX9ibYNh9forvOHw87pGMvi7Np8r7vqacz/11GcDsZzTjlqxDPaz6+1fOrYVzLiEnElL6ZgeN6ysP2MQcsXWTs2xzdSe/BB9Uhb9vUNTbOOT80hMyycpfSL5mtX0mds5070HBUoy9ZNoNA6tSJM3eSWGnFJUGj26hPRh9b5b8SfKysp47733ePHFFzl9+jSXXHIJc+fO5e6773bbvjFh6Tpj7GrAkOK+w7oGtstJOuOpLnun9EV83eGpA8pplxzxcC7j7nzBnDIsZaUGM1xMCk+CK/XbjJbYjubEHOdzBmK1+nu81WIi0ZDPxFlXkJA6lEg8mO411zZ6QqkZOiYpbSIlC7+COi4BlUYHQEfTCbbXP4uIjYKkecxIW4dKoUYU11PbW0ZiySqUqqHIJ1EUMQ/00dNejansIFOnPonJZGLNmjXcf//9rFixAtXng/YhEd1ItXbNg0ZMExORmpAXbDeAt04kd8BESvR9XUjRk1hL7ZOyesMxCy1Yy7HLtWzDFWMtR3BDFUoWjEE7X9ORxulTmXnBtwM6rzOerpnNaqG77QxJaROxWgYQBAXdbRV0t1XQWnOAwf4uDHHjaKz4lNzi5ag0Ono7ajn40UOOOtbm345COCuHA1YjRksnTftfIT1vFsd3/RO9KoUBax/Jmiy++9NbWLNmDYsWLZJcgt0TAVm6kSK8zp3VMtjvMURMLv6+crp2zkBHqt092b35oeWIdaALXDpbxL5Yx84CaxddKfENlS833JNbonFWmpS7wd91AZ3xdyDN+dy73vsNlkHjiDLjEqYzO2ktCelpqBVx9Fu6sSSk01J7kJO7nwNgfuYVpGnzRwhnWft/aemvIK5PT2vtQTQKPbPS15OkyUIQBH75y1+OOJ9cxpx7obe3noyc8bLLu150qc+uQibn9d25nDfXg6+uA285dl0/exow20rRMOs6VLkcpJByHwRjooQcRlt4XUMVfVkANZLwx4AIVHA9DWh3tZzGMmhEpdE7hDe3eDlFfRPRq1OGldWpkmBfHQd7tgFw3kW/QKNLQnC5FjU9R2jpr2BOxiVk6Yf6tyiKDmEONAlYwKIbbmvXubNaLYOIZjOqhCTZx/siNJ7EV+4rPPi+3Lq3Or0Nenj7jvY6fRXcQKffehs08ye+1xdCKbze4rjP9TwMcvqa1LWxmE201R+mfN8r6JOyGDT1MmnOVeiTstBoE4nTpyIICqxAz+fHuP7O0zWLOZQ0gEY3XCdEUeR4x8e09lexJOc6kjRnl//x1YXgiTFl6VrM/SjitJI/kLuwGE9B8nbcDbTJtXzljLa6W4RSKpG4vU65yXbsbfTnVc6b28AuvOHM1xuI8IdCeF372rmYltET/gpu/xsfcLzsPcfnlMzJZE44j/hkz9FKrg+4JE0GxsZGdr9zF/NTL8OmSaOy+zOa+ysQEFiYfQ0apXSfCEaq26CIbritXTtNlbuJ12d6LyiBnJhVAIpH3uTeOpH9xva0qrDUEj7uVnhwdlt48u96eijY6/MkOHLFLBQJZ0Llz40RHbg+DOP3tHC87D2K5n6J+JQ84pNzfbI+Xd03S8030tB3nP0t/8FqM6NSxDE5ZQk58VMQBOlkO8HKLT5mLF2b1UzNsfdJnDqU21JuqIvdAnQVWGfrzXWQx9tNLmVVugqv1ICbvZwrUoLrC3JcC+7OJ5dArNxwLtMTCvdCMMIHxyr+vG2Zuzo4tOURJs76AlkFC4LTkMUzKd5rwKAdz4C1F4N2PDpVYnDq9kLQRDcc1q7z06uj6QSJhgnkXnmTrGOdXQva41q24v712JMoyFnw0S7S9g7XhryIAXf1SbkVXIU87ajVcR45lvi5SKjcCsHIrRBjCNFqpf7fz5NVuJCcSUuDWrcgCKTr8mWVDeYKOmPG0m2q2EVWwQKfHd7Olt7ut6af3fG5ELkTPilr1llcnfePGOH3kgDf0zm9WfD2c9aVgPa40nGc8z5fiLR11dzhj4shmGvfyZ3p5YyvK0qcS9ivR2/5UUSLmfgr1sGJMDcqSAR1NeBwrqdms1mw2SxucyE4rzlmFy1Ps7JccV4nTM4xwbYeXcPPPGFvZzDyCvjrK70t68OIXmgy2AQy2yrU4WJSbQvWOoC+4im0UWqb1diHqaGG/trKUDfNLcHWtahdgt21o2YVLKCr+ZRkWVOJyfFPap8Unvyc7o7xJlB2P7EvSC1aKYXdXeK6fLrzd/dnpF57XOtzm8e62EaL68BZWKVyMIT6e3gTduf4cNeZmABWUz+dn+1EodWh1MePMJwiKaOaLwTdvRCuSIa+Qi3WqkHZ5eUIsKtFKyU+W0+dnVyATOvWl2Qs7pAzUOgu3Mx5n7e2OLtI/JlxNloEEr3ga74N+/+RNGDmy4olcpLhhLIt7n5v1+29p45iqq8macZ8NOmZ2CxmBIUCQXG2/lBfg1C8vY8Zn+64/kL2d79GW/0R0nKHfLOOzoX0RfY2COZNhKVelXxJWBOM13/XsDA54ui8P9rDrUZLbO04D5JFgtj6Q6gsxGDXmzxjHvGTplDx/+6jdev79Bw9wGBbM5Nu/znqFMOIc0fL9QiJeyEcvl2VWkvOxMW0VO93W8adi8EbUq/4vuQssA+ABeOVTs4EDme8uTRighvDjr+iGexXfdFmpfaVp+jYsw2lLp7E0tm079yCfsLQ9W764M2AzyFn8DJUOhaVlq67gQfl+bPoeGYrVsugIx0bSL+Ka49r2X18us+JYbxFAniydp3b4O9AWzDEwptYufuOwcxC5im8bjQId9KbaCbUvlRjbSW9J4/Qe/II3UcP0l8ztMRP5/4dAMRPnOK2XdFg7YZMdEfbt9tWqiSObNTJBvq66klKKxhRxjVxRlup0q3FahdPZwFyFVRnv6i38CN//VxSDwznjiUlkHYXQ6is2GBM+XXXvnCkl/SFaLipQ4m/gms19qHUxzs+2yxmjFXlgIDuZBcWs4kBYzs2m4Xu1jOo1Dri9Kmk6CZQxxmUKi3ZV36ZuleeJr5wsqz2BXKtQvm2HpWWrhRpR620FFsZbG9Bn5Qtud+O8+uQHB+tc2iZnAEoKStK7lRhd9/NuTPJ6fj+REnY2+htBYtgCaKnB4PzdqnzjbZbJCa2gVm3p/70CwyLLiBhygy6Duyi58RhdHGpqDR6mnuasVrNpGQWE5+Ug+GiS4mfOBlBqcLc1YEhXSCpdA5dh3YDYO3vA9I9nzCCCanojra1m3SwB6WgQqV2/9ro7RXfneBITRV2ZdjAHPJFwdNrrrsZZ57q8CbwnnI6eBs8DCZyhdOfBwj4v8acFJEUqTCaBNOV0L5zCx37tqNKSEYVn8jAQC99LfVnz1V3iLhFc0gqLnVsUyenkrnqEgC02VeQtfYKn9rubU06KUI9JjVmLN22UiUnfnefZAJzbykPvYWwyJmq620qsdR57bhLJu68TU7nl3KZuNbpLWLDk/AG85U/VJZqKATXjuvbUiQIcLBEMVTfRRRFGt99zfFZm5OPQqVC6DDSZ2oHQKmPJ3vD1SQUTUVQjhlJckvIv+FoWLttpUqspn5Eq4UJX/0OdAxfIdgZZ0FLO2qFo2rgrBXpyS3gGpMrlalLziCVtwQ0/uJ6Azp/V0/Tll3Lumu3ffpyMJbYkRsn7CuhEFuI3HSNga6BZq/DlWCJubHmDF2ffTpsmzolDXVBEemTriEuMyeouWq9EW4rF6LQ0nUVUkfnEAQEtYb0Oi3oh5fz1jHdiZUrrj5eT0gNxHmqWwp3uX7dnVtOOJncgSo5PulgJBkPtvi6Dnza/47hHm+TJgJBl5uPYckq2nd8SNH3f40qfnQyeUHkPBhdGRXRDZW169w5lHFaMi64iIMf/YWpS24mMXX4kj2uwmsqMX1u5Y7E2yQIqSgH1+Pd7XOHL/G3cgTYm5h5y6fr6RzhTMUoF+c2ByM0LFij4pFMMMRWFG0Yq84QXzDUvxQqNfEFxbTv+HBY9EI4iAQrF6Is94K3xCCGhSvIvPxqju14igFjp9f6pAZH7H/brVlvN6w7UfQ314EUcvMvBCMvrrc6I52YVRs+6v/zIifu+RE1zz02bHvNC48DuE0OHgrCsdy9XEbtVwj2U8TdUzlx8nTiZ8yiuXqfT/XImS3mzV8olbRjwaVHZLXDU7YldzPpPAmisyXr7gHizoKX+8AJJpEalxsNjKbl7Zo8ytRcT8vHGzHWVNB9aA+p5y2j+Ef3OMp3HhwK80qZs2jU2ujP7zGas2ijzqcrB6VWj61bOvmNazyqlLD4suqEM/74bV2PdVe33GP9OU7qYSFFqNdBi+Q8EJHuUgjGgJov2AbNVD//D4wVJwEQzUP3myrZgEITB0B/fQ2Nb71E4tTZZF181ai0y58QsdFmVN0LgTxNXF0Lnm4CpT6egf4OYPhECHfHyOmsvqSFdMYuVL6Kldw8EXIH2KTqk+s+ibZX9lC0ORrSCI7Gg8HeX6pv/Y1DcBOnziJt6Wo06Vm0bP4Ppx+9h5aP3qXmxcfJXHM5eV+8aVRcC/4K7mjniokqn65cEoqn0dlyGlG0ObZJJeXwljpRTid2J3rOr+jOSWcCGeDxduyIVz+nv725FmBkp3U35TlURKqVaycahNcdzgmXAkm8ZGqoQXtci+AkHSnfuwzzHAXpK9ajycgmcfI02rZvJnX+MgwLV8ir141hIJdIfxNxZtRF15+niq+Z9bPrk9AlZHB85z9JOdTv2O46SOZMIDeUnCxjwQhfcufXlfo+3nzAUh08GAl5AiFYrotA/dHRdAN7wpPA+tPfe08fo/Lvf8Y6YEI3Z2T+g8SSmcRl5NBz/BDjrv0G6cvXyqo32t6kAiWqLV13N4cgKJi27BsoVRrKtv+NgdYmwHtHk9ovFc3gaaDJnQA6z+iScje4rvjg+rfrOaTqd1fWH+s6FJEQ3jgXoy2CTahy/VqNfdS++DcAeo4dwrj77ABx+9NDqRYFQSB+4mQsPd3o8gpkuRSCIbjuvmskuhYgTANpozFLTaFQUTzvGmpPbqHi8T8w8Vt3ojF4T5IhNb3T0yCFfWDO3Wu81PRbKeyC48tCmK71y0nE7qkud+VCHVkgN9WkO7zNpvOVaHYjgHfB9ef7Wfp6HH83vv2S429FvI6+nYdILVwEk0qwDQwAQ1nEvJ0l3IIbLsZk9IIdQaFk3JRVVB/dyJnH7mXcNV+n9uW/DyuTNPM8ci/7ssd6EvfV0VY6fKlmXzuMq7BsPVU0bHkfV/eDq9B5O5+zz1iuoHrDXo+3jF+hwl+rPIZ7/H2gxGVkM/7Lt1Lz4uPDtqdedzFtf/sXbZ98gCYtk+YP/g2AOjE50KaGnHAtpBsV7oVAnlqCIJBVsABgmOCqDRko9fFoUtNGHCPVMQuerR722Z17wXVAQKqMO4vUblV6citIDYi52yblBnH+29OAir9haL4S6YNncoh2yxjkRXzET5qCYdEFACh0egDU9Uqy1l9Jf00FVU8/CMCEm78n63xS+NLXotHKhTBauqOZ9rFo7peYOOsL7Hr7V9isZgSlkonf+j+fw1h8jYX0dUZaKK3IQIQzHAMd4bBeXbOISW13xR/BDVeGMn8T6TuTufoylDo9LR+969hmNfZ9/pfAhG/8AF32OI91RILghsvKhSh3L/jS4RVKNeNLVlNV9h6i1Urf6RMkFE2VFFJTiYk2tMMurD2Cwl7eNfGNuwQr9u1yErC4Wn3a41o3S2oOR8qqlpNtS+5qwqFcgSJceHuA2vd7E0cpoZazUnCwhdddVj3XtrlDjuCZO9sRbVaMVacd2/rOnKDn2EEA0patjgrBDTdhFd1QJsKRujAZ4+dQd+pj5lz4Q7onpQwr60l4pSZm1JVIr53mitT+UFlxUnG1djwNsjn7gbXHtY7v54r2uHZYcnbX3LzBtNSDlbdCLlL9xVPWLVdhdbeyh9wk9IHia1ilP5x+5LcjttlX5TXcsAH9+vMw4f66RcrAWTitXIhyS9dXNLoUDDnTOPrpU6TXzEKzbjFKXfwI4XWOSKgrOStEzhc2jzzJtcuct/kjGp6sVk/C7klwnbe5rufmXL+doUToStmrX4TKYg218Dpfr0AsTzliGk6/b7DOXfLzP2Hu6qBp4+v0nioDoK/tDPFLZ5OwaoGjnNR4gJw3L094ujbRYuHaEURRdLtz/vz54t69e0PeCDnWrtST3NvroRSiaKO19hAt1fvQ6JJIveFayXJSU2ZdhbfyxuERDcHKR+rvK5inGWTurFApl4ZUfXKjItytTCHH8pZDoCLsbaKA64PTnYXqy7WW66oINqEU+pa6j2l7+k30C6bTf+AE+U/8EkHl3/m8XdNgC+5oWLqCIOwTRXG+1L6osXS9+azkIggKMsbP/lx0U4bt8/QKaHc3QJ7XC+1cj5S15s+rpqeOKSWG/r7uexNWb8v9yFloMpy4ip/rdfD2GXybeOA6WWG0xDeUgmsqMZFYspDECxdi7eqlZve9dL37CSmXyZvy68xYFFxvRETImL8/hLsLIqdTD5i6SEgZuQqFp2OdBalnXp5kHl6p8qYSU8iWe3Enks7iaM//4IxzPohoIRiDLYH+7nKODzS/QaCEWnBtg2ZE69B36z86NKjW+comn+sabcGNFKLG0gVpa9efjt1WqsQQdykn33yeKYIC26ppDvGyD55JDSTZjwVpi9vdiLX2uJa2UvmvpHL8YL4MSkiVlZrB5smFIDVRQu6xvhIsP66Uq0CqjB1v18Xb8XLbFEpGw3dc94MHUKYmMVjdCNaz38fS0Y0qNSko5wiF4EaClQsRYumGg4RJJYy7+uucOPASpqazy0C7JvL2hichda3D1xhf18+hjJd1jrKQk5HMeZ+3CI5wIPVb2zPNudvniWhIghNqwbX3P82kcQxW1A0TXIDab/9edl3eMvy5I5otXDsRI7pyn0LB+NHtF1U3roDM1ZdT+8rf4eN2x/66K81+1euwlmXcwN5cE3bsM9Ts7gLXwTKpJDpyfLmuFqs/g3fuEv74S6DZweQSadapL4RLdJz7R+a3v0zeH3/IuIf+D93sKcPKWVo7AjpPqAQ3UqxciDL3QjCxv3omz5yPzTxA3T//SkbLFWiuKh1Wzttrs7sQMXeDZd6WB3I9l5wBMU+Zyeyf7Sv3evLjyknU7ksEhSdGQ1zdPfx8sQgjTXCd3WujFYZmv6aiKFJ1488c21OuXotuehH9B05guGED7c+9g+lUNQnpqR7rc3ftx7qFaydiLF0YXWsXzgpfkW4RRfOuoeHNF1DsNHqNFPC1s9s7rb8rV4R60MtTqkpXAvUzj5Y1O9YIVvROIAiCQMqXzubI7X5nGz1b95P8hVVYe40AaEsK/ar7XBFciGJL11sn9HahXI9NzZpCzqSlVD/xJxZe8uth+9qQtnadw8KcGRo4k7ZyfRm4cYfzoJazJSuFfbl0+/+uA2jBHACTIpoFNlxWrhyBHc3JFvY+YjOa6Hz1fcf21C9fRNyUAurv+BPxS2eTcvU6rwNpUv0h1IIbSa4FiGLRhZEXxN5R/blQgqBgfMkaGit20loiICgUw2Msj6qBs53dnn9BDlKDae46muuUXOd8tu7iYt25H+xC6wlP30HqvL6K6GhP54123PXdcFq5phITlo5uFNo4kq9YRdcbQw/4tidfR0jUU/DcvbLq8VVwxyoRJ7qB5GPwRWydy9o7tFIdh81qRr+niZ5p8bRM1pFxcrgHxl20gq/C4slS8Uek7ALrbPHa/3YVX7uY7n5rumNboDl4g5HIJJIIhxhE8mu0PTJBmZKIMi0Fa1snAOMf+r+QnXMsWrkQgaIbDs6+zgmkLVvDoY8eRrFVidVqoT4hA6VaS87EJRhypyMIwrAcDZ4IVpKTraeKHBnH3M0+gyExdRZYKZeDs2XsSSi9xeU6M5YENxIsr0hyLYhWC10PD018SL1hA0nrliAIApbOHpQJOgSVPAkJh1shUomI3AtSjFauXTvOHd3Y3UScPgVTbyt93Y3EnWnjlGkv6XkzySpcRM+MRDLOqEfk45VyHYz2CLNUzgU5y+F4y7Pgq5UbE1zfsIuMHDfCaPUpS0835Q/dRfzEKeR84Xosc/2z0cIpuOGydMdE7oXRwmaz8tnmB4ChHLw2qxkBBctyb+T4wBEOb/0r1g8HqNWlMGvld+icqRt2vLvpvoHcKHIsZrtQui71LjdTmDvhlSOeUseGCl/PI2e9sGiwbkcLm8VCz9HPGGhppK/iFAA5l1+HSp+ABd8GXMMdGhaJrgWIie4IFAolS6/8A83V+1Ao1egS0tHGp2NTaZjMdEen2Fr3DLve/iV5Jy5g/NTVKBRDP6W7gbJAXA2BJMeRmlLsbP16Srzui3hGo2ULkeFOiBSMNRU0vPkCGkM6+gmTSCyZSeq8JajiE32q51yebSaHiBXd0VzOR4rM/Hlu9/UMtpGpn0iyJpOa7mr2v38/xfOvJTl9osc63SVLD5Zgya3HVXCjiWAkvfEXX1wAoSCUboW+ynLqX/8n2ZdcTeLk6d4PcCIYbxzBFtxItXIhwiZHhJPEfXWSF95mNXNq78sYe5qBoRuuo1BDZfc+mo1nWKBZw8RZV3Bi17M0Ve72eh7XmWjBmuHlDV8sV9dy4VgjLRCCneHLtW+EwyILpeD211VR89xjpK9Y71FwpXJyBCMB+bli4dqJWEsXwmPtuoaSCYICq9WM2dQDiZkAqLWJTFlwA7n1Onrm5WEgjxkrbuPwx48yI60QXWKGx3PIvYFCIXZSg2Sj5Y8NFWPdRRDSVI1N9bTtGIpySZ61wEtp+f1DzjU518TWTkSLbrixd4qShTcM297ReJzqoxtJy7kO4fNtuoQM0vJm0tZQxrjEC/w6n6e0iq6+10ATy7h+DnQ5lXONSJiWGwhWUz+dn31K+84tJE2fR/EPf4tCZviXJ+Q+AEMpuJHsWoCY6PrFuCkrUdW2oRLicO5imrhEqsreI6/ofASFf/kZpJASvWD6goMh5KNBtLQzUhGtVgbbW7D09WD8aDttdYfIu/pmn324rvjypnGuWrfORLzohntATQpBUJCz+mpcu1ru5AuoP/0Jfd2NJKQMX7jSW8f0tDaZt/Ku+JoZLFpwZ5GHOuyrZ97IJZpCbeUG26Uw2N7KmceGpusqlBpyi85n/vqf0js53e86I1FsA7Vyt27dSl9fHxdddFGQWjSS2EBaEFEq1YybsopTe1/CanWfk9e5s0olJ3eXFFyOj9ddekZvlnS0CHNs/r5/WE1DWcAQFMy/6GdMmLaeOL3nFIye8MWNEC3WrdFoZMWKFWzZsiWk54l4SxfCb+364r+Lu2wlPHeMupMfMV6xGkFQSFotup319M1Nc3yW8ts648+rtWt9Y/H1PByCG21Wrrmni5b/vk18USnjvvQ1upWB1T/a4V+jhV6vZ2BgAI1GE9LzRIXoRgttpUoEIHvD1Zx57Hdo9QYyJ8yX7KQHP/oLts2DJK5bQupVa1Do4uif0o/uhG5EWamJC+5wFlp/JjnEkGa0Bs2CLbiiKFL9z0ewDQww8X9/ijCGBTcYA2ihFlyIia5sfLF21anpgIAo2hzbXP2ORXOvovzQ61BjpP6HD6HNzqP31DHSlq0mXVyHIAgSNcu3UqNJaAMJWYuUabyRSn9NBeaONiZ955co4/zvE7EQsOARNaIbbhcDeJ+R5JxnYdw1X6firZewzM0nrz17RKc1djeRN3E5+SVrqF1QQ1/FKQyLVlL32tNYjX3EFxYTP3EKCk1cVAmorwTLlRLthCoWt+HNFxBUKtRJKX7XEQ0PtUgPE3MmNpDmB3IGBxKKS8m88FKaNr0+Yl9bqRKr2eSwZsd1j2dK2irG905gxuJvktATR9N/XqX6ucd8apecFYNDvapwIHhakdgdcgTB3QrAY53++mrMXe3kXnGjX8f7MrMvZuXKJya6AeCuo9k7anzRVAaa6mmp+WyEUMetXEDD6U+w2SzDjo1PzmXC9IsYN2UliCNF0t9JDO7ENtwi7Mv6bIESqcIbqnZ17N0OQOKUGT4d5+s06pjg+kZUiW4kvkJ46nCq+ATyb7qd8v2v0jHTMKy8Nmc8Gm0y4s4DkscqlBrUKQbHZ2dxlCO83sRUSswj1QJ2h1xhCPfrcWfTSarK3nN8FkWRwf5uRNEW0geBzdSPUqcPWf0QGYIbibrgiajx6UYbZ2/0bJr0Bsw794J+Ij3z8hzZxhLTCmhvr2XivpwRfuKMcbOp+/AjrO9UET9pis/5EvxdcDKa1jTzZRAtFMI72N9Nf28LiWkFKBRKRNGGqbeNvq56zAN9xCfnYLUOUnn4LYzdTWROmI+pr52G8k/oaDpO/KQS1DVpdH5ukU76zi/d+l798XOnn7+W+jef9+k7xdwJoSfqRDcSBtRc8RbZkF+6loP7XmHueXdiD0hJO2qlV6XFJvZLHqPS6Jg46wpOf/BvCif+2C8RlcrXEG3WrDfCFb1gGTRy4KOHUCo1WMxGElPzEUUrvR21JBomgCBQeeRt1HEJDBg7ANj//v3D6rCZzcOiVKqefghLTxcpc5eQffFVI84piiKizSp7inlcZjbmznaspn6U2pGhiK7EBHd0iDrRjVQ8CW963kyqj76PsbsRjfbsSg7xKbnUV2xkQuIsyeNSs0sQDryOflcD/YtyQ2KBylnOZyxis1lortpLRv48lEq1z8f397VhNnUTlzIOU18rHU3HHfvGTV3Nyd0vkFu0nNSsKbTWH0ITl4g+KQfLkokjXvmz1l2JaLPRV3GC2hf/Ruf+HaSvWOdIHi5arfQ+uY32nVvAZiNl/lLSllzoVUgFpQp9YTF9p4+TNG2Ox7LRKrjR5lqAmOiOGurCfGq695OcMclh3aTlzaSrpZzdtf9ixqdrSdJkDBNuQRBIHz+H8vrN5HGTx/o9rUwxlleAkCMWFrOJzuaTqFRaUrImI4oip/e/RnP1PjqbT6HRJqFQKEnJnExyZhGCoEC0Wak6upH+nmayChZhyJk6rM7E1PGcd/Gv6O9pwmoZ4NinTwOQW7Qcc383A8Z2JkxbD0BSeiEwZJW7s1EFhYKESVOZ8tMHOPPovZT/+VeMu+YbdB3ZR3/VaUREcq+4gbj0bFo+epeaFx8n/4bbUKg9B/PHFxTTV3HSq+h6I9LENpqJ2IUpvRFpLgY77qzdpglGal/6G4I6juyLvkhuy1DOXVG00Vixk5pj7zMzeTVWmxmzbYBxCUMrD1ttZra1v0zKygvRXjtkEXtae0rOwMxYSuMoR3RrT26h6sg7ACSlT8Ri7sfY1QBAfun6z9fCG6St7jBW6yAJKXlYBo10Ng+tEabWxJM+bjbxqXlkTTjPUa9+bxUV6U3ok7I4tOURUjKLmbzgBmpPfERfZy3Tz78V8D06oe/0CWpefNzxOWnWAnIvvdbxWRRFGv79HKLNRu6VN45YINUZc1cHFU/+kUm3/9zt5IhInmXmjUi1dD0tTBkT3RAgJbxtpUpEm43O/Tto2fIeunGF5FxyDVnVQ6+aXa1nOLnnBTQWNb3mVqakLKMweeiatZvq2Nf9DgW3/AjrAun1z+SKpTefrlSSnUgUYl/8uF2tZziy9a+Oz9OW3UpyxsQRYpWwt5buwWb6LB0MWvsRUNJoPEF+4kxOW4/Q21HDeRf9ArU2kc6mE1gsJk7uPjtQJSiUZBcuJq94BXH6FMd2fyMUTA21KOMTJAfXbBYzNc8/jiohkexLrvE426zuX8+gy5uAYdEFkvujdd2ySBVciIluSDjR8QkD1j5mpq+T3O9OeGEogXTzqQ8wbj3C+OtuJachBQBTXzsHP/oLlsE+AJbl3kSC2kBNzxHK2jejHVdA9u+/6bZNcsRRzkCau3oiYaaYr4NmVquZPe/cjdVionDmF8gtWjqijBxh6Rpo5tPGF1Ao1SSkjmewvxurxYR5oJeSRV+hre4QLTWfOcpP+ekffM6p7Cs2i5mmTW/QX32GvKu+SlxGtmQ5U1MdNS8+waT//dkId0S0Ci5Er+hGrU833FEMFtsA9X3H3IqupynDSq2OnBmX0WFKperlRxGuvZ3s+mS08QbGl6ym+uhG4tAjikM3RJ2qClVSCgONtZhOVqGdPIHlxeUAwwbBfLFKpY73hC9Jd4JFoFEJxp5mKg7+mwTDeApnXEZ2uRX8FBLFkjksEWfT39tCb0cNqVlTMPV1YDK2kZY7nbTc6ZgSbPQcO4hhyaqQCy6AQqUmZ8PVdB7YTd0z/4/Zq3+IWqMfYVlrs/LQ50+iddv7ZK66xLE93PHLgRDJguuNqBXdcDMldRmZes+r/8LwqAbn3AwAqeedjznTTNXLj2K74hvktmaSVbCAAWMH9eXbUAgqeublkVCWT+eJU6QuWEHzfc+Qe9+3oXjkueQIoWvYmGv0grtVKkYLf4XgzME3Uari0CWkY7NZ6euso7XuIOMmr2Ry90QU5f7V6zqwqU/MRG9fKy8ugUTDeMf+3CtvBPFGBMXozjlKmb2Avh2f0tdVT0pGkeRvOKgu5MyON0icPB3duIKo9uNGOzHR9ROVIo4MXaGsss7CK4oihjILgiDQVqokc8IqNIsSqH71cdTf/CkZpzQkXfsFeFngcNenTOzNYsK09YybsgqFUk15SipN9/+DD269ijXrB4adR+r131XopXAW3nD6cP0V3IH+Lpoqd2OzDjq25RWvYO6aOzAc6QLphG2S+JvC0Z7W05dzBRPRZkWpdB/JoNElA1D1j7+w9Mo/eKwrJrihJaqmAbsSTa8Y9o7cXn+EHW/8GNF2dn57yuersFp6uoAh8Yn/wjqsBi0VB/8DgFKlQRAEihOWYpi5nKY//pNNb3g/r5TYao9r2XqqyK1rwZ1l625Fi3DS39vKZx884BDcorlfYtFl91Aw4xLUcQmOcj3z8mT984dw53SwGvvo72tF97kFLkVq1mRZdUWD4EbTfS9FzNIdZVrrDgEwaOohTp9C2lErNQlVKEQV6tR02gxDppISJZlrv8CZh+/BZrWgUJ69VEX6RaTPMHDikZfJuON6NIV5w15pXS1e+6wtKXHw5tOVSrYTrAG1QH2K3a0VHN56NhPb9OXfIjn9rMvHWyrOsUJf5SnUmnhEm/vfU6FUg6AA0eZ2Vls0CO5YICa6o0zhzMtJNORz5uC/idOnotEmUv15MpSBqWfdBdrjWrJrE6nTp9DTXkVyxqRh9aRkTWbi1Es588DzWHq6ASi45UewwoArcuN3nfF1PTZfBThQwW2rP0L5/tcAUGuTWXDxz4ftP1cEF6C/tgJTXxu737lrxIPHmXlr/499m35HxZF3mDjzsmH7YoI7ekS1ewGi71VDo03AkDud9oYyGk5/woCx07FPtAwJ0fLickwlJtpKlRjWXsyJ3c8NC0eykzF+NgtW/4xFl90DgDopWfKco/H6K0ek7SkDgzFqrtElO0Lr5q2707HdnkIzEHdBtJGtPDuqqtEmDdvnnFJUG29ArU2ioXzbiDLRQrTd71JEvehGI1p9KgsvuZvU7FK62ytZdNk9FP/oHgTVWXG0h3QlTZtD7nVfp+LQf+jrahxRlyAoUKo0aBMyMDXUDttnT9c4GtEHnizdYAktgM1qZvvrd9B4ZicAJYu+gvJz10s4xCPc/ty0o1ZSs0sAGFeyGl3C0JLq7hLtF835EjD0pmAvF2N0GROiGy1PP+cOrtLomLr4q6g18ZTvf5X+usoR5e1iqcubQNL8RVQcepNBU69k3RNnXkb9v5+n662PGawemuK6vLjcId7ulnkPBqM5sKZQqlFp4mmu2g1AWu50YPhvey5YuM4PMkEQmLzgenrOHCZhb61HITXkDAn0mQP/jglumBgTohutCILA5POuR5icR+M7r9Hyq3+S9mwPW08VsfElM9auHmDoBitOW4lCoeL0gX9hHhh6rbbfNIn76sivS2TawpsRT/dS/4vHsLR1Sp7T+WYNlli6s6YDsXA9LbEzf/1PcY7NCqfghmOCgdQ50/NmIiJysPVd9ja9gcVmduxz/U1mLL8Nm82MzWnhVClE0cZnLe+wt+kNWvorg9L2QIgW48obsYG0UcZ1gEejTWCy9gLM58+nvvyTodH4/ToGWxrJuPAS4hevAoZCxgzXXM3pv9zN7vojTD//f2DepGGCk3tGQS7nsX/cIO3PvcNbuyaQcP5clEMZAjGVmKgriayQL29I5ctVKFUICiWizeLTKs1jAXciLwgKZqavZ3v9c4jYaDdVk6mfJPnbJKUXolJrMVo6SVCfHXhtNp5GrdRR1vYhGboCugeb6BxowioO0mqqYkLiHKYaVoTsu50rjBlLN9qegq6vduq4BPJL1zFx1heYUvpFcq+8ic5PPkb18UlHmezaxKEE2UDFobeGLfHuzAzrPLTqLAbP1NLw679iMw5Zoc7uBmdrNxgiPJpCLggKktIKAGgtTfBceBQIps86kHMkqA2ka4cm7OxveYvuubluyypVOqxO1rD9mF2Nr9BrbqWiey9tphoM2jzWT/ge+YmzqOr5jM9a3gn8y5zjxCzdCEIQBIePMglImafmxO7niZ82g5R5S8hryyL3W9/G0tNN1dMP0dtRS5JEhla1UsuCupn0zMtj+847aH3iNTK/d8OIcs6z1ZynBzsLaCCDcHb3QMDhYRL1GHKm0dVSTn35VgpnXBpQ/cHCn9A8ufXKRadOApOAQlAOW5XCGVEUGexpQ5s4/IG1Nv87VPccQKdKIk07AZXibHL3UsNKqnsO0mQ85d+XCJBoM6o8MWYs3WjE21LuqdklzL7weyT2J1D7zGMYu5sQBAXqpBRS5iyi+tgHtJlq3B4/+MnQCH/cpPEeE9wEw+oNdoSElNDY/bxtpUqyChciCArJULpwEkyr158l0EsNK1iX/x2UggpjT7NkWcugEVG0oVEMX8FCISgoSJpLlr5omOACdA00ATAx6TxiBMaYEt1ofRp6El6NNon8qWvIL13H6QOvO7anD+bQ2XScPU3/cntsq6mKtPPXknzpWT/c8uJyFlx6hLorzcOsslC+HnsaFPOXzhlaCmdeSpwuNaj1BgvnmGRffttgxDILgkBRymLKtj1B9bH3sbnMVOtuPYNFHHRrCUvRPTgk4GnafL/bFWOImHshQvA2IJRdsJDKw+9g7e9DqYvHsngiyr1xzDFscHtMgtpAw+nj9N1bzj+PnuHGnV9z3GjLi8vZfXx60L9HIDi7EaREeqT7Q4EuMWO0mhcQoXqouXtgT0icRYomh8+qN6HRJpFduMixr758G7nxUyWPc8f4xBmkafPRqZIw20yAgFoRF0jTZROtxpQ7xpSlG+1I3UB2F0TSZ40kqdIwVp3GVGLCMk9F5oar2Nf6Hwat0q/2OfElZCSXYDp6BoC3/lCDaDubtN4+680bzi4HuwvC+Z9zfaHCtW7RaqH+9DYyx88N2TkjGW+uKYDkuEzmJV1E5eG3HSsSW8wmettrmJ622udz6tXJCILA9vrn+W/NX2nsO4WnRRBiSDPmRHcsPRWdbyqTpZfOgQaU+gRHFELcl6aRMGUaDX3HJY9XCAqmGUtZeOndFHzjhxg/O87L173PxletfHzsbFpKZwvT1SILVlSCN0vPWfw9PQjs7en8bBdavYEUmdmzxgpyxHZYeU06VssAezfeS92pjzn07h+wiRbAs1gazV2Ud+5k0NoPfL78++cCOzN9PQJKDrS+w+aav8aE10fGnOhGO+5uqEbjUOiYuDjFsW15cTkJ1y3jZO9ONlY9yIDVKHls6qF2Sj6wMn/Wd0hXTqL5gWdofvA5nwXVU3lvg3DOft1g+Hg79+9g3JRVAdURTfgqts4UJy9Go9AzcPoU/dZu4lUGFIJnz6LR0kl5104OtW6iw1THpuqH2Fr/DwAM2jzWTfg2C7KuwioO0jU4cnp6sBhLRpSdmOhGIFI3V2X30Ch996btw1wEcYV5JC9YikKpZtDa57HepP31zBicgy4hA1vj8LLOQuirxeuLeMu1aD1h7u7E0tPtWNp8rOOv2NpEG3W9R0mJyyVenYpNtJKuLWBZ7o1ej20z1ZCpm0ivuY0Bq5EZaeuYl3n5sDKNxqGIGJ0qmV5zu9cZbjGGiIluhOJ6o01PWwNAz7s7eOWm/zLY47RqhNVKXvEKEjXeB5UEQWBp6pewtndQ++pTdB89gGg7e7PIXcJ9NELE3GGsKkc/YZLHpcejHbtl66/gVnUf5P3qv3C47X1q+8roNbdSnLqU+Vlf8Bq10NZfQ0X3XuKU8UxOXUp1z0HyEqYOm70GkKUfSjeqVmj5pP6f7G16Xaq6GC5E7WrA3gj3asGhoKr7AMc6tgCgzs1g9d0LOaZdQs/H+7D9p4ysggUMmnrIaFSjVSWgUyW5rcto7mJr/dMA5F55E0mlsx37nAXQXRRBKHy9ngTfOXKh4e2X0WblUhS/JChtcMWb0NmjTII9BTmYCWha+isRRZF0bT4KHxfJtIkWDrZupMlYjlLQkBqXw7xMz2Jd2b2f4x1bmZd5uexlrLwRza6FMbkEuzfGoujC0A3RUCjSZuimacubaKcXYbh2Pe0//Qd9XfUACChJ0WSzMOdLHuuyihb2Nb9Jl6KDrPVXOoRXrhAGgjvLVup8x34z/FquXLmSn//851x44YUjyi774gOAdBJzV5Ec7Sxb4Ty3Pxxv30Zlzz4ACpPmMSX1fI/lN1Y9CMD6Cd8DhmJ7W/ormZS8wK/zj1XRjcXpRhkKQUVeJeRVptGx7A7OHHiDlnueJe+mr9NXcZKE4mk0P/sMan2K17qUgooFWV+kc6CBfVveY6CxjoxVG4YlmfE2tdVVEEeDmpoaxo8fL7nvk3/9KOTn9/eBHg1C68zk1CW0D9TQM9jqmJHmiYVZV6N0mslW03OENlO136I7Vhmzli6MXWvXGVEU2Vd4ko6yT0n90lrSDAux9RupfPT3zEi6UNYy8QCD1n52dL7OE4/9kRtv9D7QEk4SExMpLy8nKysr3E2RTbT2xZqeI5S1byY3fioz09eN2G+y9FLXW0Zh8nwEBOr7TpAbX4IgCJ+Hpg3lgfCVaLZyIWbpjmkEQWB+5RRqSsZx6rW3UM63UaRbiHL2Fynb+y8UgpJ03QRZnfjw4f/lwgsvJDMzk3XrRt5gkYDVasVkMmEwjFwLLpKR8/tHojCPS5hGum4CcUq95P5+Szenuj7lVNenLM/9KofbNmG0dFKcsthrWNq5SuxXGSMc2fYEBw4cYNWqVWyufxatVstNN4ksXbqUW2+9VVYdM2bM4NVXX+XLX/4yBw8eJCMj8qbYdnZ2kpCQgFqt9l44yvAmzN5E2Wqz0D3YjNnWT5wygXh1KqJoQxAUqBQav9okCAI6VaLb/anaXFaNu5XjHVtRKeKIU8ZzumsXxSmL/TrfucCYdi9AZFoP/iLHWlqyZAk/+9nP+Pjjj3nhhRf49NNP3fo/3fHrX/+aZ599lg8//JD8/MhKcNLb20tWVhZ9fZ5jks8FysrK2LZtG1VVVfzjL/+izVSDXpWMRqnHZOmlz9KOTbQioCBOGY9VNJOhKyQ/cRZKQUmCOt2npDdyON21m1OdOwBYl/9dv+qPdtcCxNwLUUcgne5///d/ueSSS1i3bh2HDh3y6zX8V7/6FUlJSaxdu5a3336boqKR6SDDRVtbG0qlElEUgy4Y0UBnZycvvfQSTz31FPX19axfv56CggLuefhOLrjgAiZMmIBCoUAQBMf03CW5X0aj0GOcmUt9+VbKqjZjtg1gsvaQoSskJ34KufElQWnfpOQFpGnHs7PxZayimQ5TPU3GcsYnzEKnSkSjjJ5VS0JFTHTDSCie6Ndffz1XXXUVcXGBZYD6/ve/j1qt5oILLuDAgQOkp6cHqYWB8frrr3Pttdeek4L7yCOP8Itf/II1a9Zw9913s2bNGoxGI7/61a+46667EAQBq9WKKIokJSXxxBNPcM011/Bpw0sj6lo34bvU9R7lWPtHtPRXoFHoSNdNwGjuos1UTZwygUy9f/G2KXE5rMm/HaWgQkCgtvcItb1HyNIXMyfDfVa8c4UxL7rvVf457C6G0X5dClRw7dx+++0cOnSIhx56iN/85jdBqTNQEhMT6e/vD3czRp2///3vPPzww+zfv5/CwiExPHr0KNdddx2lpaU0NDSQkZHBwMAAs2fP5vjx45w5M5RdzmKxoFINv9U3VT0EwI4dO7jqqqtILOlElwK73t1Nx8DZ0LalOdfLmunoivLzQbR03QSmG9ZwpP2Dz6MZYMDax0e1fwNgXuYXyNAV+Fx/NDPmRXe0GQv+KGeuvfZa7rrrrnA3w8HixYu5//77w92MUUUURd5++23uuOMO9u/fzw9/+EMqKiqoqanh17/+NbfddpvD8tdoNBQUFFBRUcHGjRs5fPgwL774Itu2bSM7O5uysjJWrlxJUtLQbMUlS5Zw6tQpNm/ejMViYfny5bz66qtUVlZSX1/P9obnmZ95BWnafL/eLroGGjnS/gHJmizmZAwtq6RSxFGUvIhG40n6Ld3Dyo+1+0eKMT+QBqEdTBvrnWTbtm3cdtttHDx4EIUi/LkOBgYG0Gq1PPfcc1x//fXhbk7I+ctf/sJ3v/tdAObMmUNlZSX33nsvpaWlFBUV0dHRwZNPPsnChQu59tprHccdO3aM5uZmNm3axM6dO6mpqUEURU6fPk1paSlbtmzxGJ2ybds2vvWtb1FWVgbAU089xde+9jXHfrn3lNVmwWTtIV7tfoUPq2jBJlo41fkpdz/8fa6++mp0Op2s+iOVc3IasCvBEN6xLrAAr7zyCrm5uSxbtgwYejVdtmwZ//M//8NXv/rV8Dbuc+wW1+nTp5k4Ud7kj3DS0dFBfHw8Gs1Q2JbVaqW3t5fk5GSvx06ePJlTp4YWg7zuuuu46qqruOKKKxz7L7roIhQKBQcOHODpp59m7dq1Huuz2Wz8+Mc/5sknn+TCCy/kkksuYf369eTk5GA2m+np6XEMvppMJu68804eemjIFXHzzTfzwAMPkJo6UkADub+qeg5yrP0jktSZdJubWbt2LZs2bfK7vkggJrr41ynOBZF1ZcGCBXR2dnLy5Nml37dv386GDRt46KGH+MpXvhLG1g3xi1/8gt/+9rcAEZdAWxRFmpqa2LlzJw8++CAnT56kqakJnU7HzJkz6e3t5fDhw8CQn/bmm2+WVWdnZ+cIsautrXWEA86cOZPly5fz8MMPy2pnR0cH//73v3nvvffYvHkzSUlJNDc3IwgCd911F3fccYejbH9/P7/97W+59957Adi/fz9z5szxeg5P99zJju1YbAOUpq3CYhukz9zBJ3XP0djYiCAI5OTkyPoekUpMdJEnuueiyLrS0dFBZ2enY7DGziOPPMIHH3zAm2++GaaWncVoNBIfHw+EX3Tb29t55pln2LNnDydPnuTkyZPExcUxY8YMvvnNb7Js2TJSU1Ox2Wzs2bOHd955h+eee44//elP3Hjjjbz11lts2OD/iP4DDzzA5s2bWb16Nbfccoss69kVi8XCmTNnGDduHJ2dnUydOpXa2loSE89OihBFkccff5zDhw/z2GOPsWHDBpYuXcqll17K9Ony1tq7qOD7mCw9gMCRts2kxGVTlHJ27baxdP/FRBdp0R1LFznUPPjgg+zZs4fnn38+3E0B4M477+Tpp5+mqcl7IpZQIIoiv/vd7/jDH/7ApZdeypo1a5g8eTLFxcWyY6Pz8/O59NJLefTRR0PcWt9YsmQJP//5z7n44otH7LPZbJSXl7Np0ybOnDnDSy+9RHJyMmlpaTz44IOcd577Jdqfeuopvv71rzs+z0xfPyw+eCzdj36LriAIY0NxY8SIEWN0aRVFcb3UDo+iGyNGjBgxgkv4Y4BixIgR4xwiJroxYsSIMYrERDdGjBgxRpGY6MaIESPGKBIT3RgxYsQYRf4/6STveJuaWQYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.axes(projection=projection)\n",
    "ax.contourf(std_dev.longitude, std_dev.latitude, std_dev)\n",
    "ax.coastlines()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3c15304-fb95-4c3b-b516-ffa7ea855c5d",
   "metadata": {},
   "source": [
    "**Troubleshoot:** In my Jupyter notebook setup, the pyproj package gives the following error when it tries looking up the projection info:\n",
    "\n",
    "```\n",
    "CRSError: Invalid datum string: urn:ogc:def:datum:EPSG::6326: (Internal Proj Error: proj_create: SQLite error on SELECT name, ellipsoid_auth_name, ellipsoid_code, prime_meridian_auth_name, prime_meridian_code, publication_date, frame_reference_epoch, deprecated FROM geodetic_datum WHERE auth_name = ? AND code = ?: no such column: publication_date)\n",
    "```\n",
    "\n",
    "I believe this is because I run Jupyter from a conda environment that's different than the kernel Jupyter is using. In any case, there's an easy fix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "32a25d6e-2645-405c-b7d2-621645e7191c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyproj\n",
    "pyproj.datadir.set_data_dir(\"/Users/<me>/.conda/envs/<this notebook's kernel env>/share/proj\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "gefs",
   "language": "python",
   "name": "gefs"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}