diff --git a/tutorials/gfm-explorer.ipynb b/tutorials/gfm-explorer.ipynb new file mode 100644 index 0000000..bce1a1c --- /dev/null +++ b/tutorials/gfm-explorer.ipynb @@ -0,0 +1,1414 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exploring the Global Flood Monitoring data with stac and xarray\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import pyproj\n", + "from pystac_client import Client\n", + "from odc.stac import stac_load\n", + "import xarray as xr" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The pystac_client allows us to explore the data of interest with (or without) a given area and time.\n", + "\n", + "First, we select the dataset, bounding box and datetime. " + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "collection = \"GFM\"\n", + "\n", + "bbox = [15.6, 48.25, 16.2, 48.45]\n", + "\n", + "date = ['2024-09-12T00:00:00Z', '2024-09-28T00:00:00Z']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The client creates a connection to EODC's stac api and finds data that matches the requested extents." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ,\n", + " ,\n", + " ,\n", + " ]" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "client = Client.open('https://stac.eodc.eu/api/v1')\n", + "results = client.search(collections=collection,\n", + " max_items=100,\n", + " bbox=bbox,\n", + " datetime=date,\n", + ")\n", + "result_items = list(results.items())\n", + "result_items[:5]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The stac items contain all the metadata provided." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'gsd': 20,\n", + " 'created': '2024-11-04T16:22:06.314008+00:00',\n", + " 'datetime': '2024-09-27T05:02:14Z',\n", + " 'Equi7Tile': 'EU020M_E051N015T3',\n", + " 'blocksize': {'x': 512, 'y': 512},\n", + " 'proj:bbox': [5100000, 1500000, 5400000, 1800000],\n", + " 'proj:wkt2': 'PROJCS[\"Azimuthal_Equidistant\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Azimuthal_Equidistant\"],PARAMETER[\"latitude_of_center\",53],PARAMETER[\"longitude_of_center\",24],PARAMETER[\"false_easting\",5837287.81977],PARAMETER[\"false_northing\",2121415.69617],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]',\n", + " 'proj:shape': [15000, 15000],\n", + " 'constellation': 'sentinel-1',\n", + " 'flood_members': {'DLR': True, 'TUW': True, 'LIST': True},\n", + " 'proj:geometry': {'type': 'Polygon',\n", + " 'coordinates': [[[5100000.0, 1500000.0],\n", + " [5100000.0, 1800000.0],\n", + " [5400000.0, 1800000.0],\n", + " [5400000.0, 1500000.0],\n", + " [5100000.0, 1500000.0]]]},\n", + " 'flooded_pixels': 19277,\n", + " 'proj:transform': [20, 0, 5100000, 0, -20, 1800000],\n", + " 'anomaly_detected': False,\n", + " 'floodable_pixels': 39537670,\n", + " 'processing:level': 'L3',\n", + " 'processing:lineage': 'Global Flood Monitoring',\n", + " 'processing:version': 'GFM v3.1.0',\n", + " 'processing:datetime': '2024-09-27T08:37:54.957650Z',\n", + " 'processing:facility': 'EODC GmbH, Vienna',\n", + " 'processing:software': {'GFM': 'NRT v3.1.0'},\n", + " 'ratio_after_blob_removal': 0.22337500648441147}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result_items[0].properties" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We utilize odc's stac_load functionality to load the data as an xarray.DataArray. \n", + "\n", + "In the stac_load function, we need to define the coordinate reference system as well as the resolution. Both of these are specified in the item properties." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 101MB\n",
+       "Dimensions:                (y: 1353, x: 2334, time: 16)\n",
+       "Coordinates:\n",
+       "  * y                      (y) float64 11kB 1.651e+06 1.651e+06 ... 1.624e+06\n",
+       "  * x                      (x) float64 19kB 5.214e+06 5.214e+06 ... 5.261e+06\n",
+       "    spatial_ref            int32 4B 0\n",
+       "  * time                   (time) datetime64[ns] 128B 2024-09-13T05:18:08 ......\n",
+       "Data variables:\n",
+       "    ensemble_flood_extent  (time, y, x) uint8 51MB dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>\n",
+       "    reference_water_mask   (time, y, x) uint8 51MB dask.array<chunksize=(1, 1000, 1000), meta=np.ndarray>
" + ], + "text/plain": [ + " Size: 101MB\n", + "Dimensions: (y: 1353, x: 2334, time: 16)\n", + "Coordinates:\n", + " * y (y) float64 11kB 1.651e+06 1.651e+06 ... 1.624e+06\n", + " * x (x) float64 19kB 5.214e+06 5.214e+06 ... 5.261e+06\n", + " spatial_ref int32 4B 0\n", + " * time (time) datetime64[ns] 128B 2024-09-13T05:18:08 ......\n", + "Data variables:\n", + " ensemble_flood_extent (time, y, x) uint8 51MB dask.array\n", + " reference_water_mask (time, y, x) uint8 51MB dask.array" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "crs = pyproj.CRS.from_wkt(result_items[0].properties[\"proj:wkt2\"])\n", + "resolution = result_items[0].properties[\"gsd\"]\n", + "\n", + "data = stac_load(\n", + " result_items, \n", + " chunks={\"x\": 1000, \"y\": 1000}, \n", + " bands=[\"ensemble_flood_extent\", \"reference_water_mask\"],\n", + " crs=crs,\n", + " resolution=resolution,\n", + " dtype=\"uint8\",\n", + " bbox=bbox\n", + " )\n", + "\n", + "data" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The nodata value is set to 255. For further calculation, we will set all 255 values to 0. " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "data_0 = xr.where(data == 1, 1, 0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can calculate the flood frequency by taking the temporal mean. " + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "frequency = data_0.mean(dim=\"time\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'ensemble_flood_extent' (y: 1353, x: 2334)> Size: 25MB\n",
+       "dask.array<mean_agg-aggregate, shape=(1353, 2334), dtype=float64, chunksize=(1000, 1000), chunktype=numpy.ndarray>\n",
+       "Coordinates:\n",
+       "  * y            (y) float64 11kB 1.651e+06 1.651e+06 ... 1.624e+06 1.624e+06\n",
+       "  * x            (x) float64 19kB 5.214e+06 5.214e+06 ... 5.261e+06 5.261e+06\n",
+       "    spatial_ref  int32 4B 0
" + ], + "text/plain": [ + " Size: 25MB\n", + "dask.array\n", + "Coordinates:\n", + " * y (y) float64 11kB 1.651e+06 1.651e+06 ... 1.624e+06 1.624e+06\n", + " * x (x) float64 19kB 5.214e+06 5.214e+06 ... 5.261e+06 5.261e+06\n", + " spatial_ref int32 4B 0" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "flood_frequency = frequency.ensemble_flood_extent\n", + "\n", + "flood_frequency" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The reference water mask represents permanent or seasonal water bodies, which are clearly distinct from flood events.\n", + "\n", + "For the selected area, the reference water mask covers parts of the Danube. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "ref_water = frequency.reference_water_mask\n", + "ref_water = xr.where(ref_water == 1, 1, np.nan)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABHkAAAKGCAYAAADeeEcuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbL0lEQVR4nO3df5hcdX0o/vdsQrJByALSJoKRoKIYwURDEoJK0JsaLffG6L0lUmsg8BURg9BQKiI/qsGCSjReAnK5FfGKFB6s0IgUy40itkSBRB4fGrW05ddFNwGrWQkJIdn9/hFmmd2d3Z3ZmTm/5vXa5zzJzp6Z82POOXM+73l/3p9SX19fXwAAAACQax1prwAAAAAAjRPkAQAAACgAQR4AAACAAhDkAQAAACgAQR4AAACAAhDkAQAAACgAQR4AAACAAhDkAQAAACgAQR4AAACAAhDkAaCtTZ8+PU499dTEl/vYY49FqVSKG264YdR5H3nkkXjXu94VXV1dUSqV4vbbb2/5+gEAkD/j014BAGiFG264IZYvX171b5/4xCfiiiuuSHiNxu6UU06JRx99ND772c/GAQccEMccc0zaqwQA0HQ7d+6MXbt2pbLsCRMmRGdnZyrLbiZBHgAK7TOf+UwcfvjhAx476qijUlqb+u3YsSM2bNgQn/rUp2LFihVprw4AQEvs3LkzJu3/8ojdz6Wy/KlTp8ajjz6a+0CPIA8Ahfae97wn15kvTz/9dEREHHDAAaPOu3379njZy17W4jUCAGi+Xbt2Rex+LibOOCVi3IRkF75nV3Rv/nrs2rUr90EeNXkAYJD/+I//iD/5kz+Jgw46KPbdd9849thj47vf/e6Q+bZu3Rqnn356TJkyJTo7O2PmzJnx9a9/fch8v/vd7+LUU0+Nrq6uOOCAA+KUU06J3/3ud6Oux1/91V/FYYcdFhER559/fpRKpZg+fXr/30qlUmzevDn+9E//NA488MB429ve1v/cG2+8MWbPnh2TJk2Kgw46KD7wgQ/Ek08+OWQZ1113XbzmNa+JSZMmxdy5c+NHP/pRnHDCCXHCCSf0z3PDDTdEqVSKxx57bMBz77nnniiVSnHPPfcMePwnP/lJvPvd746urq7Yd999Y8GCBfHP//zPQ7atVCrFv/3bv8Wpp54aBxxwQHR1dcXy5cvjueeGfoN34403xty5c2PfffeNAw88MI4//vj4x3/8x4jY253t4IMPjhdeeGHI8971rnfF61//+mH3MQCQMeM7o5TwFOPzHdipJMgDQKFt27YtnnnmmQHTSLZs2RLHHXdcfO9734uzzjorPvvZz8bOnTtj8eLFcdttt/XPt2PHjjjhhBPiG9/4Rnzwgx+ML3zhC9HV1RWnnnpqfPnLX+6fr6+vL9773vfGN77xjfizP/uzuOyyy+L//b//F6eccsqo6/7+978/vvSlL0VExMknnxzf+MY3Ys2aNQPm+ZM/+ZN47rnn4q//+q/jwx/+cEREfPazn41ly5bFEUccEV/84hfj3HPPjfXr18fxxx8/ILj01a9+NT7ykY/E1KlT4/Of/3y89a1vjcWLF1cNBtXq+9//fhx//PHR09MTl156afz1X/91/O53v4t3vvOdcf/99w+Z/6STTorf//73cfnll8dJJ50UN9xwQ3z6058eMM+nP/3p+NCHPhT77LNPfOYzn4lPf/rTMW3atPj+978fEREf+tCH4je/+U1873vfG/C87u7u+P73vx9/9md/NubtAQDIE921ACi0hQsXDnmsr69v2PmvuOKK2LJlS/zoRz/qz4z58Ic/HG9605ti5cqV8d73vjc6Ojriuuuui5///Odx4403xgc/+MGIiDjzzDNjwYIFcdFFF8Vpp50W+++/f6xbty7uvffe+PznPx/nn39+RER89KMfjXe84x2jrvub3vSmmDx5cvz5n/95vOUtb6karJg5c2bcdNNN/b8//vjjcemll8Zll10WF154Yf/j73//++PNb35zXHPNNXHhhRfGCy+8EBdeeGHMmjUrfvCDH8SECXvTomfMmBFnnHFGTJs2bdT1G6yvry/OPPPMeMc73hH/8A//EKVSKSIiPvKRj8Qb3/jGuOiii/qzb8re/OY3x1e/+tX+33/zm9/EV7/61fjc5z4XERH/9m//Fp/5zGfife97X3zrW9+Kjo6OAcuLiHjnO98Zr3zlK+PGG2+M//pf/2v/3//2b/82ent7BXkAIE9KEfHiPUSiyywImTwAFNrVV18dd99994BpJHfeeWfMnTt3QNen/fbbL84444x47LHHYvPmzf3zTZ06NU4++eT++fbZZ5/4+Mc/Hs8++2z88Ic/7J9v/Pjx8dGPfrR/vnHjxsXZZ5/dlO0788wzB/z+7W9/O3p7e+Okk04akL00derUOOKII+IHP/hBREQ8+OCDsXXr1jjzzDP7AzwR0d+tbCweeuiheOSRR+JP//RP4ze/+U3/srdv3x7/5b/8l7j33nujt7d3xPV/+9vfHr/5zW+ip6cnIiJuv/326O3tjUsuuWRAgCci+oNIHR0d8cEPfjDWrVsXv//97/v//s1vfjOOO+64IYW3AQCKSiYPAIU2d+7cugovP/744zFv3rwhj7/hDW/o//tRRx0Vjz/+eBxxxBFDAg+V85X/fcUrXhH77bffgPmaVSdmcADjkUceib6+vjjiiCOqzr/PPvsMWL/B8+2zzz7x6le/ekzr8sgjj0REjNgVbdu2bXHggQf2//6qV71qwN/Lf/vtb38bkydPjn//93+Pjo6OmDFjxojLXrZsWXzuc5+L2267LZYtWxa//OUvY+PGjXHttdeOaVsAAPJIkAcAcmzSpEkDfu/t7Y1SqRT/8A//EOPGjRsy/+BgUy1Kw6RM79mzZ8iyIyK+8IUvxKxZs6o+Z/Dyq61jxMhd6qqZMWNGzJ49O2688cZYtmxZ3HjjjTFhwoQ46aST6nodACBlpY69U9LLLAhBHgCocNhhh8Uvf/nLIY//4he/6P97+d+f/exn0dvbOyCbp9p869evj2effXZAgKPaMprhNa95TfT19cXhhx8er3vd64adr7x+jzzySLzzne/sf/yFF16IRx99NGbOnNn/WDm7ZvCIYOVsoMplR0RMnjy5ai2ksXjNa14Tvb29sXnz5mEDR2XLli2LlStXxq9//eu46aab4sQTTxyQNQQAUHTFCVcBQBP88R//cdx///2xYcOG/se2b98e1113XUyfPr2/29Af//EfR3d3d9xyyy398+3evTuuuuqq2G+//WLBggX98+3evTu+8pWv9M+3Z8+euOqqq1qy/u9///tj3Lhx8elPf3pINkxfX1/85je/iYiIY445Jv7gD/4grr322ti1a1f/PDfccMOQYE45eHPvvfcO2IbrrrtuwHyzZ8+O17zmNXHllVfGs88+O2Tdnn766bq3Z8mSJdHR0RGf+cxnhtTzGbx9J598cpRKpTjnnHPiP/7jPxRcBoA8KpXSmQpCJg8AVLjgggvib//2b+M973lPfPzjH4+DDjoovv71r8ejjz4af/d3f9eftXPGGWfE//pf/ytOPfXU2LhxY0yfPj2+9a1vxT//8z/HmjVrYv/994+IiP/23/5bvPWtb40LLrggHnvssZgxY0Z8+9vfjm3btrVk/V/zmtfEZZddFp/85CfjscceiyVLlsT+++8fjz76aNx2221xxhlnxF/8xV/EPvvsE5dddll85CMfiXe+852xdOnSePTRR+NrX/vakJo8b3zjG+PYY4+NT37yk/Gf//mfcdBBB8XNN98cu3fvHjBfR0dH/M3f/E285z3viTe+8Y2xfPnyOPTQQ+Opp56KH/zgBzF58uT4zne+U9f2vPa1r41PfepTsWrVqnj7298e73//+2PixInxwAMPxCGHHBKXX355/7x/8Ad/EO9+97vj1ltvjQMOOCBOPPHEse9IAIAcEuQBgApTpkyJ++67Lz7xiU/EVVddFTt37ow3velN8Z3vfGdA0GDSpElxzz33xAUXXBBf//rXo6enJ17/+tfH1772tTj11FP75+vo6Ih169bFueeeGzfeeGOUSqVYvHhxrF69Ot785je3ZBsuuOCCeN3rXhdf+tKX4tOf/nREREybNi3e9a53xeLFi/vnO+OMM2LPnj3xhS98Ic4///w4+uijY926dXHxxRcPec1vfvOb8ZGPfCSuuOKKOOCAA+L000+Pd7zjHfFHf/RHA+Y74YQTYsOGDbFq1apYu3ZtPPvsszF16tSYN29efOQjHxnT9nzmM5+Jww8/PK666qr41Kc+Ffvuu2+86U1vig996END5l22bFnccccdcdJJJ8XEiRPHtDwAgLwq9dVb2RAAKLQTTjghIiLuueeeVNdjLP7+7/8+lixZEvfee2+8/e1vT3t1AIAa9fT0RFdXV0x8y4oojUv2i5q+Pc/H85vWxrZt22Ly5MmJLrvZ1OQBAArjf//v/x2vfvWr421ve1vaqwIAkDjdtQCA3Lv55pvjZz/7WXz3u9+NL3/5y8MO+w4AZFwahZALdN8gyAMA5N7JJ58c++23X5x++ulx1llnpb06AACpEOQBAAbIYy0eJQYBoCg6IkpJV5YpTiWb4mwJAAAAQBsT5AEAAAAoAN21AAAAgGxQeLkhNQd5du5u5WoAAAAAZZ1SMjLt6quvji984QvR3d0dM2fOjKuuuirmzp1bdd5/+Zd/iUsuuSQ2btwYjz/+eHzpS1+Kc889d8A8l19+eXz729+OX/ziFzFp0qQ47rjj4nOf+1y8/vWvr2u9dNcCAAAAsqHUkc5Uh1tuuSVWrlwZl156aWzatClmzpwZixYtiq1bt1ad/7nnnotXv/rVccUVV8TUqVOrzvPDH/4wPvaxj8WPf/zjuPvuu+OFF16Id73rXbF9+/b6dl9fjcNRyOQBAACAZLRbJk9PT090dXXFxDkrozR+YqLL7tv9fDz/wBdj27ZtMXny5FHnnzdvXsyZMyfWrl0bERG9vb0xbdq0OPvss+OCCy4Y8bnTp0+Pc889d0gmz2BPP/10/OEf/mH88Ic/jOOPP77mbZHJAwAAALS9np6eAdPzzz8/ZJ5du3bFxo0bY+HChf2PdXR0xMKFC2PDhg1NW5dt27ZFRMRBBx1U1/MEeQAAAIBsKBdeTnqKiGnTpkVXV1f/dPnllw9ZvWeeeSb27NkTU6ZMGfD4lClToru7uym7oLe3N84999x461vfGkcddVRdz22zBDAAAACAoZ588skB3bUmTky221jZxz72sXj44Yfjn/7pn+p+riAPAAAAkA1jKITclGVGxOTJk0etyXPwwQfHuHHjYsuWLQMe37Jly7BFleuxYsWKuOOOO+Lee++NV77ylXU/X3ctAAAAgBpMmDAhZs+eHevXr+9/rLe3N9avXx/z588f8+v29fXFihUr4rbbbovvf//7cfjhh4/pdWTyAAAAANlQUSMn0WXWYeXKlXHKKafEMcccE3Pnzo01a9bE9u3bY/ny5RERsWzZsjj00EP7a/rs2rUrNm/e3P//p556Kh566KHYb7/94rWvfW1E7O2iddNNN8Xf//3fx/77799f36erqysmTZpU87oJ8gAAAADUaOnSpfH000/HJZdcEt3d3TFr1qy46667+osxP/HEE9HR8VLHqV/96lfx5je/uf/3K6+8Mq688spYsGBB3HPPPRER8ZWvfCUiIk444YQBy/ra174Wp556as3rVurr6+urZcadu2t+TQAAAKABnW2WktHT0xNdXV0x8di/jNL4ZAse9+1+Pp7/8edj27Zto9bkybo2O2wAAACAzEqx8HIRFGdLAAAAANqYTB4AAAAgG0qlFDJ5Ei703EIyeQAAAAAKQJAHAAAAoAB01wIAAACyoaO0d0p6mQUhkwcAAACgAGTyAAAAANlgCPWGFGdLAAAAANqYTB4AAAAgG0ql5Ic0N4Q6AAAAAFkiyAMAAABQALprAQAAANmg8HJDirMlAAAAAG1MJg8AAACQDQovN0QmDwAAAEABCPIAAAAAFIDuWgAAAEA2KLzckOJsCQAAAEAbk8kDAAAAZIPCyw2RyQMAAABQADJ5AAAAgGxQk6chxdkSAAAAgDYmyAMAAABQALprAQAAANmg8HJDZPIAAAAAFIBMHgAAmmbS+Igdu9NeCwDyK4XCywXKfynOlgAAkAmTfI0IAKkQ5AEAoGkGZ/FMGi/oA7x0LXA9gNZyigEA0BIacwDUTeHlhsjkAQCg6SoDPGr0QHtzPYDk+H4FAICmqmzEyeYBBHaoS6mUfOFlmTwAADAyAR4ASJaPXgAAACAbSikMoZ74kO2tI8gDAEBL6KIBAMkqTrgKAAAAoI3J5AEAAACywRDqDZHJAwAAAFAAMnkAAACAbFB4uSHF2RIAAACANibIAwAAAFAAumsBAAAA2aDwckNk8gAAAAAUgEweAAAAIBsUXm5IcbYEAAAAoI0J8gAAAAAUgO5aAAAAQDYovNwQmTwAAAAABSCTBwAAAMiEUqkUJZk8YyaTBwAAAKAAZPIAAAAAmSCTpzEyeQAAAAAKQJAHAAAAoAB01wIAAACyofTilPQyC0ImDwAAAEAByOQBAAAAMkHh5cbI5AEAAAAoAEEeAAAAgALQXQsAAADIBN21GiOTh1yYJBwJAAAAI9J0JjfKgZ4du9NdDwAAAFpDJk9jZPKQeYOzeGT1AAAAwFCay2Tejt0COwAAAO1AJk9jZPKQC7poAQAAwMjkR5BJMncAAACgPprSAAAAQDaUXpySXmZB6K5FLsn0AQAAgIE0lcklNXoAAACKR+HlxsjkoS3JBAIAAKBoBHnIhR27X5qaRaAHAACAItHMJXOGC74IygAAABRbqRQpdNdKdnGtpNlM5pSzdSqDOs0O8KjpAwAAQNEI8pBZO3bL3gEAAGgnpUih8HKBUnnU5AEAAAAoAHkSZFplt6pyVo+uVsBYTRr/UpagawkAQPYYQr0xgjzkRhYaZJXdx7KwPkDtyudv5b/OYwAAikR3LRijSeMHTuSH94syxwIAAEXi9hbqUO1bf43EfBj8Pg3+XUZHe/K+AwBkTCmSr4NcnN5aMnmgUTt2aygWQT0ZWQJ7AABAFmmqAFQYrcC3ei4AANBCKRRe7lN4GdpHrVkbGvzZVh5RqVa1FNk24hsAAJAlumvBCEYKCpS7aWng50ej79dwx4PuW/nk3AUAoGg0TWAE5UZgtUb84O46uu+0N1k9+eOcBQDInlIK3bWSXl4rCfLAKEbL0pDFkW8a+e3Dew0AQNFpnsIYaTC2n3rr+gAAAPWRydMYNXlgGPUMqV35HLKtloLKtTwXAAAgazRZYBiyNvKrVbVWRqrRJLMLAACaoPTilPQyC0ImD4xgpIa7AFC2lTOxBk8RYx9l67x1mwc8X2AHAADIEkEeGEVlg75a415DP3tGe09GCgCNphzoqVyWYwAAANrL1VdfHdOnT4/Ozs6YN29e3H///cPO+y//8i/x3//7f4/p06dHqVSKNWvWNPyawxHkgTEarmEvwycbxhJ4GSngMzi4AwAANF+58HLSUz1uueWWWLlyZVx66aWxadOmmDlzZixatCi2bt1adf7nnnsuXv3qV8cVV1wRU6dObcprDkeQB6CKkYJ1Aj4AANC+vvjFL8aHP/zhWL58ecyYMSOuvfba2HfffeP666+vOv+cOXPiC1/4QnzgAx+IiRMnNuU1hyPIA00wXFeuvChq9lEz35PVi2dU/T8AANA8aWby9PT0DJief/75Ieu3a9eu2LhxYyxcuLD/sY6Ojli4cGFs2LBhTNvczNcU5AHaQrVaSiMFf/IYrAMAAMZu2rRp0dXV1T9dfvnlQ+Z55plnYs+ePTFlypQBj0+ZMiW6u7vHtNxmvmZBv78HalXULJ7hDA7eVP4+2r6QwQMAAMX15JNPxuTJk/t/H65rVZa1WfMOGEzGykvyvC9GClDlebsAAGgvYymE3IxlRkRMnjx5QJCnmoMPPjjGjRsXW7ZsGfD4li1bhi2qPJpmvqbuWgAFIJADAACtN2HChJg9e3asX7++/7He3t5Yv359zJ8/P/XXLHwmT7Vvt5vVGCq/tsYVkLTKa5trEAAARZFmJk+tVq5cGaecckocc8wxMXfu3FizZk1s3749li9fHhERy5Yti0MPPbS/ps+uXbti8+bN/f9/6qmn4qGHHor99tsvXvva19b0mrXKbZBn0vjaGjbNDOgMV8tDNwmgmUYLIA++5gx3PXT9AQCA5lu6dGk8/fTTcckll0R3d3fMmjUr7rrrrv7CyU888UR0dLzUcepXv/pVvPnNb+7//corr4wrr7wyFixYEPfcc09Nr1mrUl9fX18tM+5MobFQb0HYwQ2arBSUrbZeGl/AcAZfu0a6hlQGhIa75rneAADkT2dG2rNJ6enpia6urpiy/BvRMWHfRJfdu+u52PK1D8W2bdtGrcmTdYkdNlkJuCSt1m/iNcKAsagloxAAAGgPLWkWaGzUT9AHiKh+/Ryc/VdrV1VZgwAA0F6aFo4R2KlNrftJUWdoTyN1u6rlefUGgwAAIEvyUHg5y8YcmhHUGZt6G3CCPVAs9ZzT9V4vBs/vugEAAO2lrlBNqwM7lQ2SkYY+r3dkrSwFpMa6LrpdwMiqFSDO2jlTef6PtI5jzcapddQtAADIKpk8jUkt/DFaw2Okv9fbaKll/iwFgoaT1YYrZEkezuVWGCkwDgAAtIe6mkNjrRVRfm6WNbp+STawdMeAofJwLmQxuxAAACiOljc18tDwaoa0tlN3DMifVox8NVp3VwAAyAPdtRrT1KaAYMPoWtEQk9kD+dPMc3XwNcB1AAAA2lPdYYZWfAPdrgbvw2YEfdTtgfYm6AsAQK6VXpySXmZBdIzlSRoOrVH+Br4Z+1d3DQAAAGgvQgEZ1YwCrTKuAAAAyBM1eRozpkyevMlzVkujmT2Txg/d/tH2R7XnANlVvkZUZgMK8AIAQPtpiyBPI0O/Z0Uzgj21zJP3/QTtSlAHAABoqyZ9EbovNdKNq3L7R9sXed9PAAAA5I/uWo1pi0yeiOIFLca6PSMFh3TxAAAAgPxqq0yerAQwqmXRDB76fHAwptq6N6M483CKkPUEAABAvpQihUyeAo2h3lZBniwZKYhTb9BmrDWHBHIAAACgONqmu1aWVBv9phykGWtWzli7WlVbnsAPAAAA5I9MnpyoDMYkFYSR6QMAAECSFF5ujEyelFVm4NSaxdPs4smGTQcAAID807zPmFrr6wwu1DzW1xlpPcrLkc0DAABAIkovTkkvsyBk8mREvRk4tczbaHCm0TpBAAAAQHIEeXJkcMHmVgd6an3upPECQQAAAJA2TfMcGWvXqVq6bo322sP9vfJ1de0CAACgEQovN0YmT44kkZXT6OvK6AEAAIB0aJJnUFrZMNUycSqzgGrN1JHRAwAAwFjI5GmMTJ42MtbAy2g1gAR0AAAAIH0yeRhguCyc0QI5jQ7ZDgAAAKXS3inpZRaFTJ42U2/WTXnkrFoCOJWvLeADAAAAyRLkaUNj7V5Vb6AHAAAASI58C4ao7LJV2Q2r1qLLAOSX4vkAQJr2dtdKuvByootrKZk8bWq0G/jKYM2O3QI8AAAAkHWa5SRmcBDIN8UA2ePaDACkKoXCyyGThyKoJ5un3tfTSAAAAIBkCfK0uWYHY4br2lUtYKR7FwAAADSPZjYjUoATAACApJRKpRQKLxenv5ZMHhIJ4gy3DNk8AAAA0Bya2IyqWdk8tXbjAgAAoD2VUii8XKBEHpk87JVWl6xGlztp/EsTAAAAtDNNY2qSh9o81QI9WV9nAAAAXtLRUYqOjmRTa/oSXl4rCfLQb8fu/GXEDA7i5G39AQAAoFk0iRlgpEBPHrJ5sr5+AAAA0CqCPNQlD4EeAFrLZwEA0CoKLzdG4WUAoGa6xQIAZJdbNYbIY20eAJIjiwcAaJVSqRSlhFNrkl5eK8nkoaqRbuAFgADalwAPAEB2CfIwJgI9AAAAkC2a6gxLty1aoXxMyQYAAAAGU3i5MZrwjFmRR1epN7hV1P3QbJX7tcjHD7CXoC4AQLIEeRhRu2XzjHVbqz2vct9p4FSX9/1z3rrNA35fvXhGSmsC6ag1WJv3cx0ASI7Cy41po+Y7rVCEbIxWBbEGZ61Uyvs+a7a8HUeDgzuVj7cq0FN5DOVpX1F8eTt/AQCKTJCHUY2WzZPHG/y0s5M02IfKw3E0XHCnlRk81Y7VPOwriq9aILvyuEz7OgtAfskAbW8yeRrjFoymyFOjM2sNjzztO/ZKo1tW+RgpH7+OG7Ku3br7AtBc7nNgbAyhTk3qqbmQVZPGZ3cds7pezZT3bSxn8SQZ4Bmu1lN5gqyrPE4dswDUymcGjF3Om10kqZZvZbOYXZCX4IIuXNmWdlHlLJ5btK/hApDVOG4BgHoYQr0xOWn+khV5C/Q0I8CTRhZTlvZhs+i6Ub+iHQMAAEBraXLREmkHKRoNJtS77iPN3+iw7EVq6Av0QDFU1ogq0jUKAEhfKVIovBzFSeXR3KJutTbUkw5SJB3YGevr1rueRQv2VBt9pyjbBu2m1eeuawQAQH0EeRiTejIyWn2TntXgTi3La/esFg03YDiDh2h3vQAAGF2bNzFpRL1db5oZ7Emq1k6r1Rssy8I6A6TBNRAA2oPCy40R5KEhY6mxMtz8w928NzPbJYsNhHoye3RdANqV6x4AwOgEeWhYs4rptrLrUl4aB7Xsy7xsC0AjKos7AwDto1RKofBygVJ5OtJeAYohq4GHHbuzu27DGWl987YtAI3K43UcACAtvh+jabL0rWu9DYKsdYOqltFjVCoAAKDo1ORpTAaa4xRNs7pvNbL8WlRbx8GFPYcLtFQ+3sqh16sFc4w4A803UlAVAADyQpCHlkhjmPB6GmX1Fn8e6XnVAj7NCr4Mfo0sZElBO5AtBwBAHmky0nKNBHxamSVTVk+h41pHvyr/v5yN08ph42XzQGs4rwAAkqfwcmMEeUhUFhtNtWbvVBb/HCnYU9nNKonMm+G6dQH1y1JtMQAAqJfbWFqmCBkmw61/GtuVdq0jaAfl61ber10AAHml8HJjNBlpqZEKFBchCDScJIoxVyOjBxpX5GsTAADF1pH2ClBc1YI6ld2YNKLGppb9JuMH6qPIOQAAReA2lpbKUnenIhlc0Llaho9AGtTH+UIrybQEgNoovNwYmTyQQbVkEdQy6ldl9hQA6RDgAQCSovkHGVTP0Ou1BHFk9QAAALmQQuHlKE4ijyAPZFWtQZlaR93yTXJ19gvQKq4vAEDSBHkg42rJwjG8+tjYZ0CrCfAAQH3U5GmMmjyQA/XW6Gn0tdqB/QC0mgAPAJA0QR7IuHoaCQI9tam2/e2+T4DmKRe9r5wAAJLgtgPalGLM1MsxA6MbabTDCOcQAIymlELh5QL11hLkgTyoZ7St8jy+OR6eRhYAAFBEmoGQE60ITPhmOTmDg2552+eChlCbvJ3bAJA1Ci83xm075FAtwRkjbqWj1n0+0nwaiQAAwFgovAw5UxkcGC2goBBzMppdXDULBVsrl13PMQcAAKTH7TrkzOAMndFq9dSa0aOobn2SDHZUW1Yr36usBXVqWQfHLnnl2gsAAym83BiZPJBDO3a/NEU0N4OEkWVlOORWrkPlsZUFtazLSO9LFt4vqKZathwAQCPcVkDO1ZKpoz5P47K4/yrXKYmgTJqBn3oy0kZ6PEvBK6gkowcA9lJ4uTEyeSCHxjJSUzvV52lmPZusZO6MphU1fCozxrLQ+By8PmNZp7y8n+RfLbW1Ko/hLJxjAEDtrr766pg+fXp0dnbGvHnz4v777x9x/ltvvTWOPPLI6OzsjKOPPjruvPPOAX9/9tlnY8WKFfHKV74yJk2aFDNmzIhrr7227vUS5IE20g6BnmZlt+Q5GJDnda9XtcBPVoJSAADUr5zJk/RUj1tuuSVWrlwZl156aWzatClmzpwZixYtiq1bt1ad/7777ouTTz45Tj/99PjpT38aS5YsiSVLlsTDDz/cP8/KlSvjrrvuihtvvDF+/vOfx7nnnhsrVqyIdevW1bf/+vr6+mqZcacbZsiMRtL6a238562R3KzixEUPjjTjfc3jMTSW7DdoBt2wABirzoLflw7W09MTXV1dMf+z34vxnS9LdNm7d26PDZ9aFNu2bYvJkyePOv+8efNizpw5sXbt2oiI6O3tjWnTpsXZZ58dF1xwwZD5ly5dGtu3b4877rij/7Fjjz02Zs2a1Z+tc9RRR8XSpUvj4osv7p9n9uzZ8Z73vCcuu+yymrdFJg/kUCM1doqazTPWTI5mDVc+UkZJvVMrDd7esUzNXlYSBu/bvB3f5JcADwDkR09Pz4Dp+eefHzLPrl27YuPGjbFw4cL+xzo6OmLhwoWxYcOGqq+7YcOGAfNHRCxatGjA/Mcdd1ysW7cunnrqqejr64sf/OAH8a//+q/xrne9q65tEOSBNqTR0bwAQysCM+3W5SjJoI9ADwBAtpWHUE96ioiYNm1adHV19U+XX375kPV75plnYs+ePTFlypQBj0+ZMiW6u7urblN3d/eo81911VUxY8aMeOUrXxkTJkyId7/73XH11VfH8ccfX9f+c4sLOZVEAKAoIxI1uzGf5P4YvKx2CEy0umtVZSZcUY5xAAAa9+STTw7orjVx4sTEln3VVVfFj3/841i3bl0cdthhce+998bHPvaxOOSQQ4ZkAY2kDZoLQDX1dPnKY02JVgRDsrAPyuvQDsGeslYEYgbvxzwe4wAARZTmEOqTJ08etSbPwQcfHOPGjYstW7YMeHzLli0xderUqs+ZOnXqiPPv2LEjLrzwwrjtttvixBNPjIiIN73pTfHQQw/FlVdeWVeQR3ctaGP1NGqzHlRoZZefLHadyuI6tVqru3Nl/RgHACB9EyZMiNmzZ8f69ev7H+vt7Y3169fH/Pnzqz5n/vz5A+aPiLj77rv753/hhRfihRdeiI6OgSGacePGRW9vb13r55YWCmKsXVzqzeip9/WbKclGeF4CKO1YY6aZx+Hg41/3LQAARrNy5co45ZRT4phjjom5c+fGmjVrYvv27bF8+fKIiFi2bFkceuih/TV9zjnnnFiwYEGsXr06TjzxxLj55pvjwQcfjOuuuy4i9mYQLViwIM4///yYNGlSHHbYYfHDH/4w/s//+T/xxS9+sa51a5MmARTfcI3V8t/qeW4tmh3wyUqAIs+N+yQCPlkbna0ZQZlGRqsDAKC5KgshJ7nMeixdujSefvrpuOSSS6K7uztmzZoVd911V39x5SeeeGJAVs5xxx0XN910U1x00UVx4YUXxhFHHBG33357HHXUUf3z3HzzzfHJT34yPvjBD8Z//ud/xmGHHRaf/exn48wzz6xvW/r6+vpqmXFnjhs+0I7GktnTzg3dPAd3alGuOTPccZFkTZq0A1DDSTtTDQCgUmeb3Zv39PREV1dXvO2Kf4zxnS9LdNm7d26Pf7rgXbFt27ZRa/JkXZsdNtA+xtJIbZeMhrQa8GkW9y0vd7jlpzliWFmjx14zu1rptgUAkI40Cy8XgcLLwABFbdSWCxW3Q4Anr4G6Zr1HY93+asvN674EAKA9uX0FhihCRk9Rg1XtotH6QmPNxKl27BteHQCAvMh5Mw5olcpaLXmQ9UZ4FrpD5VUjx+JYAjRFCHICtTtv3eaIiFi9eMaA38vKjwOQjFKkUHg52cW1lNtYYERZC/bkMYCRRiZIEbNPxnIsjnUfVButrmj7ExhqcIAHAPImI802IOtaEewZ3CWn3LDWmG5ckfdhGpk2Rd6fAABZ0lEqRUfCqTxJL6+VBHmAurSqsTva6E95VsRtSlutgZ5GgobD1QXyfkJx6JoFQNEI8gAtIyuHVqono6eRYzErXRWB5hPUAcieUimFmjzFSeQxhDrQfJPGaxiTjFoCN40ei4OXUT6+HeMAAGSNIA/QdDt2vzQNRyOZZqk1Q2e4462W43C4ZTiOAQDIEremQGp05aJZGinGXPnckY5JQ6sDALReqVSKUsL9p5JeXivJ5AESp1YPrTDaMVXt74ODO6Nl5oyU0QMAAGlzWwokSoCHVqoM1oyk8u+VgZ7y/0fK7Bkuo8exDQDQuI7S3inpZRaFTB4gMbIdSMpoNaGqGS5wM9zr1/oaAACQFEEeAAqvkQwbgR4AAPLCrSjQcpWNXt1ZSFKzjrfhum/pugWNOW/d5oiIWL14RsprAkBmlFIohKy7FkBtBHgoFzOuNuVNtXWW0QNjd/2qa+L6VdfEgXNW9Ad8AICxE+QBmq5aI16Ap32V6+NUm/KongCVQA+M7LSLz+r///WrrklxTQDIilIpnakoBHmAlstrY57iasYxKYgJzVUZ8AEAxkaQB2iqwZkLGr9kVTMyimoJ9MjmgeGpxQPAYKWUfopCkAdomsENXgEe8qKRbmR5rjEEWVDO4BHwAYDGCfIATaGBCwAAkC7NMqDpZPBQJMMNk87I1CyiVgouA1Cpo7R3SnqZReG2FWgqjTmKYPBxXPn7WAI+k8YX+9wo75PhAmJF334a89sH1qa9CpAI10IgCYI8AFCH8g267J6hhtsn5eCPxg3Q7oqe5VgZ9IexKpVKUUp4TPOkl9dKblGBpvBhTruptxtXkW98BwdxBjdiBMQAhl4Pq10bi/gZASTLbRcAjJEAxktqGUJeNg9QafD1sx2uD6NlgxY90wdoPbemANCAsWT0FP3GfbhGTNG3GxhK0GJkgz9D8t4l2JcfNEOptHdKeplF4RQEgAa5qa3OfoH2Ntr5L+jTPtfIIndZhqxpk8sKALSWgEZ1buihfTn/h1dt3xQl6ynP6042dJRK0ZFwak3Sy2slt6MAAEDTFSVokZSi7SNddiEdHWmvAABkVb2ZObXewMr4aZx9CPmhcU+E6za1K9fkSXoqCqcaAAzSyI2oblvJqNzPGpC0m/PWbe7//+rFM1Jck5E5N9uTz0BIl0weAGgyDZtk7NhtX9M+8hLYASBd4qwAMAwBBCBLyoGe1YtnCPQAhVUqlaKUcP+ppJfXSoI8ABDN7/qj2xYA7Wi4kcN8cQLJcPsJABUEZoAsKmfuVHbbgrwQ4KEeaRRCLlAij5o8ADBYs2q9uKkFmk03LQBG4vtKAIjWBWR02wIAICluOwEgBeoTAAAM1VEqRUfC/aeSXl4r6a4FAAAAUAAyeQCgxXTZApJSOcw6QB6VXpySXmZRyOQBAIACuX7VNWmvAgApEeQBgARUq78juwdotnIGz4FzVsSBc1akvDYA9SuVSqlMRSHIAwApEugBmuW8dZuHBHYEegDaiyAPANSgGcEYo2kBrVStDs9pF5+VwpoAkBZBHgCoQSuLJ8vmAZrltw+sTXsVABrSUUpnKgpBHgBIkGweoNUqAz1G2QJoL747BIAalbN5WhGoadXrAu3ptw+s7R9OHSBP0iiErPAyALSpRgMxI3XN0m0r/yrfw0njX5qg1Q6csyLOW7d5QGBHFg9A+xHkAYAEjRYkEhDIv/J7KDOLpJQDO9evuiauX3VNymsDQJoEeQAgYRr/xbVj98D3d/Dv0AqDM3Z00wLyrlRKdioS3xcCQApGGq1LfZ762We0u3Kx5QPnrBjyt2qPnXbxWbpzARSQIA8AZJCgRX1aWRQb8mTwEOqyeoC8UXi5MbprAUAMLJKbVF0cAYnmKgd61DWCl8jWAWgvboMAIKoHXJLIDNFtq7nK+8u+g5cMzu4BfE5kWUdp75T0MotCJg8ADGOkAExS0l5+XmXhvQMg22R/UkSCPJCApLuAAM2TRLDAN4mt47oLQDXlz16fwRSNWx9osWoNjPJjPlQgH5Io6qvbVvNV1uix/wDaV+Xna+Xngc+GbFJ4uTEyeaDFBn94VH7I+IYZ8sONYD553wDaV7VMevffFJ0gDyRgx+6XJoCxcFM6durzALSnwffe7sXzoZTSVBSCPJCQat8k+KABKo12TRgtUCGQMTL7B6D9DP6y1f03RSfIAy1UblBoWAC1qjXQU1nQvTKI7HoDANC+3ApCi43U4FIMFBiL8rXDNaR2ijADQD50lErRkXAh5KSX10qCPJACDQxgJLXUkBGsAABgMEEeSEBlQ0xXCqAWAj0AQDsqlfZOSS+zKNTkgRaqNqKWBhlQq1quFwLHtXP9BbKkWl01gEa5nABAhtWa0VOel5GpZQRkRfk65JoEA5VKpSglnFqT9PJaSSYPAGScm3+A4sr6NV62EeSLIA8A5ECtXbfchI/MUPMA9alWfgDILrc4AJATtXTdipD6P5Ja9yEAA/lcISkKLzdGJg8A5EitN9myegAA2o/bPwDImcpinaNRlBkAyJOOUik6Ek6tSXp5rSSTB5rMt+dAUuoJ3LguDWWfAABFI8gDLSLYAyShXBBTYebalfeV7CYAoGjc6kGLaDwASau87owUzNGFa692LlBdeXy06z6oVEvw036iHbXTdbLys9E1Ml0KLzdGkAcACqiWuj3tHOxpx21meKOdL2kcLxqZpK18DLbDsSjLlSJxOEOTFfXDD8inWoYMb6dvamEk9RQ1B4qn8rPQ52J6SqVSlBJOrUl6ea2kJg8AFFwtNXvU62kvtXbta1dZaOjVU28LWqGyflmRj8VqWa1F3VbagyAPALSJWoM9tBeNmeqK3KgdK9eH9lP0c8AxTRE5rAGgzYyWxTH4saLf5Lcr7yu1qrwm6N7ZPtqlblvRty+POiL5bJQiZb8UaVsAgDrVM/S6bzwB2lORPwOKul20L4c0ALS5eorNtss3u8BL1HBqjrxmQRU9k6to21MECi83xmUaAIiIsQV7qj0fKC7n+VBFDHxEVN+mogX5ivi+QcFOUwCgUbUMu16NLJ/2VdRGbl55P5JRT4ZLed68vzd5Xnfyo1SK6Eg4saZAiTyCPADAUGMN9ETI8qlFUQJiRdkOqFfRMlqA4nB5AgCqqqf71mjy/u11MxUpMFKEbYAkVAbOi349FOiHdAnyAAAjaiSrp1LRGza1KFKApx3k8ZjN4zrnUaP7uMjvU7M+M+rh2losHSl010p6ea1kCHUAYFQ7djfn5rnIw/DWoln7kdbK63Gax3VuJ859KJarr746pk+fHp2dnTFv3ry4//77R5z/1ltvjSOPPDI6Ozvj6KOPjjvvvHPIPD//+c9j8eLF0dXVFS972ctizpw58cQTT9S1XoI8AJCSPDbIykGKymks8tqIptgqj8tqx3f578NNaRNEyL5Gr515Ubl9rd5WWTzFUx5CPempHrfcckusXLkyLr300ti0aVPMnDkzFi1aFFu3bq06/3333Rcnn3xynH766fHTn/40lixZEkuWLImHH364f55///d/j7e97W1x5JFHxj333BM/+9nP4uKLL47Ozs769l9fX19fLTPudNIAwLDqTb2vbBAW6cZ0LA3dIm0/+VXrOTnSMZ6FY1mDl3ZT5GO+MwPB4yT19PREV1dXfOzmB2Pivvsluuznn3s2rv7AMbFt27aYPHnyqPPPmzcv5syZE2vXro2IiN7e3pg2bVqcffbZccEFFwyZf+nSpbF9+/a44447+h879thjY9asWXHttddGRMQHPvCB2GeffeIb3/hGQ9sikwcAGjDWb/CLeDMaMbZvqbOQATFYFteJ5IyUnTPc8V3UcxqyrpmDBEBPT8+A6fnnnx8yz65du2Ljxo2xcOHC/sc6Ojpi4cKFsWHDhqqvu2HDhgHzR0QsWrSof/7e3t747ne/G6973eti0aJF8Yd/+Icxb968uP322+veBkEeAGhAow27IjcM6wn4ZKW7S1kahUNJz1gDtSMd41nqxgVF1w5d4NpJufBy0lNExLRp06Krq6t/uvzyy4es3zPPPBN79uyJKVOmDHh8ypQp0d3dXXWburu7R5x/69at8eyzz8YVV1wR7373u+Mf//Ef433ve1+8//3vjx/+8Id17T8fOwAFV9RuQVmhAVebWoMmWUq7r1znLKwPrTP4+GzGyEnVfnccAWTbk08+OaC71sSJExNZbm9vb0REvPe9740///M/j4iIWbNmxX333RfXXnttLFiwoObXcmsKUGDDNTTKNDjS047BoTym1Jcb/0Ue7pi9mvX+jnR8O44ARlcq7Z2SXmZExOTJk0etyXPwwQfHuHHjYsuWLQMe37JlS0ydOrXqc6ZOnTri/AcffHCMHz8+ZsyYMWCeN7zhDfFP//RP9WyK7loARVCtW0CeGtJ5psFWv1q6cWWpm0t5XbOyPryk0ZGusjY6FgDZN2HChJg9e3asX7++/7He3t5Yv359zJ8/v+pz5s+fP2D+iIi77767f/4JEybEnDlz4pe//OWAef71X/81DjvssLrWz8cZQIHUE3AQnEiWBmR1owVPspT5kJX1oDZjOXaG61Y1+PHhunbV2iWxlceS4aqBvOsolaIj4VSeepe3cuXKOOWUU+KYY46JuXPnxpo1a2L79u2xfPnyiIhYtmxZHHroof01fc4555xYsGBBrF69Ok488cS4+eab48EHH4zrrruu/zXPP//8WLp0aRx//PHxjne8I+666674zne+E/fcc09d6+aWE8YgS40OiNClJE21NHiGa/h5r/aqJdBTng8qDVfwuBZjma/ac+oZUr1IwV6fN0A7W7p0aTz99NNxySWXRHd3d8yaNSvuuuuu/uLKTzzxRHR0vNRx6rjjjoubbropLrroorjwwgvjiCOOiNtvvz2OOuqo/nne9773xbXXXhuXX355fPzjH4/Xv/718Xd/93fxtre9ra51K/X19fXVMuNOF3GICI0N0jPWejppBBja6TyppVirIE9tamkA22eMptYCyo2cl5UBjnY8v4cbWj4r6imirVYdWdZZoMBwLXp6eqKrqyvOvXVjTNx3v0SX/fxzz8aaP5kd27ZtG7UmT9a12WEDzeMbLPIij8Vui6bat/muH0PVcqzad4ym1uOjkePIMThU3s/NPK87FE1HJF88uEjFit3yQx2qjVTkpoAsq8yoydOxmqdvVseyblneniwYLdhTeVy7DpO2dizKPVptorSNJdPVtQQoijb7SALIr8qGRK3dCSr/n9TNaysaO+3U/YuX1FqrR+OMpLkm7VW07fe+QjakOYR6ERQpKwlabvCHvpsAklZPRk7ax2fay6cYHEekbfAQ69UKMTtO86cdaykB7UGQB2ow+Aav8nHIqqym0Y/GecVgtQQ3HTc0W72f+3m6zrY71wugyFziYBRuBMirtBocGjq0Si3dtxx/NIPPfoD0dEQpOhLuP9URxemvJZMHAMgNGT2kpZYAYjnrTLAx2yrfn+H+D5BXboVgFJUf+Hnt/gK0hmtCOuoZfQuaoXwsOaaKQ3AHskvh5cYI8kAd3ARA67XjcMSMTa2jb7l2Uy/HDAB55TYagMwZLVMjK8pBBg3C9NQSFJRxBQD50VHaOyW9zKLI+O0zpE8DjjQY2nWvPGxvHtax6OrJ/nJuAQBFpvAy1CDr2QS0j2pD+gKNB2nK59Zo55hzEADIMrcpUCPp/iSpshuQBiXUppnnSz1dwMrLBgAaVypF4kOoK7wMQMtVjuYi0AO1Set8yco5Wm+waXCX5GZtR9GDXr74ASCrMnJLAvnhZo6kaUxAffJSuLsVxrLNrdhP1V7T9QuAWhhCvTFtePsDtascfrcdGwukzxDQMHbDnTeu5+koykAGjh8AsszHFNRAQ5u0OOZotXa8vo22rZXBiNEa9L4EqE/ejzfvNQBZ56MKRuDmHSiyvDe4W6VyfwzeN5UF0e23scvr/ht8X5DHbQDIuo7S3inpZRaF5iuMovLb3LzelALQHJUF0Qc/5kuB+uQ1yJi39QWgvbgdgRq5qQNgJK36nEg7eFTrdo11PfMa7AGgNUov/iS9zKIQ5AEAyLC8BD8aHYpdtiwANE6QB4Bc0RAky0YbOrxd6rmMNeDj/AZATZ7GCPIAkBtpd1spEvuyNUYLULRjAKNa8eqR6L4FAGPXkfYKAEC9ih6gSHL70mpIHzhnRToLJnU7dgvgAECrCPIAkCtGMmqetAM8B85ZEQfOWRHnrduczoqQqtGCPeVRLbMs6+sHkEfl7lpJT0UhyAPhJg3yonL4apkAxXH9qmvSXgXGoFmfnaOdy1n8jC4HoLJyHSrvo8p/s7jfAGg9l3/a3uAbo4js3LQBIytq7Y4kticL++y0i8+KCEGePBr82dno8bRj98hBiSwFVCKytS5l1e5nAPKoVCpFqZTwEOoJL6+VfAzQ9qrdWAr4QPZVnqdZawA2oqiBq0rl4M7qxTNe/HdtmquTOVk/nqsFEZqxznkL9EQMvV9IYx3zts8AaC1BHoiRh3p1gwTZV7RztGjbM1g5uMPw2iHYV03eAj2V61v5b1bWMSvrAUBy1OSBQcq1PiprfwDZUuTuCEW/5iiyXJ9ybZUsHfOVn5Ot+KzMW42etGuEVVt22usE0AiFlxsjyAMjcIME2VR5bhbpPM1aY74Vrl91zYBRtQylPtRwx3RWj49WnIO1BHqyti8GB75q0aztqBZ4A6A9ZezjEQDaWzs10BRcHl75OBguAJClLkGtMlrXrYhidWur5T2ttj+KsO2Mrh3O+XY23LWuL9nVyIxSae+U9DKLQpAHgMQ1c0SeIina9oykXHyZkQ13TLTLsVJLoCcivw3gwcG8erYjj9tLfdSJLKasZSFSPA4x2pLRsyA9zRoVy7mbb4ovJy+vWS/1BHoqn5NXtWxHnreP2g03Aqz3P5sEb5qno1SKjoRTa5JeXis5FGl7eb3pJd8G3wi08/HnhrV9KLqcrsqskbydc6N1XxssyYDPcOtUz3KH2748vldFkZXP6cpAj2Oh9SrPuVqOAYEdsshhSVsa7puRyr8Dredcg+SVPwPTPv+Sari2urGe5MhiJKNawz3Nc8Zx0Xq1ZBkL6JAXDlXa0nCFC128ITluWtvP6sUzYvXitTJ62tBIn6/1FBOufLzRz+xWfOa3+rrmugnNN1xQR9sgPWkMaV6kIdQdtrSdkUYqiXADBc02+ObJOYZ6POlL6jwcnCU7XFeIsaxPFhtgo62P618+6D7XXga/p1m7rkC9HMIAADRdtS9Pmh3wrbdWT9rqXU8BhXTZ/+0hq9ePHbujfVvrKQyhHjJ5IJ98wwbJa8fzyre9tKM0C9U2sxtXloy1qHMzikFDkWTxuuB8pFUyeLhD62QxtZv2lLdvn6ld+T0V6CFvau22PNJ1KwvHfLtdX+vZTl920S6ydv47t0hSxg5/aL3h+t26+JIWx15xCfSQJ+UvQkb7XKyWNTOW43zwc5t9vqizUb96MoDcP5ElWTu/nReN6YhSdCTcfyrp5bVSxk4HSJ6LMGly/BVPu2URUAxjbbA3cg0bHFRqtcHrOjhYVdQuX81Qz+hoPtdotSydn453sihDpwhA+/FNaPHJ5oHhZaUb9UgBoNFkYf2zxOcarZKFc81xnYxSCoWXEy/03EIZOFUA2lMWblYAIjRcGpH0cPR5ea9k+NCItO+RHK/kmSYGQMrcSADtKu2GXJ7k/bOi8r3O+7bQGmleDxyT2dJR2jslvcyi8NFK28jbN2AUn2OxGIbrjqU2D4yusrvW4K5brpG1GWk/ZfX6433Op2a+b1k4Nh17FFUGTi9IlhsLoFkGD5deLZjsOgMjc460zlj3bZIN8HpG9CI9g9+nLARp6uWYol3k8PSE5lEQFWiGakM15+Xakqd1BZIx0pDpSfGlXDryGLwZiWMnnzpKpehIuBJy0strpYKdxgCQjMFZO3m8Mc7jOlNsWRlti6HSDPxUW047Nd5rLTnQrudO5X7JQnkGX56Qtja9FNCOqjXEXIDJEjcF+ZT3m2rHHFkzuMHmGE3HSN1Qy4Z7byrft1ZdI5v5unkJLuZhHVtttOtB2sGdtNehKAyh3hiXCtqOCy9Z5MatWPJwndF4Jg8co+mq/GyqJ5tmcF2yrH/GZX392l2ergN5WleKyyUNAEYxXEHlPDcM3IgCw2n2tS3P3VpJTuVxkrfPqDyuM8XlUgsAIxjum+wdu3UDBYqvHNAefL0by7UuCwWdSVctx03ePkcdw83XESkUXo7i9NdySAJkQN5uaNpJPQUdvY9AEQy+lg3OYmzlsjSY8yvNkSbT/rLF5z9Z4jIKADVyE0ee6U5AHgj6NEe953orRjBL8npT2YU66eCS62rzKbzcGJdNAGgCN3lkWRqNH2iGNIYNb8dzpB23uVGup2SVIA8ANMhN3ujcDKdHJgRF5rpCxNDrXKs/c1xXyTKHJwAAUDhJB5cbqQuTdk0Z6uM9aq2OF6ekl1kUgjwAwIia0VByQ5yePA9LDI0anHFRWbul/HuzFPn8yvr1Y/D7GpH9dYZWEeQBAPpVjiDWrHR0N9rZ4D2g3eRpJMTBQ9RnQfnanYdr+HD7bbRRMcmmUqkUpYQrISe9vFbK0GUEAEjKcCOplG/ofSMKkLxmXWdHC9jXUtC6PE+Wr/1ZCopBVjgtAKCNjHRDXL6pr7yhb8YNdJYbCABFMzgoP5ZrcF6u28OtpwyefCu9OCW9zKIQ5AGANjJSN6xqN8P1dNs6cM6KOO3is/p/X714xhjWEACaQxYq7UiQBwDa3Gg3wJVduIab98A5KyIi4vpV1/Q/dv2qiNMuPkuwB4BENbOuHORNkUYKA4DCKgdZhptqfW7ES12yBnfNGkk980bEgIwesq+W4wjIB5kre9kP+dVRKqUyFYUgDwBkVK1BnNGkUTxT9k5+VCuyDVAEAj20Ix/jAJARlY3rem5MR5q33I1qxwNrx7hWtfnti69frstz3rrNESHYkzcCPABkQXHyapLnoxwAMqhc/6aRRnc50BKxN/jy2xYHeiIGBnvIl+GON4VLASA/BHkAIKNqCfSM1PhevXhGXL+q+etVCzV58qHasVXtmDIcMQDkg5o8AJAR1YobV2tc11o0eXA2TWVmTyNqfR1dtfIvr/V6yuuZl/UF4CWlUjpTUQjyAEDG1BrsGU1lNk2zMmvOW7c5rl91zaiBHgGefChiZo4ADwDtzMcfAGRUta5a9dZHGVyH57x1m8ccgCkHeCIirl91TVy/aujrk0/DHVPlx/MQMMnDOgIwulKpFKWEU2uSXl4rlfr6+vpqmXFnAb/pAYC8GOvIW800XDFlgR7SNlyAZ7hzRY0hIA862yx43dPTE11dXfE39/489t1v/0SX/dyzv4//7/g3xLZt22Ly5MmJLrvZdNcCgBwYqQbPgXNWVO0+dd66zU0b5Wqk7llG0iKrJo0fOJUJ7gBkV0dKU1G0WWwQAPKtWuP0tIvP6u8+FbE3s6bZgZe9Xbz21vUpd9mqdOCcFXHaxWepxUNm1NLd0fDwABSNIA8A5Fw5AFMt+NL85USsXry3e1Y5u6eyTk+EQA/JG62mEAC0C0EeAMixwRk7p118VtOGSh9NOZhz/aqXsokEeACARii83JgidT0DgLYzuOjx3m5b1wz79yTWAfJAlg8ARSTIAwA599sH1sZpF5/V/28aVi+eIdhDLhl6HSBbSilNRSHIAwAFsnrxjP5Aj6ALAEB78d0FABREZS2etDJ6IE902QKgaGTyAEABVBY8VgAZAMircuHlpKeikMkDAAXxUmBHFg8AQDsS5AGAgpHFAwDkVUck3+WoSF2cirQtAAAAAG1LJg8AAACQCWnUyClSTR6ZPAAAAAB1uPrqq2P69OnR2dkZ8+bNi/vvv3/E+W+99dY48sgjo7OzM44++ui48847h533zDPPjFKpFGvWrKl7vQR5AAAAAGp0yy23xMqVK+PSSy+NTZs2xcyZM2PRokWxdevWqvPfd999cfLJJ8fpp58eP/3pT2PJkiWxZMmSePjhh4fMe9ttt8WPf/zjOOSQQ8a0boI8AAAAQCaUUprq8cUvfjE+/OEPx/Lly2PGjBlx7bXXxr777hvXX3991fm//OUvx7vf/e44//zz4w1veEOsWrUq3vKWt8TatWsHzPfUU0/F2WefHd/85jdjn332qXOt9hLkAQAAyLBJKqlCInp6egZMzz///JB5du3aFRs3boyFCxf2P9bR0RELFy6MDRs2VH3dDRs2DJg/ImLRokUD5u/t7Y0PfehDcf7558cb3/jGMW+DIA8AAEBGTBr/0lS2Y3d66wNJK5XSmSIipk2bFl1dXf3T5ZdfPmT9nnnmmdizZ09MmTJlwONTpkyJ7u7uqtvU3d096vyf+9znYvz48fHxj3+8of0nJgwAbea8dZtj9eIZaa8GAFXs2L03wCOwA8l78sknY/Lkyf2/T5w4MZHlbty4Mb785S/Hpk2bGh7pS5AHANrM9auuietXvfT7bx9YO/zMZNp56zZHRAjaQcEI8EA6Jk+ePCDIU83BBx8c48aNiy1btgx4fMuWLTF16tSqz5k6deqI8//oRz+KrVu3xqte9ar+v+/ZsyfOO++8WLNmTTz22GM1b4PuWgAAObV68YxYvXhGnLducxw4Z0Wct25zf+AHAPKoI0qpTLWaMGFCzJ49O9avX9//WG9vb6xfvz7mz59f9Tnz588fMH9ExN13390//4c+9KH42c9+Fg899FD/dMghh8T5558f3/ve9+rafzJ5AKDNnHbxWXH9qmv6/09xlN/XiJfeV1k+ANBcK1eujFNOOSWOOeaYmDt3bqxZsya2b98ey5cvj4iIZcuWxaGHHtpf0+ecc86JBQsWxOrVq+PEE0+Mm2++OR588MG47rrrIiLi5S9/ebz85S8fsIx99tknpk6dGq9//evrWjeZPADQZiob/QIA+XbgnBVx4JwVQx7f2yVvb8BHZg9FZtQpKJ40Cy/XaunSpXHllVfGJZdcErNmzYqHHnoo7rrrrv7iyk888UT8+te/7p//uOOOi5tuuimuu+66mDlzZnzrW9+K22+/PY466qhm7rqIiCj19fX11TLjTv1CAaAwzlu3Oa5fdY16PDlXfh9H430mzwaPMlX+Xd0aiq6zzYKYPT090dXVFbdseCT23W//RJf93LO/j6Xzj4ht27aNWpMn62TyAEAbWr14hoZ/AQyXiVXZDU+XPPJux+6XJpk7A1Ubbh3yrpTST1G4HAAA5NhwwbrViwXxKJ5q2TuGG8+OwVlXQPIEeQAA2kxlN6/KTB81moB6VcsiEuCB9OiuBQDQxsrBHgEeyIfhumel1WVLQIdmy0Ph5SyTyQMA0OYEeMizdggy1BLAycp+yMp6QLsS5AEAaGMKM0P+qE1EkZWiFB0JF0JWeBkAgNzam7lzlgweCqPow6qXRxYbvH1ZGXFsuPUDkqcmDwBAGxLgoSgqgxxZCHi0StYDKFlfP2gXBb4MAgAAFJvgCkWTRiHkIhVelskDAADkliAHwEtk8gAAALk0uA6MgA/kn0yexsjkAQBoI+et25z2KkDTDA7wTBpffWJ4B85ZkfYqAE0kyAMA0AbOW7c5DpyzIq5fdc2ARl35cciiegI1lfMMzugR6BmZ4C9ZUkrppygEeQAACqwyuDMSgR6yZHBgpxndsAR6gHYgyAMAUFDnrds8ILhz2sVnDZmncih1gR6yYsfugVMt8zdjnnZy4JwVznkoIPFsAICCGa7hNlI2T7UAEORJOYhTztgR1IF86ijtnZJeZlHI5AEAMkVtiMbU+8185f6+ftU1cd66zd4Dck1wZ3TOcSgumTwAQGZoeKSnnOVT/nf14rVprg6QkN8+4FwnW9IohFykwsuCPABAJgjwNG4s+3C0gsxA8axePEMgFwpKdy0AgAKoLLJ82sVnxW8fWDviN/Tlbl1q8UC2CYAD9RDkAQAyoXKUJ+q3evGM/sCOfQnFItBDOymV0pmKQpAHAMiMvV0IBCjGqtwQPG/d5pqHRx5uf2tUAkD+qMkDAFAgld2wag2YnXbxWWrzQMadt26zIDhtoRTJF0IuUCKPIA8AQBGUM2+MlAPFJMAD1EKQBwCgABppAO597lkDCjdrUEJ2OB9pJx2lvVPSyyyKUl9fX18tM+7c3epVAQAAACIiOtssJaOnpye6urrizo2Pxsv2m5zosrc/2xN/PPvw2LZtW0yenOyym03hZQCAjDpwzgoFkAGAmrVZbBAAIF+uX3VNXL9q7//V2wGg6Eov/iS9zKKQyQMAkFG/fWBtnHbxWWmvBgCQEzJ5AICaTBofsUONvlRVdt1SiBWAIiqV9k5JL7MoBHkAgBFNGj/w/wI96RHYAQBGorsWADCsSb4OSl15WHOgGBRUB1pJkAcAqEqAJ1vU5oF8O3DOijhwzoqIELyFkZRSmopCkAcAqMmO3bpqpUHxZdrZees2FzbrpajbBaSr1NfX11fLjDvd1AFAWyln8gjsAGkqQsHxcgbPYL99YG3Ca0KedLZZRm1PT090dXXF3Zsej5ftPznRZW//fU/80VsOi23btsXkyckuu9lk8gAAVcncAdI0OIsnrwGeiKHdLU+7+CwBHqAlBHkAAIBMW714Rq67N1UGqH77wNr+7cnzNgHZJMgDAABkVt4DPGWV2Tvl7clzdhK0isLLjWmzXn4AAEAeFC0AUi1Ydd66zYXbTiBdMnkAAIBMKwdC8prRM1zXLAEeqEIqT0MEeQAAAFokr4EpIJ8EeQAAgExph6yXom0PNEsppZ+iEOQBAAAyoR2yXlYvniHAA7SMIA8AAJC6coBn8L+V8hgcyeM6A/lldC0AACBVB85ZMeD30y4+K6U1aQ2BHqhDKaKUdO+p4vTWEuQBAADSU5SMHYAs0F0LAABIjYAOUMkI6o0R5AEAAAAoAEEeAAAgVUWrwQOQFkEeAAAgVbpsAf3012qIwssAAEDqfvvA2qqPlwszCwQBjE6QBwAAyLTrV10T16966ffhAkJA/pVe/El6mUWhuxYAAJBJ563bHNevuqbq4wAMJcgDAABkki5a0H5KpXSmohDkAQAAACgANXkAAIDMUn8HoHaCPAAAAEAmpDGieYF6a+muBQAAAFAEMnkAAACAbJDK0xCZPAAAAAAFIMgDAAAAUAC6awEAAACZUHrxJ+llFoVMHgAAAIACkMkDAAAAZEKptHdKeplFIZMHAAAAoABk8gAAAACZYAT1xsjkAQAAACgAQR4AAACAAtBdCwAAAMgG/bUaIpMHAAAAoABk8gAAAACZUHrxJ+llFoVMHgAAAIACEOQBAAAAKADdtQAAAIBMKJX2Tkkvsyhk8gAAAAAUgEweAAAAIBOMoN4YmTwAAAAABSCTBwAAAMgGqTwNkckDAAAAUACCPAAAAAAFoLsWAAAAkAmlF3+SXmZRyOQBAAAAKACZPAAAAEAmlEp7p6SXWRQyeQAAAAAKQCYPuTSp4sjdsTu99QAAAICsEOQhl8qBnUmOYAAAgMIovTglvcyi0EQmlwR3AAAAYCBNZXJJFy0AAIACksrTEIWXAQAAAOpw9dVXx/Tp06OzszPmzZsX999//4jz33rrrXHkkUdGZ2dnHH300XHnnXf2/+2FF16IT3ziE3H00UfHy172sjjkkENi2bJl8atf/aru9RLkAQAAADKhlNJPPW655ZZYuXJlXHrppbFp06aYOXNmLFq0KLZu3Vp1/vvuuy9OPvnkOP300+OnP/1pLFmyJJYsWRIPP/xwREQ899xzsWnTprj44otj06ZN8e1vfzt++ctfxuLFi+vff319fX21zLhT9xgAAABIRGebFVfp6emJrq6ueOCXv4799p+c6LKf/X1PzHn9K2Lbtm0xefLoy543b17MmTMn1q5dGxERvb29MW3atDj77LPjggsuGDL/0qVLY/v27XHHHXf0P3bsscfGrFmz4tprr626jAceeCDmzp0bjz/+eLzqVa+qeVtk8gAAAABtr6enZ8D0/PPPD5ln165dsXHjxli4cGH/Yx0dHbFw4cLYsGFD1dfdsGHDgPkjIhYtWjTs/BER27Zti1KpFAcccEBd2yDIAwAAAGRCqZTOFBExbdq06Orq6p8uv/zyIev3zDPPxJ49e2LKlCkDHp8yZUp0d3dX3abu7u665t+5c2d84hOfiJNPPrmmzKJKbZYABgAAADDUk08+OSCoMnHixMTX4YUXXoiTTjop+vr64itf+UrdzxfkAQAAADIhzRHUJ0+ePGrmzMEHHxzjxo2LLVu2DHh8y5YtMXXq1KrPmTp1ak3zlwM8jz/+eHz/+9+vO4snQnctAAAAgJpMmDAhZs+eHevXr+9/rLe3N9avXx/z58+v+pz58+cPmD8i4u677x4wfznA88gjj8T//b//N17+8pePaf1k8gAAAADUaOXKlXHKKafEMcccE3Pnzo01a9bE9u3bY/ny5RERsWzZsjj00EP7a/qcc845sWDBgli9enWceOKJcfPNN8eDDz4Y1113XUTsDfD8j//xP2LTpk1xxx13xJ49e/rr9Rx00EExYcKEmtdNkAcAAADIhjT7a9Vo6dKl8fTTT8cll1wS3d3dMWvWrLjrrrv6iys/8cQT0dHxUsep4447Lm666aa46KKL4sILL4wjjjgibr/99jjqqKMiIuKpp56KdevWRUTErFmzBizrBz/4QZxwwgm1b0pfX19fLTPu3F3zawIAAAAN6GyzlIyenp7o6uqKjY/8Ovbbv/5aNI149vc9MfuIV8S2bdvGVAcnS9rssAEAAACyqvTiT9LLLAqFlwEAAAAKQCYPTTVpfMQOXfsAAAAYi1JEKeM1ebJMJg8NmTR+71T+f+W/AAAAQHIEeRizymBOZQaPTB4AAABInpwLxqxaMEeABwAAgLHKwQjqmSaTBwAAAKAAZPIAAAAA2SCVpyEyeQAAAAAKQJAHAAAAoAB01wIAAAAyofTiT9LLLAqZPAAAAAAFIJMHAAAAyIRSae+U9DKLQiYPAAAAQAEI8gAAAAAUgO5aAAAAQCaUXpySXmZRyOQBAAAAKACZPAAAAEA2SOVpiEweAAAAgAKQyQMAAABkQunFn6SXWRQyeQAAAAAKQJAHAAAAoAB01wIAAAAyoRQRpYR7TxWns5ZMHgAAAIBCkMkDAAAAZIIR1BsjkwcAAACgAAR5AAAAAApAdy0AAAAgE0qlFAovF6i/lkyeFpokhAYAAAAkRBiiScoBnR27B/4++HEAAABgOEovN0KQp8lk7wAAAABpEJIYo8pgzo7de6dqAZ6RMngGvwYAAAC0MzV5GqMmT0oEeAAAAIBmEuQZg8EZO8N10RoueCPAAwAAADSb7lp1Gi6gU2uBZTV7AAAAoDpllxsjkydBgwM8sngAAACAZpFXUodGsnB00QIAAICRKbzcGJk8NaolwFNLDR4AAACAVhDkqUEzAzyyeAAAAIBWkGMyipECPDt27/27UbQAAACgcaUXf5JeZlHI5BlBIxk8lX8T4AEAAABaTSbPGNUauBHgAQAAgBoZQ70hMnlGUC1As2O3wA0AAACQPTJ5RlGuuwMAAAC0lkSexghf1EDmDgAAAJB1umsBAAAAFIBMHgAAACATSqW9U9LLLAqZPAAAAAAFIJMHAAAAyITSiz9JL7MoZPIAAAAAFIAgDwAAAEAB6K4FAAAAZEPpxSnpZRaETB4AAACAApDJAwAAAGSCRJ7GyOQBAAAAKACZPAAAAEAmlEp7p6SXWRQyeQAAAAAKQJAHAAAAoAB01wIAAAAyohQlpZfHTCYPAAAAQAHI5AEAAAAyQeHlxsjkAQAAACgAQR4AAACAAhDkAQAAACgAQR4AAICcmzR+7wS0N5cBAACAgpg0PmLH7rTXAsZO4eXGyOQBAADIOYEdIEImDwAAQCEI9FAEpRd/kl5mUcjkAQAAACgAQR4AAACAAtBdCwAAAMgEhZcbU3OQp1M4CAAAACCzhG4AAACATCi9OCW9zKJQkwcAAACgAAR5AAAAAApAdy0AAAAgG/TXaohMHgAAAIACkMkDAAAAZELpxZ+kl1kUMnkAAAAACkAmDwAAAJAJpdLeKellFoVMHgAAAIACEOQBAAAAKADdtQAAAIBMMIJ6Y2TyAAAAABSATB4AAAAgG6TyNEQmDwAAAEABCPIAAAAAFIAgDwAAAJAJpZR+6nX11VfH9OnTo7OzM+bNmxf333//iPPfeuutceSRR0ZnZ2ccffTRceeddw74e19fX1xyySXxile8IiZNmhQLFy6MRx55pO71EuQBAAAAqNEtt9wSK1eujEsvvTQ2bdoUM2fOjEWLFsXWrVurzn/ffffFySefHKeffnr89Kc/jSVLlsSSJUvi4Ycf7p/n85//fPzP//k/49prr42f/OQn8bKXvSwWLVoUO3furGvdSn19fX0NbR0AAABAA3p6eqKrqyu2/GZbTJ48OfFlT3l5V2zbVtuy582bF3PmzIm1a9dGRERvb29MmzYtzj777LjggguGzL906dLYvn173HHHHf2PHXvssTFr1qy49tpro6+vLw455JA477zz4i/+4i8iImLbtm0xZcqUuOGGG+IDH/hAzdsikwcAAACgBrt27YqNGzfGwoUL+x/r6OiIhQsXxoYNG6o+Z8OGDQPmj4hYtGhR//yPPvpodHd3D5inq6sr5s2bN+xrDscQ6gAAAEAm9PT0pLbMwcueOHFiTJw4ccBjzzzzTOzZsyemTJky4PEpU6bEL37xi6qv393dXXX+7u7u/r+XHxtunloJ8gAAAACpmjBhQkydOjWOOHxaKsvfb7/9Ytq0gcu+9NJL46/+6q9SWZ+xEuQBAAAAUtXZ2RmPPvpo7Nq1K5Xl9/X1Rak0cJStwVk8EREHH3xwjBs3LrZs2TLg8S1btsTUqVOrvvbUqVNHnL/875YtW+IVr3jFgHlmzZpV13YI8gAAAACp6+zsjM7OzrRXY0QTJkyI2bNnx/r162PJkiURsbfw8vr162PFihVVnzN//vxYv359nHvuuf2P3X333TF//vyIiDj88MNj6tSpsX79+v6gTk9PT/zkJz+Jj370o3WtnyAPAAAAQI1WrlwZp5xyShxzzDExd+7cWLNmTWzfvj2WL18eERHLli2LQw89NC6//PKIiDjnnHNiwYIFsXr16jjxxBPj5ptvjgcffDCuu+66iIgolUpx7rnnxmWXXRZHHHFEHH744XHxxRfHIYcc0h9IqpUgDwAAAECNli5dGk8//XRccskl0d3dHbNmzYq77rqrv3DyE088ER0dLw1mftxxx8VNN90UF110UVx44YVxxBFHxO233x5HHXVU/zx/+Zd/Gdu3b48zzjgjfve738Xb3va2uOuuu+rObCr19fX1NWczAQAAAEhLx+izAAAAAJB1gjwAAAAABSDIAwAAAFAAgjwAAAAABSDIAwAAAFAAgjwAAAAABSDIAwAAAFAAgjwAAAAABSDIAwAAAFAAgjwAAAAABSDIAwAAAFAAgjwAAAAABfD/A2rCbK2yg54BAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(16,8))\n", + "plt.imshow(flood_frequency, cmap=\"Blues\", interpolation=\"nearest\")\n", + "plt.title(\"Flood frequency\")\n", + "plt.axis(\"off\")\n", + "plt.colorbar()\n", + "plt.imshow(ref_water, cmap=\"cool\", interpolation=\"nearest\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example, we calculated the flood frequency and displayed the data by using stac. There was no need to download any of the files provided. \n", + "\n", + "We could save the data as zarr:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "data.to_zarr(\"./gfm.zarr\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "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.11.0" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}