{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Vector data cubes\n", "\n", "Exploration of possible implementation of vector data cubes in Python based on xarray and geopandas. The goal is to mimic what stars package in R can do - https://r-spatial.github.io/stars/articles/stars1.html#vector-data-cube-example\n", "\n", "![vector data cube illustration](https://raw.githubusercontent.com/r-spatial/stars/master/images/cube3.png)\n", "\n", "Because we need to use geometries as an index, we need to ensure that geopandas is using shapely 2.0 (beta) as a geometry engine. Shapely 1.8 geometries are not hashable, while shapely 2.0 are." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import geopandas\n", "import pandas\n", "import xarray\n", "import numpy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to recreate [this example from stars documentation](https://r-spatial.github.io/stars/articles/stars1.html#vector-data-cube-example).\n", "\n", "Load geometries using geopandas:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "nc = geopandas.read_file(\"https://github.com/r-spatial/sf/raw/main/inst/gpkg/nc.gpkg\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Get the GeometryArray. Note that this also contains CRS information at this point." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "origin = destination = nc.geometry.array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create dimensions and dummy data." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "mode = [\"car\", \"bike\", \"foot\"]\n", "day = pandas.date_range(\"2015-01-01\", periods=100)\n", "hours = range(24)\n", "data = numpy.random.randint(1, 100, size=(3, 100, 24, 100, 100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create `xarray.DataArray`." ] }, { "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.DataArray (mode: 3, day: 100, time: 24, origin: 100, destination: 100)>\n",
       "array([[[[[35, 56, 46, ...,  6, 79, 35],\n",
       "          [89,  5, 45, ..., 10, 23, 35],\n",
       "          [64, 19, 74, ..., 73, 85,  4],\n",
       "          ...,\n",
       "          [30, 23, 52, ...,  3, 54, 69],\n",
       "          [96, 63, 12, ..., 55, 93, 64],\n",
       "          [64,  6, 80, ..., 27, 32, 53]],\n",
       "\n",
       "         [[89, 11, 99, ..., 15, 38, 91],\n",
       "          [ 3, 97, 36, ..., 95, 98, 82],\n",
       "          [65, 71, 35, ...,  5, 80, 12],\n",
       "          ...,\n",
       "          [28, 33, 17, ..., 49, 19, 90],\n",
       "          [85, 37, 11, ..., 86, 91, 52],\n",
       "          [10, 32, 98, ...,  2, 22, 26]],\n",
       "\n",
       "         [[71, 20, 71, ..., 20, 94, 15],\n",
       "          [38, 67, 60, ..., 94, 85, 30],\n",
       "          [66, 79, 37, ..., 47, 80,  5],\n",
       "          ...,\n",
       "...\n",
       "          ...,\n",
       "          [18, 31, 36, ..., 89, 74, 25],\n",
       "          [33, 61, 48, ..., 66, 29, 65],\n",
       "          [94, 22, 18, ...,  8, 17,  3]],\n",
       "\n",
       "         [[10, 15, 91, ..., 54, 12, 28],\n",
       "          [37,  3, 21, ..., 11, 17, 53],\n",
       "          [48, 25, 50, ..., 77, 14, 16],\n",
       "          ...,\n",
       "          [35, 76, 83, ..., 15, 77, 61],\n",
       "          [ 1, 93,  4, ..., 71, 29,  6],\n",
       "          [40, 40, 24, ...,  7, 80, 56]],\n",
       "\n",
       "         [[38, 36, 12, ..., 22, 58, 60],\n",
       "          [71, 72, 88, ..., 22, 93, 56],\n",
       "          [68, 88, 61, ..., 46, 32, 54],\n",
       "          ...,\n",
       "          [57, 92, 77, ..., 96, 19,  3],\n",
       "          [77, 26, 76, ..., 45, 13, 72],\n",
       "          [ 4, 66, 93, ..., 64, 24, 65]]]]])\n",
       "Coordinates:\n",
       "  * mode         (mode) <U4 'car' 'bike' 'foot'\n",
       "  * day          (day) datetime64[ns] 2015-01-01 2015-01-02 ... 2015-04-10\n",
       "  * time         (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23\n",
       "  * origin       (origin) object MULTIPOLYGON (((-81.4727554321289 36.2343559...\n",
       "  * destination  (destination) object MULTIPOLYGON (((-81.4727554321289 36.23...
" ], "text/plain": [ "\n", "array([[[[[35, 56, 46, ..., 6, 79, 35],\n", " [89, 5, 45, ..., 10, 23, 35],\n", " [64, 19, 74, ..., 73, 85, 4],\n", " ...,\n", " [30, 23, 52, ..., 3, 54, 69],\n", " [96, 63, 12, ..., 55, 93, 64],\n", " [64, 6, 80, ..., 27, 32, 53]],\n", "\n", " [[89, 11, 99, ..., 15, 38, 91],\n", " [ 3, 97, 36, ..., 95, 98, 82],\n", " [65, 71, 35, ..., 5, 80, 12],\n", " ...,\n", " [28, 33, 17, ..., 49, 19, 90],\n", " [85, 37, 11, ..., 86, 91, 52],\n", " [10, 32, 98, ..., 2, 22, 26]],\n", "\n", " [[71, 20, 71, ..., 20, 94, 15],\n", " [38, 67, 60, ..., 94, 85, 30],\n", " [66, 79, 37, ..., 47, 80, 5],\n", " ...,\n", "...\n", " ...,\n", " [18, 31, 36, ..., 89, 74, 25],\n", " [33, 61, 48, ..., 66, 29, 65],\n", " [94, 22, 18, ..., 8, 17, 3]],\n", "\n", " [[10, 15, 91, ..., 54, 12, 28],\n", " [37, 3, 21, ..., 11, 17, 53],\n", " [48, 25, 50, ..., 77, 14, 16],\n", " ...,\n", " [35, 76, 83, ..., 15, 77, 61],\n", " [ 1, 93, 4, ..., 71, 29, 6],\n", " [40, 40, 24, ..., 7, 80, 56]],\n", "\n", " [[38, 36, 12, ..., 22, 58, 60],\n", " [71, 72, 88, ..., 22, 93, 56],\n", " [68, 88, 61, ..., 46, 32, 54],\n", " ...,\n", " [57, 92, 77, ..., 96, 19, 3],\n", " [77, 26, 76, ..., 45, 13, 72],\n", " [ 4, 66, 93, ..., 64, 24, 65]]]]])\n", "Coordinates:\n", " * mode (mode) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (destination: 100)>\n",
       "array([85, 66, 30, 37, 86,  3, 86, 54, 45, 89, 52, 61,  6, 38, 89, 66, 88,\n",
       "       63, 17, 95, 87, 99, 35, 67, 97, 55, 38, 18, 79, 59, 46, 21, 10, 42,\n",
       "       28, 35, 73, 10, 56, 23, 94, 87, 14, 66, 97, 69, 16, 95, 54, 54, 69,\n",
       "       10, 53, 78, 48, 10, 32, 54, 36, 40, 61, 35, 50, 66, 32, 31, 31, 13,\n",
       "       94, 73,  2, 64,  2, 19, 54, 30, 54, 12, 12, 70, 15, 39, 41, 81, 46,\n",
       "       62, 50, 24, 29, 59, 95, 34, 17, 37, 64, 40, 88, 47, 63, 27])\n",
       "Coordinates:\n",
       "    mode         <U4 'car'\n",
       "    day          datetime64[ns] 2015-01-01\n",
       "    time         int64 12\n",
       "    origin       object MULTIPOLYGON (((-81.4727554321289 36.23435592651367, ...\n",
       "  * destination  (destination) object MULTIPOLYGON (((-81.4727554321289 36.23...
" ], "text/plain": [ "\n", "array([85, 66, 30, 37, 86, 3, 86, 54, 45, 89, 52, 61, 6, 38, 89, 66, 88,\n", " 63, 17, 95, 87, 99, 35, 67, 97, 55, 38, 18, 79, 59, 46, 21, 10, 42,\n", " 28, 35, 73, 10, 56, 23, 94, 87, 14, 66, 97, 69, 16, 95, 54, 54, 69,\n", " 10, 53, 78, 48, 10, 32, 54, 36, 40, 61, 35, 50, 66, 32, 31, 31, 13,\n", " 94, 73, 2, 64, 2, 19, 54, 30, 54, 12, 12, 70, 15, 39, 41, 81, 46,\n", " 62, 50, 24, 29, 59, 95, 34, 17, 37, 64, 40, 88, 47, 63, 27])\n", "Coordinates:\n", " mode " ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAADDCAYAAAAiPnOsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACFHklEQVR4nO2dd5jc1L33P5Kmz+zO9l7d1t3GBVwwHdMxJQESWgKkk16Be0nIDWn3vukJaYQACZAQQgtgqk11771v72V2epPO+8faa6+3TV0bW5/nmZ0Z6ejozGhW+up3fkUSQgh0dHR0dHR0dMYI+UQPQEdHR0dHR+f0QhcfOjo6Ojo6OmOKLj50dHR0dHR0xhRdfOjo6Ojo6OiMKbr40NHR0dHR0RlTdPGho6Ojo6OjM6bo4kNHR0dHR0dnTNHFh46Ojo6Ojs6YYjjRAzgeTdNobm4mIyMDSZJO9HB0dHR0dHR0YkAIgcfjoaSkBFke2bZx0omP5uZmysvLT/QwdHR0dHR0dBKgoaGBsrKyEducdOIjIyMD6Bt8ZmbmCR6Njo6Ojo6OTiy43W7Ky8v7r+MjcdKJjyNTLZmZmbr40NHR0dHR+ZARi8uE7nCqo6Ojo6OjM6bo4kNHR0dHR0dnTNHFh46Ojo6Ojs6YctL5fOjonAhULURrYBNG2Y5Rth1+9L2WJf3fREdHRyeV6GdVHR3ggOcVVrX/dMh1imTCKNswyDaM0rHiZKBI6QmPpzVkRpFkDLKCQZL7XkvK4WV97xVJwSDLGI5Zd2SbvvZH2x5Zb9QMmAwKiiLr+W90dHQ+9OjiQ+e0RxMq23ueGHa9KsKoahhU14j97PJ8lNfbOlM8OpjqqWLvix4AJAlMBgWj0YBRkTEalf73fc/KwGeDgslkwGw0YDH3PZtNfQ+LyYDZZKQgV6W6ugWQkZCRJAWQkFCQJPnwcgUJ6fA6GUnqW8bhZRJy/3IZI4qcRVSY0IRAIFCFQBMCTWh9zxx5f8yDo23Ecd/B8XJLOmaJFpX69nO4HwH9fYrj9iWOW281GLCbTDiMJmyKEQUZWZKQZJAlqe/14WdZkkAauHw4opEoQX+YaERFMcgYDAqKQUExyMgJCkghNDQRJhDRCEYEYVUlomqEo9Gjr1WViKoSjqpHXx9+rwoBhz+/EH3f8JHX2XaVSWX7D+9HHFmDGHQkxHHvjnkvjm87RJtB22sIoaKh9j0LFUEUTah4duWy/R8gyRKyLCHJ8uHnge9lRT68TEacXUyXc+CQBKJ/aEc+0bGf/8h30r9egMUAi8t3Y1T8h5dqg4Z/7OfSog5WrJ5EONzXRNMO70nQ/z/Q9/Uf3ebIb+C/PnoBL63fzZ7mDqoKchhXmMP4olwKsxy8vGE3kiRx+dzJQ363H2Z08aFz2lPvfRtPpDHpfrJNvhSMZjBilxnoEx9CQCiiEoqoKet/2dIg8+2Pp6w/gP2+K/nZLufoDZOkJrOAHfs8KelrsaGULXtb4trmWCEiHRYmVS/X42nsGXk7RcZg7BMjinL42aj0v/70s01ElQY0EUYTITQRQRABYEvtzfzpPXvCn3MoFo5zstTyh5T2mSyiGBrqzmPra3Uxb1P21YW8NSGckv2/fDCDj8/YhM20L6b2k2dO4pGn51DXFt9+Lpo1kWsXTOf/Pf8O//fc2/3Lz58+nqvPnMrPX3iXy+bUnHIWT93hVOe0p867MiX9FJlf5qLCvJT0dYSqUCH79qfemnIsDkfqhMxRxuZEaZaMKetLU4e+cx9xGyGIalq/tSEYiRJ0+WPYl0Y4GCHgDeHt9dPb5aW7tZfyM0Lc9lgjPtYTUtuIaD2owt8vPABKs1Mjto7F7pBQhSPl/SaDJMMZv1jJx583sPSbVVTNLBp1m65/bEdJ0UX6YK+Hn6+ZR6v7opjaB9jLbTe8ynlztJjaO20WzppUjhACo0HhO9efz8fPmd2/vrowh7OnVuHyB9jf0pXIRzip0cWHzmmPWUndHfoE+/qU9WXUDFi3ZKWsv+Gw26Jp6HVsxIdBSt0pTI3GdtEYjYg/sTvvipkWzv/+a6gZG0Zsl23vTqj/4agstrPWdIBf7V/IC603sMNzA92R84hqqRXSiSBJYJq0j5I7VnDu39cwYX7piO2zzqseNFGUDL5IhKd2ZlFunYJJHl2chUUPCxY9w13X9E3XDIfFaOCui8/kt5+5lotmTexfbreY+1/XlOZjMhhYOnsSb26NzfryYUKfdtE57YnlpBJ7X4e4sHAebybp+2HSDEzYMJ69abZ6AFitkdEbxYk2RuJDFqkTH5EUiA+DIoOW2OXvgs8pMbWzGNuAqQntYyicmUa6Dl8o93m66ZvFMgCzKbc5megwkmv2E1DNeCIGDLIgyxiiyLwdi5z+i6IWMhLeVYNk0MguN8G6odtln1nGxvPsaJHUTLscodQhcIU24pAzMZum0RLczkjiWqCSW/YC9961mJ8/XkYgNLhtMBLl/z3/Ds+u3s63rjuPhTWVAJw7bRwOiwmXL8Ds6hIArpg3hR8+/RafvXRhSj/XiUYXHzqnNaoW5qDn1ZT2OdXxNBblOl5qHnnefyQMKBw4mNo73OGwWFJ7sgZAjNH8dAr3E4kmP/1kUmITEEORP6WLWOSPSiN9d9Wp+ewW6/ACrsHfS0P/LNLA6aQZzrksLUi/+Gj44zms+N3Bw+8ODduu97OT8QQS/58bjgxT33NEcxMJr6fYPBGfJnBHmkbcTrW8z7c/PZ7HnpnPweah2xxs6+azD/2b82eM5xvLzmFGZREzKgdOL82uKsEXClPX0UNlfnYqPtJJgT7tonNas6f3WXzR1pT2KUkw3vZvLhzF/8MkGSmz5g0ZJOCXQ4wfl5PScQ2H2RRKeZ9ijCwfWmpmSgAIRZKffjIZEhcfUozfmSYCjM9P3eSCwZRYX9t6O9nhuYGwVpKysQxF94HRRaHiMLEv6ErL/t9tBJM8of+9N7IPoR6g1DKBXNM4RppeCYgDfPz6V7h4/sg/1BXbDnDtjx/j2dXbB62TZYkr5k7mra37E/4MJyO6+NA5bQmrPrb2PJq2/gvM7gHvs4x2ptjHUWYqxSnloUat7O1xU2AoYmpG5aDtbcVj8+9pNAVT3udYiQ81ASfR4QiFUyA+krB8iDisONNLUuckLBkT/w6Xt3bzu4PTWev6GH41dVNBx9LTNPrv05xrQ06h/8+xCCSWHzpz0FJ3eBvh6A6KzCUUWaYwnAiJiF7mLXiOSSNXmGd2dQnqYTV9qK17QFju5XMn88YWXXzo6JwSbGt5gdAouTuSIce4gUJLX2npalsRwbCRDV2N7HV30BJw44n0nVTrfN2s72hhnK24f9vpGdWESzSQUuk+NzSKIfXiY6yIRFJn+giEk/d9McpJnFLjONTVeakL61bl5ESXKjTe7ejgz4fKUzSio4R3T6J+++ixq/46F4vXRJjuyMUkJy4Ah+PZfQoGadaQ63yRg/jCGymzTsEgWYZso4kw+dkjH+ANBxrJy+wLof7He1sIHzMNOK4oFyEErT2pj3Q6UejiQ+e05a2HrWz/9c1EauekpX9F6uSKohamZ4xjT4+L7tDIIZi9IY1SSy55pkw2drawOnyICTflYLKm/mR6LLIyemhovMRzF58MgWiKInUEBFJg+UhGfAgt9m2LMt2jN4qRsJQah2NDMsJrGA48Gbugaf7rJsK3v8Ylj/fwHes0vjH3bC6sGE+WeWhBEB+CA52TGMnPpje0Cbvko8w6Y8j155/3CjddHGY4lalqgm/+9SXe3nGQ/6zfhS800Bfr8rmTeWvbqWP90B1OdU5b5i2cwE/u38aW1QXUzPw4M5bVYap5n0TSBES0MvZ6FtPglSmwQr7Zh93Qy8buUtZ3xZbArME32FluY6ierKutTHEVs+vVjvgHFgOSnHrxMVahtv5QaqYfbAYDxOTuOTIGWSLhS3kcgi3D2glUJLqnAfi11Fi+FuWlNqmc1p3DhmdjTzB2hD0f7GPPB/vILnRyxacu4sd3fZxeG6xpaeDXm1bR6vPG1I+E4JrqCNdVNTLOvh5JNBEyLGG7b9fwYyZMb2gdJeaJRLDQETrqkBsRbipr/s39E87gmeU1SBIU5mqs2ynjPmzICkdVvvSn5wHwBsLkOGz92186p4Z7//YKHz/njLi/k5MRXXzonBbUHmintdmFzW7CZjdjtZmZPbeKm+88h8f/9DZ7tnazZ2sGVZNuYuHdq5Bz4zvp1fvO5MUGV9+bfg1hAZJPDuSKBFhrP8jUaWU07OhNur/jEZIvLpN/LIxVqO3xd4eJYpWNhEje8TaZvCPDZCcfej+G+DKxDkd2hommYGr6muzYlJJ+jtD63CyikdqEt+9p6+VvP3iGJ3/8HGdfdybLPncJL193O99+51VerxvegmAzCu6b7+P8nBdAHJ7yOXxszNH3KTbPoSU00rgEnsheTHIORtlORBs4RRZSNnHlFUe/q8Vz57Fq4zheXzvwt/Pyht189tIF/e/zMu0YFYVur3+AKPmwoosPnVOa9rZeHv39St54eUtMJ/favS56vz+XS+41oRTGHkZYYKkHchMf6CgIwDbdADtS37dGbHeCcZF+VxUAPMHUROrYZENqxEcyoisOy0dEa8FiEASjyYm8SeMd7EjRwfKrhVjkPQBowkBIq8EVrcQghck3vRFXX1pYYd3jqQk1V6Mqb/9zFW//cxXjZlZyy3eu4Zyzq/jBqpWE1Cg2A1w3KYpJ0bAaNKbkraLYnHdUeAwcGaVyA92yjZA2ssUwrHVjk8zk2OfTHNhPSBv6xiFkWM+iBY1owTt4e0ctUbXPAvfYig3ctGQWWXZrf9tL59SwctsBrls49NTOhwldfOickvi9AZ545H2e/ccaIuH4TPM9nX5eemAal99rxlA2OPRtKCxKesUHgKyl3ppgUASaCKS837GYdrEqRjwpEjlmOTWnwmS8c4SINdgWQGNqsWBjQ3Lfs8gMHp++I2EerTUzxXkTqiY46HMR1lSgzzJ1S+UlFJpiz6fjfXs+Pa1xFkmJgYNb6/jRrb/mF+98nxeuvYXX6h6lOus1wtrRKc2wBmY5C4Y5bQitjenWeWzwHWK033llxvVMyfkqIc3NqrYHaAsMnSEtqLXy9Rtz+W+xlJc37Oah5avwBEI88uY6vnr1Of3tLpw5gfufeO2UEB+6w6nOKYcQEUzBz9FwcFXcwuMInt4gLzwwDv+WC1C7h59b39xzE2s7P86bLRcmOtyYkdQ4LzQSFORnICsS+XkZTJ9WQrbTOqBJTlZ6TBRjMe3iMJpHbxQjJilF4iOZrzMOy4dQF7ApSeGR6zRzIDhM9qsE0BDs6O1kt6frsPA4ypP1gpA2YZgtB7P1sfQ5WWuqxk8/+VsqzHaqs/4zQHgcwSJGySwcXc8k2+ihxa3+lUQ1LzZDPheU/IrZuXcfrgY9mH29z5DtsHLRrAkEQn2eQ0++u5n23qOWSZvZhNVsxO3/8EaoHUG3fOiccgjvz5GjH3DFZfmsej8j4X4C/jDP/tQETOXMi5Yw/toVyFl9J+uwWsz7HeeztrNrrGYYkIbQUZMmFGCxmwhrKiZJQWiCQ4c6yHLawK5wsK0byQFNYS9NDV5mVhZRHMlC0wThcBS7sY2Jtpn9Zew1jpT/7nvVVw788Gu0w+3E4dLnfc99r/vWafS9VyQJSUiINIYK2xUTkBqrjSlFOSKUJESXELHZi3qi81jZsIS50yV6XFFaOgP4g/GL7AmzDOzwp6Ouz2BUodEQmMV42/5RHbpDuyazf33qRNFQmMxG2nu3Eh5iKqTAMgWbtnLUPjLUdeSaJtEVHn6sgWgza9ruZlHxwxhkG1Ozb6XQOpf32/4b73GVtJv9q/BGmqnrUCnPy+JQezehiMqfX1/LvR+5oL/dJbMn8e7OQ1wxb0rsH/gkRBcfOqcUIrwFfH8B4I03M0mV88HaN7rY/O4cllx7NtIFLl5ry6TJP7aVJo3CwLSpJchGmWA0SpfXz/b2wXdtZouBXr8b1dv32Y/1ddlaNzCb6+VnmMlR1xBzpId03PMw7AophNwmZEnCJCuYFAWjrGBUZAyyjFFWMMgSiiSjyBIGSUaWj5anlw+XqJcQSJKEOHwcBQKkPlGUETFj8QUwKDIGg4yiyCiKhCTLSDJ9JchlEJJAk0CTBCoCk8OPzeknqkFUhYgKhUoUe3RghVzpmBfHflxX0EI0KqNpAk3T0FSBpoGqaRiScfmIYVqtPbKEP9cXEVD3gQ2wgbkE8g1WcoyZ2LCiREwEfRJrt3Qz3IGaPSWbHf6xLVb2fHMvFbbrubZ4BQZ5eH+O1o1nMuNGCS0YJuILEuz1oUU1ZIOErChIiowsy0iKjBCClq21RMNRFKOCbFBQDAqyQT783PdaVhRkRUZSFCrH5/PNP36G+ujfh9SupVJjjKeNIFXSNqrss+jSbNQG9jDU990b2sHG9u8wv/DnSJJCrmUql5U/xvqO/8chz0vHtBTs732WMyd+gWfvuY01e+v5zEP/5pkPtnHbeXMoy8sCYEFNJT9+ZoUuPnR0ThaECCPc99J3/y6xcYNGKn0PwqEobz7lpsJUQVNpQ8r6jXn/PbC5cfTIhHjShL+8SeXWJRdSkvF6MkMbhC/aNyWiCUFQjRJUU3+HPcNcQF1b/LU8rrlUIVL4zqDlZaNkoDzCvgNLaAsMneSrKpjEKXUY8aEKCy913oJfFez1uIiIwcG87mgAd/SYK6kZFp9bRf1uCIU1fIEovmDfMbBZFHqyW4+4Y4wpnmgERRooPLRAFprfiQg78NSO561Xjz+m5qNXKpVjfDEOK4QJVXGN4dN3X4Etw0p706oh12uSOY57FgHRzeQCufa5bPG3ERWDp0Ta/CvZ0fW/TM/7DgBG2c7Cwvsptp3Fuo6f9EfEHHC/wIzcT6FIJhbUVDK7upjNh1r4/aur+cHNlwJ9xQsLD/lpa+ymsGxsSjCkA93nQ+eUQXgfgmjf3ZyE4PNfSM+8aOANd1qcP0clTTMYgUgqkjANxK8aR2+UJHKiwjKJqSBJcwwrPACSSRY6XGK2Hf7rWdfTyQ53FxER+/TKzkAt3spaIhPrMc1sZtY5gnMWZTF/bhZd4ROTKdMoyXRHL0TTTP3Ltj++jMe/OpO/fXsczz+U/knM6nH5RDQfXcGtQ64Pk2DxtugGFuR/HbuxasjVh9xPcLD38QHLqjIu4bLyv5Fn6XMgDWkumnzv9q+/8ezZAPxn/S72txz1QxlfU8KTv34tsXGeJOjiQ+eUQER2ge8PA5adt/AZnn92NX/8YzMVg0unJExWVQbaGKQ9P5548kDEg9mQ+ltgXzT9RtVExUcyJz1ZjPxDkqKJHyR3ixmhlaOJo99dQCvlxZbUHJ/6YAe+SJiNm71M9kyjuHkSxc2TmEhVSvqPhfaQl7/WSvy57mL2+j6CJgw4clJXp2Y0bDYTBUVOOgLrEAytFEMk7shs0JpYUvp3Cm3nDbl+R9f/0eZ/d8Ayh7GEi0p/z/TsO8i3zKbYdjS3x8WzJpKbYUMI+N0rRy01M84cx+v/WkfTofQkHhwL9GkXnQ89QkQRvffBECcTs7KfiqL9PPSrYr7yzYvYtze5K3jV9ELeXdI8Vgk8B5LCCq7HYlJSk177WHzR9KaEB1ASSUULSR27UCQXRsiLIkUS/3395bN+8p6+jFqtkyJLBiVWI+6IRlBLTb6LDOwc3C1wB0Os2380jHWiyMFSaiI4hvMwnmiIF5tDTM64npIaM9XTCzi0vT3t+60aV4AkSXQGNgzbJpjEKSIcXk+G/DXmF/6CRu9L9IZ24AnvxxM5QEjtIts8m2zz4BoxsmRgZu5nEKLPx+kIRoPC9Qtn8MfX1vDm1v1sr29lekUR+SXZ5BU5+dsvlvPtX96a+IBPILr40Dkp2bK1nn89vQ5FkSkpyWLRoolMnz7MpLzvLxAdOR+HIrWwZEmEfXsT/8mXTchly2XdROQ0qYBR0NJk+jAaUpOo61h80fQbVaVE68ckYbVyBS2kS3xU3TKeLWrfBbgp0EtTCtOvCAHOjirq/IPvlPe1d2PtyWFyaQ5miwRGlaghhDFiwe+WiKoaiiJhzAmS025Aer0H9UwTGya0oMrJ/SZ3e7rYbelCWSYzyVVAW2PqM/gei8nc9//vNE8ato1f9SdsHguHNyKEhiTJlGdcRXnGVf3rQmoPBsmKIg8/zSkNIag/smgmD7+xFlUT/OalD/j9564DYMZZ43nr2Q3c+PmLqKopHrTdyY4uPnROKEIIdu9uYc3aAyiyjNFkYNWqfWzbNjAM7fkXNvGbX99KdXX+0W3VdoTv9+B/KqZ9FRdHOP4nr5gV7Hk2rLlWjBkmtLBK1B8l4g0T8oYJ9obIzrXjuMTJ+8UNiBM4USlSrHmcNo1bz/HiNB9KbceAL/XGlEEkPO2ShOWj1TfKQUiwyq4AvMvMwya2SpYp0Sl80DC8iT4QibKp9njLw3F+IfWwdLvCgQ/q4X2YVJqF+pkM9uYnH/Wlyho5MzLSLj42b6iloa6LkrLz2IQRbYhKPIFoN5iG2DgGhHATje7DaKwZtM6sJOZLUpjl4MKZE3ht8z5W7alj/f5G5k0oY/qZ43nz3+t5/Gev8N9/uCOxAZ9AdPGhc0IIhSI89/xGlr+ylbr60U9egUCY+/7rX3z9a5fhdgeYUv0+BZZfQhwpsRfNf4+b7ruRFw9m4I1EcAVD+MJHTM2BvocdjvqbGQADWk4mWqmMAYVIuq4OMZBK8fG9j3Zy/qSngfQ45frGIH1E4g6niW0mBNR6RnbUVHMSc7Stun0CW9T0zN8XUsCGna6k+ykXFg6uOhqe29XkwvADD4vuHEdttZsWuycxcS6gKOgknB/FYJCJRtNrWXz2n2v50jcvo9C2iBb/24PWhzQXLnkJGeo6FGn4/w9JysBgGI/BMO6Y53EYDCl0MDvMFy5bRHO3h+31razf33BYfIwD4INXt7FvWwMTZ8ReAfhkIC7x8dBDD/HQQw9RW1sLwLRp07j//vu57LLL+tvs2rWLb3/727z99ttomsa0adP45z//SUVFaiow6nz42b2/hcee+IB9u1vpbPPEfC1obe3lm9/qs3J8+na4YWl80wWK1MYnzvoVE8uX8oMVM48RHiPT2O2msRsU2UFlTia52VY6nF00iLHN8yFSOO1y7qQ3SZfwAPAlaAGIh4SnXRIMGzJQjD86skmnI09FEJ++EYD7SlNarB4GoRCocxKKJm9RmNYgse+432A0rLL/oT5BUmxSKBqfS91dKl220eeMZE1i4YYJNHzQidcd5ACto26TCl5/ZSuf+Mx5lDkuGVJ8AOz278UglZBlKsMqm7BKKrmmMjItsw4LjfHIcu6Q0yTpoKowh8e/chMvrtuJ93A9o9LqfLLzMujp9PDY/73M/zz6mTEZS6qIS3yUlZXx4x//mAkT+tLkPvrooyxbtoxNmzYxbdo0Dhw4wNlnn82dd97JAw88gNPpZNeuXVgsqQ/l0zl50DQNNaphNI3+czrU2MXdP/gX/kDfhX/qmWUcXBtbyflj2bbLyQ1L494MgMUlr/Hsje/y0zWf45+bYp9YVzVBfWcv9Z29GBWZeZMm4bZ72CtSUxW0H9FXHVUVfblFj5zfUik+gmoBNiV9WSR90fRbiBI97Sd6wdCipYzm9bsz1MVZv12Ixa1RmBFFyfGhmVTkkIz7pQD7nqwfpH3yzypgf5qsHuN8U1jdlbwjpxMD9W8fHLFNNKzSuKudHF/ZqOJDFjJnNkxn31sHiSZYAiFRgsEILz+3ketvOQdZMqOJEGYlh5A60LE3Knx0hvYcXeDfSZapgYlZt1FuzhvTMQPIssSys6b1v5ckielnjePdl7aw/u3dbF93kOnzx435uBJFEkme0XJycvjf//1f7rzzTm666SaMRiOPP/746BsOg9vtxul00tvbS2ZmZjJD00kzHc09vPjYeyx/ajWSJPGZ717LBdfMHba92xvkrvv+TmOrq39ZZXE2bdviPzlmZJj51+/+gJKg86eqydz5n2+yvt6d0PZHkCWJsrkm6kTyFw+TrDHFNo69HW66QwNP3rLUlw3UoiiYZANmRcEkHc4aKikokoQRmayACcu2wydz6fCF9vispBKU5gSoKmwmFDZy/dLfJj3241nw2tfwjWIlSJYlxgq2bIlfQN14jcCb9UHc2wV8F/LiodiTml1X46AutGvAsiwlk8LGfNqe6qHhjb6Ik+JLy9j1idQ7/Y4TVazfmJqqcRd3Z3DwX7GVVK65bgIfXNCCOkyhG1nITN4/me172qkpysN4yE/r7rENGc3Lz+DxZ+7GE91NZ3AT/mgzEjL7e/8+7DYWJZ+a7DuoyrgWRU5dXaFk8Lh8eHoDbHh7N+++tJmfPPWFMbPGDEU81++EfT5UVeXpp5/G5/OxcOFCNE3jpZde4lvf+haXXHIJmzZtorq6mnvuuYdrrrkm0d3onGQIIdi9sY7n/vI27y3fiqYevfj/+QfPc+5VZ6Aogyd+o6rG/b/8zwDhAdDR4yUn34Ez1040HKXhYGxTGR5PiBffuZlrzktM6P5526eSFh7QF4HSvCnCvEmT2GtvwH1MhVgDgjuqI0zP7CGkGQiqCiFNocjso9Kyg5DIZZ+vEqcxQImpDqPYzZ8P3sHq0GBLoSb66qhEjqkUOhRzIwUc3DK6mOu7jPR51YUjX+YjS3+P0RjbBbBTPpMI5v4beAmwCA8O6jGKbjRB2oUHkHDSNZMlsbF1BuIPHy6UZyAbu2gJ9Ykkl+rGVeyGr8KErxfj3OPE1yARj+9SLNiFlYO7U2MpMwiJ7nfqY26/59/7yX9eJrfESWa5g+hZZjaMb0ZTBgoPgD2tnWCF6VdUY3Op1K9vIhpKvyWks8PDO2/t4oKl08m29FkThOira3Sg94kBbU9G0XGEjCw7GVl2im9dzDv/2cTm9/dyxtmDnV1PRuIWH9u2bWPhwoUEg0EcDgfPPvssU6dOpbW1Fa/Xy49//GN+8IMf8JOf/ITly5dz3XXXsWLFCs4999wh+wuFQoRCR//x3O7kLwg66WH169t56jevs2fL0Ceiy29eNKTwAPj9k++ydmvdoOX+YAQ/0NEbxWRQKK3Kpam2i5o5ZSgGhb2bGolGhj4Z/fVJO8vOZdhCVZoG8hDD2dq5kN+9r5CqlKFRVWPDrlaKs3IwTnBTmOHhM1XtTLS8hxDDCAENrOxjpm113/vDQ7mpeitPHDgz4bHYogrxeqI8+W+Vhuav8a1P/Sim9vtDQbzRA0OskTArE8g0xJinPEkSvb9TjInls2jyxrfds3s8CCQurMgFebCFpk3rpG1iJ0yEhdZSZJHJdlc3vZHkhUixZwJr/akpSb8okEFT++D/3ZHQVI2Ohh46GnrgAzj7IxPYvLibsoZx/cLjWLY39y2zn5HJdKuTAytqUzH0EXnmqTWcf/G0fkuBJEnMzP0GlRlXo4oAUS2IEBHyrWeedKLjeCRJ4vPfv57f/Ne/mL140gm1fsRK3OKjpqaGzZs343K5eOaZZ7j99tt5++23ycrKAmDZsmV89atfBWD27Nl88MEH/P73vx9WfPzoRz/igQceSPwT6IwJzXWdfP/TfxnR7yDgH3zSdHsD/PLRlbzyzs4R+49EVCIRlTabxJQFlWyvayMSUSnJstLdMXRehZxs85DCQ1UVfvWbu3jvLS9XXGfmE7f+pV+EeCNZfPvVxaja8CmyE6XF5eGx6dvIsbwKiIQykjqkD5idO5fNXYkl6TKHE9tu/aYAOw5cxLTxbyCEhEeezOZgCAkJSVKQJQMSMrKk4IsO7+MSUntwCQmYkdA44iJR7WgIxb+tMNPgje/GSByWR7EkPW0MNAFNZFtlZuZUEFItbOnuJKTFHzZUSjHr9qVGeCAEYk3yPiO7/7UfyzMSkU+M/Hl8oTCHzH4mXDyO/a+P7GOSLGUVg+uiSJJElnlky4EvGiKsRsk220dsJ4RAaB1IshNJSr94qZ5cwuTZlax5cwcLLpqe9v0lS9ziw2Qy9Tuczps3j3Xr1vHLX/6SX//61xgMBqZOnTqg/ZQpU3jvvfeG7e+ee+7ha1/7Wv97t9tNefmHK2TodKCtsXtUh8fN7+3l3Zc3UzGhEFenl5UvbKT2YDublNhPoF5/iE37j94lWmzDB9yrmiAUtmI2DfSPePLpW3j9P33OY//6W5iN6z7NpMlQVe1mq7Gchp7hk0QlhyDXuiZJx1ATrYHEk4kYwond8QRDUe754WSuvGYmNUv/TURL/IIjYU1427hI8GtWZX/ckSWKqEhwdwKnvQdvjIFFqtCo9dcCMN5pxxPKocEfRx0WIRFoyEQQf8G9oZijZtK8LzZfj1ERAvF6LZkL8nEHhrfutLt9VJWk199vyrRSvn7PVUNaCHzREC0BF82BHloCPTQHXLT4e2g+/Lo34sdptPHows9Sbs8dsK3QvERDHxANrSAaehuh9jnTS3IeklKMrJQcfi7HZLsJSc5I6ef6+Jcv4cdffIwzL5iKPJTZ9yQi6TwfQghCoRAmk4n58+ezZ8+eAev37t1LZeXwcc9msxmz+eQ2aZ1O/POhN9m3rQGbw8I1nzyH6ikl+L1BHvnxi6Nue3BXMz/8/KODls+5aAqrR8mPMBzOAgfNdUOnl25q8vJ/j3yGez/9i34LyD/dN7IpNxc46mh3cE8PBw//LA2GBs65fDzvSKlJWX0s1bkCIZI76bvURbT6kzCZBpILcf3Pc2EmL5kE5uHTT4+GxBhFtyUoPqJD1VEfhUikgCHrr4/CZeOyaQpui3s7AE/UR47FSCBqozMcm+PoZLWGDzpS99vO3e4llWm/3I0uJrsLWDtKOhRJG/rgCvqiPsQw62MhI9PKrXeew8Zdh6jL7hxSXIxGb8TPVzc8zl8WfBq7VEs0uJJo6G3U8DqGKvMgtE6E1okWOfpbCPv+gtX5EwyWcxL+LMdjz7Bw4XXzeO/lLZxz5Rkp6zcdxCU+7r33Xi677DLKy8vxeDw89dRTrFy5kuXLlwPwzW9+kxtvvJFzzjmH888/n+XLl/Piiy+ycuXKdIxdJ8Xs2VzHIz/5T//7N/61los/eibtzT3s2xZ/OGw/Scw/bjrUwqS5JRgCGgd3Ds4DsOIdN5Vln+KWK/4EQKWxnUfLNDJKLXiaBt9uRqMq0V0uzpyezWbRSziFmbvKspPPrPV8Q3LOYmqS4gPA78rBUJjEGCJW7JoBnxRJ6tiPSoLXn5CI3/LlDtmIRXyUOSxIQJs/xKXjnLSpu0bdZiS6Iy6KHBrTTaVs6eqhJzK8CcWOje17UjedWK1ZObh2z+gN46T+tb2c+cl5rG0f3rq2tbuLkkvKsRgUjJKMAQlFSHREg9R39TIvt4D9rw7ldzQ6HneA+775JJXfKWRjS22CnwIO+TrwdH0CxOqEthdqE/7uWzDabsOa9YMR2/rcAQ7tacHV6WHRJTNGtGqce9UZ/O7+Z1h86UwUQ/prLCVKXOKjra2NW2+9lZaWFpxOJzNnzmT58uVcfPHFAFx77bX8/ve/50c/+hFf+tKXqKmp4ZlnnuHss89Oy+B1UovDaRvwXtMEr/5jTdL9Cn8EkywT1hK7MO6t6yA70zps8qa/PqHS6PoqxqW9tIQkAloHk75Qxp57hxZMjfvaYV87Z8+vYEWJr39uPlk2NRjoKwqRoACQivnbvuSsgOEU5DXXkiwK17WzgPI/uJAkMFtNWGwmzFYjRosBg8WAYlGQzDKSSUaYZDQTqEaJqEEQMQhCikZQ0QgqKn4pileO4iFC9Lg6IonNbgmkBI53e4wRq3PLQjQG6ykAmlKU5bU34qY34qYyMweD10ZHaOjBFPeOpyGYIl8PoKZOsG/0ZgnR8vQ2pPMKhz2GgXCEA+3DW3AiiSWSBUBSwP49GxvDtYl3ApRZJDJE8udHLXpU4KmqRktdJ4d2N3NoVzOHdrVwaHczbY1Hv4uFF0/n6//v49gzh57elCSJq25fwjsvbeb8ZcOnPjjRxCU+Hn744VHb3HHHHdxxx4cvz7wO5JdkI0lSSpNZAex+fx9FuQ4K5leyvicxI643EB6x1tPq9QHcszv73+9xNDBaPMT+dfVceN543sh0peQO3RuSQSoDEXtY4rHs9S1CTargOwTijMg4litvjDJ10Q5Cyp6k4oCiwb5pFyEg6A8T9Cc2Jom+bPd2oBAwGJV+IWOyGjHPTERoSYTq50HZOzFvIQTUuke3KORbTTQF0+ck2RXupsSRzXinjVDUwabuo5aDMooHVKpNlhxhpHaUpGLJkFGSmaB47MPg7nPakWSJ8jklKEVWmgM+ijwydWub+ttlFtjJq8mjdVs7flef5Wr8zcWsDSf/2RqDGh5pLhlifVL9GCx9N+/tTT18+qIfEwqM/P+y6vXtfHnZz/mv339y2IJyFRMK2bZ6P5FwNKbkjyeCk9sjRWdMaTrUjpRM1a0R6O3ysm/5DiozR/YQHw6H1TSklJAUmHxtLtmfG+hBKGK8fO5deYCLXE7ylNT4HQmRuNXgT7tLk96/1514yvRt641EDQ2IJPN8R4MJVuUard+Iirc3QFerm5ZDXQnng/jPG1GusFdwY14GV+Zkc05mDjPsWZSYrBiO+5UZkMhV8ii2BDgjW7AwT+P8Qo3LS1TOKzIyOz+DWYcfi8usMf/uEqU73ENToJne6EGmZR3Osikk/A2ZKd3zvA4zkXD6CvRkVOeO3mgEgj1BSqYXUnH5ODaqLtY1tdDU7WZTtJeyWUX97fLOLGatt5Pg7AyqlhQz8fPFbKiqTXL0R5D4rwNTSfYyekR8NBxoG1V4HKHpUAdfueYX9HYPP4V4wbXz2PRe6qfNUsXJKYl0xpxwKMr/fuXvA5KGpYO87hB1CfzqbBYTAfpExfRb85EKVQK2AO30sDO6h+OLU1ZQQAOdQ/Z1PHvfOUiW2cDsJVVsyPDRoyU2dSEhkEgsZXmYOWxIMLz2CLlmK5aPZJEXNiC7NLyHfHTuccW07eQZUFSqoggn0eOrmcZJJGAinXVjjqAZNazq0UuuxFEDlnT8ZNqRbK+HCfeqFFpe7dvIePhx2E9Wkpx9PQg/RxK63TJEQdJ/u2/k7Z6jn7NhDPKqHSEiohhoJMeURYG/mg/ak6szJBCMK8zGbDLg8YfZUwVF087AuMNNx+v1RFNcpljOs4Ir8eyrWwyHf6NNAy++mhBYCvqmjy2ZZrZ19GVObe/1YrkeGkNdqUrvA8BHigIkPM0KyIYJKIZqABoPxhdhFgqE6WzpxZnjGHK91W4mw2knFAxjtqTnhiAZdPGhA8DjP3uFQ7vTV+vjCPs31DLhksnsd8Xn9GcyKDiyLFTfZmOLY3ffwuNuzDIUK2VqHuY6M3t+3RBX/+FQlL1v7CffamTukirW2r244xQhOXbBSJkqNVXG3VWNq62cqhkr+5dLUiY/3bY4rn0NxfjMXNY2HzY5W4GpsHB+GeHVHjr39w6yHF16XZSC0hAFlc1grSMsOlKSZzMSMDIW4uOBTz6BwZCiMNBjEGL0qUFNnsTbPQEST3WWPAE1yExnPm9sSjx03Gk1M748h9aohz3uzj6fWgnw0ufvMRVK5o4j63/2IcLJ3Zjkzsoif1kerrd6CVjS970Jk4wl04zj7CIaW/puQKYvyGFnKLUeLIuzJZbYnk+qD4Plov7XTYfiTzHf3e5m/LThLaaT51RSv6+VyklDT8+cSHTxocO2NQd45o8r0ta/rTwD5xmFUG3Hmy1TJ4fwRkxkm6yUd1s5+MHo1SwjmSEKvyyxJbJ32DZZf8ng0KbkEiIFAxH2vraP0gwLRUuL2BuN/cTe5ZPRqCHk6aCrqYbOpjLaG7JprTfTUhulpdbTb8r+r0fnUTx+PWDlF3tu5a2W5P8VO32D7yRX+ZpgBuTOszGjK4va148RmEIiq+aFvnv7FN4Nhv1j42Evy+kXOMOxKTCTQea2MUUww1mKnR7cQdvozY9BkmBqRT7CAtu721jdPbJQbw/6yAxrCcksg11BCJj47Wo2lR2iEQ98EopNEou8BdRuDtPcmNrjuLmjg+hEI2rLUcunVBpKaQZ7qwzfrfoAtMR9rAAM5ov7XzcdTEB8dIyc+E6SJHILs4hGVAzGkyvyRRcfpxG+cJi9nZ0UOhwUOBwYZJnNzS08umJ1yp1M7VVOMpdVszcnyC6/h75bqsCAaEV3OER2fmy+FqaLetkXHvmf03yHgZzfZdC9K7lpAwCfJ4hzZQfF52bTogZBCAySjF02YEKiQw0NclIdJ1v5/k3n09HoOmbp0CeHl/+ykE9+fyeP1d/O83XJm0TL7U4OuobPMdIVCrDSEeDMW0tRV3npPuRm11YDUy6yoIrUnvzDSSRJi4cTJj6kLJ5uT58/xGhMcBTgNAZoC+2gByh0LqKtd3T/l5KcDIoKHOz1drLJ1wIxRuVGhUbOnEJ6Nsbn0Fp1fSkHr+zCIClsUg8NWNcS7qLF1AVnwowlhdhduax+xUUqLEmhyPHHRhCUUls75zvjgxi1/Un1Ick5KKY5/e8TsXz0jCI+ABxOK2qap9MTQRcfpzChaJR1jU20ejy8sf8A79bWEjpc6lyWJAodDjp8PjDCjKosvLWupPYnmWSKr5lI93QrG/2daKLj2FxfQ9KrhRCjhD8680w0hkcvW79fbcb6OROT95Sz67cNCYVUDhhblw/Hy2FqFJlQMDLAH+b8OeXsqYJmtU9NTZEcSM8dosMdWyKqjW934P3253j7QgUMyd9Bl9gyaXCNfiJa622GGWCcJVMhlbDhhWXMvuofSe//WEKBsZqKSE3F1ngJSRMIJJHkKlHKrNmUWGWag3sIHnMtPWuyxAvDRHzaTEYmVeTiIsheVye1nYklwVMuL4EYxYcAJt5Ryd4lbYS0CKFRLET1oTawtlFWOZHGOh+pnsqadXYu24LDW0wTwR5H1ubhMJgvRJL6rBHhYIT2pviPTXd7bCn/FUVGjaonVd4PXXycorR6PHzh+RfZ0jL0lIYmBC2Hs45mWyx46xPPY5h9RiHGS0rYanKzM+QFX+xTFYc8Ls65tIz9y5uGXJ+Ra6TgphB1MVo3A1qYzRMPMOmXZbh/5cd1ILmkS6Hg0CfOAxsbsGxTWFiZQ26Fk32v7447pHTv6iam1TnpuqWA3bYkUr4LqO+N7/hFNI23aWRpdQkRNRujkpp03ADDpKFIA/FnHE0FPuEc0/3lmRxMzHDQGNhNc3Cw6MkpaKEqv5jajr4LokBQU5qHJcPA9p421rqSSBB4mA1yNxOzLYR7RrY2FZyVg/kzDraJxrj9MD1z95E3z0CGwUpeRznr3kpNblVPUU/qXZBSoD3FMRbHlvrOhKzPPR2xW3k9vQGycod2Tj0R6OLjFKLT52NXewe7Ozr4y/qNfVaNGOgJBrF9YwoVWPE+sptgx+jbmXIs5F0/gcZyiQ983SDaE55TfSfYyNRl+WTWS9RvOmp6zCszY7/GQ104fk/+vaIR21fM1GwrZfcfG5O2ggxFNKLStL8DW1RLOJdFd0svJc+b2XOjjEgwzHliVh77OhOLdnitrY13u67i4ioL84pX4jAlf4cY9KffxCsrgrFwah2KFb15Y7Jvh8HMDGcezcHdNASGv9PuUA8y88w6ZnQuxBvMYqvLy3ZfW0p9HPzRCMavTSb83S0whNXHUmCm8p5ytmbWocXgsDscERGlO+KhO2snS+bNp+VAAJvDzLaGloSSAdrsBupTmHjtCKlITBgNvowWbUQ2lNGYgL8HxG75gL609K5OD1l5qa0nkyi6+DhFWL5nL3e/8J/RGw5Dk89DEx7O+VgNLb/aOGQbARQtrSa4IIvN0W62qD0xzxuPxs7eDjLyTJTlG6k+y4EoC1AnNdMSTvwM6ldDbJl6kPFLSmh6N7lQxOGYMLGQ/RsPjd5wBFr2tnP5G5V4Z9nZkOvGa4gvf0WWMblaKoFolBf2e3lx/zzOKr4QmxH8EYnzKzeTb/8g7v5CYyA+MjJPTJSJkCt4J81RLibZwBlZRXSG99MQGN0ZG0AVKuS+h1MyMd4wn5bUF21mfaCdJffOpuMHm45+elmi5qvVHJzRwWY1uf+DYxnnGcfubT109/aZ0UpzM8nPc7C1vhU1VguBEEytzifisaNJAoGGkPssQ0ISA541BBGriiaDQENFQxVHXmmoou85KjQ0tKH0VwKohH1/weK8n6Y4w2yPEI/lIzPbzo71B3XxoZM6Gj29PH1gO3PHFdPjCXKwI3ET+ka5l6ocK1nzC1HnZmPuVhH+KKFyMweNfnYHvBBKvsT2UOTarASu3cNOSCZ0fhAZF9sgDeIjLz+DnsbU9HtwVR2sgnEmA5bPTWBtVox3jwL2dadmDAJY3XJ0v/OL88hPICdcwJ9ckrJYcGadGPFxKDKbdAqPcfZ8zEp7wsXoIiJMnj19/ijvam1c8JVZtPxiC+VXFRP8iGBztC7uKsGj0WPpRoijichauty0dLnJy7JTVpTNtsYWwtGhTxK5DivVeTk0tbnYtD729AFTFpfyvvvY6V+J4S6Rn97jxGH8Jk6LEadZwWmGLLPAaY7ykao3mJAxMOvpqy23cMB3Fp+a5SDD2IlQW9HUVoTWhRBRGhNwNoU+y4cQYsjqvENhMhvpaO4hv2SIxDVjjC4+PuRs72zjk688Q0eg71ZHAs6cUM6G/Ynl7PBGwmR+ezYru+tBbQcnfY8oQxVrTCnlORbiy84RGw057QhJIInUXTRy8xwIXyAus2csRMJR5D8eYNbnxrElI7a+0zGlBPDs3kK+MLcAsyE+sRlIcUKqochwjr3DJ0i80GkiPSG2gnk55bQHdxKOJnclL3S6gfRVCt9R4GHeX8azQ9Sm7ZzQbXBRmFVKz3EO3J0uH50uH06HheqqXHa1dOAPRVAkmFpWiIgK9tS2s6kzfj+XPauaqFqYRa3bFVN7bySCNxKh6RjjgwR8atJAa5Um4A/bKjnoOsBTO8x8Zt585hQvojcYojcUpPfQFuonGyn93GyIahDRICIQYZXW5/eijZBfJRQIEw5FYk4iNmF6GS8+9h5X374kpvbpRBcfHzJUTeOxnZt4q/4gvaEge7s7CapHzwACCIvkzggHunuQBKTwWh0TwbDoyzSZYrqibiZeVEbD64ndXRxPbp4DKRCiqy21wuMIoUAE+5/rqPlMGXtso9jPJahwOulpT73zZaPXxxM7rufWGY9giDGsVQgIJOj/Eg8ZmWnfxSDC8hkcGsYBORlMsoE52Tk0BhKzdhxPp7qL8sz5NLhTnJUUicWGMpTszj7hkUYkCQzG4UO2e71BNu9uwm41cUZZMQ1tPezcF9sU1XBoGuQbbdTiSriPWYUZ5JoHCp93O67joKvv/9gdCvG/77839MZmBmnGC26aStNj24fd38UfOTOu7KWSJLFncx1tF06jsCwn5u3SgV7b5UPE5vYWrnv+CR744C3ebaxla0frAOFxBE8kRKYl8TufNq+Ps/LKAJDG8AZzV7uLLDk95kD7eam5E8zOcSAHw3S1psYTfzh87iDZf2mlKoYEUhZD+u4htnX28odNt7Ot7RYi6uhRHiJsRxuDMFR7Rur3EY5a8Ybyhl2/yT8u5fsssjiZkqnQGEhdDQ4NlYtqUmuSqDHnMqMjj22bWtj/lkRxqCSl/Q+FNrudsqKRVaYvEGbr/mZ6PKkR3wZPcr+rC8oHVuIVAh7ZNS3h/sSMLGacNX7IdVa7mU986/K4+ywsy+Gp376e8JhShS4+TnJ6Q0Ee37GJK555lGue+xtbOmLId9HbTXFWciFV6+uamWzOIytqYULm2ChkbyRMlpp8cbWhqHW2Dfi1C4A4o0synVaMkTCdLa5UDm1Yeru8VLlGF02ucHpDTg+4PDyyXePfe25G1UbOE6CFxsaZzZGRWieDTt94vvn8Pdz48F38bOV32dV+8cCqq5Kd5ztTe0Gf4SzBbmilK5zcHftQBJVdWA3Jmy4dspFztQraN/RS19rnS+YNhml620JhuDDp/kei1dyGbVF6rIvDsWdDC4sdpRTZEzt/tvgcbHOdS/Tw/8ma7ivY3JZY0sMsi4X/uXQpP/3H3fzs319m4dIZA9bf+PmLyCmIP+y7tDqfrauSS5CWCiSR6tSWSeJ2u3E6nfT29pKZeQJsqycJQgj+tnMzD65eOaR1YzQmGHOp7XSlZCzZVgvGDIVmf/KZQ0diTnEBXY6tqahuPyQ1r5dT91I75lwzzs+Xs05uocDkYNJ+Ow2PjF5ie8qkAnZvqE3P4IbA7rSy7Ys5hJSRvW+NsgwRiYiW/iiTaybYOafyT8OuD7dP5JG7p6R9HB+7M8LSZX9NSV+72i/key8uots3cGrpxjMN3LXgR/3vhVxJuzaZDV4nK7pDhBM8c0rAgtxSGgLDm9NTQY6xFC1czKEuM1pTHqomsVProEcdZgpNgElWsEoGQqrKTEsBnXu8dHuGTtySZbdQdK6PVlPqxdMRJE1iUvcUGjdEaetKIhdOnAigYnwO3mLBgd74HfgdRiPVWRa2xRGNcjx/uvoaLhw30NpWv7+Nf/3+LXZtrOW3L38DkyX+eerdm+r450Nvcv8f70h4bMMRz/Vb9/kYQ6KaRqvPQ7PXQ4vPTZPXg8No5OLKiRQ7jt4xBqMR/uu9N/jX3sRPTtlZFmpjK+o6Kj2BIBXm9CdW2tjSzrnjZ9Ki7CGcyiQFgE22EbzQSWCRxIGQG2gCAc0hN9aJBrKnZ9Gz3dXfXkBfJdTDF5hxE/PHVHgAFJ9XwXpl9Du/iKYx0ZnHvp70hBMfy3P7fYzLXkpZ5mtDrleDCYTHJIDNnpzQEgJCUQd7Os7nO/+eSlQbfEH+13qVSybPpTxrAwCSVkchdVxmh7r6+9jlyiAvO0LE0kyv1ByTaM40WKnJNKVdeAB0R5pAaiLfmcmqV2biC0aQJYnZRflISESiKpGISjgSJRSJEoxED1t7VM6dUs6a7SO7f7t8QUJvmJl6ziQOOFKbQfQIQhbsydvJ+EXjaXsxLbsYEgloONCNsV7mKx9bwOqORlY3xu7E6o1E2NaRuM/NXXPmDhIeABUTCvna/30sqVotpePyqZ584gvN6eIjTQghWH5oL6/W7qPJ66bZ66HV5xkyRv3+999kel4hNqORnmCANp8XdxL5LQDWdTQxPjeHuq7U+Ca0erx9lVLTzNsH2rEZCzijJIeQtZVObfQEQQoKBslASAz/nRUZK3l3GA/4A8FuZnzWgSmaTX7IgccUok1zIwElUhYuKUBDQENUTyYQjGA0KBgVGZOiYFBkDLKEQZJRZAkFCUUCWUjI9IkXWRN9Lu+aQFIFqBoievgR0dAiKmpERYSjqCGVaCSKGoyya1LsFq8cyxgcnMO83ziBG6cOJz7iK3KWKFZ7/NZAd7AIs8FNMJrFz1d8BlWT2NsaGVJ4AKia4Fcrl3HfJa1kWY+GYP679m7+viMEhKAJwE6BfRYTCuxkZHrxGGsJyYNF45Ew2pZg+qtHH0t0zwJ8wb7zgCYEdS2j38mHI7FNawXCUTa8EeWss2azu2hzyq2WZtVCpbsSOWpg4fQCNu5pJBTj2FJBJKKxYvlu/nb/LfSEAzz47tu8uj+9Uxazi4r45uKzR2yTTJE4vydIlS4+Tk3q3S7uf/8NVjbEnnRne2dqs/BNzylk16GOlIVhGmW5zxwwBhEw/kiE9+vaAIn5ZTPJzA7i1rrpjHQj0JCQcRoycSgOTFIGu7vdSEhUZ0u0hVuxK3asig0TZhTJiKYZ2N418gk3qIVpC/fSJPUMiKTcQ59J2SQbUPwGguE0xBYaDz9sAMrhh5lJ1hz2x1i/JJRkhFM8rG1xcdXEKmzG2kHrosHkEp7FisUa++cVAtY1Xs99z45jblUGhzoidMdYAmBzvYdb/3ob180VXDfzMbb3ns0v1g/+jO0+P+2HjhyrQsbn1lCRq2B0dONSDnJGbhEdoZ14kgyjjRdjMJ931sc/XSHiNCytWdPBtAmzMI73oKgKBxzJl6+XhEz+xkls3NdnwrVbfPz0C1fj8gR4de1u1uysT2nBNKNBoSzfSUVhNhVF2VQUZFFRlE1lYTY2ixG71YRJSW9tlEyzmV9ddgXGNO6nfl8b1ZPT7zA8Grr4SBFv1O1nRf1B9rm62NLeSigBP41kKFRsTAnnIUcESBLh3pGLtcXLpMI81rmGrr+SDooz7Jizo2z218NhH1ujnI3TaKE75KcRAXgPP/pobwE4Yvb3E0/hMX905PBQFUFBto2mNIXXDkVej8y5Ig9NAiEJNEkcft03Hk0S5JV7yMv1kWFsYFqmSjAq4Y/K+KMS/qiEL9qXrdQbAV9EIqBCUgpSCGyqzIFD5zBlkny4JKB0ON20RNSaTfmC46wwx+1uUEKkw2+FBqrbimKUkY0KikFBMkhIBhlJkfschmUZIcG7B0Ns6a5GkQWypCHLAkXWkCUN6fBrRdKQZI1afzFPb+r7nWyojX8OPhiJ8sRq2Nz+FbaHuxAxZNQ60OXiwOFZsArnTLomt6Fax1Z4APh3zSUUid9nQYtXfQA79nfCflBkiTPOr2F/RnIRPJO6J7Nh39G5Y18wzLd++wL/+4Wr+dVXrsPlCfDmhn28unY3m/Y2koz3YpbDyqs//wyKPHwMhqpprG5Ivk7OSPz04ksoc6Z3irutsZs559SkdR+xoIuPFPB2wyE+89pzsaf9TSUCFiulNK3rYmfwqDiQJDh3dhnvhhtTkiw0IKVDTAlqCrIozDbh8WtsaujkyJWo0Glnh79+QOuIptIZSkPeaAS+6MjTXKqsYpofQnqJpE5y8aDUBml/c2QT7x1/aqeg8qU4epUAE5JkBkyIYx4aRjRhRMOIqhlRMbLqpalseSZCuCeAv9OLp9NDOBBmvcNC/d/H4c497q46swlujveT9qFIMi2PTyTg10bPmDnoGiBx1Gp0zFIJCiclnzxmXHkOOyM9hNT4BUR9r4f6NTaumX0OIuudpMcSK2Z/Ce9uTGzaVUviR65qgo4dCixIuIu+MbQPzl8Riqh87dfP84NPX85F8yZx/Xkzuf68mbR1e3h93V5eW7ubnbXxW5Fd3gCqqg0rPpo9Hr6+/BU6/Ok4//TxidlnsHTChLT1fwRJ6qtye6LRxUcSBKIRntq9lZ+te++ECI8sycysnjx27RrsbS4EbNvUyML5JbzvT26OOdtqZWdPalKqmxWF6rwM9rW7WDSxkM2+AzQfvn7NmJSHOZRBiytI0Ogfs7phdoMJLYadHbQ0M2f+VDatTX2hqqFwGUaXje/8p4CPfMFM7FXEBBBCHOMfIx1+DDgdydBVN4GXv2Ug6Bvcd9AbhO9GGHdNKd6aAG1lXTGneB4OVWhMOsPKlvdTd4KfMaWUVV3J/f6ry3I4JLkJRJJL2vXc5iBXTL8QJffNtEV0HUvP1lmoWvfoDYcg2fNZY7uHimg2PYbESz20NAwdQh5VNe79/UsEPhHhqrP7cmgU5mRwyyVzueWSudS39fDa2j28umY3h1pi//ztLi9l+VmDli/fv497Xn+N3lBqneCPZXpBId8+e2yyjiqG9E4dxYouPhJk+WMreHXLDjwiytm7eujZ2UbOBePoqbTgs4LbpNFtUnEbVFJ9pjk/WkF3rYe2Dg+7xMhhbs27ujl3YildhhDbfZ0JWdw9oRCl+Zk0+I7eRdlNBmaUZmM0StS2+2lwDZ5XlhBYjAYUSaIiJ4OcbJl9vmaaVBeFFRbqIwPHXuvrZHaWg3ZfB+2esUuv6jCaYi7Ovr1sH4UHSsYk7K8tPLogWvE8XPzRJTgL3kjpvqNhIw9/aTZB3/C/r/rNTdRv7rO2nfHALA5dmPy0XG61Bu8n3U0/O3Y3c9a0Eta0JyZAKkuyqVc8+MOpyRb60nYPSycvxVo0tMNuqrB4xrFyW2LCA0g6UZyqaeQ2VdFTmZj4yIxm0tw5/BSZJgQPPPIqvmCImy6aM2BdRWE2d121gDuvPIv9jZ28unY3r63dQ3PnyFOmrV2eQeLjN2tW87NV8RdXjAeHycRvrrgCcxqTBR5BCHFSWD1AFx9xo0ZV/vCNx3j2Vy8PWtf16NFQRyNQCJQYFJyFmeRMK0ab6KS1xMDurBDjPCbGbfARKrWxtxzq7bGd3CQB9Ts68QxxNzoUvd4g2zb1XRTOnVnK29GmuAVIVNPQfIIZWQXscnewaGIeO/317Aj3QBiww7T8XGyKGU0IzJKRzm6Nwhwjmz19TrcN9NBwzP++Jzr0hbXW3cuYeLUeg91oJBDj3FRYjpC5QKMtnpmOBKlRzeyOod27/ynnyhSH7L/yi49wYFPslXa2/3AnlTVVdJYlfqcLEMrsJpVhVaom2LG9CWuJkUA0vqnD8uIsmow+vOHUpot/bbeLj2bNIGxJTTr1oWjfUgMkLj7UFOSM2bqti6KyTNwxhIsfT6G3mGZG98/5vydXsuNQGxfOm8hZUyuxmo9OsUmSxMTyfCaW5/OF685m24EWXl27mzfW7aXLfdQfrDDbwUXza6gozBrU//VTp/Hs7l0c6knudz0SP754KRXOwftOB93tbvJKxmZfo6EnGYuDoD/Ed6/9KRtf35pUPwajghB9QgZgwiXTeOXa2PIj5Ms2lPfDCc/JTp9WwjuiOaHru9mgcNbUYtYFd6Yt7/qc7Co+aB07x1aAqTm5NIbjcyRb1DyLUJsAARoCTRNoQkPTBFFNQ9U0oqogHFWJqioRVevPp6CqIqbvf2HIyv6NowuAGz6rcu51j8U1/pHY+945/L+Px3/xKZ5USPQhiJgTd6xUJJnmxycS9Kc2YdrMqaXUhnpp8cY2pVNW5KTdEkibqX1aQTbjp7yatumXjX8/B7c/8bFPLM2j1eXBHYi/jyX2AiIKaJIgI0dlw6x4RZbgDzVN7NxXxU+ezCeqxlix1aAwb0o558waz9kzqynKHfr6EVU1Nu5pZNPeRs6aVsnM8SXII2Q77vD5+PvWLaxvbmZTS3PcInbY8SoKX1mwkM/OPzMl/cXCpvf3Uj25hKzc5DJgD0c8129dfMRBwBdkWeZtpPorm3jlDF6+Mva7vQvCFWxfn/gFevz8Itb5481KKJg2roAtXa3MKS7ioDiEmuo62oDTYKdMKgUkhBBoAoQmkGQJIUfZ6mlEpFj4zMkvYH+gLq5tzJqR/DdLaU0gg6EkgclowGRU+nKGHH4YFBnl8MNmNOBZUUtPDNE1X/2plwmzn457HEPh6Sjify6dTW9HYtNKM26aSsPdyfkHZW2fzZYPUu/YZ7OYyKi0c9DlGrFdSX4m3Y4QPcH0Oh1dOcOBkvtWyvtVIpm8+dfks8zOmVDKhoPxnWfOcObS+l4T0cPl7m12M8GvdRMyxGY9cihwzzgP1UpfRjFNjOeJty7hn2/Hnz9mUnk+Z88cxzmzxzG1qmhEgRErEVVlZ0cH65ubWNfUxIbmJroC8Zc3uHJSDd8+ewmlY3yNe/flzSy5fHba+tfFR5oQQnBD8adwtae2qJj007PZmxm7aVcREnMa8qhrTMysKkkwfXop66V2fOro0z0CwdwJJaztOGodmJ6fT6uhkZBIbeXMKcpE1jYNX7/mrLJSNvpGT4UeD2cVFbPLG3+f44LFND0XQU1DIbX5OTm0bWrC0+mLyUh1w2dVzrvuX4g4wouH4tEv38wHz8YnxI5nyl+n0DIh8QrC4/0zeetv6alXM21iMWvdwwvvorwMejMjdCdwQUmEa88woTmHqXKaIBb3JJb/I/kCjQPEhxBkmy3kmyzs87oOh1aDQZKoycgmRzMQ7Q5xaHfboNwb0z5WxObJu4bZi+C6IoVnW1UEEo9N3YJB2zGolTe4mJ/+cz6b9yfmKZCTaWPxjGqWzBrHWVMrsVtjrwQ7EkIIDrlcrGtqZH1zM+ubm6gbQdzOLCzkv849j3kl6alhNRofvLqNRZfMGL1hgujiI02sW76Jey//Ycr6UwwKZXcv4I3J8V/Az9PK2bk6OS/+nEwb0QlGakMji6m540pY0zXY/D8xOwevtQOvlpoT9SRLORtqR59bPaM8n+3e1MXbLy4uZpsnMUEz3V1N+7tRunuTu+gPx/zSAg68sDO2tudKfPK+5QgS+10cWLOYn340ubvDjFwHticzCdkS95MoM5Tx/u/Sk7FVkSXMFTY6AoOPV2GOA2+2Sqc/PcdyKGRJ4tp5gqh9bcr6NNeey6uvJ2e1qXBkYN3rx+8LYbEa8XlDhEN90w0TzqtiVW87szJziezupWMU61xuvoOca0zUFjThNx0dl0WG/5vUQg4r8EqL2Oit4Bz7swwfuWWktv1KHni0mi5P4k6TRoPC3Joyzp41jmuWTMdiSj4U+1jafd4+IdLUxPrmJnZ2dFBgt/OtxWdz9eQpyGMR6jQMa97cwVkXJl5ldzR08ZEGhBB89Zz/Zsf7qSl9XTylBM8t49mUm9i87AWRCravS943oqosh83O7mFzgcwsL2STu5nhfiQVmU6kTDfd0eSKztllK5GeDLqGuCgMams0Mr04n23eeiIi2akfwfQCB/X+xAvhWDQT0+snsvmD1BTYkiSYX1SAe38Xbbvb4nLPqZ4i8Y1fPgbENy+tRgz86IqraNiduMUCYM73ZnHwouR+l+ny+zhCfraDgjInGztbiB62WuXnOAjmaLT70pfHYTiMsszViw4RNSQewm1tn0fjzjxCIcHeelfSY5rnd1C7b/jps8ISJ23N8VmAL7pEJXeRiZWWPovtg+O2YNQScLqVcnln29X8/F/ZqFpykRvzJpfzsy8uw2ZJjSVkKLzhMIokYTWmVuQkwroVu5h/fvoKP+qF5dJAJBxNifAon1OJ/+oK3ikMgjS68FA0OLc7F1uPhhQWENYQYQ1zbmo88Gsbuzm7oJR3QgMvGAWZdsrzMtnkahlWeADUu3spjNopyjXQmkAmRejLSVIiStkYiO3i7YtEWFPfTLE9G3OmSlMwcU/0ic5c6v3JXSyDcphAbuouWlNK89n/fJ+1I957pEO7BIHeBVjjNOWvevJqGnYn56tRMqWIg+clL4hVoVEz25oWvw+Ajh4vHT1eSpx2iiuzqI948J8g4QHgVA3YWgpwl8cnPqyuKfjqKxAClq9zIUTqIjKi4ZFFfTzCw2jS+P6PGykp7fPjuLxvKWgJTtmKLs6Z/giLpk3mkVcv4sUPEk/pv353A1/8+b/55ZevxWEzA7DNdYB8czYFlixkKfmwVIcpfcImHnq7vWTnZ4zecIzQxUeM9HZ7UUwG1CRre/ReW8GG/BDDXVbMmkyBZiEnYqDIZcT3VistLYNFT3ZBBoy3piSHyPaNTZxVnU+0WoEMmZagm3pPL81dsYXItfl9hFQLefmZdEbjDKsTElMNE0b08xiOFp+HjIiZifmF7PMndtdY6DDSlgIXnoghdRlgLR3JmcxXPj+By26LXXx4Owt55ifJp43P/UIefkNqrD8541RIb3oFunp9hPZHEdPNtMZY6yXVnO3LIvhCA3tkieof5xB0xubHZWufz5svGwlFXGkZVySSut/znZ/19QuPY/aQdL8GaTefunQ3Hz3nfH70xCx21Sd2Oduyv5nP/+xf/Por13Eo0sg3t/wOgAyDjT/O/yZ55qykx3oy0Hyok/HTToyvyVDo0y4xEAyE+fWDL9K0v40d/1yRVF9Trp+L94wCzH4w+DTwRIj0hgm4gni7/fg9sV94Ci+dxE536pxfqy4rZG1v4n4kUyosNIbjm76oNBeyvS65O06jrDCzNJet3tjzUgDYFAMZjiCBGJxuR0XAjM2T2bUruWmLgmwH0TfrksrhLkmC3y5/AyHF5hfz5LdvZuWTyTmZAkz450Q6S1Jz911gzGfj77NIIJt5zBgNMjnzstkR7Bq9cRo425dFx6NHbywqp+Zhv3crmmno36MhYCNj3wSCWzLYIow0h9LnFDupSaa7MzWCrLxc8IP/+yeQ+NTm6JjZ23w133+8HLcvMWvFlOp8zEubqTsmEnBJ/izun/aJFI1x7NE0jUO7Wtj47m7q9rbyjZ8lWPsgRvRplzh47bmN1B3s4JqPLyC/aGBBn94eHy/+Yy0v/mMNvT19vgjTr1rA9hdXj9inyWbG7LDgd/kGWEqyynIJuhVa/poav5GM9iCYRcoyqLpFclM52XI2jXGeYCySGUhOfEQ0lQ0N7Swor2aDN/ZKwrML8tmeoKPpICTYdsZuJk+qYP/ziX+eKouN/UneDwgh8dB/L+WWrx7Ckfs2jFDdp37L/JQIDwCTN3Wnk/ZIB+fcUMyKJ9N3gZWQyI9aMSIRGXFycRQElNkyKDdmIkX6QsOh71kAsiyhmGUCUpTOSIDmgJviqAXvM7UDuqnb2Un1T2eTN8OAoSCMyPeiWoMoewvpWgf7NnegRnuBXgorcvBOs+FOUc6J44mksHpzTw8E/BOw2tIpPkJMKnmav327kKdX38DjL5uId9IyZ4GfncelIHi3YwsbuvcwN+fEF2KLle72Xja+u5eN7+5h03t7cB0WkedcMfvEDuw4TlvxEQ5HeegnL/PKM+sBePbvq5g6q5xIOIrPG8LnCeLu9aNGB5643ZqC2W4h5DtqoZBliXGLp2IryaehxYPb5ccPWKqNVFbmIEsS7Z0+Av4wDYdS9w+4f2M9i5ZO5gNvDwhBjtlMd4LZGAWCrkhyJ/q9bb0QR0FGIaC1K3UnuZ1tXUeL2o5CvsXKoUDqk5k1WtqJeRBDkSJH+G1rBd/+WBVT51Zz69fqycxfwfHV2jRV4qn/rgRSM1WiPhGB76WkKwDqMrczftp0DuxIjwAJR1W2r29kSoETbxUcDMdvRaxx5BBqDtPR5KWH2CwFdkWmujZMq3ewz9ehre0cGpTDcPCUYlt9N5PtZjaXy4QTqEA7GqFQ4v+XRVPzsVZlEGjwUO5o5/NfehVZPpDC0Q1Ps5jDm/ntLLx9PFv+peL3DTadGU0SOfkmnLlGbFky5gwJQ2aUnYaB+YQtsomPV17MDOe4MRl7MnQ09/DcX95h47t7qN0z9BR29dSSMR7VyJyW0y7dnR4e+MoT7Nme2MXHmW2jvKDvAqPKCo1NLjzu0adLSityad6dXHjsUBRW5BD0hvC4/CjnV9IaZ44C2SBRdVERq93Jh6/mFAeJiNhPXDn+UupSNHU0LiubBim2C+nZpUVs7Y3dShIzAhzPFBJM8M5xSlk+bc8PlxMhcWpmS9z2tUayit7iSCTM6n8s45FvpnbKYfpD02mckRoxA1BoLGTjHzOJpjadzCAyHRYs0+zsDA3zfQgotNopNjmwCyNqSKWj20erK/4or0VZBdQ/llyW5COUTy7COCGLvSJAZzhF2ViFIG9r4oIvf0IOW4r7zoczS1V+cPW7WJV1qRnbsJh4xvURXug4KjYKjNnkBQrRTFHCSoiAFMCtevGpI5+rJSSWFs3nk9WXk2tOb3n7VLD5g3386O5HcXf3WVwtNhOVE4uorCmmalIRFZOKqKopJqcgM+nij6OhT7uMwvtv7kpYeAD09vj7p2FipaQ8h87a5PwBhqOt/qiTmkOJr2JhRrYZ41m2lAgPgFxDxqhRL5eVTufq8tl8YfXfmZibQ727F4Msc0ZBCWUZmYRUlZCqElajRDWNqKYR0VQimka0/7nvEVbVw69Vsi1WGmI8//rUNJnzJcjJttEcQ2bSoWju8ZCOmpN7NgueeGISNZ+xM9keIrKpm7/9V2xjLJ6XQ/CWRWT0hKl/YHBBNAE4qnLx1XahvRyBFOYwaou0seTGAlb8LXXZRnNLbOQstqCoCv49Yep2uHB7g7jXBFlQnU+kSMJv0MiSzMhR8HrDtHS58Yb87EsyiZvJoBA54ErNBwEadrfC7lYMwNyJBdirs9ms+JKajjEnWfW040A31opMApEoW5sUrn7oXO65dDoXTHwRRHqmXvZELx8gPADaIz20G3r6Zh1jNA5Ny6zm8xOvZVJGeeoHmUICkQgGWWbbB/tZ8dwGrr/rvMNio5j80ixk+eQoHjcSp6X4aKpL59zj0IRCkf4kPenEGGdomHmBnZ29qRNF+UruIPFhlBWmOUvoCvm4oWoet49fhCRJ/PPczzIxs4DeUAiLwYDVkHwcvCo0/NEw3kgYXzSMLxoiqmkIBJrom4MXQoAUxacG8UfD+NXQwOdoCL86/HNwBAfVTNVGV5zC9FiyHJYYymnFi+Di/8ogNGst7dEo7b1gHm/nqv+Us/zjjfg7Qn3TPQLO+nIllmwDbZt8lC5x4DzfR3fUwqOvd2Mwy8w8q4quNbV9vUpQcduZSOOy2emIUvbITty1sY8+I5CH3Z9D1BQkaPLiN/aioeGIZGEJO9AklS57Mwdt27jgk1N552+RpCwgdqeR8gucbBQNHAwdvlCNh8KpGVRH8wjVqXhdAYqsmeze2kJ94rsaFlmS6NiYeusnQNO+dtjXzvyPnkHAAFFZ0BMNIwF5spmeNi/u3gAWi4Gsqiw2tB02zx93N5xtsaAqPoSaoFFcQFVWFrs6jpxnJX60PI9fr7idb13SwcLKV0EkntNkKP7YaOnbcYJYZBNfn3wT5+bPTrt1IBGEEOzu7OSdulravF6unFTDnJIS5iypYc6SD48/yrGcpuJjbL3bFUXGYTbQPQYzXLki9ruWzBxLSoUHwM4WF8XFWeRYbJTbcji3qIZzCyeRYRwci1/jLAIg25K6bJaKJJNhtAy5v1ShCo3AcaIlEI0gSxLGiIHQeEEgFCEQjBAIRQiGIviDh59D4QHvj7Tz+EP0egPkB6WUig9btsLSn0fptA+MWw0JH6GiHVz4Vj4SFrx0YZEc+EWfM3TR9X0eIt1AxHc20FfduOvjkzHdPBUjErkbO3h3ggHwQAAsd81gXNCFl2OceDWJrEAh5ogNJWxBDloItZhpOOinodtH3y2pCcgBcpCkgYE+BXmzqZpiIZjvY9H8XNZv6MIfik+BGE0yEy/OZZu5mTXRwVa5tpCHNjxQBpSBzZ6+ZFDBSJSiRZU0vpk+H4j2Lh91TQNDdo+dCPN4oaPTS1W2Hb8/TFlFDm6bRqbRBL1RDu5pp2hSHg4hU7s7fpEgG2RMQ9wEeUMy979QiNV4C99Y2sWSca8jieT9rrq4kM5IcufWoBbmkLeFc/NnJz2eVOIKBnh6xw6e2bmDcdk5XDmphjvnzMXwIbBsjMZpKT4a69MnPmSzQt6ZeRinOvDka9TSiwZkqk4KwtlEjRA2aDgboel3qYl6OZa9r+0+6oQ6Cr3dAXItNrqCyaeTLrI7uHr8FJZNmMLU3IKT8u4hVSiSjMNowTGcwCmMv89///U9/vTTl0lR7A0A4xZamfz1ejrF8DlUPBwVn37hGrJNwO8E+qY96jzHTNVMGHj6aPJ66a0uo7wuG6GBCBg4sCXItu4j0T/hw4/hOV6ft3d6aX/3iCNngCnjithR3xpzhFfF9Cxap7hYFToUc9JXUyi9mSgjrvTWjfHGWM22u6fvuOzf14bRoNChqv3ff3Nrnx/WtDMr6TjQiasrtgguR54N/8JMtrQPn7AuEJH5n5fyMRlu4qsXurhg4hvIJB5xpUhR+qweyZ1znqh/nWJrLpcWn5VUP6lge3sbf9uyBW84zEXjx/OvGz920iQrSxWnnfiIRKK0NbvS0nfZXeN4L6+N3VpfONyx59k1+EGh70YvDNmlVqYuKabz3fiTa43GgRX7KD23gqZRalSMm13AlkjiU1B3TJ9LaUYmU3ILOKuoDOUUUOMniouvncuLT66mtSGxYoHHs+TuDJTz1uMWyTshuv1mjoiP0WhtD9O6drADcXlJNlabkb37k8uiuutgK7Mml7LlwOhTF5Wzstk3rpVAnJaSqD99Fkqr0UD75tRHWR1BAB5v/P4xkejQCVV27GuhsiwnJvFhshtRz86msSu2XC/hqMxPXs3h/73+Ub54fi+XTnkLmfgtQlnibe6rvoYHD5nj3vZ4TrSD6bt1tfxy9Sq6AwF+eOHFLCg/uX1PkuG0u1q0NrnQ1NSHppV9opq3cpoJa7FlReqJBNiwNEzx9VUpH0s0omJY18psZ9awbSYuKGKDvZ2gGr8fiixJFNoc3DrtDO6cMY9FJRW68EgSq81EaWXugGVGk4HLbzwTQxwOgEaLxNW/MyLOfZ9oCoQHQE8cRbz2+XqoObOUqVOLyc2yU1Wey+T5pWzNcrHa3EnNxATMQsdhiMHqUT07m73jWhNKIGewpc9qN8eanbgvRQwoBhktxdO7dY3d5JdkjdhGNsg4LinhQIzC41iimsTP38zi0t9cy3PbP43KpBFaW5DlYiTJMWDpJMNznJ+T3L20IslMd1Yn1UcyrGqo53P/eZEFZeW8csttp7TwgDgtHw899BAPPfQQtbW1AEybNo3777+fyy67DIBPfOITPProowO2Oeuss1i9euSkXGNJOpxNS2+u4q2C+C0YQTXKOzNdnJc/iabf703pmLy9AbIP9kLu4HWTzi3hvUhDTP5Zk7LzWFRSwYKSciZm55FjseI0W05oZcZTDVXV+N/v/JMN7+3rX3bOZTO442uXUliazfgpJfz6e8+N2k/xFDPzHminW0pN0rAjtHtiv4BHEKzytfRZwEsBAuDn6DRJdvLG1tEuruPm5LCrojkhYQ0QMaXHMbwo00Ht37alKpXLkGhRjWmVBWw5kLpwZwCzefjjJoCiZVWsb09unwKJ36508tuVV/P0p9aRZX77mLVGioo2ISs5fW2FSjR6gEh4E+HIZiLhTZQkafiYnFGJVUneepIIezo7+f26dTx9w01Myc8/IWMYa+I6E5SVlfHjH/+YCRMmAPDoo4+ybNkyNm3axLRpfWV6L730Uh555JH+bUwn2TxVKp1NC84pxHuehRXBxP/pNCF4q6iN8741hZafpja/gzTOCb2uAcuKqpx9wmMYKjKcLCypYFFpJQtLyimwOYZtq5MaPL1+dm7qi62YOK2UT3/nCqbPrepff/kNZ3JwdwsvPbVm2D7m35xB5rVb6BGprVEiNIVWd+pKzK/yNrNwUiF79iYe7bBldxNVxdnYrSZkWcJgUKjvcNHl9jN+fg7bS5titkAOhZv0+GQEo1EyDDIimp5qvUeItKU+XmqojKcCqFhYSnuBSFp4DETis0/M46+3eTHIfjY1zeGNXYV894rPoBhKUZRSDEoZilKG2XIRNvuNALR2/Q3YkNAe7YqFmysvTuFniI+Jubn89drrTmlfueOJS3xcddVVA94/+OCDPPTQQ6xevbpffJjNZoqKilI3whTTmCLxkb+4gPcW9EIwNQmy3rG3cNbcPLo3pM4y07u2ieLZebQck3SspdbFvCklbHS3oAlBntXGopJKFpf2CY7yjJM/qc6pRlaOg+///nb272jiwmVnDBmj/7n7rsLnCbDypWOTUwmq5tuY/jGBq/wDQmmIppKjRUS1FF4sJYluZ/IFW+pbBpr3HTYTC86sZEXx7qSEB0BLqBeTyUo4nFqR4PIHmXxeNY1vpDfbZ3CIzKnJEj5OfJTMLsRTbWJNZ1dfSFSK6fLJXPXQBccsEdyzdCuE3x/QzmAYT17+K8iyA4uc2I1uTUYF9029jWLrEGbiMeJ0tCQnbANVVZWnn34an8/HwoUL+5evXLmSgoICsrKyOPfcc3nwwQcpKChIyWBTQaqmXQzzM4HkK4EeQROCjcsEi8+YSNOf942+QQz0dnnJ26kQnuKkK9h3QspzOpjU5eTjc2YyfVIxE7NzTyu1fbJSPamI6knDi3ZFkfn6Dz9KwB9mzYpdnPfVTExn78CvuXClcVxqKPU3EoWahVQFeOdk2yiryOVgRzfvba1nZqCS/VWtuKTErTU+Ncy8JcVsfTP1V9VQlS3lfR5PoNcPzlRbnPvOEQU1uWjTM1jf3p7eOnFD7F8VxSjHBaJHowfodX2LrOzfcm3ZObijft7r2IoYZU4535zFzKzxzMqawEWF8zDKsV0KWwO9FFn1G7RUELf42LZtGwsXLiQYDOJwOHj22WeZOnUqAJdddhkf/ehHqays5NChQ/z3f/83F1xwARs2bMBsHnouLRQKEQodVepud+ou6EORimmXspureMuUmoygxxLSoqyr6mV6lZNobxSDUcZgUFAMCgajgsGgYDDKfe/7l8koxmPfKyjHbGc0KjgKM7CNz+GMmjLKC7N0sfEhxWBUuPdnH+P5vQ+xV36NNFvvAQgFskkmedNQiBT8/IoKnWTk2djV0E7rwaPRI7v2dZDd6qDkTCc7lcQjydYH6pg9u4qdm1Mblr+9q4uqHBuB7tRNZR2P35V68WEvzaBsQQmv9TYg2tMbKjwc0jCx0oHAc5jMC6m038r90z5Bra+FJ+reYGX7pn4RUmLNY6ZzPDOzxjPDOZ4ia05c+97YVcdDe1ay193Giku+gRxnMkedwcRd2yUcDlNfX4/L5eKZZ57hz3/+M2+//Xa/ADmWlpYWKisreeqpp7juuuuG7O973/seDzzwwKDl6art0tXhASEQoq+YGqIve5w4/AxH38uyhCRJSBJIstRnGpMkVJuEokhISDEXlJVicDOTJQmH0XxamuB0Yieoenns0LfpCKXWsXQoOuuW8dLm1PqRyMBsdxa1CYYVV1bk0hzw4QsOnzPEYTXB2WHapMT9HzINFmw7MunuSu3FdkFOAQ1/3ZpWx9Po9ELCkeSnt45nyqxSNkrdWJCZoWYihzSiNpleRaVZ9dMRCaSsyvbxvPbFJ5HEcDd9JvLyX8Bkmtm/pMHfzn5PIzOyxpFnzhq1f1VohNQIITVKSIsSUqM0+V38Zf97rOk8moHn70s+xczssgHb+qMhbIYT46x6MhFPbZekC8tddNFFjB8/nj/84Q9Drp84cSJ33XUX3/72t4dcP5Tlo7y8PK2F5XR0Puy4I5389eA38ETTa/s+sPMjvLcv/vDJ0Zhtz6d+bWJOp9NnlbHx4Oh5PsZX5LCjph5VSvwUV2MrpO5NP1qK3WkWOfOpezx9kS+WsyrodqXHulJSkIl3Xzeh4OAoKIvVRF5JJvZcG+EcA5vUHnxaaqKHXv/ioyCGzxOjKJXkF7yKLMd33fjxtpf5R+06ojFWB/7UxHP40pQL+9+v6jjAA1te4JnzPo/9NBcgY1pYTggxQDwcS1dXFw0NDRQXFw+7vdlsHnZKRkdHZ2gyjXl8pOJeHjn4tbTup9Odnrmdg6HehE8+kRjz9Byo72Z+RTWrbYnljbVgIONgBuVOE2aLwt7W1E3BfNDbwYJPzqTliR2oodRbKGwWE91JFsEbjuZ2N3NnlrJ9be2gdcFAmMYDnRzJFVZgMlA9u4TuPIktoS7EsXJLiDitJCMnT1PVOlw9XyEj8zsoSgmyHFuk3oXFU/n7oeEjyY7nX3Xr2edpw2EwY5YNvNC4hQJLJoFo5LQXH/EQ1///vffey2WXXUZ5eTkej4ennnqKlStXsnz5crxeL9/73ve4/vrrKS4upra2lnvvvZe8vDyuvfbadI1fR+e0xSil/0TX6k595ASAOxpman4mrR0J+HjFEYUi9SqQgI9nARlk781ie+tR68wZ1cXUdrvo8aVmGmZ1VzuVN00gf1MvHVtTm5fDYkpv8uodrV1YrCaCgZHT5YfDUfas7Qsjn5bnoGB6Hs22CCVeI4fWNGBzmMnMtSPG21kd6hhZjMSQNC8YXE4wuBwAScpEUUoGPwxHXhcjSVbm5lZQbHXSEogtcrEn7Gdl69HSGAWWDP688HbyLHpagniI6xfa1tbGrbfeSktLC06nk5kzZ7J8+XIuvvhiAoEA27Zt47HHHsPlclFcXMz555/PP/7xDzIyMtI1fh2d0xZ3EqnxY0J14AqkrpT9sRiRMBgSdNqLYwrE1RmE4Q2vQ1IjiujdpFHndQ1YvuVQCzazkVnlRWxpSI1YqOvppWm8zLzcahpXHEpJnwAGOb1+Y2pUIxyOL3tsT6eXnpV9/kNHLt3hULQvdfteOG9OGVvzgnSrQ4kMAcQnhIVwE426iUZ3D9tGlnNQlBI+XzWPH+61xZ0RN9tk408Lb6fMnh3Xdjpxio+HH3542HVWq5VXX3016QHp6OjERm8ktRWJBxEuSUu3khCcqRWyqyX+aJQJ4wrY0xn79EdzuxuDkIlKsVlL5kfGsWttF9Fhpnb8oQg769qZWJzLvrbUTMNEVY3teSrVE3Lp2Z98nwazgbpWV/IDG4GinAzam2MrNhcr+zY2UpBpxbook6bIwCkja5pq/UXVbnb5p7PNk41FicYlPjIMFv6w8DbGZRzNSKoJTY+EiZHTrrCcjs6pgjvN4iMSTEOaZyFYYihh+9b4i6uNr85nv6uH8DBF0IYiqmqUk8uhUTKLyELiLPcENm0Z3aKhalrKK4x6Q2Fci3IxtXkIe0aeyhiNyrmVbOtKTfLD4ci2m0muRODQeNwBJoWLaTouT0uGOXUev0LAwfBC3u09kxUdEq1BHxBfRJdVMfG7BbcwxXnUrLa9p4lKRy4Zw1W71hmALj50dD6kpFt8BANO4jV1j8ZsRwHb1yZW1VU1SQmFjxYHnXRbvURQiQqVMOqA6usOzExoLGXTodinUg60dCFLfUWqy3Oc5GTaRi2uKEmHZ4wOV3+XkPrD+yUgpKoEb5yA6ekDhHuT+N6zrJBm8WFMRbKWYXDv68FUIxM+JvrEMcz1PKKZCIhsAlo2PpFJQM0koNnxqTZ8mpWAasanGg8/ZHyqRJ0vQlPACwk65JpkA78682PMzukr/BaIhvnbwdVcUTZTFx5xoIsPHZ0PKaqI4DAcLrSFgCN5a468PuZvX06bY5cNXHssR94HgjaMcoRICtOrW13xtZckidLiLCRZwmFJzNqgtsuE6gyAARmwAIosYVAUjIpMVNXYHaeQ8wbDVBVkY3YY2dHUTp0rdRf7CR8Zh+XftQR74ndsNZgN7G9KbWK0oQj5krPOjERrQw9FzTKBi3LpiPb5HDW64L7/fAUhyWQXO9lo7MIXjRAaNYw3eviRGidhRZL5f/NuYEH+OAD2u9v52c7XuHfGFZTYslKyj9MFXXzo6HxIubZ86Nw5KWMKcLickyYEUVUlqgmiWt+zqmlEDz9UTSOqakTFMa8HrBdE1CiKW6BdJtA0gXZ4+ZHXmuh7rR55rwlavV7+vmITqiYGFUlMBlUTqFqUUHz+hQOobe+hSk69o+H+7h7GXVeJ5bl6gl3x3Z2PxZQLQFd7ejNRF5RlsSEa4IiJSlKM+B15dOZ1sCHUCOnTPsMiIfHDM67jvKIaAN5t28f/7ljOb8+6WXc4TQBdfOjo6IyKLEmYDAb6bA9p8v4bhiWzx/OdP71EV4LVdT3u9ETsANS29nDGxGI21Seeyn0oDna7qFpWju0/Tfjb4/BHyLKkfcrFajbSXZ9e8dFS1835xZX0Xilo0DrpDLvYSVeqZwFjptKey92TL+DS0ukAvNCwmd/tWcEfF9xGuT2+VO06fehuuTo6Oic1cyeW8fd7bmZGdZwxs4exWYzpyvgNQFObC4sx9fdxtT29+K4owZZvj6m9YlLY35yGErPHUZwzNvks9qyuw7I6gkE5MZcpCYnzimr4w4JbeeGCu/uFx2MHPuCXu97kDwtuo8Jx4irhftjRxYeOjs5JT0G2gz9//aN89NxZcW03oTyPvR1dJFdEYmS63QFmlBSmpe8Gl5u8S8fF1DanNBurOdXVbAeTmaDvTSLU/qeF30++mY9Uzh2zfWaZbNwx4WxeuejL/PrMj7OoYEJ/+Ozypu38df8H/HnR7VTqwiMp9GkXHR2dDwVGg8I9H7uAaZWF/PCJN0cNua0oyqLB3RtXaG6ibDvQQnVxNoc6Ul8HZ78UwGKU0SIjO/52HOqE2i6mL6imLhwm0ymYv9hDRlYYe0YYqz2I0Rzk4M5SVr1twu6QUAywa/swHUoCh13C7gC7A6w2sNkEtuDYzn0c2NLKd6+8mouLp/G9Lc/HnIl0NGQk7AYzNoMJm8FEntnBsvIzuLR0GmZl8NTiDlczv9j1Bn9edDvVjryUjOF0JunCcqkmnsI0Ojo6pye76tv4xu//Q0v30L4HRbkZ+Ijg8qXP3+N4CpwOPETwpOHivDCngKaHt8bc3p5tY/KCKGd9+flR25rFuD6/TimMRghBCE0E0Ybx6jQF5/D/PhGbNSYVXHLVbL7+31cD4I0E+fnO19nU3YDNYMJuMGE3mLEefm1TTMcICnPfMoMJu3JUZBxZb5YNSDHOx3UEPXxhzd/50ZzrGJ9RkM6P+6FmTKvaphpdfOjo6MSCyxvg3odfZvWu+gHLc502MEu0u1ObgTMWZk0oYUPD6BV3E2FOXgHig1a6doyc3iszT+L6/4rimPIBqogveVasvPXzG9ixJj1FB4/njPnV/OS3t47JvoZCCMFX1j3FXROXMCO77ISN48NAPNdv3edDR0fnQ0mWw8qvv3gtn7x0fv+yDJsZo914QoQHgCql715uY2c7oUUj33WbbXDnXw9gnfxa2oQHwIW3p7YQ3nDMmlvJfT/8yJjsazhebd7B9KxSXXikGF186OjofGhRZJkvXnM2//eZq8jJtJGTb6epO/15LobCabOwvalt9IZJsLOzk5Jzq4ZZK/jUn3oJGfYMsz6F5KyipCp9IURlFbn8148+wk9/dxuZTmva9jMaveEAr7fs5I6JZ5+wMZyq6OJDR0fnQ88FZ0zgp5+5goNt6Q81HY5xZbnDFqRLJdL4oc3ZH/+BjJb7Qdr3DyCIcOXnUp/rI68gg6/eeyV/eupznHPh1Jh9MtLFn/a9w9enLkXRi8WlHP0b1dHROSWYM76M//vEFWnJuTEqEtT1uMZkVzsCvTirBmbUPOdjJrLnjW1VcXvVKuzO1PSVkWnhri9exCP/upvLrpmDYjjxl6YfbP0P07NK9bTpaeLEH2EdHR2dFHHx7Ek8/MWPkpdpG9P9Tq0sHDM/E3cwRO3CTErPrwZg3Gwjsz+xAkH6Q4qPJSrcXPO55PKcm80GbvrEYh599kvccOsizJaxzZ47EreMW9CfWEwn9ejiQ0dH55RiekURf//qx6kpzR+7nRrGdnrAEwyxOT9C0cxslv1wK1GR3nTnw1E0ex2yHL+TraxIXHndXP767y9yx+cvxJFx8lWDrdJzeaQVXXzo6OicchRlZ/DXL93AedPHJh+FfAJ8EwKRKBd9x01Irh3zfR8hTAuXfzK+bc69eBoP/+PzfOk7V5Cbn5Gegemc9OjiQ0dH55TEZjbxszuu4hMXpD81t1E+MadSZ2bnCdnvsUy+YBswuvVj7lnj+O1jn+K+B6+ntEJPTX66o4sPHR2dUxZFlvnq1efw3ZsuxpBOgXCCUjVu3FuNLJ3YKYuQspdFVyjDrq+ZVsJPfnsrP/r1LUycnFhxQJ1TDz3DqY6OzmnBun0NfO2RF3H7U5v+3G42IdsVenyBlPYbKxV5gjsvbcJiX4EgekLGYPTP42d3VA1YVlaZyx2fv4DF500+4SGzOmODnl5dR0dHZwjqOnq4+4/PUd/hSrovo0Fm8qQi3HIEq6aw7UArx15iBWAzGbBbzdgsJlRNpakjfY6hk0oEn7xkP5LpvbTtY3gkWlcv4/W/GTHIDm791HksvWLWSREyqzN26OJDR0dHZxh6fUG+/tf/sG5fQ0Lbz5heSlTR2Ofqxh06akXJtdnItlgIqlF84TCeUJiodjTpmMVgYGZ+IaiCvbUd+EORpD/LsYwvEnxh2VrC7Expv/FQ6biFmsyvnlQhszpjhy4+dHR0dEYgoqr88Om3+Pfq4erJD01VWQ47o8lnUbUYDEzLL8DfG2J/UxfJTkp8dHGIs2b854SE3GabZ6KKEBOzPkWJY+mY71/n5CGe6/cJSAWoo6Ojc2IxKgr333gRVYU5/PyFd4j1FkzJMcLIRWVjIhiNsqGlr/ptRbWTYrODfXUdeAPxJe0yGwX33tCE2fE60RN0G5ltmc203G+cmJ3rfGjRxYeOjs5piSRJ3H7+XCrzs/jOY68QCA8/DVJc4MReZGFre+oLx9W7eqmnF5NdYXp1MUFPhH0NnSNbQySYVl1EVWkA2fp0SscT0TIJi1zsyqER2ymSjRl591GecVVK969zeqBPu+jo6Jz27Gnq4It/eo4219Bl6GfMLmVNa9OYjacsM5NSewZCA03ViERVwhGVYChCfo6DhoCHFo+nr63Tzs1nbqcweyUAqiqjBvMxmHuQDYMtKZqQ2NR7MxYlykT7q+z1LiXD6MWudLG1dxpvt/nINVv43PinMcpDp4x3mqYyp/AnOIyVafsOdD586D4fOjo6OnHS0evlyw+/wI76gdaNmTNKWd0xdsIjUS6dbCff4aF1dSFbNreSYTcz9wwTE2fuJrv8bWSlLwxXEwZ+svsmesIBMoxmPJGhQ48X5uVwdcnjKFJwwPLxztuZnPMlZEl3KtUZiC4+dHR0dBIgEI5w/xOv8trmfQBkZljozBgYtXIys8hcxP73BwulzAwz8+YYKZzbyVbhZF2Xm6A6ek6QWdk53FD2T4xyLyYlhzPyH6TAtjgdQ9c5BdDFh46Ojk6CaJrgt698wAf19fQoYQ50Jx/dMhaUmzKQNvgIjhDCm+W00nh5D2E59gq4kzKz+FpNI2cVfRuLQS+2pjM88Vy/9QwwOjo6OscgyxJfvGIxNy6ZRX1v74keTkwYkChskEcUHgCu3gBVkfjqqkzMLGNxyU904aGTUnTxoaOjozME10yfyqM3XUe29eQr9348S6Ri6htis9A4t1iRYpxFuqp8Jj+ee2166+LonJbovygdHR2dYZhfXsbTt32M6pzsEz2UYZltyWPHutiztR7Y18nivRNGbXdNxSx+NHeZLjx00oL+q9LR0dEZgcrsLJ6+7SYWVpaf6KEMIlsx498af1bTnRtaOMtVPez66ypn84M5V6NI+iVCJz3ovywdHR2dUXBaLDx8w7V8dOb0EzYGGSgy2ZlmyWGBpYhzjcVMbLXicidWTdcUVIZc/tGqOfzPGbrw0EkveoZTHR0dnRgwKgoPXnYR1bnZ/O+KdxnLMMECoxX7lgg+v5sOoCMFfcqBwTlULyudxvdmX4ksJVttRkdnZHRpq6OjoxMjkiTxqbPm8Ztrr8JiGJt7N0XAtGg2Pv/QycASJVMe7Eh7beVsXXjojAm6+NDR0dGJk6U1E3ji5o9S4LCnpL9JlmzOF8WcFylikbmIHIMZgAKjjXld2XE5lMZKpsGK03hUgMhInJFz8vm16Jya6OJDR0dHJwFmFBfxr9s+xuSC/KT6cShG2Opjx4ZGdm5pYv/7TRjX+DmzI4uMbREO1XamaMQDEZrgrPyjTqdTsopwGM1p2ZeOzvHo4kNHR0cnQYozM3jqlhs4f/zwkSOjMcWUM8hpVAiorevC4w0Os1XyqJpgUcG4/vfz8/QicTpjhy4+dHR0dJLAbjLx0PVX84l5ZyS0vcN/Yk7DmqaxMP+o+JiXq4sPnbEjrl/9Qw89xMyZM8nMzCQzM5OFCxfyyiuvDNn2M5/5DJIk8Ytf/CIV49TR0dE5aVFkmfsuOo/vLb0AJU6HTVeLN02jGhlNE1Q4ciizZWFRDMzJrTgh49A5PYnLXbusrIwf//jHTJjQlx3v0UcfZdmyZWzatIlp06b1t3vuuedYs2YNJSUlqR2tjo6OzknMzXNmUZHl5EvPv4Q3FB61/QRzFk3NXWMwssFoh2uK/nXJ7RRYMjDKQ+f90NFJB3FZPq666iouv/xyJk2axKRJk3jwwQdxOBysXr26v01TUxN33303f//73zEajSkfsI6Ojs7JzJJxVfzzlhspc45c1bPYaMOwI7EEYalA0/rER6ktSxceOmNOwpONqqry1FNP4fP5WLhwIdA3h3jrrbfyzW9+c4AlZCRCoRBut3vAQ0dHR+fDzMT8PP5128eYXVI8aF2B0cr5WjHyej/dLv8JGF0fmhZjdTkdnTQQt/jYtm0bDocDs9nMZz/7WZ599lmmTp0KwE9+8hMMBgNf+tKXYu7vRz/6EU6ns/9RXq7Hmevo6Hz4ybXbePxjH+GKKZMGLJ8WzWLHxkYiEfUEjawPWS8Yp3MCiTtFX01NDZs3b8blcvHMM89w++238/bbbxMIBPjlL3/Jxo0bkeJwuLrnnnv42te+1v/e7XbrAkRHR+eUwGI08LOrL6cqO5vffrAGhCDYM7ovSLopzMvgy588/0QPQ+c0RhJCJFWi4KKLLmL8+PFMmTKFr33tawPUtKqqyLJMeXk5tbW1MfXndrtxOp309vaSmTnynKmOjo7Oh4Xntu/kxTe3sfvd1GcrjYe8HAe//f5NlBVlndBx6Jx6xHP9Tro4gRCCUCjErbfeykUXXTRg3SWXXMKtt97KJz/5yWR3o6Ojo/Oh5prpU6nAwb0b2/H4UlunJVZysmz86rs36MJD54QTl/i49957ueyyyygvL8fj8fDUU0+xcuVKli9fTm5uLrm5uQPaG41GioqKqKmpSemgdXR0dD6MzJlewe8f/Dhff/AZWjvS71xvtRjJsFuw20xkOqx889MXUVmak/b96uiMRlzio62tjVtvvZWWlhacTiczZ85k+fLlXHzxxekan46Ojs4pRVVZLn/44cf51o+eZc/BtmHbWcwG7DYzGXYzdpsZh+3os8N+5L1pmOVm7FYTiqI7leqcnCTt85FqdJ8PHR2d0wF/IMzba/Zht5kGiwerCYNBz72h8+FiTH0+dHR0dHTix2Y1cdl5seVD0tE51dBtcjo6Ojo6Ojpjii4+dHR0dHR0dMYUXXzo6Ojo6OjojCm6+NDR0dHR0dEZU046h9MjwTd6gTkdHR0dHZ0PD0eu27EE0Z504sPj8QDo9V10dHR0dHQ+hHg8HpxO54htTro8H5qm0dzcTEZGRlwF6saCI0XvGhoa9BwkJxH6cTk50Y/LyYl+XE4+TpVjIoTA4/FQUlIyatXkk87yIcsyZWVlJ3oYI5KZmfmh/oGcqujH5eREPy4nJ/pxOfk4FY7JaBaPI+gOpzo6Ojo6Ojpjii4+dHR0dHR0dMYUXXzEgdls5rvf/S5ms/lED0XnGPTjcnKiH5eTE/24nHycjsfkpHM41dHR0dHR0Tm10S0fOjo6Ojo6OmOKLj50dHR0dHR0xhRdfOjo6Ojo6OiMKbr40NHR0dHR0RlTdPERI3v37mXZsmXk5eWRmZnJ4sWLWbFixZBtu7q6KCsrQ5IkXC7X2A70NGO047JlyxY+9rGPUV5ejtVqZcqUKfzyl788gSM+PYjl/6W+vp6rrroKu91OXl4eX/rSlwiHwydoxKc2K1euRJKkIR/r1q3rb7du3TouvPBCsrKyyM7OZunSpWzevPnEDfwUJ9bjAvDXv/6VmTNnYrFYKCoq4u677z5Bo04NuviIkSuuuIJoNMpbb73Fhg0bmD17NldeeSWtra2D2t55553MnDnzBIzy9GO047Jhwwby8/P529/+xo4dO7jvvvu45557+M1vfnOCR35qM9pxUVWVK664Ap/Px3vvvcdTTz3FM888w9e//vUTPPJTk0WLFtHS0jLgcdddd1FVVcW8efOAvnocl1xyCRUVFaxZs4b33nuPzMxMLrnkEiKRyAn+BKcmsRwXgJ/97Gfcd999fOc732HHjh28+eabXHLJJSdw5ClA6IxKR0eHAMQ777zTv8ztdgtAvPHGGwPa/u53vxPnnnuuePPNNwUgenp6xni0pw/xHJdj+fznPy/OP//8sRjiaUksx+Xll18WsiyLpqam/jZPPvmkMJvNore3d8zHfLoRDodFQUGB+P73v9+/bN26dQIQ9fX1/cu2bt0qALF///4TMczTjqGOS3d3t7BarSOe0z6M6JaPGMjNzWXKlCk89thj+Hw+otEof/jDHygsLGTu3Ln97Xbu3Mn3v/99HnvssVGL6ugkT6zH5Xh6e3vJyckZw5GeXsRyXFatWsX06dMpKSnp3+6SSy4hFAqxYcOGEzX004YXXniBzs5OPvGJT/Qvq6mpIS8vj4cffphwOEwgEODhhx9m2rRpVFZWnrjBnkYMdVxef/11NE2jqamJKVOmUFZWxg033EBDQ8OJG2gqONHq58NCY2OjmDt3rpAkSSiKIkpKSsSmTZv61weDQTFz5kzx+OOPCyGEWLFihW75GANGOy7H88EHHwij0Shee+21sRvkachox+VTn/qUuPjiiwdtZzKZxBNPPDGGIz09ueyyy8Rll102aPn27dvF+PHjhSzLQpZlMXnyZFFXV3cCRnh6MtRx+dGPfiSMRqOoqakRy5cvF6tWrRIXXnihqKmpEaFQ6ASNNHlO69vz733ve8M6+xx5rF+/HiEEn//85ykoKODdd99l7dq1LFu2jCuvvJKWlhYA7rnnHqZMmcItt9xygj/Vh59UHpdj2bFjB8uWLeP+++/n4osvPgGf7MNNqo+LJEmD9iGEGHK5ztDEekyOpbGxkVdffZU777xzwPJAIMAdd9zB4sWLWb16Ne+//z7Tpk3j8ssvJxAIjOXH+tCTyuOiaRqRSIRf/epXXHLJJSxYsIAnn3ySffv2DRv08GHgtE6v3tnZSWdn54htqqqqeP/991m6dCk9PT0Dyh1PnDiRO++8k+985zvMnj2bbdu29Z84hRBomoaiKNx333088MADaf0spxKpPC5H2LlzJ+effz533XUXDz74YNrGfiqTyuNy//338/zzz7Nly5b+9T09PeTk5PDWW29x/vnnp+1znErEekwsFkv/+//5n//h17/+NU1NTRiNxv7lDz/8MPfeey8tLS3908bhcJjs7GwefvhhbrrppvR8iFOQVB6XRx55hDvuuIOGhgbKysr6lxcWFvKDH/yAT33qU6n/AGOA4UQP4ESSl5dHXl7eqO38fj/AID8OWZbRNA2AZ555ZsDdwbp167jjjjt49913GT9+fApHfeqTyuMCfRaPCy64gNtvv10XHkmQyuOycOFCHnzwQVpaWiguLgbgtddew2w2j+ivozOQWI/JEYQQPPLII9x2220DLnDQd9xkWR5geTry/tj/J53RSeVxWbx4MQB79uzpFx/d3d10dnZ+uH1xTtiEz4eIjo4OkZubK6677jqxefNmsWfPHvGNb3xDGI1GsXnz5iG30X0+0k8sx2X79u0iPz9f3HzzzaKlpaX/0d7efoJHf+oSy3GJRqNi+vTp4sILLxQbN24Ub7zxhigrKxN33333CR79qc0bb7whALFz585B63bt2iXMZrP43Oc+J3bu3Cm2b98ubrnlFuF0OkVzc/MJGO3pw0jHRQghli1bJqZNmybef/99sW3bNnHllVeKqVOninA4PMYjTR26+IiRdevWiaVLl4qcnByRkZEhFixYIF5++eVh2+viY2wY7bh897vfFcCgR2Vl5Ykb9GlALP8vdXV14oorrhBWq1Xk5OSIu+++WwSDwRM04tODj33sY2LRokXDrn/ttdfE4sWLhdPpFNnZ2eKCCy4Qq1atGsMRnp6Mdlx6e3vFHXfcIbKyskROTo649tprB4REfxg5rX0+dHR0dHR0dMae0zraRUdHR0dHR2fs0cWHjo6Ojo6Ozpiiiw8dHR0dHR2dMUUXHzo6Ojo6Ojpjii4+dHR0dHR0dMYUXXzo6Ojo6OjojCm6+NDR0dHR0dEZU3TxoaOjo6OjozOm6OJDR0dHR0dHZ0zRxYeOjo6Ojo7OmKKLDx0dHR0dHZ0xRRcfOjo6Ojo6OmPK/wfnFkJUb7WCPAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "gdf = s.to_pandas().reset_index().set_geometry('destination')\n", "gdf.plot(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The issue is that we have lost the CRS along the way. CRS is stored on a GeometryArray level but when we create DataArray, xarray uses only the underlying numpy array and ignores our custom attribute. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gdf.crs is None" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What needs to be done\n", "\n", "It seems that the it will not be very complex to implement convenient vector data cubes. We need to figure out:\n", "\n", "- CRS handling\n", " - pyproj.CRS class storing the info on the GeometryArray level should be stored on a DataArray level as well (per dimension) and retained all the way back to a geopandas object.\n", " - It may require some custom logic during I/O and figuring out in which case we lose the info. Hopefully, most of it should be doable using xarray attrs. But some custom `to_geopandas` method may be needed as `to_pandas` still uses geometry as an index and that cannot contain CRS. So we need to convert to DataFrame/Series, reset index and create GeoDataFrame to which we pass a CRS stored on a DataArray level.\n", "- Figure out which geometric ops need to be available on an xarray level\n", " - I have to check stars docs to see how the actual data cube interacts with the vector data apart from using them as an index.\n", "- ... and probably more stuff I missed now." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wrap `shapely.STRtree` as an `xarray.indexes.Index`\n", "\n", "Create a new `ShapelySTRTreeIndex` for using with Xarray Dataset and DataArray objects (note: this works only with the last release of Xarray 2022.11.0):\n", "\n", "- It holds the CRS information, which must be provided explicitly. \n", "- It implements label-based data selection (`sel`) with two different modes:\n", " - \"nearest\" mode (default): selects the nearest geometry for each of the given input geometries\n", " - \"query\" mode: selects one or more geometries given a single input geometry and a predicate\n", "\n", "Next steps:\n", "\n", "- Figure out if we can leverage `shapely.STRtree`'s query capabilities for implementing custom alignment (spatial join) of Xarray Dataset / DataArray objects.\n", "- `ShapelySTRTreeIndex` only supports 1-dimensional geometry coordinates, but we could probably make it work with n-d coordinates.\n", "- It would be nice to automatically extract the CRS from the coordinate data. Unfortunately Xarray currently converts any geopandas.array.GeometryArray as a numpy array so we loose access to the CRS. I think this could be pretty easily fixed in Xarray." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "import shapely\n", "from xarray.indexes import Index\n", "from xarray.core.indexes import IndexSelResult\n", "\n", "\n", "class ShapelySTRTreeIndex(Index):\n", " \n", " def __init__(self, array, dim, crs):\n", " assert numpy.all(shapely.is_geometry(array))\n", " \n", " # only support 1-d coordinate for now\n", " assert len(array.shape) == 1\n", " \n", " self._tree = shapely.STRtree(numpy.ravel(array))\n", " self.dim = dim\n", " self.crs = crs\n", " \n", " @classmethod\n", " def from_variables(cls, variables, *, options):\n", " # only supports one coordinate of shapely geometries\n", " assert len(variables) == 1\n", " var = next(iter(variables.values()))\n", "\n", " return cls(var._data, var.dims[0], options[\"crs\"])\n", " \n", " def sel(self, labels, method=None, tolerance=None):\n", " # We reuse here `method` and `tolerance` options of\n", " # `xarray.indexes.PandasIndex` as `predicate` and `distance`\n", " # options when `labels` is a single geometry.\n", " # Xarray currently doesn't support custom options\n", " # (see https://github.com/pydata/xarray/issues/7099)\n", " \n", " # only one coordinate supported\n", " assert len(labels) == 1\n", " label = next(iter(labels.values()))\n", " \n", " if method is not None and method != \"nearest\":\n", " if not isinstance(label, shapely.Geometry):\n", " raise ValueError(\n", " \"selection with another method than nearest only supports \"\n", " \"a single geometry as input label.\"\n", " )\n", "\n", " if isinstance(label, xarray.DataArray):\n", " label_array = label._variable._data\n", " elif isinstance(label, xarray.Variable):\n", " label_array = label._data\n", " elif isinstance(label, shapely.Geometry):\n", " label_array = numpy.array([label])\n", " else:\n", " label_array = numpy.array(label)\n", " \n", " # check for possible CRS of geometry labels\n", " # (by default assume same CRS than the index)\n", " if hasattr(label_array, \"crs\") and label_array.crs != crs:\n", " raise ValueError(\"conflicting CRS for input geometries\")\n", " \n", " assert numpy.all(shapely.is_geometry(label_array))\n", " \n", " if method is None or method == \"nearest\":\n", " indices = self._tree.nearest(label_array)\n", " else:\n", " indices = self._tree.query(label, predicate=method, distance=tolerance)\n", "\n", " # attach dimension names and/or coordinates to positional indexer\n", " if isinstance(label, xarray.Variable):\n", " indices = xarray.Variable(label.dims, indices)\n", " elif isinstance(label, xarray.DataArray):\n", " indices = xarray.DataArray(indices, coords=label._coords, dims=label.dims)\n", "\n", " return IndexSelResult({self.dim: indices})\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Set two `ShapelySTRTreeIndex` instances for the `origin` and `destination` coordinates, respectively (first drop the two `pandas.Index` objects that were set by default for these two dimension coordinates)." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (mode: 3, day: 100, time: 24, origin: 100, destination: 100)>\n",
       "array([[[[[35, 56, 46, ...,  6, 79, 35],\n",
       "          [89,  5, 45, ..., 10, 23, 35],\n",
       "          [64, 19, 74, ..., 73, 85,  4],\n",
       "          ...,\n",
       "          [30, 23, 52, ...,  3, 54, 69],\n",
       "          [96, 63, 12, ..., 55, 93, 64],\n",
       "          [64,  6, 80, ..., 27, 32, 53]],\n",
       "\n",
       "         [[89, 11, 99, ..., 15, 38, 91],\n",
       "          [ 3, 97, 36, ..., 95, 98, 82],\n",
       "          [65, 71, 35, ...,  5, 80, 12],\n",
       "          ...,\n",
       "          [28, 33, 17, ..., 49, 19, 90],\n",
       "          [85, 37, 11, ..., 86, 91, 52],\n",
       "          [10, 32, 98, ...,  2, 22, 26]],\n",
       "\n",
       "         [[71, 20, 71, ..., 20, 94, 15],\n",
       "          [38, 67, 60, ..., 94, 85, 30],\n",
       "          [66, 79, 37, ..., 47, 80,  5],\n",
       "          ...,\n",
       "...\n",
       "          ...,\n",
       "          [18, 31, 36, ..., 89, 74, 25],\n",
       "          [33, 61, 48, ..., 66, 29, 65],\n",
       "          [94, 22, 18, ...,  8, 17,  3]],\n",
       "\n",
       "         [[10, 15, 91, ..., 54, 12, 28],\n",
       "          [37,  3, 21, ..., 11, 17, 53],\n",
       "          [48, 25, 50, ..., 77, 14, 16],\n",
       "          ...,\n",
       "          [35, 76, 83, ..., 15, 77, 61],\n",
       "          [ 1, 93,  4, ..., 71, 29,  6],\n",
       "          [40, 40, 24, ...,  7, 80, 56]],\n",
       "\n",
       "         [[38, 36, 12, ..., 22, 58, 60],\n",
       "          [71, 72, 88, ..., 22, 93, 56],\n",
       "          [68, 88, 61, ..., 46, 32, 54],\n",
       "          ...,\n",
       "          [57, 92, 77, ..., 96, 19,  3],\n",
       "          [77, 26, 76, ..., 45, 13, 72],\n",
       "          [ 4, 66, 93, ..., 64, 24, 65]]]]])\n",
       "Coordinates:\n",
       "  * mode         (mode) <U4 'car' 'bike' 'foot'\n",
       "  * day          (day) datetime64[ns] 2015-01-01 2015-01-02 ... 2015-04-10\n",
       "  * time         (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23\n",
       "  * origin       (origin) object MULTIPOLYGON (((-81.4727554321289 36.2343559...\n",
       "  * destination  (destination) object MULTIPOLYGON (((-81.4727554321289 36.23...\n",
       "Indexes:\n",
       "    origin       ShapelySTRTreeIndex\n",
       "    destination  ShapelySTRTreeIndex
" ], "text/plain": [ "\n", "array([[[[[35, 56, 46, ..., 6, 79, 35],\n", " [89, 5, 45, ..., 10, 23, 35],\n", " [64, 19, 74, ..., 73, 85, 4],\n", " ...,\n", " [30, 23, 52, ..., 3, 54, 69],\n", " [96, 63, 12, ..., 55, 93, 64],\n", " [64, 6, 80, ..., 27, 32, 53]],\n", "\n", " [[89, 11, 99, ..., 15, 38, 91],\n", " [ 3, 97, 36, ..., 95, 98, 82],\n", " [65, 71, 35, ..., 5, 80, 12],\n", " ...,\n", " [28, 33, 17, ..., 49, 19, 90],\n", " [85, 37, 11, ..., 86, 91, 52],\n", " [10, 32, 98, ..., 2, 22, 26]],\n", "\n", " [[71, 20, 71, ..., 20, 94, 15],\n", " [38, 67, 60, ..., 94, 85, 30],\n", " [66, 79, 37, ..., 47, 80, 5],\n", " ...,\n", "...\n", " ...,\n", " [18, 31, 36, ..., 89, 74, 25],\n", " [33, 61, 48, ..., 66, 29, 65],\n", " [94, 22, 18, ..., 8, 17, 3]],\n", "\n", " [[10, 15, 91, ..., 54, 12, 28],\n", " [37, 3, 21, ..., 11, 17, 53],\n", " [48, 25, 50, ..., 77, 14, 16],\n", " ...,\n", " [35, 76, 83, ..., 15, 77, 61],\n", " [ 1, 93, 4, ..., 71, 29, 6],\n", " [40, 40, 24, ..., 7, 80, 56]],\n", "\n", " [[38, 36, 12, ..., 22, 58, 60],\n", " [71, 72, 88, ..., 22, 93, 56],\n", " [68, 88, 61, ..., 46, 32, 54],\n", " ...,\n", " [57, 92, 77, ..., 96, 19, 3],\n", " [77, 26, 76, ..., 45, 13, 72],\n", " [ 4, 66, 93, ..., 64, 24, 65]]]]])\n", "Coordinates:\n", " * mode (mode) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (mode: 3, origin: 2, destination: 2)>\n",
       "array([[[43, 51],\n",
       "        [83, 38]],\n",
       "\n",
       "       [[37, 61],\n",
       "        [16, 22]],\n",
       "\n",
       "       [[87, 19],\n",
       "        [37, 36]]])\n",
       "Coordinates:\n",
       "  * mode         (mode) <U4 'car' 'bike' 'foot'\n",
       "    day          datetime64[ns] 2015-01-01\n",
       "    time         int64 0\n",
       "    origin       (origin) object MULTIPOLYGON (((-79.91995239257812 34.807918...\n",
       "    destination  (destination) object MULTIPOLYGON (((-80.72651672363281 35.5...
" ], "text/plain": [ "\n", "array([[[43, 51],\n", " [83, 38]],\n", "\n", " [[37, 61],\n", " [16, 22]],\n", "\n", " [[87, 19],\n", " [37, 36]]])\n", "Coordinates:\n", " * mode (mode) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (mode: 3, day: 100, time: 24, origin: 41, destination: 100)>\n",
       "array([[[[[88, 26, 19, ..., 19, 63, 33],\n",
       "          [77, 97, 29, ..., 36,  8, 63],\n",
       "          [21,  8, 10, ..., 53, 53, 89],\n",
       "          ...,\n",
       "          [76, 53, 77, ..., 42, 10, 18],\n",
       "          [55, 16, 93, ..., 78,  8, 37],\n",
       "          [84, 84, 88, ..., 84, 65,  4]],\n",
       "\n",
       "         [[99, 13, 41, ..., 20, 94, 92],\n",
       "          [67, 23, 68, ..., 95, 57, 69],\n",
       "          [28, 53, 99, ..., 22,  4, 81],\n",
       "          ...,\n",
       "          [94, 48, 16, ..., 22, 16, 92],\n",
       "          [57, 60, 49, ..., 12, 48, 24],\n",
       "          [20,  7, 88, ..., 89, 21,  1]],\n",
       "\n",
       "         [[91, 53, 73, ..., 66, 59, 25],\n",
       "          [52, 96, 39, ..., 14, 97, 48],\n",
       "          [ 3, 90, 26, ..., 39, 20, 56],\n",
       "          ...,\n",
       "...\n",
       "          ...,\n",
       "          [81, 12, 92, ..., 88, 64, 90],\n",
       "          [97, 77, 62, ..., 35, 72, 48],\n",
       "          [69, 19, 94, ..., 44, 23, 43]],\n",
       "\n",
       "         [[89, 58, 43, ..., 94, 28, 94],\n",
       "          [67, 78, 14, ..., 29, 24, 37],\n",
       "          [ 1, 33, 15, ..., 96, 76, 40],\n",
       "          ...,\n",
       "          [20, 81, 87, ..., 27, 62, 51],\n",
       "          [29, 11, 93, ..., 41, 34, 76],\n",
       "          [20, 24, 18, ..., 56, 11, 52]],\n",
       "\n",
       "         [[62, 50, 43, ..., 92,  5, 23],\n",
       "          [74, 22, 77, ..., 12, 52, 86],\n",
       "          [96,  4, 47, ..., 87, 74, 25],\n",
       "          ...,\n",
       "          [19, 26, 15, ..., 93, 84, 94],\n",
       "          [29, 93,  7, ..., 41, 89, 47],\n",
       "          [43, 52, 79, ..., 18, 37, 35]]]]])\n",
       "Coordinates:\n",
       "  * mode         (mode) <U4 'car' 'bike' 'foot'\n",
       "  * day          (day) datetime64[ns] 2015-01-01 2015-01-02 ... 2015-04-10\n",
       "  * time         (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23\n",
       "    origin       (origin) object MULTIPOLYGON (((-78.11376953125 34.720985412...\n",
       "  * destination  (destination) object MULTIPOLYGON (((-81.4727554321289 36.23...\n",
       "Indexes:\n",
       "    destination  ShapelySTRTreeIndex
" ], "text/plain": [ "\n", "array([[[[[88, 26, 19, ..., 19, 63, 33],\n", " [77, 97, 29, ..., 36, 8, 63],\n", " [21, 8, 10, ..., 53, 53, 89],\n", " ...,\n", " [76, 53, 77, ..., 42, 10, 18],\n", " [55, 16, 93, ..., 78, 8, 37],\n", " [84, 84, 88, ..., 84, 65, 4]],\n", "\n", " [[99, 13, 41, ..., 20, 94, 92],\n", " [67, 23, 68, ..., 95, 57, 69],\n", " [28, 53, 99, ..., 22, 4, 81],\n", " ...,\n", " [94, 48, 16, ..., 22, 16, 92],\n", " [57, 60, 49, ..., 12, 48, 24],\n", " [20, 7, 88, ..., 89, 21, 1]],\n", "\n", " [[91, 53, 73, ..., 66, 59, 25],\n", " [52, 96, 39, ..., 14, 97, 48],\n", " [ 3, 90, 26, ..., 39, 20, 56],\n", " ...,\n", "...\n", " ...,\n", " [81, 12, 92, ..., 88, 64, 90],\n", " [97, 77, 62, ..., 35, 72, 48],\n", " [69, 19, 94, ..., 44, 23, 43]],\n", "\n", " [[89, 58, 43, ..., 94, 28, 94],\n", " [67, 78, 14, ..., 29, 24, 37],\n", " [ 1, 33, 15, ..., 96, 76, 40],\n", " ...,\n", " [20, 81, 87, ..., 27, 62, 51],\n", " [29, 11, 93, ..., 41, 34, 76],\n", " [20, 24, 18, ..., 56, 11, 52]],\n", "\n", " [[62, 50, 43, ..., 92, 5, 23],\n", " [74, 22, 77, ..., 12, 52, 86],\n", " [96, 4, 47, ..., 87, 74, 25],\n", " ...,\n", " [19, 26, 15, ..., 93, 84, 94],\n", " [29, 93, 7, ..., 41, 89, 47],\n", " [43, 52, 79, ..., 18, 37, 35]]]]])\n", "Coordinates:\n", " * mode (mode) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (mode: 3, day: 100, time: 24, points: 2, destination: 100)>\n",
       "array([[[[[13, 61,  7, ..., 45, 83, 86],\n",
       "          [21, 56, 21, ..., 10, 86, 94]],\n",
       "\n",
       "         [[79, 71, 90, ..., 61, 14, 64],\n",
       "          [95, 63, 57, ..., 78,  3, 59]],\n",
       "\n",
       "         [[24,  1, 15, ..., 83, 43, 97],\n",
       "          [48, 21, 84, ..., 47, 41, 42]],\n",
       "\n",
       "         ...,\n",
       "\n",
       "         [[64, 84,  4, ..., 26, 97, 49],\n",
       "          [32, 25, 17, ..., 96, 16, 74]],\n",
       "\n",
       "         [[66, 51, 55, ..., 89, 40, 48],\n",
       "          [17, 51, 19, ..., 14, 94, 58]],\n",
       "\n",
       "         [[83, 97, 86, ..., 64, 74, 88],\n",
       "          [29, 63, 63, ..., 85, 86, 30]]],\n",
       "\n",
       "...\n",
       "\n",
       "        [[[32, 43, 89, ..., 99, 21, 52],\n",
       "          [76, 18, 26, ..., 10,  7,  8]],\n",
       "\n",
       "         [[56, 48, 34, ..., 22, 64, 84],\n",
       "          [62, 31, 59, ..., 75,  1, 42]],\n",
       "\n",
       "         [[76, 37, 32, ..., 50, 77, 57],\n",
       "          [11,  3, 85, ..., 86, 19,  3]],\n",
       "\n",
       "         ...,\n",
       "\n",
       "         [[76, 15, 85, ..., 25, 69, 56],\n",
       "          [33, 39, 23, ..., 77, 74,  2]],\n",
       "\n",
       "         [[26, 82, 76, ..., 96, 17, 22],\n",
       "          [87, 55, 94, ..., 96, 19, 28]],\n",
       "\n",
       "         [[59, 41,  5, ..., 29, 53, 63],\n",
       "          [50, 38, 17, ..., 99, 43, 54]]]]])\n",
       "Coordinates:\n",
       "  * mode         (mode) <U4 'car' 'bike' 'foot'\n",
       "  * day          (day) datetime64[ns] 2015-01-01 2015-01-02 ... 2015-04-10\n",
       "  * time         (time) int64 0 1 2 3 4 5 6 7 8 9 ... 15 16 17 18 19 20 21 22 23\n",
       "    origin       (points) object MULTIPOLYGON (((-79.91995239257812 34.807918...\n",
       "  * destination  (destination) object MULTIPOLYGON (((-81.4727554321289 36.23...\n",
       "Dimensions without coordinates: points\n",
       "Indexes:\n",
       "    destination  ShapelySTRTreeIndex
" ], "text/plain": [ "\n", "array([[[[[13, 61, 7, ..., 45, 83, 86],\n", " [21, 56, 21, ..., 10, 86, 94]],\n", "\n", " [[79, 71, 90, ..., 61, 14, 64],\n", " [95, 63, 57, ..., 78, 3, 59]],\n", "\n", " [[24, 1, 15, ..., 83, 43, 97],\n", " [48, 21, 84, ..., 47, 41, 42]],\n", "\n", " ...,\n", "\n", " [[64, 84, 4, ..., 26, 97, 49],\n", " [32, 25, 17, ..., 96, 16, 74]],\n", "\n", " [[66, 51, 55, ..., 89, 40, 48],\n", " [17, 51, 19, ..., 14, 94, 58]],\n", "\n", " [[83, 97, 86, ..., 64, 74, 88],\n", " [29, 63, 63, ..., 85, 86, 30]]],\n", "\n", "...\n", "\n", " [[[32, 43, 89, ..., 99, 21, 52],\n", " [76, 18, 26, ..., 10, 7, 8]],\n", "\n", " [[56, 48, 34, ..., 22, 64, 84],\n", " [62, 31, 59, ..., 75, 1, 42]],\n", "\n", " [[76, 37, 32, ..., 50, 77, 57],\n", " [11, 3, 85, ..., 86, 19, 3]],\n", "\n", " ...,\n", "\n", " [[76, 15, 85, ..., 25, 69, 56],\n", " [33, 39, 23, ..., 77, 74, 2]],\n", "\n", " [[26, 82, 76, ..., 96, 17, 22],\n", " [87, 55, 94, ..., 96, 19, 28]],\n", "\n", " [[59, 41, 5, ..., 29, 53, 63],\n", " [50, 38, 17, ..., 99, 43, 54]]]]])\n", "Coordinates:\n", " * mode (mode) 1\u001b[0m \u001b[43marr\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msel\u001b[49m\u001b[43m(\u001b[49m\u001b[43morigin\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m[\u001b[49m\u001b[43mshapely\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mPoint\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m80\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m35\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mshapely\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mPoint\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[38;5;241;43m76\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m35.6\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mintersects\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n", "File \u001b[0;32m~/miniconda3/envs/xarray-vector/lib/python3.11/site-packages/xarray/core/dataarray.py:1526\u001b[0m, in \u001b[0;36mDataArray.sel\u001b[0;34m(self, indexers, method, tolerance, drop, **indexers_kwargs)\u001b[0m\n\u001b[1;32m 1416\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21msel\u001b[39m(\n\u001b[1;32m 1417\u001b[0m \u001b[38;5;28mself\u001b[39m: T_DataArray,\n\u001b[1;32m 1418\u001b[0m indexers: Mapping[Any, Any] \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m,\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1422\u001b[0m \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mindexers_kwargs: Any,\n\u001b[1;32m 1423\u001b[0m ) \u001b[38;5;241m-\u001b[39m\u001b[38;5;241m>\u001b[39m T_DataArray:\n\u001b[1;32m 1424\u001b[0m \u001b[38;5;124;03m\"\"\"Return a new DataArray whose data is given by selecting index\u001b[39;00m\n\u001b[1;32m 1425\u001b[0m \u001b[38;5;124;03m labels along the specified dimension(s).\u001b[39;00m\n\u001b[1;32m 1426\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1524\u001b[0m \u001b[38;5;124;03m Dimensions without coordinates: points\u001b[39;00m\n\u001b[1;32m 1525\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1526\u001b[0m ds \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_to_temp_dataset\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msel\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 1527\u001b[0m \u001b[43m \u001b[49m\u001b[43mindexers\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mindexers\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 1528\u001b[0m \u001b[43m \u001b[49m\u001b[43mdrop\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mdrop\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 1529\u001b[0m \u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 1530\u001b[0m \u001b[43m \u001b[49m\u001b[43mtolerance\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtolerance\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 1531\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mindexers_kwargs\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m 1532\u001b[0m \u001b[43m \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1533\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_from_temp_dataset(ds)\n", "File \u001b[0;32m~/miniconda3/envs/xarray-vector/lib/python3.11/site-packages/xarray/core/dataset.py:2554\u001b[0m, in \u001b[0;36mDataset.sel\u001b[0;34m(self, indexers, method, tolerance, drop, **indexers_kwargs)\u001b[0m\n\u001b[1;32m 2493\u001b[0m \u001b[38;5;124;03m\"\"\"Returns a new dataset with each array indexed by tick labels\u001b[39;00m\n\u001b[1;32m 2494\u001b[0m \u001b[38;5;124;03malong the specified dimension(s).\u001b[39;00m\n\u001b[1;32m 2495\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 2551\u001b[0m \u001b[38;5;124;03mDataArray.sel\u001b[39;00m\n\u001b[1;32m 2552\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 2553\u001b[0m indexers \u001b[38;5;241m=\u001b[39m either_dict_or_kwargs(indexers, indexers_kwargs, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msel\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m-> 2554\u001b[0m query_results \u001b[38;5;241m=\u001b[39m \u001b[43mmap_index_queries\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m 2555\u001b[0m \u001b[43m \u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mindexers\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mindexers\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmethod\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mmethod\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mtolerance\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mtolerance\u001b[49m\n\u001b[1;32m 2556\u001b[0m \u001b[43m\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 2558\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m drop:\n\u001b[1;32m 2559\u001b[0m no_scalar_variables \u001b[38;5;241m=\u001b[39m {}\n", "File \u001b[0;32m~/miniconda3/envs/xarray-vector/lib/python3.11/site-packages/xarray/core/indexing.py:183\u001b[0m, in \u001b[0;36mmap_index_queries\u001b[0;34m(obj, indexers, method, tolerance, **indexers_kwargs)\u001b[0m\n\u001b[1;32m 181\u001b[0m results\u001b[38;5;241m.\u001b[39mappend(IndexSelResult(labels))\n\u001b[1;32m 182\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 183\u001b[0m results\u001b[38;5;241m.\u001b[39mappend(\u001b[43mindex\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43msel\u001b[49m\u001b[43m(\u001b[49m\u001b[43mlabels\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43moptions\u001b[49m\u001b[43m)\u001b[49m) \u001b[38;5;66;03m# type: ignore[call-arg]\u001b[39;00m\n\u001b[1;32m 185\u001b[0m merged \u001b[38;5;241m=\u001b[39m merge_sel_results(results)\n\u001b[1;32m 187\u001b[0m \u001b[38;5;66;03m# drop dimension coordinates found in dimension indexers\u001b[39;00m\n\u001b[1;32m 188\u001b[0m \u001b[38;5;66;03m# (also drop multi-index if any)\u001b[39;00m\n\u001b[1;32m 189\u001b[0m \u001b[38;5;66;03m# (.sel() already ensures alignment)\u001b[39;00m\n", "Cell \u001b[0;32mIn [10], line 39\u001b[0m, in \u001b[0;36mShapelySTRTreeIndex.sel\u001b[0;34m(self, labels, method, tolerance)\u001b[0m\n\u001b[1;32m 37\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m method \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m \u001b[38;5;129;01mand\u001b[39;00m method \u001b[38;5;241m!=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mnearest\u001b[39m\u001b[38;5;124m\"\u001b[39m:\n\u001b[1;32m 38\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(label, shapely\u001b[38;5;241m.\u001b[39mGeometry):\n\u001b[0;32m---> 39\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m 40\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mselection with another method than nearest only supports \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 41\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124ma single geometry as input label.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 42\u001b[0m )\n\u001b[1;32m 44\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(label, xarray\u001b[38;5;241m.\u001b[39mDataArray):\n\u001b[1;32m 45\u001b[0m label_array \u001b[38;5;241m=\u001b[39m label\u001b[38;5;241m.\u001b[39m_variable\u001b[38;5;241m.\u001b[39m_data\n", "\u001b[0;31mValueError\u001b[0m: selection with another method than nearest only supports a single geometry as input label." ] } ], "source": [ "arr.sel(origin=[shapely.Point(-80, 35), shapely.Point(-76, 35.6)], method=\"intersects\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python [conda env:xarray-vector]", "language": "python", "name": "conda-env-xarray-vector-py" }, "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" }, "vscode": { "interpreter": { "hash": "716119b7942b628d25d105f893f626cc168d483792084e8427e38a31c13f35ec" } } }, "nbformat": 4, "nbformat_minor": 4 }