diff --git a/python/examples/kerchunk_hrrr_subhourly.ipynb b/python/examples/kerchunk_hrrr_subhourly.ipynb index c8bab0c..b0ccfd9 100644 --- a/python/examples/kerchunk_hrrr_subhourly.ipynb +++ b/python/examples/kerchunk_hrrr_subhourly.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -14,11 +14,19 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Read 19 HRRR files\n" + ] + } + ], "source": [ - "hrr_subhourly_member_files = fs_read.glob('s3://noaa-hrrr-bdp-pds/hrrr.20230720/conus/hrrr.t23z.wrfsubhf*.grib2')\n", + "hrr_subhourly_member_files = fs_read.glob('s3://noaa-hrrr-bdp-pds/hrrr.20230721/conus/hrrr.t11z.wrfsubhf*.grib2')\n", "\n", "files = sorted(['s3://'+f for f in hrr_subhourly_member_files])\n", "print(f'Read {len(files)} HRRR files')" @@ -26,7 +34,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ @@ -42,7 +50,9 @@ " return f'{json_dir}{date}_{name[0]}_{name[1]}_{name[2]}_message{message_number}.json'\n", "\n", "def gen_json(file_url):\n", - " out = scan_gribberish(file_url, storage_options=so, only_variables=['prate', 'ugrd', 'vgrd', 'tmp'], perserve_dims=['hag'], filter_by_attrs={'statistical_process': ''}) \n", + " out_precip = scan_gribberish(file_url, storage_options=so, only_variables=['prate'], skip=1)\n", + " out_wind = scan_gribberish(file_url, storage_options=so, only_variables=['ugrd', 'vgrd'], filter_by_attrs={'statistical_process': '', 'fixed_surface_value': '10'}, skip=8) \n", + " out = out_precip + out_wind\n", " for i, message in enumerate(out):\n", " out_file_name = make_json_name(file_url, i) # get name\n", " with fs_write.open(out_file_name, \"w\") as f: \n", @@ -51,9 +61,20 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/matthewiannucci/Developer/gribberish/python/examples/env/lib/python3.9/site-packages/distributed/node.py:182: UserWarning: Port 8787 is already in use.\n", + "Perhaps you already have a cluster running?\n", + "Hosting the HTTP server on port 49897 instead\n", + " warnings.warn(\n" + ] + } + ], "source": [ "from dask.distributed import Client, progress\n", "\n", @@ -62,9 +83,24 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "d0c3742a7f0747e9ac961966ca103fb9", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "VBox()" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "futures = client.map(gen_json, files[1:], retries=1)\n", "progress(futures)" @@ -72,7 +108,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ @@ -81,9 +117,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 162 reference files\n" + ] + } + ], "source": [ "from kerchunk.combine import MultiZarrToZarr\n", "\n", @@ -94,7 +138,7 @@ "# combine individual references into single consolidated reference\n", "mzz = MultiZarrToZarr(reference_jsons,\n", " concat_dims = ['time'],\n", - " identical_dims=['x', 'y', 'latitude', 'longitude', 'hag'])\n", + " identical_dims=['x', 'y', 'latitude', 'longitude'])\n", "\n", "d = mzz.translate()\n", "\n", @@ -104,9 +148,944 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:    (y: 1059, x: 1799, time: 72)\n",
+       "Coordinates:\n",
+       "    latitude   (y, x) float64 dask.array<chunksize=(1059, 1799), meta=np.ndarray>\n",
+       "    longitude  (y, x) float64 dask.array<chunksize=(1059, 1799), meta=np.ndarray>\n",
+       "  * time       (time) datetime64[s] 2023-07-21T11:15:00 ... 2023-07-22T05:00:00\n",
+       "  * x          (x) float64 -2.701e+06 -2.698e+06 ... 2.69e+06 2.693e+06\n",
+       "  * y          (y) float64 -1.581e+06 -1.578e+06 ... 1.59e+06 1.593e+06\n",
+       "Data variables:\n",
+       "    prate      (time, y, x) float64 dask.array<chunksize=(1, 1059, 1799), meta=np.ndarray>\n",
+       "    ugrd       (time, y, x) float64 dask.array<chunksize=(1, 1059, 1799), meta=np.ndarray>\n",
+       "    vgrd       (time, y, x) float64 dask.array<chunksize=(1, 1059, 1799), meta=np.ndarray>\n",
+       "Attributes:\n",
+       "    meta:     Generated with gribberishpy
" + ], + "text/plain": [ + "\n", + "Dimensions: (y: 1059, x: 1799, time: 72)\n", + "Coordinates:\n", + " latitude (y, x) float64 dask.array\n", + " longitude (y, x) float64 dask.array\n", + " * time (time) datetime64[s] 2023-07-21T11:15:00 ... 2023-07-22T05:00:00\n", + " * x (x) float64 -2.701e+06 -2.698e+06 ... 2.69e+06 2.693e+06\n", + " * y (y) float64 -1.581e+06 -1.578e+06 ... 1.59e+06 1.593e+06\n", + "Data variables:\n", + " prate (time, y, x) float64 dask.array\n", + " ugrd (time, y, x) float64 dask.array\n", + " vgrd (time, y, x) float64 dask.array\n", + "Attributes:\n", + " meta: Generated with gribberishpy" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "import xarray as xr\n", "\n", @@ -119,17 +1098,17 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "import pyproj\n", - "to_xy = pyproj.Transformer.from_crs('epsg:4326', ds.apcp.crs, always_xy=True).transform" + "to_xy = pyproj.Transformer.from_crs('epsg:4326', ds.prate.crs, always_xy=True).transform" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], "source": [ @@ -139,9 +1118,30 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA04AAAG2CAYAAABI/tbcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxIElEQVR4nO3de3zO9f/H8ee1za6NMQnbMOaQEjlEDp1UVksiHengVKEDylJSMqpv6leJpINUOomUdCCHFinkPDlEFCZsSIzZwXa9f3/ouuyy83XYdY3H/Xbbrfb5vD+f6/157bou1+t6vz+vt8UYYwQAAAAAKFSArzsAAAAAAP6OxAkAAAAAikHiBAAAAADFIHECAAAAgGKQOAEAAABAMUicAAAAAKAYJE4AAAAAUAwSJwAAAAAoBokTAAAAABSDxAk4A02dOlUWi0U7d+702Dl37twpi8WiqVOneuycRenbt69iYmLK5LEgxcTEyGKxyGKxaNCgQcW2X7x4sSwWixYvXuz9zqHUYmJidOONNxbbbvbs2Y6/u8Vi0erVq0v9WGX93lAcb7z/FSUmJkZ9+/Ytk8cC4FskTgCcTJs2TePHjy9R27lz52r06NFe7U958L///U8Wi0XNmjXLt2/BggW677771KxZMwUGBrqVDP75558KCQkp8APukiVL1K1bN0VHRyskJESRkZG6/vrrtXTp0hKf/4orrtDHH3+sPn36uNxHlJ3Nmzdr9OjRbiUIbdq00ccff6wBAwZ4rmN+6M033/SbxA5A+UXiBMBJYYlTvXr1lJGRoV69ejm2zZ07V2PGjCnD3vmfv//+Wy+88IIqVapU4P5p06Zp2rRpCg8PV61atdx6rKFDhyooKKjAfX/88YcCAgL0wAMPaNKkSRo2bJhSUlJ05ZVXat68eSU6f4MGDXTPPffokksucaufKBubN2/WmDFj3Eqc6tSpo3vuuUcdOnTwXMd8rFevXsrIyFC9evUc20icAHgCiROAErFYLAoJCVFgYKCvu+JXhg0bpvbt26tNmzYF7n/hhReUlpampUuXqkWLFi4/zvz58zV//nwNHTq0wP3333+/Zs+eraefflr33Xefhg0bpmXLlqlGjRolHkEsD2w2mzIzM33dDZ/KzMyUzWbzdTf8VmBgoGNkFgA8icQJOEt8/fXX6tKli2rVqiWr1aqGDRvqueeeU25urqPNVVddpTlz5mjXrl2Oex7sU8tOv4+hb9++mjRpkiQ53SMhFX7/S2H3QsyePVvNmjVTSEiImjVrpq+++qrAa7DZbBo/fryaNm2qkJAQRUREaODAgfr333/dD5ALlixZoi+++KLIxKRWrVqqUKGCW49z4sQJPfLII3rkkUfUsGHDEh9XsWJF1ahRQ4cPH3br8f/++291795dlSpVUs2aNTV06FBlZWUV2HbFihW6/vrrFR4erooVK6pjx44FThdcvHix2rRpo5CQEDVs2FDvvPOORo8ene/Drv2eq08//VRNmzaV1Wp1jKDt2bNH9957ryIiImS1WtW0aVO9//77+R4rKytLCQkJatSokaxWq6Kjo/XEE0/ku4aFCxfq8ssvV9WqVRUWFqbzzz9fTz31VKliZb+36JdfflHbtm0VEhKiBg0a6KOPPsrX9q+//tLtt9+uatWqqWLFimrfvr3mzJmTL04Wi0XTp0/XyJEjVbt2bVWsWFGvv/66br/9dknS1Vdf7Xj9nf6aK0k/vOXHH3/UFVdcoUqVKqlq1aq66aab9Pvvvzu1sf/Nt2/frr59+6pq1aoKDw9Xv379dPz4cae2GRkZGjJkiKpXr67KlSurW7du2rNnjywWi9OU4dPvcYqJidGmTZv0008/OeJ01VVXOT3+6Qq6T8oYo+eff1516tRRxYoVdfXVV2vTpk0FXvvhw4f16KOPKjo6WlarVY0aNdJLL71EwguUcwXP+QBwxpk6darCwsIUHx+vsLAw/fjjjxo1apTS0tL08ssvS5KefvppHTlyRH///bdee+01SVJYWFiB5xs4cKD27t2rhQsX6uOPP3a5XwsWLNCtt96qCy+8UGPHjtU///yjfv36qU6dOgU+5tSpU9WvXz8NGTJEO3bs0BtvvKF169Zp6dKlRSYoWVlZOnr0aIn6VL169WLb5ObmavDgwbr//vt10UUXlei8rho/frz+/fdfjRw5UrNmzSqybVpamrKzs3Xw4EF99NFH2rhxY6k//OeVkZGhTp06KTk5WUOGDFGtWrX08ccf68cff8zX9scff1Tnzp3VunVrJSQkKCAgQB988IGuueYa/fzzz2rbtq0kad26dbr++usVFRWlMWPGKDc3V88++6xq1KhRYB9+/PFHff755xo0aJCqV6+umJgYpaamqn379o7EqkaNGvr+++913333KS0tTY8++qikk8l2t27d9Msvv2jAgAFq0qSJNmzYoNdee01//PGHZs+eLUnatGmTbrzxRjVv3lzPPvusrFartm/fXqp7xOy2b9+u2267Tffdd5/69Omj999/X3379lXr1q3VtGlTSVJqaqouvfRSHT9+XEOGDNG5556rDz/8UN26ddMXX3yhm2++2emczz33nIKDgzVs2DBlZWXpuuuu05AhQ/T666/rqaeeUpMmTSTJ8d+S9sNbfvjhB3Xu3FkNGjTQ6NGjlZGRoYkTJ+qyyy7T2rVr893rd8cdd6h+/foaO3as1q5dqylTpqhmzZp66aWXHG369u2rzz//XL169VL79u31008/qUuXLsX2Zfz48Ro8eLDCwsL09NNPS5IiIiJKfU2jRo3S888/rxtuuEE33HCD1q5dq+uuu07Z2dlO7Y4fP66OHTtqz549GjhwoOrWratly5ZpxIgR2rdvX7EjwMeOHSvRqGqFChUUHh5e6usA4AYD4IzzwQcfGElmx44djm3Hjx/P127gwIGmYsWKJjMz07GtS5cupl69evna7tixw0gyH3zwgWPbww8/bAp6G1m0aJGRZBYtWlTsOVq2bGmioqLM4cOHHdsWLFhgJDn14+effzaSzKeffup0znnz5hW4/XT2mJTkpyTeeOMNEx4ebvbv32+MMaZjx46madOmRR5TWGyLsm/fPlO5cmXzzjvvOF3HqlWrCmwfFxfnuI7g4GAzcOBAk5GRUezj1KtXz/Tp0yff9vHjxxtJ5vPPP3dsS09PN40aNXL6G9tsNnPeeeeZuLg4Y7PZHG2PHz9u6tevb6699lrHtq5du5qKFSuaPXv2OLZt27bNBAUF5Yu/JBMQEGA2bdrktP2+++4zUVFR5uDBg07be/bsacLDwx3P948//tgEBASYn3/+2and22+/bSSZpUuXGmOMee2114wkc+DAgeJCVaR69eoZSWbJkiWObfv37zdWq9U89thjjm2PPvqokeTUr6NHj5r69eubmJgYk5uba4w59Vpq0KBBvtfwzJkzC3ydlaYfdsU9r4pS2Ou6Zs2a5p9//nFsW79+vQkICDC9e/d2bEtISDCSzL333ut0zptvvtmce+65jt/XrFljJJlHH33UqV3fvn2NJJOQkJDvWvK+/zVt2tR07NgxX9/tj3+608+xf/9+ExwcbLp06eL0/H7qqaeMJKfXznPPPWcqVapk/vjjD6dzPvnkkyYwMNAkJyfne7y8+vTpU6L3qYKuB4B3MVUPOEuEhoY6/v/o0aM6ePCgrrjiCh0/flxbtmzxSZ/27dunpKQk9enTx+mb02uvvVYXXnihU9uZM2cqPDxc1157rQ4ePOj4ad26tcLCwrRo0aIiHysuLk4LFy4s0U9x/vnnH40aNUrPPPNMoaMknjJ8+HA1aNBA999/f4nav/jii1qwYIHee+89tW/fXtnZ2crJyXH58efOnauoqCjddtttjm0VK1bMV4UtKSlJ27Zt01133aV//vnH8fdJT09Xp06dtGTJEtlsNuXm5uqHH35Q9+7dnYplNGrUSJ07dy6wDx07dnR6Phhj9OWXX6pr164yxjg9H+Li4nTkyBGtXbtW0snnTZMmTXTBBRc4tbvmmmskyfG8qVq1qqSTU1rdnU514YUX6oorrnD8XqNGDZ1//vn666+/HNvmzp2rtm3b6vLLL3dsCwsL04ABA7Rz505t3rzZ6Zx9+vRxeg17qh/eYH9d9+3bV9WqVXNsb968ua699lrNnTs33zEPPPCA0+9XXHGF/vnnH6WlpUmSY3rmQw895NRu8ODBnu5+gX744QdlZ2dr8ODBTlP77CObec2cOVNXXHGFzjnnHKfnXGxsrHJzc7VkyZIiH+uJJ54o0fvUq6++6unLBFAMpuoBZ4lNmzZp5MiR+vHHHx0fRuyOHDnikz7t2rVLknTeeefl23f++ec7PvxK0rZt23TkyBHVrFmzwHPt37+/yMeKiopSVFSUG709ZeTIkapWrZrXP7T9+uuv+vjjj5WYmKiAgJJ9z9WyZUvH/99zzz26+OKL1bdvX33xxRcu9WHXrl1q1KhRvvtAzj//fKfft23bJklFljI/cuSIMjMzlZGRoUaNGuXbX9A2Sapfv77T7wcOHNDhw4c1efJkTZ48ucBj7M+Hbdu26ffffy80wbW369Gjh6ZMmaL7779fTz75pDp16qRbbrlFt912W4ljb1e3bt1828455xyne/F27dqldu3a5Wtnn2q3a9cup/L2p8fAU/3wBvvr+vTniHTy+ubPn6/09HSnSpSn9/Wcc86RJP3777+qUqWKdu3apYCAgHxxKOw542mFvVfVqFHD0Ve7bdu26bfffiv2OVeYCy+8MN8XR644duyYjh075vg9MDDQ61/0AGc6EifgLHD48GF17NhRVapU0bPPPquGDRsqJCREa9eu1fDhwz1+w3Jh1azyFqIoLZvNppo1a+rTTz8tcH9xHwgyMjJKnCBGRkYWum/btm2aPHmyxo8fr7179zq2Z2Zm6sSJE9q5c6eqVKni9E27q5544gldccUVql+/vuMm9YMHD0o6+a1+cnJygR+O7YKDg9WtWze9+OKLysjIKPWIRWnYn0Mvv/yyU/KWV1hYmEsV8U7vt/2x7rnnnkITtebNmzvaXnTRRRo3blyB7aKjox2PsWTJEi1atEhz5szRvHnzNGPGDF1zzTVasGBBqapJFtbWGFPic5zOlb+dN/rhLb7qq7feq6699lo98cQTBe5v3LhxkccfOXJEGRkZxT5OcHBwke8zr7zyitNyEfXq1SuzRYGBMxWJE3AWWLx4sf755x/NmjVLV155pWP7jh078rUtTQnfwtrav4E9vZqb/VtbO/s6K/bRiry2bt3q9HvDhg31ww8/6LLLLnPpQ+SMGTPUr1+/ErUt6sPanj17ZLPZNGTIEA0ZMiTf/vr16+uRRx7xSAnw5ORk7dq1q8DRhm7duik8PLzYinkZGRkyxujo0aMuxa1evXrauHGjjDFOf++C/j6SVKVKFcXGxhZ6vpo1ayokJETbt2/Pt6+gbQWpUaOGKleurNzc3CIfy96v9evXq1OnTsU+twMCAtSpUyd16tRJ48aN0wsvvKCnn35aixYtKvZxSqtevXr5YijJMW027xpEhfHXctv2vhd2fdWrVy903bOizmmz2bRjxw6nUZ+SPmdK8l5ln64pFf1e1aBBA8f2AwcO5BvBa9iwoY4dO+byc+aRRx7Rhx9+WGy7jh075quimFfv3r2dpoJ684sT4GxB4gScBezf5uZNCLKzs/Xmm2/ma1upUqUSj8zYP/yc/qGjXr16CgwM1JIlS9S9e3fH9tMfLyoqSi1bttSHH36oJ5980nGf08KFC7V582anD4933HGH3nzzTT333HN64YUXnM6Tk5OjY8eOOfXhdPZ7nNxVWLn0kSNH6ujRo5owYUKpSobbnThxQn/++afCw8MdUwonT56cryTzjz/+qIkTJ+qVV17RBRdc4Ni+f//+fNMYDx8+rC+//FLR0dGFTnEszg033KAFCxboiy++cJS/Pn78eL4pcq1bt1bDhg31yiuv6K677spXjfHAgQOqUaOGAgMDFRsbq9mzZ2vv3r2O+5y2b9+u77//vkR9CgwM1K233qpp06Zp48aNTlPa8j6WdPJ5M3fuXL377rv57svKyMiQzWZTpUqVdOjQoXzf3ttHzgorve6OG264QePHj9fy5csdi8+mp6dr8uTJiomJKdFUrbyvP3+S93U9YsQIx+ty48aNWrBgge65555SnzMuLk5PP/203nzzTUfFT0maOHFiiY6vVKlSgXGyv1aXLFmibt26STr5dzg9cYmNjVWFChU0ceJEXXfddY5ErKAvSO644w6NHj1a8+fPV1xcnNO+w4cPKywsrNCFrKWTI80lidHpUwRP16BBA6ckL699+/bpyJEjatiwoaMa6ZEjR7Rv3z5FRUU53osLel8CzmYkTsBZ4NJLL9U555yjPn36aMiQIbJYLPr4448LHFlp3bq1ZsyYofj4eF1yySUKCwtT165dCzxv69atJUlDhgxRXFycAgMD1bNnT4WHh+v222/XxIkTZbFY1LBhQ3333XcFzu0fO3asunTpossvv1z33nuvDh06pIkTJ6pp06ZO8/M7duyogQMHauzYsUpKStJ1112nChUqaNu2bZo5c6YmTJjgVMDgdJ66x6l69epOyaCd/QPU6ft+++03ffPNN5JOJgdHjhzR888/L0lq0aKFI7Z79uxRkyZN1KdPH8c6V9ddd12+x7F/+OvYsaPTorudO3dWnTp11K5dO9WsWVPJycn64IMPtHfvXs2YMcPl6+3fv7/eeOMN9e7dW2vWrFFUVJQ+/vhjVaxY0aldQECApkyZos6dO6tp06bq16+fateurT179mjRokWqUqWKvv32W0kn185ZsGCBLrvsMj344IPKzc3VG2+8oWbNmikpKalE/XrxxRe1aNEitWvXTv3799eFF16oQ4cOae3atfrhhx906NAhSVKvXr30+eef64EHHtCiRYt02WWXKTc3V1u2bNHnn3+u+fPnq02bNnr22We1ZMkSdenSRfXq1dP+/fv15ptvqk6dOk7f2nvKk08+qc8++0ydO3fWkCFDVK1aNX344YfasWOHvvzyyxLdV9WyZUsFBgbqpZde0pEjR2S1WnXNNde4nCSfzl76/4MPPlDfvn1LdezLL7+szp07q0OHDrrvvvsc5cjDw8Od1lwqqdatW+vWW2/V+PHj9c8//zjKkf/xxx+Sih99a926td566y09//zzatSokWrWrKlrrrlG1113nerWrav77rtPjz/+uAIDA/X++++rRo0aSk5Odhxfo0YNDRs2TGPHjtWNN96oG264QevWrdP333+fb/mCxx9/XN98841uvPFGR/n39PR0bdiwQV988YV27txZ5JIHnrrHqSgjRoxwPN/speG/+uqrfH/vgt6XgLOar8r5AfCegsrxLl261LRv396EhoaaWrVqmSeeeMLMnz8/XznjY8eOmbvuustUrVrVqSR4QSWHc3JyzODBg02NGjWMxWJxKut74MABc+utt5qKFSuac845xwwcONBs3Lgx3zmMMebLL780TZo0MVar1Vx44YVm1qxZpk+fPgWW7p48ebJp3bq1CQ0NNZUrVzYXXXSReeKJJ8zevXs9EDnXFVaOvKgy6HlLGNvjW1BJ8ILOd3rZ6DfeeMNcfvnlpnr16iYoKMjUqFHDdO3a1akcdVEKK0dujDG7du0y3bp1MxUrVjTVq1c3jzzyiKMM/OmlsNetW2duueUWc+655xqr1Wrq1atn7rjjDpOYmOjULjEx0bRq1coEBwebhg0bmilTppjHHnvMhISEOLWTZB5++OEC+5WammoefvhhEx0dbSpUqGAiIyNNp06dzOTJk53aZWdnm5deesk0bdrUWK1Wc84555jWrVubMWPGmCNHjjj6c9NNN5latWqZ4OBgU6tWLXPnnXfmKyldnHr16pkuXbrk296xY8d85aP//PNPc9ttt5mqVauakJAQ07ZtW/Pdd985tbGXI585c2aBj/fuu++aBg0amMDAQKe/R2n6YUzBz6uJEycaSWbevHlFXnNB7w3GGPPDDz+Yyy67zISGhpoqVaqYrl27ms2bNzu1sZcDP70MfEHvYenp6ebhhx821apVM2FhYaZ79+5m69atRpJ58cUXizw2JSXFdOnSxVSuXDlfKe81a9aYdu3ameDgYFO3bl0zbty4As+Rm5trxowZY6KiokxoaKi56qqrzMaNGwt87Rw9etSMGDHCNGrUyAQHB5vq1aubSy+91LzyyismOzu7yHiWBXvJ87zXZ7/mvH/Hkr4vAWcLizF+eJcoAKBMxcTEqEOHDpo4caJCQ0NLfQ+KJ3Tv3l2bNm0q8J43eEd2drbS0tI0ffp0DR48WKtWrXKMZN5xxx3auXOnVq5c6eNeFi4pKUmtWrXSJ598orvvvtvX3QFwhmMdJwCAJGn69OmqUaOGhg8f7vXHOr1q2LZt2zR37lxdddVVXn9snDJ37lzVqFEjX2l9Y4wWL17smFbqDwqqNDd+/HgFBAQ4Fb0BAG9hxAkAoKVLlzo+mEZHRxe4Bo8nRUVFqW/fvmrQoIF27dqlt956S1lZWVq3bl2B63r52oEDB4osUV1caWh/deDAAa1fv97xe7t27VS5cmUf9qhwY8aM0Zo1a3T11VcrKChI33//vb7//nsNGDBA77zzjq+7B+AsQOIEAChz/fr106JFi5SSkiKr1aoOHTrohRde0MUXX+zrrhUoJiYmX4nqvIorDQ33LVy4UGPGjNHmzZt17Ngx1a1bV7169dLTTz9dZJU6APAUEicAAIqRd0SuIOecc46jyiQA4Mzk08RpyZIlevnll7VmzRrt27dPX331VYFlfguydOlSdezYsVTlawEAAADAFT4d205PT1eLFi1077336pZbbinxcYcPH1bv3r3VqVMnpaamluoxbTab9u7dq8qVK/vtqusAAAAAvM8Yo6NHj6pWrVrFrqHnN1P1LBZLiUecevbsqfPOO0+BgYGaPXt2qUac/v77b0VHR7veUQAAAABnlN27d6tOnTpFtil3d1N+8MEH+uuvv/TJJ5+UqExqVlaWsrKyHL/b88Tdu3erSpUqXusnAAAAAP+Wlpam6OjoElUULVeJ07Zt2/Tkk0/q559/LnEFnbFjx2rMmDH5tlepUoXECQAAAECJbuEpNwvg5ubm6q677tKYMWPUuHHjEh83YsQIHTlyxPGze/duL/YSAAAAwJmo3Iw4HT16VKtXr9a6des0aNAgSScLPRhjFBQUpAULFuiaa67Jd5zVapXVai3r7gIAAAA4g5SbxKlKlSrasGGD07Y333xTP/74o7744gvVr1/fRz0DAAAAcKbzaeJ07Ngxbd++3fH7jh07lJSUpGrVqqlu3boaMWKE9uzZo48++kgBAQFq1qyZ0/E1a9ZUSEhIvu0AAAAA4Ek+TZxWr16tq6++2vF7fHy8JKlPnz6aOnWq9u3bp+TkZF91DwAAAAAk+dE6TmUlLS1N4eHhOnLkCFX1AAAAgLNYaXKDclNVDwAAAAB8hcQJAAAAAIpB4gQAAAAAxSBxAgAAAIBikDgBAACgXDLG6OX5WzRz9W5fdwVngXKzAC4AAACQV/Kh45q06E9VrVhBt7eJ9nV3cIZjxAkAAADl0tHMHEnSsf/+C3gTiRMAAADKpaycXElSjs0oJ9fm497gTEfiBAAAgHIp88SpZCkzh8QJ3kXiBAAAgHIp80Rugf8PeAOJEwAAAMolpxEnEid4GYkTAAAAyiXnESem6sG7SJwAAABQLmXmMFUPZYfECQAAAOVS3lGmrBwSJ3gXiRMAAADKJabqoSyROAEAAKBcyps4ZWQz4gTvInECAABAueQ04sRUPXgZiRMAAADKJedy5EzVg3eROAEAAKBcYgFclCUSJwAAAJRLmTksgIuyQ+IEAACAcilvspSVw1Q9eBeJEwAAAMolpuqhLJE4AQAAoFzKOsFUPZQdEicAAACUS3lLkFNVD95G4gQAAIByial6KEskTgAAACiXnNZxojgEvIzECQAAAOUSI04oSyROAAAAKJdInFCWSJwAAABQLuWdnpdFcQh4GYkTAAAAyp1cm1F2nsQpgxEneBmJEwAAAMqdrBznRImpevA2nyZOS5YsUdeuXVWrVi1ZLBbNnj27yPazZs3Stddeqxo1aqhKlSrq0KGD5s+fXzadBQAAgN84fd2mzBwSJ3iXTxOn9PR0tWjRQpMmTSpR+yVLlujaa6/V3LlztWbNGl199dXq2rWr1q1b5+WeAgAAwJ+cPsLEArjwtiBfPnjnzp3VuXPnErcfP3680+8vvPCCvv76a3377bdq1apVgcdkZWUpKyvL8XtaWppLfQUAAID/yJ84MeIE7yrX9zjZbDYdPXpU1apVK7TN2LFjFR4e7viJjo4uwx4CAADAG04fYaKqHrytXCdOr7zyio4dO6Y77rij0DYjRozQkSNHHD+7d+8uwx4CAADAG+z3NIVZT06gys61KddmfNklnOF8OlXPHdOmTdOYMWP09ddfq2bNmoW2s1qtslqtZdgzAAAAeJt9al54aAUdy8qRdLLSXsXgcvvxFn6uXI44TZ8+Xffff78+//xzxcbG+ro7AAAAKGP2qXnhoRUc2ygQAW8qd4nTZ599pn79+umzzz5Tly5dfN0dAAAA+IB9xKmSNVAVAi1O2wBv8OlY5rFjx7R9+3bH7zt27FBSUpKqVaumunXrasSIEdqzZ48++ugjSSen5/Xp00cTJkxQu3btlJKSIkkKDQ1VeHi4T64BAAAAZc9+j1NIhUCFBAXqRG4OiRO8yqcjTqtXr1arVq0cpcTj4+PVqlUrjRo1SpK0b98+JScnO9pPnjxZOTk5evjhhxUVFeX4eeSRR3zSfwAAAPiGfVqeNShQ1gqBTtsAb/DpiNNVV10lYwqvfjJ16lSn3xcvXuzdDgEAAKBcsI8uhVQIUEiFk2MB9lEowBvK3T1OAAAAgH10KaRCoEIcI04kTvAeEicAAACUOxkFjTiROMGLSJwAAABQ7mT9lySFVghUKPc4oQyQOAEAAKDcOXWPE1P1UDZInAAAAFDu5L3HyRrEiBO8j8QJAAAA5Y69gp41iHucUDZInAAAAFDuFDhVj3Lk8CISJwAAAJQ7zuXIA5y2Ad5A4gQAAIByx2kB3P/uccpiqh68iMQJAAAA5U5mzn8jTkFU1UPZIHECAABAuZPldI8TU/XgfSROAAAAKHecpupRHAJlgMQJAAAA5Y7TOk5M1UMZIHECAABAuWMfXTpZHIKpevA+EicAAACUO/bRJSvFIVBGSJwAAABQrhhjTlvHicQJ3kfiBAAAgHIlK+fUlLzQ4ECFOhInpurBe0icAAAAUK7kHVkKCQo4VY6cqnrwIhInAAAAlCv2kaWgAIuCAgOYqocyQeIEAACAciUzz+K3J/9LVT14H4kTAAAAypW8pcilk5X1JEac4F0kTgAAAChX7CNL9oTJPvKUlWOTMcZn/cKZjcQJAAAA5cqpqXoBTv+VnCvuAZ5E4gQAAIByJf89ToH59gGeRuIEAACAciXv4reSVCEwQIEBFqd9gKeROAEAAKBcyTqtOIR0cj0niREneA+JEwAAAMoVx1S9oFNT9BxrObEILryExAkAAADlyulT9fL+P1P14C0kTgAAAChX7CNO1jxT9awVmKoH7yJxAgAAQLmS8V9yFJpnxMn+/xkkTvASEicAAACUK0VN1csicYKX+DRxWrJkibp27apatWrJYrFo9uzZxR6zePFiXXzxxbJarWrUqJGmTp3q9X4CAADAf5y+AG7e/+ceJ3iLTxOn9PR0tWjRQpMmTSpR+x07dqhLly66+uqrlZSUpEcffVT333+/5s+f7+WeAgAAwF84ypHnraoXZC8OwYgTvCPIlw/euXNnde7cucTt3377bdWvX1+vvvqqJKlJkyb65Zdf9NprrykuLq7AY7KyspSVleX4PS0tzb1OAwAAwKeKrqpH4gTvKFf3OC1fvlyxsbFO2+Li4rR8+fJCjxk7dqzCw8MdP9HR0d7uJgAAALyooKl6jqp6OUzVg3eUq8QpJSVFERERTtsiIiKUlpamjIyMAo8ZMWKEjhw54vjZvXt3WXQVAAAAXnKqHDkjTig7Pp2qVxasVqusVquvuwEAAAAPKXCqXhAL4MK7ytWIU2RkpFJTU522paamqkqVKgoNDfVRrwAAAFCWMh3FIQqqqseIE7yjXCVOHTp0UGJiotO2hQsXqkOHDj7qEQAAAMpakes45ZA4wTt8mjgdO3ZMSUlJSkpKknSy3HhSUpKSk5Mlnbw/qXfv3o72DzzwgP766y898cQT2rJli9588019/vnnGjp0qC+6DwAAAB/IchSHyJs4sY4TvMunidPq1avVqlUrtWrVSpIUHx+vVq1aadSoUZKkffv2OZIoSapfv77mzJmjhQsXqkWLFnr11Vc1ZcqUQkuRAwAA4MxT8AK4FIeAd/m0OMRVV10lY0yh+6dOnVrgMevWrfNirwAAAODP7CXHCy4OQeIE7yhX9zgBAAAAjhGnoFOJk5WpevAyEicAAACUG8YYZdgTp+BTH2VD/xt9ymDECV5C4gQAAIByIzvXJvudHgVV1WOqHryFxAkAAADlRt6peHmn6p0qR85UPXgHiRMAAADKDXsp8gCLVCHQ4tjOArjwNhInAAAAlBt5F7+1WPImTkzVg3eROAEAAKDcyMzJv/itlLccOVP14B0kTgAAACg3TpUid/4Y65iql5Nb5DqhgKtInAAAAFBu5J2ql5f1v9+NOVl5D/A0EicAAACUG/YRJ+vpU/UqBORpQ+IEzyNxAgAAQLnhmKpXwfljbHBggOy1IrIoEAEvIHECAABAuZH53zpNeddwkiSLxUKBCHgViRMAAADKjcJGnPJus1feAzyJxAkAAADlRtaJgsuR593GWk7wBhInAAAAlBuFVdXLu42pevAGEicAAACUGxklGHHKYMQJXkDiBAAAgHKjRPc4kTjBC0icAAAAUG4UOVUviHuc4D0kTgAAACg37BXzTi9HLp0accriHid4AYkTAAAAyo2ip+r9N+JEOXJ4AYkTAAAAyo2sElXVI3GC55E4AQAAoNwoWXEIpurB80icAAAAUG447nEqYMTJSnEIeBGJEwAAAMoN+2iStcDiECyAC+8hcQIAAEC5UaKpehSHgBeQOAEAAKDcOJU4URwCZYvECQAAAOVG0Qvgso4TvIfECQAAAOVGVk4J1nFixAleQOIEAACAcsM+4hRawIhTaDAL4MJ7SJwAAABQbmQUcY+TvdJeRjaJEzyPxAkAAADlwolcm3JtRpIUUmA5chbAhff4PHGaNGmSYmJiFBISonbt2mnlypVFth8/frzOP/98hYaGKjo6WkOHDlVmZmYZ9RYAAAC+kvfeJWtR9zgxVQ9e4NPEacaMGYqPj1dCQoLWrl2rFi1aKC4uTvv37y+w/bRp0/Tkk08qISFBv//+u9577z3NmDFDTz31VBn3HAAAAGXNPpJksUjWoMITJ6rqwRt8mjiNGzdO/fv3V79+/XThhRfq7bffVsWKFfX+++8X2H7ZsmW67LLLdNdddykmJkbXXXed7rzzziJHqbKyspSWlub0AwAAgPLHPuJkDQqQxWLJt//UVD1GnOB5PkucsrOztWbNGsXGxp7qTECAYmNjtXz58gKPufTSS7VmzRpHovTXX39p7ty5uuGGGwp9nLFjxyo8PNzxEx0d7dkLAQAAQJk4VYo8//1N0qn7nkic4A1BvnrggwcPKjc3VxEREU7bIyIitGXLlgKPueuuu3Tw4EFdfvnlMsYoJydHDzzwQJFT9UaMGKH4+HjH72lpaSRPAAAA5ZBj8dsCCkNIee9xYqoePM/nxSFKY/HixXrhhRf05ptvau3atZo1a5bmzJmj5557rtBjrFarqlSp4vQDAACA8ifzROGL3+bdnmszOpFL8gTP8tmIU/Xq1RUYGKjU1FSn7ampqYqMjCzwmGeeeUa9evXS/fffL0m66KKLlJ6ergEDBujpp59WQEC5ygMBAABQCo4Rp8Km6uXZnnkiVxUC+WwIz/HZsyk4OFitW7dWYmKiY5vNZlNiYqI6dOhQ4DHHjx/PlxwFBp58gRhjvNdZAAAA+JyjOEQhiVPeSnus5QRP89mIkyTFx8erT58+atOmjdq2bavx48crPT1d/fr1kyT17t1btWvX1tixYyVJXbt21bhx49SqVSu1a9dO27dv1zPPPKOuXbs6EigAAACcmezrM4UUUIpckiwWi6xBAcrKsVEgAh7n08SpR48eOnDggEaNGqWUlBS1bNlS8+bNcxSMSE5OdhphGjlypCwWi0aOHKk9e/aoRo0a6tq1q/73v//56hIAAABQRoqbqmffl5Vjc1TgAzzFp4mTJA0aNEiDBg0qcN/ixYudfg8KClJCQoISEhLKoGcAAADwJ8UVh7DvO5LBVD14HnfMAQAAoFywJ06hRYw42fcxVQ+eRuIEAACAcuHUiFPRU/UkKYPECR5G4gQAAIByoST3OFkdI05M1YNnkTgBAACgXDhVjryIe5z+q7jHVD14GokTAAAAyoVT5ciLn6pH4gRPczlxOnz4sKZMmaIRI0bo0KFDkqS1a9dqz549HuscAAAAYFeycuT/jTjlMFUPnuVSOfLffvtNsbGxCg8P186dO9W/f39Vq1ZNs2bNUnJysj766CNP9xMAAABnuZKVIz+ZVGUx4gQPc2nEKT4+Xn379tW2bdsUEhLi2H7DDTdoyZIlHuscAAAAYFeiEacgpurBO1xKnFatWqWBAwfm2167dm2lpKS43SkAAADgdFk5JVsAV6KqHjzPpcTJarUqLS0t3/Y//vhDNWrUcLtTAAAAwOkcU/UoDgEfcClx6tatm5599lmdOHFCkmSxWJScnKzhw4fr1ltv9WgHAQAAAKmU6zjlkDjBs1xKnF599VUdO3ZMNWvWVEZGhjp27KhGjRqpcuXK+t///ufpPgIAAAAlW8eJqXrwEpeq6oWHh2vhwoVaunSp1q9fr2PHjuniiy9WbGysp/sHAAAASMqzjhPFIeADLiVOH330kXr06KHLLrtMl112mWN7dna2pk+frt69e3usgwAAAIB0ahQptIjEKTQ40Kkt4CkuTdXr16+fjhw5km/70aNH1a9fP7c7BQAAAJwuM7sEI06OqXqMOMGzXEqcjDGyWCz5tv/9998KDw93u1MAAADA6TJLUo6cqXrwklJN1WvVqpUsFossFos6deqkoKBTh+fm5mrHjh26/vrrPd5JAAAAnN1ybUYnco2kEpYjp6oePKxUiVP37t0lSUlJSYqLi1NYWJhjX3BwsGJiYihHDgAAAI/LO4JUdDlyqurBO0qVOCUkJEiSYmJi1KNHD4WEhHilUwAAAEBeeRMna1BR5ciZqgfvcKmqXp8+fTzdDwAAAKBQmTknR5CCgwIUEJD/Xnu7U/c4MeIEz3IpccrNzdVrr72mzz//XMnJycrOznbaf+jQIY90DgAAAJBOjSCFFDHaJJ0qHJHFiBM8zKWqemPGjNG4cePUo0cPHTlyRPHx8brlllsUEBCg0aNHe7iLAAAAONs5Eqci7m/Ku5/iEPA0lxKnTz/9VO+++64ee+wxBQUF6c4779SUKVM0atQo/frrr57uIwAAAM5y9ql3JU2cTuQa5dqM1/uFs4dLiVNKSoouuugiSVJYWJhjMdwbb7xRc+bM8VzvAAAAAJ2aelfUGk6n76dABDzJpcSpTp062rdvnySpYcOGWrBggSRp1apVslqtnusdAAAAoLyL3xYz4pRnjScSJ3iSS4nTzTffrMTEREnS4MGD9cwzz+i8885T7969de+993q0gwAAAIBjql4Ri99KUkCARcGB/63llENlPXiOS1X1XnzxRcf/9+jRQ/Xq1dOyZct03nnnqWvXrh7rHAAAACCdGj2yFjNVz94mO9fGiBM8qtQjTidOnNC9996rHTt2OLa1b99e8fHxJE0AAADwCvuIU2gxU/XytiFxgieVOnGqUKGCvvzyS2/0BQAAAChQRgnLkedtQ+IET3LpHqfu3btr9uzZHu4KAAAAULDMElbVy9vGPkoFeIJL9zidd955evbZZ7V06VK1bt1alSpVcto/ZMiQEp9r0qRJevnll5WSkqIWLVpo4sSJatu2baHtDx8+rKefflqzZs3SoUOHVK9ePY0fP1433HCDK5cCAACAciCLESf4mEuJ03vvvaeqVatqzZo1WrNmjdM+i8VS4sRpxowZio+P19tvv6127dpp/PjxiouL09atW1WzZs187bOzs3XttdeqZs2a+uKLL1S7dm3t2rVLVatWdeUyAAAAUE7YK+SVKHEKsidOjDjBc1xKnPIWhjDm5IrMFoul1OcZN26c+vfvr379+kmS3n77bc2ZM0fvv/++nnzyyXzt33//fR06dEjLli1ThQoVJEkxMTEuXAEAAADKE8dUvaCSVdXLewzgCS7d4ySdHHVq1qyZQkJCFBISombNmmnKlCklPj47O1tr1qxRbGzsqc4EBCg2NlbLly8v8JhvvvlGHTp00MMPP6yIiAg1a9ZML7zwgnJzC39RZGVlKS0tzekHAAAA5cupcuSlmKqXQ+IEz3FpxGnUqFEaN26cBg8erA4dOkiSli9frqFDhyo5OVnPPvtssec4ePCgcnNzFRER4bQ9IiJCW7ZsKfCYv/76Sz/++KPuvvtuzZ07V9u3b9dDDz2kEydOKCEhocBjxo4dqzFjxpTyCgEAAOBPHAvgluoeJ6bqwXNcSpzeeustvfvuu7rzzjsd27p166bmzZtr8ODBJUqcXGGz2VSzZk1NnjxZgYGBat26tfbs2aOXX3650MRpxIgRio+Pd/yelpam6Ohor/QPAAAA3lGqqnpBTNWD57mUOJ04cUJt2rTJt71169bKyckp0TmqV6+uwMBApaamOm1PTU1VZGRkgcdERUWpQoUKCgw89U1DkyZNlJKSouzsbAUHB+c7xmq1ymq1lqhPAAAA8E+O4hBBJR9xyiJxgge5dI9Tr1699NZbb+XbPnnyZN19990lOkdwcLBat26txMRExzabzabExETH9L/TXXbZZdq+fbtstlPDrn/88YeioqIKTJoAAABwZsgsVTny/0accpiqB89xacRJOlkcYsGCBWrfvr0kacWKFUpOTlbv3r2dpsaNGzeu0HPEx8erT58+atOmjdq2bavx48crPT3dUWWvd+/eql27tsaOHStJevDBB/XGG2/okUce0eDBg7Vt2za98MILpVo3CgAAAOVPVqkWwGUdJ3ieS4nTxo0bdfHFF0uS/vzzT0knp95Vr15dGzdudLQrrkR5jx49dODAAY0aNUopKSlq2bKl5s2b5ygYkZycrICAUy+O6OhozZ8/X0OHDlXz5s1Vu3ZtPfLIIxo+fLgrlwEAAIBywrXiECRO8ByXEqdFixZ5rAODBg3SoEGDCty3ePHifNs6dOigX3/91WOPDwAAAP9nLy1ekhEnq6M4BFP14Dkur+MEAAAAlJXS3OMUGsyIEzyPxAkAAAB+LyO7FMUh/qu8l0HiBA8icQIAAIDfc5QjL8U9TllM1YMHkTgBAADAr9lsRtmOdZxKUlXPXo6cESd4DokTAAAA/FpWnvWYqKoHXyFxAgAAgF/LmwCVagFcpurBg0icAAAA4NfsU+4qBFoUGFD0OqGSZA1ixAmeR+IEAAAAv+ZY/Dao+NEmial68A4SJwAAAPg1ewJkLcE0PSlvcQim6sFzSJwAAADg104tfluyj672EafsHJtsNuO1fuHsQuIEAAAAv+aYqlfiEadT7bIYdYKHkDgBAADAr9mLQ5R4xCnPWk/c5wRPIXECAACAX8uyT9UrYXGIoMAABf1XfY9FcOEpJE4AAADwa/apeqHBJUucJCnUUVmPqXrwDBInAAAA+DVHVb0SjjhJpyrwMVUPnkLiBAAAAL+WUcqqennbZpA4wUNInAAAAODXSltVL29bRpzgKSROAAAA8GulXccpb9ss7nGCh5A4AQAAwK85ypGX4h4ne1tGnOApJE4AAADwa1nuTNWjHDk8hMQJAAAAfs2dqXqUI4enkDgBAADAr51KnChHDt8hcQIAAIBfs48aWUszVS+IBXDhWSROAAAA8GunikO4MlWPESd4BokTAAAA/JorU/UoDgFPI3ECAACAX3NtAVzWcYJnkTgBAADAr7lUVY91nOBhJE4AAADwa1k5J0eNQksx4hQaTOIEzyJxAgAAgF9zrxw5U/XgGSROAAAA8GsZLk3VC3A6FnAXiRMAAAD8mn3EyRrkQlU9Eid4iF8kTpMmTVJMTIxCQkLUrl07rVy5skTHTZ8+XRaLRd27d/duBwEAAOATxhgXq+rZy5EzVQ+e4fPEacaMGYqPj1dCQoLWrl2rFi1aKC4uTvv37y/yuJ07d2rYsGG64ooryqinAAAAKGtZeRKfUk3Vc5QjZ8QJnuHzxGncuHHq37+/+vXrpwsvvFBvv/22KlasqPfff7/QY3Jzc3X33XdrzJgxatCgQRn2FgAAAGUp7zpMLo04kTjBQ3yaOGVnZ2vNmjWKjY11bAsICFBsbKyWL19e6HHPPvusatasqfvuu6/Yx8jKylJaWprTDwAAAMqHzJyTiU9ggEUVAl1Zx4mpevAMnyZOBw8eVG5uriIiIpy2R0REKCUlpcBjfvnlF7333nt69913S/QYY8eOVXh4uOMnOjra7X4DAACgbDhKkQeV7mOrfaqePfEC3OXzqXqlcfToUfXq1UvvvvuuqlevXqJjRowYoSNHjjh+du/e7eVeAgAAwFNcKQyRtz1T9eApQb588OrVqyswMFCpqalO21NTUxUZGZmv/Z9//qmdO3eqa9eujm0228kXU1BQkLZu3aqGDRs6HWO1WmW1Wr3QewAAAHibK4vfSpLVPuJ0wiZjjCwWi8f7hrOLT0ecgoOD1bp1ayUmJjq22Ww2JSYmqkOHDvnaX3DBBdqwYYOSkpIcP926ddPVV1+tpKQkpuEBAACcYRxrOJWiop7knGhlUZIcHuDTESdJio+PV58+fdSmTRu1bdtW48ePV3p6uvr16ydJ6t27t2rXrq2xY8cqJCREzZo1czq+atWqkpRvOwAAAMo/+zpMIaVY/Pb09lknbKUesQJO5/PEqUePHjpw4IBGjRqllJQUtWzZUvPmzXMUjEhOTlZAQLm6FQsAAAAecmqqXuk+D1YItCjAItnMyQIR4argje7hLOLzxEmSBg0apEGDBhW4b/HixUUeO3XqVM93CAAAAH7BnjiFBpduxMhisSi0QqDSs3MpEAGPYCgHAAAAfsu+AG5pp+pJeSvrcY8T3EfiBAAAAL+V4WJVvbzHZDDiBA8gcQIAAIDfcrWqXt5jmKoHTyBxAgAAgN9ydQFc6dT0PhIneAKJEwAAAPxWZs5/U/Vcusfp1CK4gLtInAAAAOC3XC1HfvKYk8lWVg4jTnAfiRMAAAD8lltT9SowVQ+eQ+IEAAAAv5Xl1ogTU/XgOSROAAAA8FuOe5woDgEfI3ECAACA38p0YwFcKwvgwoNInAAAAOC33FnHyTFVj+IQ8AASJwAAAPitU1X1KA4B3yJxAgAAgN+yT7MLdSFxCmWqHjyIxAkAAAB+y63iEP9N1ctixAkeQOIEAAAAv5XlWMfJ9QVwuccJnkDiBAAAAL+V4c49Tv9V4svIJnGC+0icAAAA4LccxSFcKkfOArjwHBInAAAA+CVjTJ6qekzVg2+ROAEAAMAvncg1spmT/291qxw5I05wH4kTAAAA/FLekSKXRpyCqKoHzyFxAgAAgF+yT9OzWKTgQDem6pE4wQNInAAAAOCXHKXIgwJlsVhKffype5yYqgf3kTgBAADAL7lTGCLvcYw4wRNInAAAAOCXMh2L35a+METe4zJP5MoY47F+4exE4gQAAAC/ZC8O4XLi9N/aTzZzskIf4A4SJwAAAPgl+xQ7a5BrH1mteab4sZYT3EXiBAAAAL9kn6oXGuzaiJM1KED2mhLc5wR3kTgBAADALzmKQwS5ljhZLBbHsVksggs3kTgBAADAL7lbVS/vsYw4wV0kTgAAAPBLpxIn10ac8h6bQeIEN5E4AQAAwC+5W44877GZTNWDm/wicZo0aZJiYmIUEhKidu3aaeXKlYW2fffdd3XFFVfonHPO0TnnnKPY2Ngi2wMAAKB88sRUPXtFPqbqwV0+T5xmzJih+Ph4JSQkaO3atWrRooXi4uK0f//+AtsvXrxYd955pxYtWqTly5crOjpa1113nfbs2VPGPQcAAIA32UuIW10sDiE5L4ILuMPnidO4cePUv39/9evXTxdeeKHefvttVaxYUe+//36B7T/99FM99NBDatmypS644AJNmTJFNptNiYmJZdxzAAAAeJNnpur9N+KUw1Q9uMeniVN2drbWrFmj2NhYx7aAgADFxsZq+fLlJTrH8ePHdeLECVWrVq3A/VlZWUpLS3P6AQAAgP/zTFU9RpzgGT5NnA4ePKjc3FxFREQ4bY+IiFBKSkqJzjF8+HDVqlXLKfnKa+zYsQoPD3f8REdHu91vAAAAeJ9HRpwc6ziROME9Pp+q544XX3xR06dP11dffaWQkJAC24wYMUJHjhxx/OzevbuMewkAAABX2O9xCgnyxDpOTNWDe4J8+eDVq1dXYGCgUlNTnbanpqYqMjKyyGNfeeUVvfjii/rhhx/UvHnzQttZrVZZrVaP9BcAAABlJ8uD6zgxVQ/u8umIU3BwsFq3bu1U2MFe6KFDhw6FHvd///d/eu655zRv3jy1adOmLLoKAACAMubRdZxySJzgHp+OOElSfHy8+vTpozZt2qht27YaP3680tPT1a9fP0lS7969Vbt2bY0dO1aS9NJLL2nUqFGaNm2aYmJiHPdChYWFKSwszGfXAQAAAM/yyDpOTNWDh/g8cerRo4cOHDigUaNGKSUlRS1bttS8efMcBSOSk5MVEHDqxfLWW28pOztbt912m9N5EhISNHr06LLsOgAAALzIcY+TGyNOoUzVg4f4PHGSpEGDBmnQoEEF7lu8eLHT7zt37vR+hwAAAOBzHp2qx4gT3FSuq+oBAADgzJXpieIQQfYFcBlxgntInAAAAOCXPLoAbjaJE9xD4gQAAAC/5JiqF0RVPfgeiRMAAAD8kkem6lFVDx5C4gQAAAC/k5NrU47NSHK3HDlV9eAZJE4AAADwO5k5p0aI3CsOQeIEzyBxAgAAgN/Jm+hYg9wpDsFUPXgGiRMAAAD8jj1xsgYFyGKxuHwe+2hVFsUh4CYSJwAAAPgdTyx+m/d4RpzgLhInAAAA+B1PrOGU93jucYK7SJwAAADgd+xT69wecfqvOESOzSgnl1EnuI7ECQAAAH7HPrUu1M3EKTT41PF5K/UBpUXiBAAAAL/jKA7hZuKUtyIf0/XgDhInAAAA+B1HcQg3SpFLksVicSRPJE5wB4kTAAAA/M6p4hDujTjlPQeV9eAOEicAAAD4nQwPVdXLew5GnOAOEicAAAD4He+MOJE4wXUkTgAAAPA7WTn2e5w8kDgFMVUP7iNxAgAAgN/x1AK4ec/BiBPcQeIEAAAAv+PJqXr2kuaZOSROcB2JEwAAAPyOfVqdu+s4SVTVg2eQOAEAAMDveHSqHus4wQNInAAAAOB3Mj1ZHIKqevAAEicAAAD4Hc+WIz/5kddeqQ9wBYkTAAAA/I5nq+ox4gT3kTgBAADA72T9V8gh1AMjTqEkTvAAEicAAAD4HXvpcI+WI6eqHtxA4gQAAAC/Yx8dsrIALvwEiRMAAAD8jn10yCPFIYLsC+Ay4gTXkTgBAADA72TYi0N4sBx5RjYjTnAdiRMAAAD8jmer6tnLkZM4wXV+kThNmjRJMTExCgkJUbt27bRy5coi28+cOVMXXHCBQkJCdNFFF2nu3Lll1FMAAACUhSxPTtWjqh48wOeJ04wZMxQfH6+EhAStXbtWLVq0UFxcnPbv319g+2XLlunOO+/Ufffdp3Xr1ql79+7q3r27Nm7cWMY9BwAAgDfk2oyycz2ZONmLQ3CPE1xnMcYYX3agXbt2uuSSS/TGG29Ikmw2m6KjozV48GA9+eST+dr36NFD6enp+u677xzb2rdvr5YtW+rtt9/O1z4rK0tZWVmO39PS0hQdHa0jR46oSpUqXriikvt89W59v2GfT/sAAADgb3KNtOSPA5Kkzc/GqWJwkFvnW7b9oO6askKVggPVtn41T3QRHjC+RyuFV6zg0z6kpaUpPDy8RLmBe89CN2VnZ2vNmjUaMWKEY1tAQIBiY2O1fPnyAo9Zvny54uPjnbbFxcVp9uzZBbYfO3asxowZ47E+e9KfB45p0dYDvu4GAACAXzqnYgWPFIeIqhoqSUrPzuWzlx/Jys2V5NvEqTR8mjgdPHhQubm5ioiIcNoeERGhLVu2FHhMSkpKge1TUlIKbD9ixAinRMs+4uQPulwUpYY1wnzdDQAAAL90cd2qCgiwuH2e+tUracaA9tp16LgHegVPqWwtP0mT5OPEqSxYrVZZrVZfd6NAzetUVfM6VX3dDQAAgDNeuwbnql2Dc33dDZRjPi0OUb16dQUGBio1NdVpe2pqqiIjIws8JjIyslTtAQAAAMBdPk2cgoOD1bp1ayUmJjq22Ww2JSYmqkOHDgUe06FDB6f2krRw4cJC2wMAAACAu3w+VS8+Pl59+vRRmzZt1LZtW40fP17p6enq16+fJKl3796qXbu2xo4dK0l65JFH1LFjR7366qvq0qWLpk+frtWrV2vy5Mm+vAwAAAAAZzCfJ049evTQgQMHNGrUKKWkpKhly5aaN2+eowBEcnKyAgJODYxdeumlmjZtmkaOHKmnnnpK5513nmbPnq1mzZr56hIAAAAAnOF8vo5TWStNrXYAAAAAZ67S5AY+vccJAAAAAMoDEicAAAAAKAaJEwAAAAAUg8QJAAAAAIpB4gQAAAAAxfB5OfKyZi8imJaW5uOeAAAAAPAle05QkkLjZ13idPToUUlSdHS0j3sCAAAAwB8cPXpU4eHhRbY569Zxstls2rt3rypXriyLxeLr7vhEWlqaoqOjtXv3btaycgHxcx2xcx2xcw/xcx2xcw/xcx2xcx2xKzljjI4ePapatWopIKDou5jOuhGngIAA1alTx9fd8AtVqlThxeQG4uc6Yuc6Yuce4uc6Yuce4uc6Yuc6YlcyxY002VEcAgAAAACKQeIEAAAAAMUgcToLWa1WJSQkyGq1+ror5RLxcx2xcx2xcw/xcx2xcw/xcx2xcx2x846zrjgEAAAAAJQWI04AAAAAUAwSJwAAAAAoBokTAAAAABSDxAkAAAAAikHiBAAAAADFIHECAABwE0WKXUfsSufEiRP6+++/JUk2m83HvTm7kDidYXjzcV1WVpaSkpIkSbm5ub7tTDmUmZmp9957T+vWrfN1V8qdvP8ISryOS4t4eQZxLJ2JEyfqrrvu0ujRo/XXX3/JYrH4ukvlxvvvv69hw4Zp4sSJSk1NJXalsGzZMkVEROi5555TRkaGAgICeO2WIRKnM8ikSZN05513atCgQfr555+VnZ3t6y6VGzt27FBYWJi6d++uw4cPKzAwkG9xSuGNN95QzZo1NWPGDB04cIDnXimMGzdOLVq00M0336zu3btr48aNslgsPP9KaPz48erWrZseeughrVixgudeKbzzzjsaMGCARo8era1bt/LhtYQOHTqkzp07a9y4capWrZo+++wzxcXF6f333/d11/zeihUr1KxZM7366qvau3evxo4dqxtuuEF79+71ddfKjcTERGVmZiolJUXz5s2TJF67ZYjE6Qywbt06XXLJJZo4caLOO+88rVixQgMHDtTMmTN93bVyY8uWLapfv77q1aunF198URJvRCU1ffp0vfPOO5o8ebIWLFig6667TsHBwb7ult87evSoevbsqbfeekvPPvus7r33Xh05ckQjRoyQJAUE8PZclPT0dN16662aMGGCWrVqpdWrV6tXr1567rnnfN01v7d69Wq1atVKEyZMUFBQkD755BN1795dy5Yt83XXyoU1a9bor7/+0o8//qg33nhDW7duVbt27TR58mTNnz/f193zW5s2bdLIkSPVuXNnrVq1StOmTdO2bdu0ZcsWJSYmSmLUsyj2L9OCg4M1dOhQHT58WN9884327dvntB/exb/M5VxKSoomTpyoVq1aacWKFXruuee0atUqVa9eXatWrZLEG1FR7LE5fPiwLrjgAl1zzTX65ptv9Ntvv/GtfzHs0xlnzZqlzp07q2fPntqzZ4+mTZum1atXKzU1VRJv5oXZsGGD1q9fr1mzZum2227Tgw8+qDZt2qh69eqSTj43ee0WbvPmzdqwYYNmzZqlZ599VitXrtSdd96pzz77TNOnT/d19/zWjh079OKLL+qKK67QypUr9eabb2r79u1KT0/XkiVLJPFvRmHs72W7d++WzWZTlSpVHPueeOIJRUVF6dVXX/VV9/xelSpVFBgYqAEDBqhixYrKzs5WpUqVdPnll2vlypWS+MKyKPYv077++mtdd911euCBB7Ru3Tp9++23TvvhXUS5nLL/wxYcHKywsDA98MADCg8PV1ZWliSpTZs2jsSJN6LC2WOzcuVKxcbGqlevXqpRo4ZeeuklSbwRFSUwMFDZ2dlasWKFbrzxRn366adq3ry53nzzTXXt2lVdu3bV8ePHiWEhjh07puTkZFWoUMGxbePGjapTp442bdoki8XCa7cA9g+vhw8f1j///KMaNWo49t1///266qqr9Mwzz/iqe34vNDRUwcHBGjhwoMLCwpSZmSlJuuKKK7RixQpJ/JuR14IFC/TTTz/pn3/+cbyXpaenq0KFCjp48KCjXfPmzdWtWzcdPHhQn376qa+661fyxs5msyk6OlozZ87UeeedJ0mOmQm7d+9W+/btfdlVv2OPXd7nWHZ2tmw2m6KiolShQgXdfvvtaty4sRYuXKjHHntM//vf/3zY47MHn2jKmQ0bNkiSYzSkWrVqevnll3XxxRdLkqxWq6STb0RXXHGFz/rpj+yxy8v+IaxixYo6evSo6tevr759+yopKUn33nuv7rvvPqc3rrPZ6fEzxignJ0fNmzfX5MmTNX36dE2dOlVz5szRzJkzdfz4cfXu3VsSo055Y2f/0qNmzZpq27at4uLiNGzYMFWrVk07duzQokWL1KVLF/Xv399X3fU7ixYt0tKlS7V//37Hh/q0tDTVrl1bv//+u6NddHS07r77bgUEBPDN/3/yxi43N1eRkZGaPHmymjZtKkkKCQmRJCUnJ+vKK6/0ZVf9yqZNm9S8eXPde++96tu3r6699lq98847kqQ+ffror7/+0g8//OB0TKdOnVStWjWtXbtWOTk5vui2Xygodu+9954kqXLlyk7/HiQnJ+vEiRNq2bLlWf/vhJQ/dtddd52mTJki6WSiabFYtHnzZoWHhysoKEht27bV999/rzfeeEPh4eE+7v1ZwqBcSEpKMu3atTO1atUyP/zwgzHGmJycHMd+m83m9P+XXnqp+fzzz8u8n/6ooNjl5uY6tbn00kvNN998Y4wx5pNPPjFVq1Y1AQEB5vXXXzfGOMf3bFPUcy8nJ8c88MADJioqylx66aUmKyvLcVxiYqKxWCxmx44dvui2XygodidOnHDs3717t/nqq69M27ZtTUJCgjlx4oQ5fPiwWbRokbFYLGbNmjXGmLP3+ffrr7+a5s2bm4YNG5oGDRqYli1bmk8++cQYY8zx48dN3bp1zZNPPmnS09Mdxxw8eND07dvX3HLLLeb48eO+6rrPFRS7adOmOfbnfQ88fPiwueiii8xPP/3ki676pYEDB5pbb73VHDx40GzcuNE89NBDJiIiwsyePdsYY8yjjz5qoqOj872/9ezZ03Tp0sUYc/a+bguL3XfffWeMOfncsz//Zs6caerVq2eOHDniOD4jI8Mn/fYHxcVu+/btplu3bmbbtm3mxhtvNFar1TRr1sxcfvnlZunSpcaYs/d5V1YYcSoHfv75Z/Xv31/h4eFq3LixPvnkE+Xm5iowMNDx7XXeqRXbt2/Xhg0b1KxZM8e2f//9t8z77Q8Ki529fKf9G67GjRvr2LFjuummm3TfffepXbt2uuCCCxQUFCTp7B0xKeq5l5OTo8DAQN16663KycnR4cOHnYpC1KlTR3Xr1tX69et9eAW+U1jsgoKCHM+nOnXqqF69etq1a5fuv/9+BQUFKTw8XM2aNVPt2rX1448/Sjo7p0599dVX6tu3rzp37qzly5dr1qxZOv/88zVt2jTt3LlToaGhGjx4sN5++22tXr3acdy5556ratWqaf/+/QoNDfXhFfhOYbH75JNPtGfPHklyKmG8ceNG7d27V40aNXKc4+jRoz7puz/4559/NH36dHXv3l3nnnuumjZtqpEjR+rGG2/UQw89JEl64YUXlJ2drTFjxiglJcVxrDFGNWvWlDHmrHzdFhW7gQMHSjr53LPfIzt//ny1a9dOVapU0aFDh3TvvfdqwoQJjtsOziZFxW7AgAGSpMjISH333Xdq3LixTpw4oU2bNumLL77QiRMnNGXKFB07duysfN6VJRKnciAmJkatWrXSm2++qS5duuj333/XtGnT8rWz/yM4b948NWzYUE2aNHG8EfXq1eusnHJWVOyMMY456z/++KPuvvtu2Ww2rVu3TtOmTVOnTp30wgsvKCUlRYGBgb68DJ8pKn72N+fY2Fj169dP+/fv16RJkxzHbt++XVWrVj1r566XJHaSdPDgQTVs2NDpw1dSUpIqV66suLi4Mu+3P8jKytLmzZt1880366mnnlL16tXVokULde/eXVu2bHHclD9s2DBFR0dr/Pjx+u233xzH5+Tk6Jxzzjkrp0sVF7tKlSo52tr/zfj+++/VvHlz1apVy/FvxvDhw5WWluary/CpgIAANWzYUPv373dsi4qK0tChQ3XixAmNHj1aoaGhmjJlin766SfdfffdmjZtml544QUlJibq1ltvPWs/vBYXO3vVywoVKsgYo7/++ku33XabJk6cqPr162vt2rXq2bOn47aDs0lxsXv++edVqVIlTZw4UV9//bXmzJmjhg0b6vzzz9c111yjGjVqcE9xWfDZWBdK5ejRo8YYY/bu3Wtuu+02c+ONN5rU1FRjTP5pZw8//LAZOXKkGT9+vKlcubK55JJLzPbt28u8z/6iqNjZp5b98ssvZvbs2U7TqObOnWuGDx9uDh06dFYPfZckfsnJyWbw4MHGYrGYO++80zzyyCOmZs2a5sknnzTZ2dlnbfyKip39ubZ+/XpzxRVXmHbt2pl3333XPPvssyYyMtIMGjTIHD9+/KyN3fLly82hQ4eMMaemnqxYscLExMSYffv2ObYtXbrUdOjQwVxwwQVmwoQJZvTo0aZatWrm008/9Vnffa242J3ujjvuMBMnTjQTJkwwlStXNhdffLH566+/yrTP/uTo0aPmhhtuMP3793e8Xo05+X6XkJBg6tSpY44dO2aMMWbhwoWmR48e5rLLLjNNmzZ1TKk6WxUVu9GjR5vo6GjHFNo1a9YYi8ViAgICTGRkpGMa5NmquOdd3tjZ2T//5f3sAu8icSpH7PeVfPLJJ6Z9+/Zm7Nix+dr8+++/pnbt2sZisZg6deqYr776qox76Z9KEjs7+weNs/UDa0FKGr/JkyebRx55xHTp0sV8++23ZdlFv1WS2C1cuNDcfPPNjgSK2J1is9kcHw7+97//mXbt2jm2223ZssU8+OCDpnPnzqZVq1bE7z8lid3OnTtNaGiosVgsJioq6qz4N+P0e4LzssdrwoQJ5vzzzzezZs1y2j9nzhzTokUL8+uvvzptz/tB90zmiditWLHCGHPyCzer1WomTZrk5V77B3dj17x5c7Ny5UrvdxRFInHyI8V9ULfvP378uOnfv7+56qqrzG+//WaMOfXh7MCBA6ZTp07mnXfe8W5n/YwnYnc2J0ruxO9s/6bLndhlZ2c7tUlJSfFiT/1TSV93NpvNdO7c2Tz//POFtrGP8J0tPBG7vXv3mvPOO8+89dZbnu6eX5owYYIZMGCASUhIKHBULW9M27dvb+644w6zfv16x7bExEQTFBRktmzZkq/9mc7TsTPm1Hvgmc4bsYNvMBnSxz7++GPNmTNHUvE3gNtLkIeGhqpnz54yxmjq1KnatWuX7r77bv3222+qXr265s2b57iR8Ezmydht3rz5rJuT7qn43XPPPdq8eXNZdNlveCp2vXr1cqzZJEkRERFe77s/KE387LKysrRr1y7HPXPHjh3TK6+8ot27dzvahIWFeb6zfsaTsduxY4eioqK0efNmPfDAA17rsz/45ZdfdOGFF+q9995TYGCgPvjgA91yyy1OhUWkkzG1Fy4YM2aMdu7cqaeeekrbt2/XkSNHtHDhQl177bWKjIx0tD/TeSt2kpzWsTsTeTN28BHf5m1nrzVr1pjLL7/cWCwWc8cdd5idO3caY0r37dWwYcNM1apVjdVqNTExMWb79u1nxbdfxM49xM91xM497sRv5cqVpmHDhiY1NdVMnjzZVK5c2bRv397s3bvX2932C96I3Z49e86K597WrVtN165dzbBhwxylrm02m6lZs6aZOHGiMSb/vcJ23333nWnTpo2pW7euadCggYmKijILFiwos777GrFzHbE7MzHi5APHjh3TF198ofr162v8+PFasWKFfvjhB9lsthJ9e3XixAnNnj1bn376qc455xxNnz5dO3bsUMOGDc/4b7+InXuIn+uInXvcjd93332nf//9V1deeaWefPJJvfXWW1q+fLmioqLKoPe+5a3Y1apV66x47gUHByssLEwDBgxQSEiIMjMzZbFY1K5dO61atUqS8lUjM/9VHOzSpYt++OEHffHFF/q///s/7d27V9dee22ZX4OvEDvXEbszlI8Tt7NSdna2WbJkiVm+fLkxxpjbbrvNXHPNNU7zWYuSnp5uLrjgAvPEE094s5t+idi5h/i5jti5x5345ebmmnvvvdeEhoaaMWPGeLurfofYlc7ChQvN0qVLzYEDBxzf6OddJNnu4osvNlOmTCnr7vk1Yuc6Ynd2IHEqA3lfTPZpEfaCBMYYs3nzZlOnTh3z/PPPO15khU2fsB+XmZnp5V77B2LnHuLnOmLnHk/Gz5iTSwYcPnzYu532E8TONYsXLzYXXHCBOf/88010dLRp166d+eKLLxz7806LSk1NNU2aNDGrV6/2RVf9DrFzHbE7u5A4eVFxLyZjTv1j+Pjjj5sWLVqYRYsW+aCn/ofYuYf4uY7YucfT8Tsb7sGxI3ausdlsZtq0aeb88883CQkJ5uDBg2bVqlXmlltuMTfddJPZv3+/U1tjjJk/f76JiIgwBw4ccOxLS0sr8777GrFzHbE7O3GPkxcYY/TZZ59p4MCB6tGjh5YuXapZs2apdu3a+vjjj3XgwAFHW/v88lGjRikrK0tffvmlDh48KEn6888/Hec7WxA79xA/1xE793grfmfDPTjEzj0ZGRn6/fffdeedd+qJJ57QOeecozZt2uj666/X1q1bValSJUdbe2y+++47tW/fXtWrV1dKSoruueceJSQk6Pjx4766DJ8gdq4jdmcnEicvKM2LKSAgQLm5uQoLC9PQoUO1YMECTZ06VV27dtUdd9yhf/7556z5x08idu4ifq4jdu4hfq4jdu6pWLGirr/+esXHx6tixYqO62/YsKFsNpsyMzMdbQMCAmSz2bRz507ddNNNeumll9SwYUPt2LFDw4YNU8WKFX11GT5B7FxH7M5SPhjlOissXbrUMfxqH6JNTEw0jRs3Nv/880+Bxxw6dMhUrlzZWCwWc91115nk5OQy668/IXbuIX6uI3buIX6uI3aeYbPZHPeUDB8+3MTGxjq2223dutVYLBZjsVhMvXr1zLx583zSV39D7FxH7M4ejDh5yaWXXqrKlSvLnLyPTJK0YMEC1a1bV9WqVcs3jWfBggWKiIhQ3bp19fPPP2v+/PmKjo72Rdd9jti5h/i5jti5h/i5jtgV7dChQ3rvvfeUkpKSb5994VB7eXb7qNzy5cvVtWtXSc7TFoOCgtSgQQN98MEH2rlzp+Li4srmInyE2LmO2CGfss7UzhT//POPmTJlitm3b1++ffabd/NWUsnJyTFXXnmlmTBhQoHn27x5s3n99de901k/Q+zcQ/xcR+zcQ/xcR+xcN2LECMc39f/++69je95Kg6c7cuSIady4sfn111+NMcb8+++/5sUXX3TEv7CFR880xM51xA4FIXFygSdfTHmrrpwNiJ17iJ/riJ17iJ/riJ1rPvnkE1OtWjXTrFkzM27cONO4ceN8VQaNMeazzz4zl1xyiVOFwcTERNOkSRNz6NAhM2HCBBMaGmquuuoqpxLvZzJi5zpih6KQOJWCN15MBw8eLMMr8B1i5x7i5zpi5x7i5zpi55rc3FzTvXt3Y7FYzFtvvWWMMebAgQMmMjLSfPbZZ8aYk0lnamqqad++vYmMjDQTJkww2dnZjnMMHTrU1KxZ0zRs2NDUrFmzwLifiYid64gdSoLEqQR4MbmO2LmH+LmO2LmH+LmO2Llv1apV5ujRo8aYU9ObWrVqZQYNGuRok56ebp577rl80x9PnDhhevToYSpXrmxeeumlsuu0nyB2riN2KA6JUwnxYnIdsXMP8XMdsXMP8XMdsSudZcuWmb179+bbbp/elJGRYXr37m26detmMjIyij3fTz/9ZNLT0z3eT39E7FxH7FBaJE6F4MXkOmLnHuLnOmLnHuLnOmLnmh9++MHUr1/f1KtXz9SpU8f079/fbN261Rhj8t0T8uCDD5p27doZYwq/yf5suo+E2LmO2MFVlCM/TWJioho0aKA777xTbdu21YABA/THH39IOrnys720ZEhIiCpVqqTU1FSFhITIZrMVeD7zXwnZK6+88oxf4IzYuYf4uY7YuYf4uY7YuW737t0aOXKk7rnnHiUmJuq1115TYmKiRowYoV27dslischmszliFRsbq82bN2vfvn0KCCj448vZsvgvsXMdsYNbfJOv+afk5GTTvn1788wzz5jt27ebmTNnmgYNGphbbrnF7Ny50xhz8tsG+zcOX375palcuXKB3zKebYide4if64ide4if64idexYsWGBCQ0PN9u3bHdu+/PJLc+WVV5qBAwfmaz979mzTuHFj89NPP5VlN/0SsXMdsYM7GHHKY8uWLVq/fr369Omjhg0b6rbbbtPLL7+sgwcPauzYsZKkgIAAxzcOgYGBioqK0rZt23zZbb9A7NxD/FxH7NxD/FxH7Nxz6NAhNWnSxLGQqCTddNNNiouL088//6yffvpJkpSTkyNJat++vZKTk5WZmSlJ+RYFPpsQO9cRO7iDxCkPXkyuI3buIX6uI3buIX6uI3buadq0qTZv3qwtW7Y4tgUGBuqGG25QdHS0vv76a0lSUFCQbDabgoODFRkZqfXr10s6u6dHETvXETu4g8QpD15MriN27iF+riN27iF+riN27mnWrJmuvvpqjRs3TseOHXNsb9mypWrWrKm//vrLca9JQECAKlasqJkzZ+rxxx/3Ya/9A7FzHbGDW3w2SdBPde7c2XTs2NFRRtauV69e5qabbnKar56ZmWlWrVrli276JWLnHuLnOmLnHuLnOmLnnqSkJBMUFGTeeustk5WV5dj+9NNPm0aNGvmwZ/6P2LmO2MFVjDidZuzYsVq6dKk++eQTZWdnO7bXrVtXmzZtcpqvbrVa1aZNG1911e8QO/cQP9cRO/cQP9cRO/e0aNFCw4cP13PPPaePP/5Y6enpOnr0qFavXq177rnH193za8TOdcQOrgrydQf8Td4XU4UKFdSzZ0/ZbDZeTCVA7NxD/FxH7NxD/FxH7Nz3/PPP699//9Uzzzyjd999VykpKapUqZLGjRvn6675PWLnOmIHV1iMOcvvTi3Eww8/rK+++kp169Z1vJhmzpypCy+80Ndd83vEzj3Ez3XEzj3Ez3XEzj2ZmZn6/ffftXbtWlmtVpLOUiB2riN2KC0Sp0LwYnIdsXMP8XMdsXMP8XMdsQOAMx+JEwAAAAAUg+IQAAAAAFAMEicAAAAAKAaJEwAAAAAUg8QJAAAAAIpB4gQAAAAAxSBxAgAAAIBikDgBAAAAQDFInAAAAACgGCROAOAFixcvlsVi0eHDh8v8sS0WiywWi6pWrVrmj+1pFotFs2fP9vh5R48erZYtW3r8vACAMxeJEwC46aqrrtKjjz7qtO3SSy/Vvn37FB4e7pM+ffDBB/rjjz988tjlwbBhw5SYmOjTPmRmZqpv37666KKLFBQUpO7duxfYbtKkSWrSpIlCQ0N1/vnn66OPPnLaf9VVVzmS5bw/Xbp0cbQxxmjUqFGKiopSaGioYmNjtW3btmL7mJycrC5duqhixYqqWbOmHn/8ceXk5Dj279u3T3fddZcaN26sgICAfK+DokyaNEkxMTEKCQlRu3bttHLlSqf9kydP1lVXXaUqVar47EsIAMiLxAkAvCA4OFiRkZGyWCw+efyqVauqZs2aPnns8iAsLEznnnuuT/uQm5ur0NBQDRkyRLGxsQW2eeuttzRixAiNHj1amzZt0pgxY/Twww/r22+/dbSZNWuW9u3b5/jZuHGjAgMDdfvttzva/N///Z9ef/11vf3221qxYoUqVaqkuLg4ZWZmFtm/Ll26KDs7W8uWLdOHH36oqVOnatSoUY42WVlZqlGjhkaOHKkWLVqU+NpnzJih+Ph4JSQkaO3atWrRooXi4uK0f/9+R5vjx4/r+uuv11NPPVXi8wKAVxkAgMv69OljJDn97NixwyxatMhIMv/++68xxpgPPvjAhIeHm2+//dY0btzYhIaGmltvvdWkp6ebqVOnmnr16pmqVauawYMHm5ycHMf5MzMzzWOPPWZq1aplKlasaNq2bWsWLVpUZJ8kma+++sppW1JSkrnqqqtMWFiYqVy5srn44ovNqlWrHPt//vlnc/nll5uQkBBTp04dM3jwYHPs2DGnfjzxxBOmTp06Jjg42DRs2NBMmTLFsX/x4sXmkksuMcHBwSYyMtIMHz7cnDhxwrG/Y8eOZvDgwebxxx8355xzjomIiDAJCQlOffzjjz/MFVdcYaxWq2nSpIlZsGCB07VkZWWZhx9+2ERGRhqr1Wrq1q1rXnjhhULjsGjRInPJJZeYihUrmvDwcHPppZeanTt3GmOMSUhIMC1atHC07dOnj7npppvMyy+/bCIjI021atXMQw89ZLKzs0scgw0bNpjrr7/eVKpUydSsWdPcc8895sCBA4X/ofKwP/7pOnToYIYNG+a0LT4+3lx22WWFnuu1114zlStXdvz9bDabiYyMNC+//LKjzeHDh43VajWfffZZoeeZO3euCQgIMCkpKY5tb731lqlSpYrJysrK175jx47mkUceKfR8ebVt29Y8/PDDjt9zc3NNrVq1zNixY/O1Pf21BAC+wogTALhhwoQJ6tChg/r37+/4xj86OrrAtsePH9frr7+u6dOna968eVq8eLFuvvlmzZ07V3PnztXHH3+sd955R1988YXjmEGDBmn58uWaPn26fvvtN91+++26/vrrSzTNKq+7775bderU0apVq7RmzRo9+eSTqlChgiTpzz//1PXXX69bb71Vv/32m2bMmKFffvlFgwYNchzfu3dvffbZZ3r99df1+++/65133lFYWJgkac+ePbrhhht0ySWXaP369Xrrrbf03nvv6fnnn3fqw4cffqhKlSppxYoV+r//+z89++yzWrhwoSTJZrPplltuUXBwsFasWKG3335bw4cPdzr+9ddf1zfffKPPP/9cW7du1aeffqqYmJgCrzcnJ0fdu3dXx44d9dtvv2n58uUaMGBAkSOAixYt0p9//qlFixY5RlemTp1aohgcPnxY11xzjVq1aqXVq1dr3rx5Sk1N1R133FGyP1AhsrKyFBIS4rQtNDRUK1eu1IkTJwo85r333lPPnj1VqVIlSdKOHTuUkpLiNKoVHh6udu3aafny5YU+9vLly3XRRRcpIiLCsS0uLk5paWnatGmTy9eUnZ2tNWvWOPUnICBAsbGxRfYHAHzO15kbAJR3BX3TXtCIkySzfft2R5uBAweaihUrmqNHjzq2xcXFmYEDBxpjjNm1a5cJDAw0e/bscTp3p06dzIgRIwrtjwoYcapcubKZOnVqge3vu+8+M2DAAKdtP//8swkICDAZGRlm69atRpJZuHBhgcc/9dRT5vzzzzc2m82xbdKkSSYsLMzk5uYaY07G6PLLL3c67pJLLjHDhw83xhgzf/58ExQU5HSt33//vdO1DB482FxzzTVOj1OYf/75x0gyixcvLnB/QSNO9erVcxrtu/32202PHj2MMabYGDz33HPmuuuuc9q2e/duI8ls3bq12P4WNuI0YsQIExkZaVavXm1sNptZtWqViYiIMJLM3r1787VfsWKFkWRWrFjh2LZ06dIC299+++3mjjvuKLRP/fv3z3dN6enpRpKZO3duvvYlHXHas2ePkWSWLVvmtP3xxx83bdu2zdeeEScA/oIRJwAoIxUrVlTDhg0dv0dERCgmJsYxamHfZr/PY8OGDcrNzVXjxo0VFhbm+Pnpp5/0559/luqx4+Pjdf/99ys2NlYvvvii0/Hr16/X1KlTnR4jLi5ONptNO3bsUFJSkgIDA9WxY8cCz/3777+rQ4cOTqM5l112mY4dO6a///7bsa158+ZOx0VFRTmu9ffff1d0dLRq1arl2N+hQwen9n379lVSUpLOP/98DRkyRAsWLCj0eqtVq6a+ffsqLi5OXbt21YQJE7Rv374iY9S0aVMFBgYW2L/iYrB+/XotWrTIKYYXXHCBJJX6b5XXM888o86dO6t9+/aqUKGCbrrpJvXp00fSyVGa07333nu66KKL1LZt21I9TufOnR39btq0qcv9Pd3PP//sFJNPP/3UY+cGgLIW5OsOAMDZwj41zs5isRS4zWazSZKOHTumwMBArVmzxukDvSSnZKskRo8erbvuuktz5szR999/r4SEBE2fPl0333yzjh07poEDB2rIkCH5jqtbt662b99eqscqTFHXWhIXX3yxduzYoe+//14//PCD7rjjDsXGxjpNbczrgw8+0JAhQzRv3jzNmDFDI0eO1MKFC9W+fftS9y80NLTIvh07dkxdu3bVSy+9lG9fVFRUSS6vQKGhoXr//ff1zjvvKDU1VVFRUZo8ebIqV66sGjVqOLVNT0/X9OnT9eyzzzptj4yMlCTH8XapqamOkuxTpkxRRkaGpFNxiIyMzFfpLjU11emcxWnTpo2SkpIcv0dERMhqtSowMNBxrrznLul5AcAXGHECADcFBwcrNzfX4+dt1aqVcnNztX//fjVq1Mjpx5UPmI0bN9bQoUO1YMEC3XLLLfrggw8knUxINm/enO8xGjVqpODgYF100UWy2Wz66aefCjxvkyZNtHz5chljHNuWLl2qypUrq06dOiXqW5MmTbR7926nUaFff/01X7sqVaqoR48eevfddzVjxgx9+eWXOnToUKHnbdWqlUaMGKFly5apWbNmmjZtWon6c7riYnDxxRdr06ZNiomJyRdD+71G7qhQoYLq1KmjwMBATZ8+XTfeeGO+EaeZM2cqKytL99xzj9P2+vXrKzIy0qn8elpamlasWOEY1atdu7ajv/Xq1ZN0csRvw4YNTpXuFi5cqCpVqujCCy8sUb9DQ0OdYlG5cmUFBwerdevWTv2x2WxKTEzMN8oIAP6ExAkA3BQTE6MVK1Zo586dOnjwYKlGUYrSuHFj3X333erdu7dmzZqlHTt2aOXKlRo7dqzmzJlT4vNkZGRo0KBBWrx4sXbt2qWlS5dq1apVatKkiSRp+PDhWrZsmQYNGqSkpCRt27ZNX3/9taM4RExMjPr06aN7771Xs2fP1o4dO7R48WJ9/vnnkqSHHnpIu3fv1uDBg7VlyxZ9/fXXSkhIUHx8fIHTyQoSGxurxo0bq0+fPlq/fr1+/vlnPf30005txo0bp88++0xbtmzRH3/8oZkzZyoyMrLAhX537NihESNGaPny5dq1a5cWLFigbdu2Oa65tIqLwcMPP6xDhw7pzjvv1KpVq/Tnn39q/vz56tevX5FJ9ebNm5WUlKRDhw7pyJEjSkpKchqh+eOPP/TJJ59o27ZtWrlypXr27KmNGzfqhRdeyHeu9957T927d89XZt1isejRRx/V888/r2+++UYbNmxQ7969VatWrULXjpKk6667ThdeeKF69eql9evXa/78+Ro5cqQefvhhWa1WRzt7n48dO6YDBw4oKSlJmzdvLjKe8fHxevfdd/Xhhx/q999/14MPPqj09HT169fP0SYlJUVJSUmOEc8NGzY4YgUAPuHrm6wAoLzbunWrad++vQkNDS22HHlepxcoMCZ/kYDs7GwzatQoExMTYypUqGCioqLMzTffbH777bdC+6PTikNkZWWZnj17mujoaBMcHGxq1aplBg0aZDIyMhxtVq5caa699loTFhZmKlWqZJo3b27+97//OfZnZGSYoUOHmqioKBMcHGwaNWpk3n//fcf+kpQjP71wwE033WT69OnjFMfLL7/cBAcHm8aNG5t58+Y5XcvkyZNNy5YtTaVKlUyVKlVMp06dzNq1awuMQUpKiunevbujv/Xq1TOjRo1yFKsorBx5Xo888ojp2LFjiWPwxx9/mJtvvtlUrVrVhIaGmgsuuMA8+uijRRazqFevXr5y9nn/ad68ebNp2bKlCQ0NNVWqVDE33XST2bJlS77zbNmyxUgyCxYsKPBxbDabeeaZZ0xERISxWq2mU6dOJSpasXPnTtO5c2cTGhpqqlevbh577DGnv6sxpsD+16tXr9hzT5w40dStW9cEBwebtm3bml9//dVpf0JCQoHn/uCDD4o9NwB4g8WYPHMrAADlnsVi0VdffVXkaAIAACgdEicAOMNYLBaFhITo3HPPdapqBwAAXEdVPQA4w9gXxz29Eh8AAHAdI04AAAAAUAyq6gEAAABAMUicAAAoQzExMbJYLLJYLDp8+LCvuwMAKCESJwCAT02aNEkxMTEKCQlRu3bttHLlSse+zMxMPfzwwzr33HMVFhamW2+9VampqcWec+bMmbrgggsUEhKiiy66SHPnznXab4zRqFGjFBUVpdDQUMXGxjruDSvK4sWLdfHFF8tqtapRo0aaOnVqqa5HklatWqUvv/yy2McCAPgXEicAgM/MmDFD8fHxSkhI0Nq1a9WiRQvFxcVp//79kqShQ4fq22+/1cyZM/XTTz9p7969uuWWW4o857Jly3TnnXfqvvvu07p169S9e3d1795dGzdudLT5v//7P73++ut6++23tWLFClWqVElxcXHKzMws9Lw7duxQly5ddPXVVyspKUmPPvqo7r//fs2fP7/E1yNJNWrUULVq1VwNGQDARygOAQDwmXbt2umSSy7RG2+8IUmy2WyKjo7W4MGD9eCDD6pGjRqaNm2abrvtNknSli1b1KRJEy1fvlzt27cv8Jw9evRQenq6vvvuO8e29u3bq2XLlnr77bdljFGtWrX02GOPadiwYZKkI0eOKCIiQlOnTlXPnj0LPO/w4cM1Z84cpwSsZ8+eOnz4sObNm1fs9Tz55JOO4xYvXqyrr75a//77r6pWrepi9AAAZYkRJwCAT2RnZ2vNmjWKjY11bAsICFBsbKyWL1+uNWvW6MSJE077L7jgAtWtW1fLly93bIuJidHo0aMdvy9fvtzpGEmKi4tzHLNjxw6lpKQ4tQkPD1e7du2cznvVVVepb9++JT5vcdcDACjfSJwAAD5x8OBB5ebmKiIiwml7RESEUlJSlJKSouDg4HwjMvb9dg0bNlT16tUdv6ekpBR6Tvt++7aizlu3bl1FRUUVe960tDRlZGQUez0AgPKNBXABAOVaYmKiV8770UcfeeW8AIDyiREnAIBPVK9eXYGBgfmq5KWmpioyMlKRkZHKzs7OV7Lbvr8wkZGRhZ7Tvt++zRPnrVKlikJDQ4u9HgBA+UbiBADwieDgYLVu3dppxMhmsykxMVEdOnRQ69atVaFCBaf9W7duVXJysjp06FDoeTt06JBvFGrhwoWOY+rXr6/IyEinNmlpaVqxYoVb5y3uegAA5ZwBAMBHpk+fbqxWq5k6darZvHmzGTBggKlatapJSUkxxhjzwAMPmLp165off/zRrF692nTo0MF06NDB6RzXXHONmThxouP3pUuXmqCgIPPKK6+Y33//3SQkJJgKFSqYDRs2ONq8+OKLpmrVqubrr782v/32m7nppptM/fr1TUZGhqNNr169zJNPPun4/a+//jIVK1Y0jz/+uPn999/NpEmTTGBgoJk3b16Jr8du0aJFRpL5999/PRJHAID3cY8TAMBnevTooQMHDmjUqFFKSUlRy5YtNW/ePEeBhddee00BAQG69dZblZWVpbi4OL355ptO5/jzzz918OBBx++XXnqppk2bppEjR+qpp57Seeedp9mzZ6tZs2aONk888YTS09M1YMAAHT58WJdffrnmzZunkJAQR5vk5GQFBJyamFG/fn3NmTNHQ4cO1YQJE1SnTh1NmTJFcXFxJb4eAED5xTpOAACUMdZxAoDyhxEnAADKUNOmTfXXX3/5uhsAgFJixAkAgDK0a9cunThxQpLUoEEDp+mAAAD/ReIEAAAAAMXgay4AAAAAKAaJEwAAAAAUg8QJAAAAAIpB4gQAAAAAxSBxAgAAAIBikDgBAAAAQDFInAAAAACgGCROAAAAAFCM/wdl3hOAgNtgfQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], "source": [ "from matplotlib import pyplot as plt\n", "\n", @@ -150,10 +1150,55 @@ "selected_ds = ds.interp(x=x_sel, y=y_sel, method='linear')\n", "\n", "selected_prate = (selected_ds.prate * 3600) / 25.4 # convert from kg/m^2/s to in/hr\n", - "selected_apcp = selected_ds.apcp / 25.4 # convert from kg/m^2/s to in/hr\n", "\n", - "#selected_prate.plot.line(ax=ax)\n", - "selected_apcp.plot.line(ax=ax)" + "selected_prate.plot.line(ax=ax)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "\n", + "ROTCON_P = 0.622515\n", + "LON_XX_P = -97.5\n", + "LAT_TAN_P = 38.5\n", + "\n", + "angle2 = ROTCON_P*(ds.longitude-LON_XX_P)*0.017453\n", + "sinx2 = np.sin(angle2)\n", + "cosx2 = np.cos(angle2)\n", + "\n", + "un = cosx2 * ds.ugrd + sinx2 * ds.vgrd\n", + "vn = -sinx2 * ds.ugrd + cosx2 * ds.vgrd\n", + "\n", + "wind_speed = np.sqrt(un**2 + vn**2)\n", + "wind_dir = (270 - np.arctan2(vn, un) * 180 / np.pi) % 360" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Time: 2023-07-22T05:00:00.000000000\n", + "Wind speed: 3.69 m/s\n", + "Wind direction: 197.23 degrees\n" + ] + } + ], + "source": [ + "wind_speed_ll = wind_speed.interp(x=x_sel, y=y_sel, method='linear').isel(time=-1)\n", + "wind_dir_ll = wind_dir.interp(x=x_sel, y=y_sel, method='linear').isel(time=-1)\n", + "\n", + "print(f'Time: {wind_speed_ll.time.values}')\n", + "print(f'Wind speed: {wind_speed_ll.values:.2f} m/s')\n", + "print(f'Wind direction: {wind_dir_ll.values:.2f} degrees')" ] }, { diff --git a/python/gribberish/kerchunk/mapper.py b/python/gribberish/kerchunk/mapper.py index 5510cc7..10e513b 100644 --- a/python/gribberish/kerchunk/mapper.py +++ b/python/gribberish/kerchunk/mapper.py @@ -146,23 +146,23 @@ def scan_gribberish( for coord_name, coord_data in dataset['coords'].items(): coord_values = coord_data["values"] - if "offsets" in coord_values: - _store_array_ref( + if isinstance(coord_values, (list, np.ndarray)): + coord_array = np.array(coord_data['values']) + _store_array_inline( store, z, - coord_data['values']['shape'], + coord_array, coord_name, - offset, - size, coord_data['attrs'] ) else: - coord_array = np.array(coord_data['values']) - _store_array_inline( + _store_array_ref( store, z, - coord_array, + coord_data['values']['shape'], coord_name, + offset, + size, coord_data['attrs'] )