diff --git a/doc/api.rst b/doc/api.rst
index 97fa1d7..9e5306e 100644
--- a/doc/api.rst
+++ b/doc/api.rst
@@ -36,6 +36,7 @@ properties listed below. Proper use of this accessor should be like:
Dataset.xoak.set_index
Dataset.xoak.sel
+ Dataset.xoak.query
DataArray.xoak
--------------
@@ -58,6 +59,7 @@ The accessor above is also registered for :py:class:`xarray.DataArray`.
DataArray.xoak.set_index
DataArray.xoak.sel
+ DataArray.xoak.query
Indexes
-------
diff --git a/doc/examples/index.rst b/doc/examples/index.rst
index bc48c62..a4035f7 100644
--- a/doc/examples/index.rst
+++ b/doc/examples/index.rst
@@ -7,3 +7,4 @@ Examples
introduction
dask_support
custom_indexes
+ query
diff --git a/doc/examples/query.ipynb b/doc/examples/query.ipynb
new file mode 100644
index 0000000..fe728c1
--- /dev/null
+++ b/doc/examples/query.ipynb
@@ -0,0 +1,382 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Querying for indices\n",
+ "\n",
+ "This notebook briefly shows how to use Xoak's `query()` method to build more complicated queries."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import numpy as np\n",
+ "import xarray as xr\n",
+ "import xoak\n",
+ "\n",
+ "xr.set_options(display_style='text');"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Load example dataset"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<xarray.Dataset>\n",
+ "Dimensions: (lat: 25, lon: 53, time: 2920)\n",
+ "Coordinates:\n",
+ " * lat (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0\n",
+ " * lon (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0\n",
+ " * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
+ " lat2d (lat, lon) float32 75.0 75.0 75.0 75.0 75.0 ... 15.0 15.0 15.0 15.0\n",
+ " lon2d (lat, lon) float32 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0\n",
+ "Data variables:\n",
+ " air (time, lat, lon) float32 241.2 242.5 243.5 ... 296.5 296.2 295.7\n",
+ "Attributes:\n",
+ " Conventions: COARDS\n",
+ " title: 4x daily NMC reanalysis (1948)\n",
+ " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n",
+ " platform: Model\n",
+ " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...
"
+ ],
+ "text/plain": [
+ "\n",
+ "Dimensions: (lat: 25, lon: 53, time: 2920)\n",
+ "Coordinates:\n",
+ " * lat (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0\n",
+ " * lon (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0\n",
+ " * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
+ " lat2d (lat, lon) float32 75.0 75.0 75.0 75.0 75.0 ... 15.0 15.0 15.0 15.0\n",
+ " lon2d (lat, lon) float32 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0\n",
+ "Data variables:\n",
+ " air (time, lat, lon) float32 241.2 242.5 243.5 ... 296.5 296.2 295.7\n",
+ "Attributes:\n",
+ " Conventions: COARDS\n",
+ " title: 4x daily NMC reanalysis (1948)\n",
+ " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n",
+ " platform: Model\n",
+ " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly..."
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "ds = xr.tutorial.load_dataset(\"air_temperature\")\n",
+ "\n",
+ "# create artificial 2d coords\n",
+ "ds.coords[\"lat2d\"], ds.coords[\"lon2d\"] = xr.broadcast(\n",
+ " ds.coords[\"lat\"], ds.coords[\"lon\"]\n",
+ ")\n",
+ "\n",
+ "display(ds)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "ds.xoak.set_index(['lat2d', 'lon2d'], 'sklearn_geo_balltree')"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Create trajectory dataset"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<xarray.Dataset>\n",
+ "Dimensions: (along_track: 23)\n",
+ "Coordinates:\n",
+ " * along_track (along_track) int64 0 1 2 3 4 5 6 7 ... 15 16 17 18 19 20 21 22\n",
+ "Data variables:\n",
+ " latitude (along_track) float64 75.0 72.27 69.55 ... 20.45 17.73 15.0\n",
+ " longitude (along_track) float64 200.0 205.9 211.8 ... 318.2 324.1 330.0
"
+ ],
+ "text/plain": [
+ "\n",
+ "Dimensions: (along_track: 23)\n",
+ "Coordinates:\n",
+ " * along_track (along_track) int64 0 1 2 3 4 5 6 7 ... 15 16 17 18 19 20 21 22\n",
+ "Data variables:\n",
+ " latitude (along_track) float64 75.0 72.27 69.55 ... 20.45 17.73 15.0\n",
+ " longitude (along_track) float64 200.0 205.9 211.8 ... 318.2 324.1 330.0"
+ ]
+ },
+ "execution_count": 4,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "ds_trajectory = xr.Dataset(\n",
+ " {\n",
+ " 'latitude': ('along_track', np.linspace(75.0, 15.0, 23)),\n",
+ " 'longitude': ('along_track', np.linspace(200, 330, 23)),\n",
+ " },\n",
+ " coords={\"along_track\": (\"along_track\", np.arange(23))},\n",
+ ")\n",
+ "\n",
+ "ds_trajectory"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Select and check for distances to grid points"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<xarray.Dataset>\n",
+ "Dimensions: (along_track: 23, time: 2920)\n",
+ "Coordinates:\n",
+ " lat (along_track) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0\n",
+ " lon (along_track) float32 200.0 205.0 212.5 217.5 ... 317.5 325.0 330.0\n",
+ " * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
+ " lat2d (along_track) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0\n",
+ " lon2d (along_track) float32 200.0 205.0 212.5 217.5 ... 317.5 325.0 330.0\n",
+ "Dimensions without coordinates: along_track\n",
+ "Data variables:\n",
+ " air (time, along_track) float32 241.2 244.7 244.4 ... 296.7 295.7 295.7\n",
+ "Attributes:\n",
+ " Conventions: COARDS\n",
+ " title: 4x daily NMC reanalysis (1948)\n",
+ " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n",
+ " platform: Model\n",
+ " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...
"
+ ],
+ "text/plain": [
+ "\n",
+ "Dimensions: (along_track: 23, time: 2920)\n",
+ "Coordinates:\n",
+ " lat (along_track) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0\n",
+ " lon (along_track) float32 200.0 205.0 212.5 217.5 ... 317.5 325.0 330.0\n",
+ " * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
+ " lat2d (along_track) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0\n",
+ " lon2d (along_track) float32 200.0 205.0 212.5 217.5 ... 317.5 325.0 330.0\n",
+ "Dimensions without coordinates: along_track\n",
+ "Data variables:\n",
+ " air (time, along_track) float32 241.2 244.7 244.4 ... 296.7 295.7 295.7\n",
+ "Attributes:\n",
+ " Conventions: COARDS\n",
+ " title: 4x daily NMC reanalysis (1948)\n",
+ " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n",
+ " platform: Model\n",
+ " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly..."
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "# ds.xoak.sel() does the following two steps\n",
+ "\n",
+ "indices = ds.xoak.query(\n",
+ " lat2d=ds_trajectory.latitude,\n",
+ " lon2d=ds_trajectory.longitude,\n",
+ ")\n",
+ "selected = ds.isel(indices)\n",
+ "\n",
+ "display(selected)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "selected.coords[\"along_track\"] = ds_trajectory.coords[\"along_track\"]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## We can now easily calculated the mismatch"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "selected[\"diff_lat\"] = ds_trajectory[\"latitude\"] - selected[\"lat\"]\n",
+ "selected[\"diff_lon\"] = ds_trajectory[\"longitude\"] - selected[\"lon\"]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "<xarray.Dataset>\n",
+ "Dimensions: (along_track: 23, time: 2920)\n",
+ "Coordinates:\n",
+ " lat (along_track) float32 75.0 72.5 70.0 67.5 ... 20.0 17.5 15.0\n",
+ " lon (along_track) float32 200.0 205.0 212.5 ... 317.5 325.0 330.0\n",
+ " * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
+ " lat2d (along_track) float32 75.0 72.5 70.0 67.5 ... 20.0 17.5 15.0\n",
+ " lon2d (along_track) float32 200.0 205.0 212.5 ... 317.5 325.0 330.0\n",
+ " * along_track (along_track) int64 0 1 2 3 4 5 6 7 ... 15 16 17 18 19 20 21 22\n",
+ "Data variables:\n",
+ " air (time, along_track) float32 241.2 244.7 244.4 ... 295.7 295.7\n",
+ " diff_lat (along_track) float64 0.0 -0.2273 -0.4545 ... 0.4545 0.2273 0.0\n",
+ " diff_lon (along_track) float64 0.0 0.9091 -0.6818 ... 0.6818 -0.9091 0.0\n",
+ "Attributes:\n",
+ " Conventions: COARDS\n",
+ " title: 4x daily NMC reanalysis (1948)\n",
+ " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n",
+ " platform: Model\n",
+ " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...
"
+ ],
+ "text/plain": [
+ "\n",
+ "Dimensions: (along_track: 23, time: 2920)\n",
+ "Coordinates:\n",
+ " lat (along_track) float32 75.0 72.5 70.0 67.5 ... 20.0 17.5 15.0\n",
+ " lon (along_track) float32 200.0 205.0 212.5 ... 317.5 325.0 330.0\n",
+ " * time (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
+ " lat2d (along_track) float32 75.0 72.5 70.0 67.5 ... 20.0 17.5 15.0\n",
+ " lon2d (along_track) float32 200.0 205.0 212.5 ... 317.5 325.0 330.0\n",
+ " * along_track (along_track) int64 0 1 2 3 4 5 6 7 ... 15 16 17 18 19 20 21 22\n",
+ "Data variables:\n",
+ " air (time, along_track) float32 241.2 244.7 244.4 ... 295.7 295.7\n",
+ " diff_lat (along_track) float64 0.0 -0.2273 -0.4545 ... 0.4545 0.2273 0.0\n",
+ " diff_lon (along_track) float64 0.0 0.9091 -0.6818 ... 0.6818 -0.9091 0.0\n",
+ "Attributes:\n",
+ " Conventions: COARDS\n",
+ " title: 4x daily NMC reanalysis (1948)\n",
+ " description: Data is from NMC initialized reanalysis\\n(4x/day). These a...\n",
+ " platform: Model\n",
+ " references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly..."
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "display(selected)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Visualize"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from matplotlib import pyplot as plt"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjQAAADQCAYAAAAQwfu+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABH5ElEQVR4nO3deXhcZfn/8fdn0iTdW0jLWkrZZC1LKUtYA0VAQAoWEEVZtaCgoILCF0UUEZSfiooKVcEii6AVWQRZAsNmWFpAyiqLRQoVaGhLS2kmmbl/f5wn7TSd5WSdmeR+Xde5Zs5+z8mZk2eeVWaGc84551wlS5Q6AOecc8657vIEjXPOOecqnidonHPOOVfxPEHjnHPOuYrnCRrnnHPOVTxP0DjnnHOu4nmCZgCTNE/SAeV+zC7E0CBpfoH1JmnzPo7pSknf6ctzOheXpD9I+oGkvSW9nLV8S0lPS1oq6auShki6XdISSX8uZczOdeQJmiyl/mdc6vOXi/aHa6nj6CpJJ0p6JHuZmZ1mZheVKibn4jCzh81sy6xF3wSSZjbCzH4BHAWsC9SZ2dElCdK5PDxB04MkVZU6Buec60EbA893mP+3mbV19kCSBvVYVM7l4AmaQNIfgfHA7ZKWSfpmWP5nSf8LWawPSdo2a58/SPqNpDslfQjsJ2lSVhbtnyXdlJ3bIOkwSc9IWizpn5K2L3T+DjGuJekOSe9JWhTej8tan5R0kaRHw/nvkTQma/3nJb0hqVnS+UWuxyGSXgjHeUvS2cU+Q45jJCSdK+m1cM6bJa2dtX6vsP9iSW+GnI3pwHHAN8N1uD1su4GkWeGz/0fSV7OOMyT8LRZJegHYpdBn6xDjKEnXhuO+IenbkhJZ678o6cVwHV6QNCksb/9c7cuPDMu3Bq4E6kP8i8Py1XKdwnFflfS+pNskbZC1ziSdJumV8Jl+JUlh3eaSHgz340JJN8X9rM61k7STpKfC/XsTMDgsX1lcK+l+YD/ginAv3whcAHw6zJ8Stjs5fEcWSbpb0sZZ5zFJp0t6BXglLMv7/FCUS322pGfDPX6TpMFZ66eGfT8I37+Dw/JRkn4vaUF4Xv1A4Qemf2cGEDPzKUzAPOCADstOBkYAtcDlwDNZ6/4ALAH2JEocjgTeAM4EqoFPASngB2H7ScC7wG5AFXBCOGdtvvN3iKUOmAYMDTH9Gfhb1vok8BrwMWBImL80rNsGWAbsEz7LT4G2fOcDFgB7h/drAZM6+xmAs4DHgHHhnFcBN4Z144GlwGfCtaoDdsy6rj/IiiUBzCF6mNYAmwKvAweF9ZcCDwNrAxsBzwHzC1xHAzYP768Fbg3XcwLwb+CUsO5o4C2iBJKAzYGNs9ZtEGL7NPAhsH5YdyLwSIdzrvxMwP7AwnAta4FfAg91iO8OYHS4Tu8BB4d1NwLnh/MOBvYq9ffGp8qawnfoDeBr4bt3FNAK/ABoyP7uED1DvpA1fyFwXdb8EcCrwNbAIODbwD+z1htwb/huDiHe8+OJ8N1aG3gROC2s25XoefvxcP9vCGwV1v2N6PkyDFgnHOPUsM6/MwNkKnkA5TRRPEExOnxBR4X5PwDXZq3fJ/wDVNayR7L+kf0GuKjDMV8G9o1z/hzx7AgsyppPAt/Omv8y8I/w/gLgT1nrhhEltvIlaP4LnAqM7LA89mcID6MpWdutT/TgHAScB9yS59x/YPUEzW7Afztscx5wTXj/OuEffpifTowETXigtgDbZK07lajOAMDdwJkx/xbPAFPD+xMpnKD5PfDjrHXDw3WZkBXfXlnrbwbODe+vBWYA40r1PfGpsqfwnHq7w3Pqn3QtQXMX4QdAmE8Ay1mV8Ddg/6z1cZ4fn8ta92PgyvD+KuBnOT7PuuF7PCRr2WeAB8J7/84MkMmLnAqQVCXp0pC1+QHRlw1gTNZmb2a93wB4y8K3KMf6jYFvhKzWxaE4YqOwX5x4hkq6KhSNfAA8BIzW6nV3/pf1fjnRP8v22FbGYmYfAs0FTjcNOAR4I2TX1nfhM2wM3JK13YtAmugBtBFRblIcGwMbdDjn/4XjrPHZiH59xjGGVb9Ws/fdMLzPG6Ok47OyzRcD27H6fVHIBtnnNLNlRH+LDbO2yfd3/CZRbtETkp6XdHLMczrXLtdzKu53pqONgZ9nfQ/eJ7o/s+/lzj4D8937+b6PGxPlNC3IOuZVRDk14N+ZAcMraa2u49DjnwWmAgcQJWZGAYuIvhy59lkAbChJWQ+L7C/hm8DFZnZxzPN39A1gS2A3M/ufpB2BpzvEk88ComxhIEocERXz5A7E7ElgqqRq4AyiXIKNYnyGbG8CJ5vZox1XSHqTKAs55+lzHOc/ZrZFnu0XhNjaKy+OjxEbRMU+rUQPxBey9n0r67ybddwp1BH4LTAFaDKztKRnWPV3KPZ3fDucs/14w4j+Fm/l3aP9wGb/A74Y9tsLuE/SQ2b2arF9nQtyPafGE/8HRrb258H1Bbbp+AMv7vMj17nW+D6G5S3AGMtRWdm/MwOH59Cs7h2i+hntRhB9UZqJ6q38sMj+TUQ5EGdIGiRpKqv/0/4tcJqk3RQZJulQSSPynL+jEcBHwGJFlWu/G/eDAX8BDlNUEbcG+D55/v6SaiQdJ2mUmbUCH4TPFeczZLsSuLi9kqCkseGaAFwPHCDpmHCt6kICDda8Dk8AH0j6lqIKwFWStpPUXvn3ZuA8RZWmxwFfiXNBzCwd9r1Y0ogQ59eB68ImvwPOlrRz+Kybh22GET2k3wuf6ySiHJp27wDjwnXO5QbgJEk7Sqoluq8eN7N5xWKWdLRWVQRfFOJIF9jFuY6aiOrPfTV89z5F/h8XxVxJ9N3bFlZWzi3UnLszz4+Ofk/0vZmiqMHBhpK2MrMFwD3ATySNDOs2k7RviMm/MwOEJ2hWdwnw7ZBteTZR2esbRL+cXyCq4JqXmaWIKgKfAiwGPkdUubMlrJ9N9EvhCqIv1qtE9S3ynb+jy4kq1i0Msfwj7gczs+eB04n+mS4I58/b+RzweWBeKNo6LXyWOJ8h28+B24B7JC0NMe8WjvNfoiKtbxBlUz8D7BD2+z2wTbgOfwsJj08S1Rn6T/j8vyPKMQP4HtHf6T9ED7Y/Fr8iK32FqELv60T1nW4Arg4x/hm4OCxbSlTxcG0zewH4CdE/hneAiUB2LtT9RLlF/5O0sOMJzawR+A4wi+hvsRlwbMx4dwEel7SM6NqeaWb/if9x3UCX9Zw6keg7/Gngr1081i3Aj4A/hWfFc8AnCmzfmedHx32fAE4CfkZUOfhBVuV0Hk9UfPxCOO5fiOrsgX9nBgytXozqepqkx4kqtV1T6licc865/spzaHqYpH0lrReyck8AtqcTOSnOOeec6zyvFNzztiSqlzGcqJLdUaGM1znnnHO9xIucnHPOOVfxvMjJOeeccxWvYoqcxowZYxMmTCh1GM4NWHPmzFloZmNLHUdv8ueM60868509aL9h1vx+/tbsc55tudvMDu6x4HpBxSRoJkyYwOzZs0sdhnMDlqSu9iZbMfw54/qTznxnF76f5vG7x+VdX73+a3F7Qi+ZiknQOOecc653GEarVXZ/g56gcc45N2C8uOB17n3lIWqqajhiu4NYb1TeEWAGFE/QOOeccxXiosaf89aI2VALmHjmpX9wQM1xHLXTQaUOreQMaCVT6jC6pV8laObOmEHzrFnUTZvGxOnTSx2O68daW1uZP38+K1asKHUoPW7w4MGMGzeO6urqUofiXI9pfOlx3ho+h8Sg9q5Kotf7UtczZXk9aw0dWbrgykSm6Li65a3fJGjmzpjBZqeeytZA6p57mAueqHG9Zv78+YwYMYIJEyYgxRnsvDKYGc3NzcyfP59NNtmk1OE412PufeMeVLdmDoRl4La5jZyw25EliKp8GNBa4f3S9Zt+aJpnzaKGKIVWHead6y0rVqygrq6uXyVmACRRV1fXL3Oe3MCWJg25vq6CNmvr83jKjZmRKjBVgn6ToKmbNo0U0BqmumnTShyR6+/6W2KmXX/9XG5g23O9vbHWNf/lKWF8YquGvg+ozBiitcBUCfpNkdPE6dOZC16Hxjnn3BoOn7gfD9+b5MPRb6BBGSwDZMT2rQexweh+3V9kLAZkKiMjJq9+k6CBUGfGEzJuAFi8eDE33HADX/7ylzu13yGHHMINN9zA6NGjeycw58pUIpHgso9/j3tfauLhtx6lRjVM3fIT7LDRlqUOrSwYkKrwQpt+laApqqmJN669lgeBLY4/nvr6+lJH5FyXLF68mF//+tdrJGjS6TRVVVV597vzzjt7OzTnylYikeCgbfbkoG32LHUoZSeqFNz1BI2kjYBrgfWADDDDzH4uaW3gJmACMA84xswWdTfeXCo7OdYZTU2k99uPDa+8kqOuvJLzGhpoamoqdVTOdcm5557La6+9xo477sguu+zCfvvtx2c/+1kmTpwIwBFHHMHOO+/Mtttuy4wZM1buN2HCBBYuXFiqsHuMpKslvSvpuTzrJekXkl6V9KykSX0do3OVxBCtVpV3iqEN+IaZbQ3sDpwuaRvgXKDRzLYAGsN8rxg4OTTJJEqlqCJKie7Z2koymfRcGtdnmpqaSCaTNDQ0dPu+u/TSS3nuued45plnSCaTHHrooTz33HMrm1pfffXVrL322nz00UfssssuTJs2jbq6ftUj6h+AK4h+EebyCWCLMO0G/Ca8OudyMCDdjcq/ZrYAWBDeL5X0IrAhMBVoCJvNBJLAt7oRal4FEzSSPhXjGCvMrPzzsRsasJoaWltaaAUera7mkoaGUkflBoimpiamTJlCKpWipqaGxsbGHk1M77rrrqv1G/OLX/yCW265BYA333yTV155pV8laMzsIUkTCmwyFbjWzAx4TNJoSeuHh65zroMoh6ZgkmCMpOyRW2eY2YxcG4bv5k7A48C67d87M1sgaZ0eCnkNxXJofgvcSu7W++32AXImaCRtSVR21m5T4AKiX1V9Uqa2Un09VQ88wPxQh+YSr0Pj+lAymSSVSpFOp0mlUj2eOzhs2LDVznXffffR1NTE0KFDaWhoGIj9ymwIvJk1Pz8sWyNBI2k6MB1g/PjxfRKcc+XGTKQKFy0tNLPJxY4jaTgwCzjLzD7oy24giiVo7jKzkwttIOm6fOvM7GVgx7BdFfAWcAurytQulXRumO+VLKjV1NezcX09x8fdvqkJkkloaABP/LhuaGhooKamZmUOTUM3cwdHjBjB0qVLc65bsmQJa621FkOHDuWll17iscce69a5KlSup2jORqnhV+YMgMmTJ+dtuLpo+RIWL1/C+LU2LFjx2rlKFI3l1L37WlI1UWLmejP7a1j8TnvuqKT1gXe7F2l+BRM0Zva5YgeIs00wBXjNzN6Q1Gdlal01d8YMtjrjDAal06i2FhobPVHjuqy+vp7GxsYeq0NTV1fHnnvuyXbbbceQIUNYd911V647+OCDufLKK9l+++3Zcsst2X333Vfbd4B0nDcf2ChrfhzwdlcOtHj5En7+xAVUrf0qINrmVbPdoOM4asejeiJO58qESHevlZOA3wMvmtlPs1bdBpwAXBpeb+1OlIXErhQsaQ+iIqKV+5hZvgp5uRwL3BjexypTK1VWcFNTE3edfjoXtLUhwFpaUDLpCRrXLfX19T1azHTDDTfkXF5bW8tdd921xvJ0Os3SpUsZOXJADMJ3G3CGpD8RVQZe0tX6Mz+b/S2q155P1SADjKpBLbzYOpOHX9uAvTfboydjdq5kombb3cqh2RP4PDBX0jNh2f8RJWRulnQK8F/g6O6cpJBYCRpJfwQ2A54B0mGxkb+FQcf9a4DDgfM6E1zcrOCelkwmuT+T4VxCHnVVFdVegdhVuG233ZYvfOEL/WIUbUk3EuXyjpE0H/gu0TBumNmVRPX6DgFeBZYDJ3XlPK+/N4/q0fOpGrT6oIZVVWnun3+DJ2hcvxGjUnDh/c0eIX992yldPnAnxI1+MrBNaDHQFZ8AnjKzd8J8n5WpdUVDQwMX1dZyYEsL+ycSHH3FFUz03BlX4V566aVSh9BjzOwzRdYbcHp3zzN/8Xwy6cQaCRolIF39fncP71zZyFC0UnDZi5ugeY6o97+uNnn8DKuKm6APy9S6omN9B0/MODcwbbP+Njz0emaN5em0GNK6SY49nKtcmQrva7dYPzS3E5W6jABekPQE0NK+3swOL3YCSUOBjwOnZi3uszK1ropV38FbQTnXr40ZvjZDl9Tz4cjHGFQdlbZnMpBpq+LYbb5Q4uic6zlm6m4dmpIrlkPz/7p7AjNbDtR1WNZMH5Wp9ZowlIJSKaymhqoHHvBEjXP90Fl7nsf1c67jleV3kqhuIbFsPEdt+WU2H+s5NK7/MCDVjTo05aBYs+0HAST9yMxWa1Yt6UfAg70YW1l749pr2bClhSqgtaWF+ddey8aeoHGu30kkEnx+l+Mhfg9WzlWc9rGcKlncArOP51j2iZ4MpNI8CKSA1jAN2JSdK3tXXnkl1167ZoPEefPmsd1225UgIudcuTEgY4m8UyUoVofmS8CXgc0kPZu1agTwaG8GVu62OP54Drn6avZsbY3GhTref7258tPW1sZpp51W6jCcc2WuP+TQFCswuwG4C7iE1Yf8XmpmA7rNYn19PZckkySTSS7J1/OrVxp2veyiiy7i+uuvZ6ONNmLMmDHsvPPO3HHHHeyxxx48+uijHH744SxdupThw4dz9tlnM2fOHE4++WSGDh3KXnvtVerwnXNlot9XCjazJZKWAhPN7I0+iqliFGwJ1dQEU6ZAKgU1NT50guvxBO7s2bOZNWsWTz/9NG1tbUyaNImdd94ZgMWLF/Pgg1FB6IUXXrhyn5NOOolf/vKX7LvvvpxzzjndjsE51z/0QE/BJVe0YMzMMsC/JPkwtJ2RTEaJmXSa9IoVvJGjDoMbQNoTuN/5TvTa1NTtQz7yyCNMnTqVIUOGMGLECD75yU+uXPfpT396je2XLFnC4sWL2XfffQH4/Oc/3+0YnHP9RTSWU76pEsRto7U+8Hzoh+bD9oVx+qEZsBoaSA8aRCadptWME66+mkuOP75Hx/JxFSQrgUsqFc13814o1HH3sGHDcm4/QAamdL3g7SVvcfsLP8Nq55JuG8qEYcdw6LYFO2x2FWRA5NAE3wMOA74P/CRrcvnU13P9SSdxocQU4JF0mmQyWeqoXKk0NERFj1VV0WsPjA221157cfvtt7NixQqWLVvG3//+94Lbjx49mlGjRvHII48AcP3113c7BjcwvPPBAv7+xnEMWftRRo5czFprv83CxC+55vELSx2a6yGGaLOqvFMliJVDY2YPSloX2CUsesLMymr8pXK0xfHHc9rMmaRSKWpqamjwAS4Hrvr6qB5VD9ah2WWXXTj88MPZYYcd2HjjjZk8eTKjRo0quM8111yzslLwQQcd1O0Y3MBwx4u/pnZUC4OqVg0DUV2dpm3UPTR/eDp1w8aWMDrXE8ygNVMZRUv5KM54k5KOAS4DkkSjae4NnGNmf+nV6LJMnjzZZs+e3Ven6zFNTU0rx4QqWNzkLaIqyosvvsjWW29d6jBYtmwZw4cPZ/ny5eyzzz7MmDGDSZMmdfu4uT6fpDlmNrnbBy9jlfqc6W2/enIqa631vzWWt7RUs2nt+ey7+YDulqxsdeY7u842dfbp6w7Ou/6KnW8o++9/3Do05wO7tOfKSBoL3Af0WYKmUsUeE8pbRLkumD59Oi+88AIrVqzghBNO6JHEjHMdWdvaZDL/I9HhB3wikWHssA1LE5TrUVEdmsrOoYmboEl0KGJqJn79G1dMMom1tKBMJnrtgQqjbmC44YYbSh2CGwAmr3ciL684j0QivXJZOi2WLRvLNltvX8LIXM9RxdSVySduouQfku6WdKKkE4G/A3f2XlgDy9y6Oj7KZGgFPspkmFtXV3QfV3pximsrUX/9XK7rdt9kX0amTqGlpZpUahBtbQmWLB7HtC1/U+rQXA+J6tBU5Z0qQdxKwedImgbsSVSHZoaZ3RJnX0mjgd8B2xHlap0MvAzcBEwA5gHHmNmiTsbeb9zR3MzfEwn2zmR4OJHg0OZmJpY6KFfQ4MGDaW5upq6url81hTYzmpubGTx4cKlDcWXmyB1OoaXtc7z0zlxGD65j4618tPH+xBAZq+xnWeyxws1sFjCrC+f4OfAPMztKUg0wFPg/oNHMLpV0LtGwCt8qdJD+rKGhgYtqa3kstIa6zFtDlb1x48Yxf/583nvvvVKH0uMGDx7MuHHjSh2GK0O1g2rZYcOyrhfqusiAtoFQh0bSp4AfAesQ5dAIMDMbWWS/kcA+wIlEO6SAlKSpQEPYbCZR66kBm6Cpr6+nsbExXmsoVxaqq6vZZBP/heqc6ydMtFVI0VI+cXNofgx80sxe7OTxNwXeA66RtAMwBzgTWNfMFgCY2QJJ6+TaWdJ0YDrA+PH9e+SFWK2hnHPOuV7QH3Jo4kb/ThcSMxAlmCYBvzGznYiGTTi38C6rmNkMM5tsZpPHjvWOm5qamrjkkkto6oFxgJyrdJIOlvSypFdD0XXH9Q2Slkh6JkwXlCJO5yqBARlT3qkSxM2hmS3pJuBvQEv7QjP7a5H95gPzzezxMP8XogTNO5LWD7kz6wPe63ARTU1NTJkyZWWvw42NjZ6j4wYsSVXAr4CPEz1nnpR0m5m90GHTh83ssD4P0LkKY4i2Cu8pOG6CZiSwHDgwa5kBBRM0ZvY/SW9K2tLMXgamAC+E6QTg0vB6a2cDH2iSySSpVIp0Ok0qlSKZTHqCxg1kuwKvmtnrAJL+BEwlerY412sWL1/Io6/9kRVtC5mw9r7stNEBJDr2OFiJrPKLnOI22z6p0HpJ55nZJXlWfwW4PrRweh04iaio62ZJpwD/BY6OH/LAdFhdHR9J3J9I8JSPC+XchsCbWfPzgd1ybFcv6V/A28DZZvZ8xw0GUl091z1z/nsf76/4KrWDMwyrSrMoNYubntqSo3b8M9WDakodXrcYdCuHRtLVRINYv2tm24VlFwJfJKpLC/B/ZtZrfdj1VHIsb4LEzJ4J9WC2N7MjzGyRmTWb2RQz2yK8vt9DcfRPTU1MPOssvpfJ8EBVFY9ffrnnzriBLlehfsceAZ8CNjazHYBfEhWZr7mT19VzMaTTaf637GwG16SorW6jKmHUVrcxdtTL3PNS5Xcw2N4PTTfq0PwByDUY1M/MbMcw9WqHvD2VoKmMGkOVKpmEVAplMlRnMkxsbi51RM6V2nxgo6z5cUS5MCuZ2Qdmtiy8vxOoljSm70J0/cnzC/5JTXXLGstrq9v4sLV/1JpIWyLvVIyZPQSUNHOipxI03ld6b2poiAatrKqKXgsUN3lLKDdAPAlsIWmTUJx9LHBb9gaS1lPoxlnSrkTPO/814Loo/7859YN/gWaQziTyTsAYSbOzpukxD32GpGclXS1prV78CPF7Ci7Cc2h6U319NAJ3MhklZvIUN3lLKDdQmFmbpDOAu4Eq4Goze17SaWH9lcBRwJcktQEfAceaD1Tlumjb9fdk/r9rGFzdutryltZBDBnUHxrSqT3hks9CM+tsN9G/AS4iSg1eBPyEaPij/FFI3zezC7Lmq4Brzey4YifrqQTNn3voOC6f+vqiI3Bnt4Sa1NJCy4UXwoUX+sjdrl8KxUh3dlh2Zdb7K4Ar+jou1z9VVVWxzpAfs6j1LBIYgwa10dpWzcIPNmXaDl8pdXjd1t4PTY8e0+yd9veSfgvcEWO38e0NjSTVEqUvnopzvlhFTpJ+LGmkpGpJjZIWSvpcVtA/jHMc17saGhqoqalhz0SCezIZ9r3vPpgyBbz4yTnnum3yhIPYaf0kyz86iXcXHcmIQZdxzE63V3wLJwAM0qa8U1eEfubaHQk8F2O3k4CJks4DbgceMLML45wvbg7NgWb2TUlHElXGOxp4ALgu5v6uD7SPCdVy4YUMue8+lMlAKhUVVXkujXPOdduY4evyye3739CDVrzIqSBJNxKN0ThG0nzgu0CDpB2JMoDmAacW2H9S1uzPgauAR4EHJU0ys6K5NHETNNXh9RDgRjN7P9S1c2Wmvr4+KmZ6+OEoMVOkErFzzjkHkMl0/f+6mX0mx+Lfd+IQP+kwvwjYJiw3YP9iB4iboLld0ktEFeu+LGkssKITgbq+FKMScVNTk4/u7ZxzDohaOVkJx2wys/26e4y4PQWfK+lHwAdmlpa0nKibcVeuClQi9tZQzjnnOkp3I4emp0j6IfBjM1sc5tcCvmFm3y62b9xKwUOB04maYAFsAHS2+ZYrE8lkkkktLZwTWkMlk8lSh+Scc66EDJHJJPJOfegT7YkZADNbRFTdpai4RU7XAHOAPcL8fKKmVHGaYLkyc1hdHWdmMtQAqUyG1+rqSh2Sc/1eS+uHzF/0InXDNmL0sHVLHY5zq7Oeb7bdRVWSas2sBUDSEKA2zo5xEzSbmdmnJX0GwMw+ktcKrlgTm5uxRAJlMlQlEj6UgnO9rPH5Cxk37BoAmlekefK1ndlrq+sYUjOixJE5l6U8up28DmiUdA1RRCcDM+PsGDdBkwqpJAOQtBmw5qAWrjI0NKDa2mh8KG8F5Vyvanr1asYPv5rBg1b1MLvJqDk8/NIJHLj9X0sYmXOr604rp55iZj+WNBeYQjQKwUVmdnecfeMmaL4L/APYSNL1wJ7AiXF2lDQPWAqkgTYzmyxpbeAmYAJR2/RjQjmZ6wsxh1JwrrskbWJm/ym2rD9Lp37D4CGrd5dfOyjNZqOfYOlHCxkxxMfLdKVnBta3dWXyMrO7gLs6u1+s6M3sXuBTRImYG4HJZpbsxHn2C0OHt1ckPhdoNLMtgMYw7/pSfT2cd17BxIwPdOl6wKwcy/7S51GU0IiaJTmXZyzB4o/eybnOuVKwTP6pr0jaXdKTkpZJSklKS/ogzr6xcmhCfZlPAJua2fcljZe0q5k90cWYpxL1KAhR2VgS6H9dL1awuTNmcNfpp3N/JsNFtbXetNt1iqStgG2BUZI+lbVqJDC4NFGVxrvLt2RU7RyqEqtXUGhNVzFh1BYlisq5jlTSfmiyXAEcS9TwaDJwPLB5nB3j5i/9GqgH2nsCXAr8Kua+BtwjaU7WcOPrmtkCgPC6Tsxjub7Q1MRWZ5zBBW1t3JPJeNNu1xVbAocBo4FPZk2TgC+WLqy+t/n636clXU1bVv2EFW3VvJs6g+qqfjAGkOsfDCyjvFOfhmL2KlBlZmkzuwaI1ele3Do0u5nZJElPh5MtkhT3m7inmb0taR3g3tDjcCwhATQdYPz48XF3c92VTDIonUaE/qYTCRq84rDrBDO7FbhVUr2ZDegyy03G7sQ8/Y1/v30RYwa/xNLU2gwb9lUatjqm1KE5t7ryyKFZHtIXz0j6MbAAGBZnx7gJmlZJVaxq5TQWiFWqZmZvh9d3Jd0C7Aq8I2l9M1sQRuN8N8++M4AZAJMnTy6PBmUDQWgFZS0tUFXF0VdcwUQvbnJd87Sk04mKn1YWNZnZyaULqe9NGLMTE8Z4iyZX5vqwrkwBnycqPToD+BqwETAtzo5xi5x+AdwCrCPpYuAR4IfFdpI0TNKI9vfAgUTDh98GnBA2OwG4NWYcri+EVlD6wQ+ofvBBJk6fXnwf53L7I7AecBDwIDCOqMjaOVdOjCiHJt/UB0LGycVmtsLMPjCz75nZ10MRVFFFc2gkJYD/AN9kVbvwI8zsxRjHXxe4JfTBNwi4wcz+IelJ4GZJpwD/BY6OE6zrQwXGgsrWPsjlYXV1UQd93gzcrW5zMzta0lQzmynpBiBWnxLOub7Vl62Zcp4/GityrKQaM0t1dv+iCRozy0j6iZnVA7Hrv4R9Xwd2yLG8mShx5CpY+yCXk1paODOTiXofrq2N+rjxRI2LtHfAsljSdsD/iPqfcs6VGZVBx3pEfdM9Kuk24MP2hWb202I7xi1yukfSNB/uwGVLJpOkUin2DuNCKZOBVCrqsM+5yIwwWu53iIqaXwB+XNqQnHNrMEGmwNR33iYaJzIBjAjT8Dg7xq0U/HWiWsZtklYQFTuZmY3sfKyuv2hoaKCmpoaHW1pIhXGhfCgFl83MfhfePghsWspYnHNFlEfTmxfM7M/ZCyTFqpYSK0FjZj6CmltDfX09jY2NJJNJXvM6NC6LpK8XWh8n+zjGOQ4Gfg5UAb8zs0s7rFdYfwiwHDjRzJ7q7nld16RSK7j9vm+x/YR/MqS6jafnbc4OW1/M+HEfK3Vorl15tHI6j6hTvWLL1hC3p+BJORYvAd4ws7Y4x3D9U319vfcg7HLp1R9BoTXEr4CPA/OBJyXdZmYvZG32CWCLMO0G/Ca8uhJINk3joO1fZ+igNABjtpnLe8s/zcL372LM2uuVODqHlbYOjaRPEP342FDSL7JWjQRipTPiFjn9mqiHz7lhfiLwL6BO0mlmdk/M47gBaO6MGTTPmkXdtGneBHyAMLPvxdlO0nlmdkkXTrEr8GpoeICkPxENqZKdoJkKXGtmBjwmaXR7/1ddOJ/rhmeeS7LXZq8zOCRmAGqqMqw1uIW7my5l2qGXly44t0ppi5zeBmYDhwNzspYvJeqPpqi4CZp5wClm9jyApG2Ac4CLgL8CnqBxOc2dMYPNTj2VrYHUPfcwFzxR47IdDXQlQbMh8GbW/HzWzH3Jtc2GRD2Puj40b/4DbDoyAaRXWz60uo21RnWq8azrRSphgsbM/gX8S9INZtaabztJs8wsZ0d7cVs5bdWemAknfgHYqf3XkXP5NM+aRQ1Ryrk6zDuXpat53Ln26/g4jrMNkqZLmi1p9nvvvdfFcFwhI4ZNQDn+W6bSCZZ+WFeCiNwajLJo5VQoMRPkbVwQN0HzsqTfSNo3TL8G/i2pllX9TDi3hrpp00gR3SStYd65LF39TTifqEv0duOIsqw7uw1mNsPMJpvZ5LFjx3YxHFfIfnscz7vLhq02QCdAWybBZhufXqKoXEfK5J/KSN5nRtwEzYnAq8BZRGVZr4dlrcQcBdMNTBOnT+e1q67i0QMP5LWrrvLiJtdRV3/6PQlsIWmTMJDdsUT93GS7DThekd2BJV5/pjQSVVVY7e95dsG6tKQTfNRWxVtLh/HP17/OdlvtUerwXDsrMFWAuM22Pwq5MneY2csdVi/r+bBcfzJx+nQolJBpaoo64/Mm3/2GpB+Z2bckHd2xT4kOijbFzMXM2iSdQTSMQhVwtZk9L+m0sP5K4E6iVhOvEjXbPqkr53I9Y/NNtwce4o03X2L58iVsuflkNtqiqtRhuUDWvZwYSVcDhwHvmtl2YdnawE1EvYPPA44xs0XdDTXfilg5NJIOB54B/hHmdwzdEjvXPU1NMGUKfOc7pPfbj2u/9CWamppKHZXrvkMkVRP1H5GXmRUd5LbAvnea2cfMbDMzuzgsuzIkZrDI6WH9RDOb3dVzuZ6z8UZbsfWWu5Go8sRM2eleHZo/AAd3WHYu0GhmWwCNYb67vpVvRdwip+8SNZNcDGBmz+DjsbiekExGwyWk02RaWnj5qquYMmWKJ2oq3z+AhcD2kj7ImpZK+qDUwTnn1tSdOjRm9hDwfofFU4GZ4f1M4IiiMUhzJT3bYXpY0s8k1RXqJiZugqbNzJbE3Na5+BoaoKaGtEQrcL8ZqVSKpI8HVem+bWajgL+b2cisaYQPmeJcGbJQ7JRnAsa0twYMU5wKkeu211sLr+vE2Ocu4O/AcWG6HXiIaGDbPxTaMW4/NM9J+ixQJWkL4KvAP2Pu61x+9fXQ2Mj8a6/lhKuv5sl0mpqaGhp8PKhK10TUGafnxjhXKQrnxCw0s8l9EMWeZrZn1vxcSY+a2Z6SPldox7gJmq8A5wMtwI1EFfEuihtd6KZ8NvCWmR3WSxWFXKWqr2fj+nouOf54kskkDQ0NPpxC5auRdAKwh6RPdVxpZn8tQUzOuQJ6oWO9d9p755a0PvBujH2GS9rNzB4HkLQrq0bbLjgEQtxWTsuJEjTnx9k+hzOBF4nGZIBVFYUulXRumM9b0ccNDLHHhfJWUZXgNKLs4tHAJzusM6Iexp0bcFpSrVx5+/08/t832GjkaL429eOsN2atUocVxnLq8aPeBpwAXBpeb42xzxeAqyUNJ2rR9AHwBUnDKNKreMEEjaTbKdAC3cwOLxaZpHHAocDFQPsIvFOBhvB+JpDEEzQuhrkzZrDVGWcwKJ1GtbXQ2OiJmjJkZo8Aj0iabWa/L3U8zpWD997/gP1+92s+HCYyIxIkWhbx5z/N4Le7HcqUXbYrdXjd6m9G0o1E/9fHSJpP1JjoUuBmSacA/yUa6qRwCGZPAhMljQJkZouzVt9caN9iOTT/L7x+ClgPuC7Mf4aoqCiOy4Fvsvrou6tVFJKUs6JQqHQ0HWD8+PExT+f6q6amJu46/XQuaGtDgLW0oGTSEzRlSNL+ZnY/sMiLnJyLnDnzJpaNEFYdtcfJ1CQgY3w1eQfPlzhBI7qXQ2Nmn8mzakqn4ohGIJhGVCVlkKT243+/2L4FEzRm9mA4wUVmtk/WqtslPRQjsPZOduZIaii2fY7zzwBmAEyePLlC+ip0vSWZTHJ/JsO5hB8SVVVUe+XhcrUPcD9RcZMRPS+zXz1B4wac2daMVXfofychPhwuXpr3NltN2KA0gUFvFTl1xa3AEqIRt1s6s2PcSsFjJW3aPhilpE2AOIOe7AkcLukQYDAwUtJ1dK2ikBvgGhoauKi2lgNbWtg/keDoK65goufOlKulkr4OPMeqhAxUTCfqzvW8vJVuJaoHlUFHg+WRoBlnZh076IslboLma0BSUvvo2hMIRUGFmNl5hJ5CQw7N2Wb2OUmX0fmKQm6Aq6+vp7GxcWVLqLyJGa80XA7aWyVsCexC9B0XUY5N0dxd5/qjPWvW5f7Wd1cWOQGQMUYuMTYbt27pAgt6oZVTV/xT0kQzm9vZHeO2cvpH6H9mq7DoJTPrVFZQB52uKOQcxGgJ1T6UQioFNTVeabhEzOx7AJLuASaZ2dIwfyFdHL/JuUr385OPZb9fXkHziAyZKpFIG1Vtxm8PParUoUV5p+WRQ7MXcKKk/xAVOYloJJPti+1YrJXTJDN7iuhoLcC/Cm1TiJkliVozYWbNdLKikHOxJJNRZeFMxisNl4fxQCprPoUPm+IGqBHDhvDEN7/Bzfc/zqOvvMb4urU47fD9GTF0SKlDA8qmDs0nurpjsRyaa0JRUaGRqX4P7NTVAJzrSXPr6tgsk6EaaM1keK2ujomlDmpg+yPwhKRbiH4DHsmqsV2cG3ASiQTHHlDPsQeU3w+tckjQmNkbknYA9g6LHjazNTJTcimWoBlFVNO4UILmvTgncq4v3NHczN8TCfbOZHg4keDQ5mZP0JSQmV0s6S5WPZxOMrOnSxmTcy4Hoyyq7Es6E/giq1pCXidphpn9sti+xZptT+h+eM71nfaWUI+lUtTU1HCZN+suuVAkXbRY2jlXOqJsKgWfAuxmZh8CSPoR0dhw3UvQOFdpOraEylmB2FtBOefcGsqhyIkobZXOmk9TuJRoJU/QuH6nYEsobwXlnHNrKp+O9a4BHg/17gCOIKqrW1Si+CbO9SPJZJSYSaexlhaSF15IU1NTqaNyrqhMJkPyXy/ywDMv0JZOF9/Buc6yAlNfhWD2U+Ak4H1gEVG9u8vj7Bsrh0bRYArHAZua2fcljQfWM7MnuhaycyXS0AA1NVhLCx9lMnz7vvt46uGHaWxsjDfSt3MlcOcT/+IbT99B6+BoftCTcMm2B3PkHjuXNjDXr5Qyh0bS2lmz88gaL1LS2mb2frFjxM2h+TVQTzQoJcBS4Fcx93WufNTXQ2MjDx5wAAcmEjyayZBKpUgmk6WOzLmc3lv8AWc+fwctI0WmJppSw8U5//4Hby0s+ox3Lp5Q5JRv6gNzgNnhtf397Kz3RcWtQ7ObmU2S9DSAmS2SVNP5eJ0rA/X11F54IU89/DBVoTVUQ6HWUF6J2JXQr+95AMvx09MEv7i7kR8d5x2tu+7r7mjb3WVmm6yMJcqt2YJoDMjY4iZoWiVVEUrSJI2lXDpJdq4LYrWGAq9EXIbCw+4moh6H5wHHmNmiHNvNI8pNTgNtZja576LsOe8s+4BMjp+PVg3vfbis7wNy/VcZNNuW9AXgTGAc8AywO/BPYowuELfI6RfALcA6ki4GHgF+2JVgnSsX9fX1nHfeeYXrzmRVIiaViuZdqZ0LNJrZFkBjmM9nPzPbsVITMwBTPrYViRx1gBOtsO+mW/R9QK5/MlDG8k596EyiAW3fMLP9iEYiWBhnx1gJGjO7HvgmcAmwADjCzHyAOdfvza2rozWRwBKJKIfGO+orB1NZNXzCTKJmnf3WkXvszHpLalFq1T8VpYwxH1Tz+f32KGFkrr8pcR2adivMbAWApFozewnYMs6OBRM0ktZun4B3gRuBG4B3OtRIdq7faWpqYrezzmK/dJrvJhLMvfxyL24qD+ua2QKA8LpOnu0MuEfSHEnT8x1M0nRJsyXNfu+98hvJJZFI0HjamXy6ektGNScY1SymVW1B4xfPJJHwnjdczymTBM18SaOBvwH3SroVeDvOjsXq0MwheiiIaNTcReH9aOC/wCZ59wQkDQYeAmrDuf5iZt+NWwbuXCklk0lSqRSPZjI8JjHEx4XqM5LuA9bLser8ThxmTzN7W9I6RA/Gl8zsoY4bmdkMYAbA5MmTy6AWwZpqa6r54WeP8nJ+13usPIY+MLMjw9sLJT1ANKbkP+LsW2wsp00AJF0J3GZmd4b5TwAHxDh+C7C/mS2TVA08Egaq+xRRGfilks4lKgP/VpyAnesrDQ0N1NTUkIrTEsr1KDPL+3yR9I6k9c1sgaT1iXKPcx3j7fD6buh1dFeiH1jOuQ5K3copFzN7sDPbx82v3KU9MRNOchewb4xgzMzaq+FXh8kYYGXgrjK1t4S66KKLvOO98nIbcEJ4fwJwa8cNJA2TNKL9PXAg8FyfRehcBSqTSsFdFrfZ9kJJ3wauI0qQfA5ojrNjaO49B9gc+JWZPS5ptTLwkCWca9/pwHSA8ePHxwzVuZ5TcFwoVyqXAjdLOoWo6PtoAEkbAL8zs0OAdYFbok7OGQTcYGaxsq2dG5AMVOEjasRN0HwG+C5R022Ism0/k3/zVcwsDewYKvncImm7uMFVQtm2c01NTcX7s3E9xsyaydEnRShiOiS8fx3YoY9Dc66yVfh/2VgJmjCGwpndOZGZLZaUBA4maiVVtAzcuXLX1NTElClTVtaz8aIp51xFCv3QVLK4g1M+QI60m5ntX2S/sUBrSMwMIapI/CNWlYFfSp4ycOcqQTKZZFJLC3tnMjzc0kIymfQEjXO95IZH/skV/06ydGiK4curOW2zfThh371LHVa/0d1KwaXunTtukdPZWe8HA9OAthj7rQ/MDPVoEsDNZnaHpCZylIE7V2kOq6vjzEyGGiCVyfBaXV2pQ3KuX7om+SA/WfgAFr5iHwxp5bL3GvnovhSnHVC0V3xXhKzHKv/uZ2axevbtaXGLnOZ0WPSopKLNqczsWaJuizsuz1kG7lylmdjcjCUSKJOhKpFgYnOsuvLOuU761RsPYaNXX2bVMOPtRznN/530iHLoh6Y74hY5ZfcKnAB2JnenV84NLA0NqLYWUinkQyM412s+GtlG1FvK6lpGpslkMt5rcg/ogX5o2nvnNuCq0LCnz8QtcsruMbgN+A9wSm8F5VzFqK+PRuBOJqPETIH6M94ayrmuG7Q8QdvwNbMQqj6SJ2Z6ggHpglk0YyTNzpqfkSPBEqt37t4SN0GzdftgUe0k1fZCPM5Vnvr6omM8eWso57rnk4O34W+p57GaVcuUggMTHytdUP1MkTo0C4tV8i1179xxk7X/zLGsqScDca4/ax8Xapd0mq+tWMEr115b6pCcqyjfP3IaB6S3ILEC1AqJFtg3tQk/PvrYUofWb8jyT0X3LYPeuQvm0EhaD9gQGCJpJ1YVYI4EhvZybM71Gw0NDexVVcWd6TQ1Zuiaa+D44330budiSiQSXH7scbSkWnn7/fdZb621GFJbU3xHF4u63w9NyXvnLlbkdBBwIjAO+GnW8qXA//VSTM71O/X19cw8+WQGX3UVCTNoa4vq3XiCxrlOqa2pZpP11i11GP2SCtehKagceucuNtr2TKJ+ZKaZ2aw+ism5fmnj44+HmTMhlQJvEeWcKydm0J97Cpb0OTO7Dpgg6esd15vZT3Ps5pzLJUaLKG8J5Zwrlf7eD82w8Do8x7oK/+jOlUCBFlHeEso5VzLWvSKnclCsyOmq8PY+M3s0e52kPXstKucGIB8XyjlXUv25yCnLL4FJMZY557rIx4VyzsVx17+eYtbLTRjGpz5Wz6E77twjx1Wm+10Fl1KxOjT1wB7A2A51aEYCVb0ZmHMDjY8L5Zwr5vS/XMW/hr0Go6P5Sxb8l9v+/QRXHfOl7h3YgMpOzxTtWK+GqP7MIGBE1vQBcFTvhubcANM+LlRVVfTqraCcc1maXnk5SsxUE/33TgDV8PyI//DQS89369jCUCaTd6oExerQPAg8KOkPZvZGZw8uaSPgWqKBLDNEYz/8PAx2eRMwAZgHHGNmizp7fOf6lZjjQnlLKOcGpr89/zgMzrEiAbe+8AT7bLVt1w9efCynshe3Ds1ySZcB25J1Oc1s/yL7tQHfMLOnQpfIcyTdS9RZX6OZXSrpXOBc4Fudjt65/qbIuFBzZ8zgrtNP5/5Mhotqa70llHMDSE1iUO72xQY1VXH/nedXKTkx+cQdy+l64CVgE+B7RLkqTxbbycwWmNlT4f1S4EWioRSmAjPDZjOBIzoTtHMDUlMTW51xBhe0tXFPJsOk0BLKOTcwfG7yvrlXGBw3Kc+62CzqXC/fVAHiJmjqzOz3QKuZPWhmJwO7d+ZEkiYAOwGPA+ua2QKIEj3AOnn2mS5ptqTZ7733XmdO51z/k0wyKJ1mEFER+v6JBA0DsJ6NpKMlPS8pIynv6L+SDpb0sqRXQ06wcxVtyw025FOD9oJWIAWkBK1wuOrZbqPx3Tt4e5FTvqkCxM2jag2vCyQdCrxNNL5TLJKGA7OAs8zsgzB4VVFmNgOYATB58uTKuKLO9ZZQadhaWqCqiqOvuIKJA7O46TngU8BV+TaQVAX8Cvg4MB94UtJtZvZC34ToXO/4+scP5+jmPfjT7EcA49OT9mL82LE9cuxKL3KKm6D5gaRRwDeI+p8ZCZwVZ0dJ1USJmevN7K9h8TuS1jezBZLWB97tXNjODUCh0rCSSaobGgZqYgYzexGgyA+jXYFXw4B5SPoTUVG3J2hcxduobgznHHREzx7UDNIDIEFjZneEt0uA/QAknVVsP0VPnN8DL3YY9+k24ATg0vB6a/yQnRvAilQabuctodgQeDNrfj6wW64NJU0HpgOMH9/NbHvnKlmF1JXJpzvVor8OXF5kmz2BzwNzJT0Tlv0fUULmZkmnAP8Fju5GHM65LNljQu1VVcXMk0+ORvquoISNpPuIunvo6Hwzi/MDKFf2Tc6ntRdtO0eoQzMAcmjyKFoRxsweKbDdlG6c2zmXRzKZJJVKsUs6zZ3pNIOvugpmzoz6uKmQRI2ZHdDNQ8wHNsqaH0dU9885l5NBJl3qILolbiunXPyXjHNlqKGhgZqaGvaXqAESZpBKRR32DRxPAltI2kRSDXAsUVG3cy6X9hyafFMFKJigkbRU0gc5pqXABn0Uo3OuE+rr62lsbGTLU09dOZQCNTX9ZigFSUdKmg/UA3+XdHdYvoGkOwHMrA04A7ibqP+rm82se33D9wPPvvk6f33qYd5s9nYYLocK74em2NAHI/oqEOdcz6mvr48qAx9/fNGhFCqNmd0C3JJj+dvAIVnzdwJ39mFoZat56Qd89YEf8dGoxZAW1y0zxi3elMsP+xqJRHcy6l2/YQbpgVvk5Jwrd/X1cN55eRMzc2fMIHnQQcydMaOPA3N96Zz7L+ej0YtQtaHBGTTImD/ydS67/0+lDs2Vk0wm/1QBuj/4g3OuIs2dMYPNTj2VrYHUPfcwF5g4fXqpw3I9bNGHS1m81v9Qh6e9aownP3oS+GxJ4nJlxgzzHBrnXCVqnjWLGlg5lELzrFkljsj1hiXLPyRfY9NMTVvfBuPKW8byTxXAEzTODVB106aRIhrXpDXMu/5nfN066KM1M+MtA6OX5hxGzw1E7XVo8k0VwIucnBugJk6fzlyinJm6adO8uKmfSiQSHL32VG5ePguqDCXA2oC2BF/b+bhSh+fKSKUXOXmCxrkBbOL06eAJmX7v2F32Z9y/xzDzxdv5ILGEDdiA03c9ms3X3bDUoblyMVDGcnLOOVfZ9vrY9uz1se1LHYYrZ9b1BI2kg4GfA1XA78zs0p4KKy5P0DjnnHMDnHWjlZOkKuBXwMeJhh15UtJtZtano9t7gsY555xz3alDsyvwqpm9DiDpT8BUoE8TNLIK6dJY0nvAGzE2HQMs7OVwuqqcY4Pyjq+cY4Pyjq+nYtvYzMb2wHHKVieeM8WU2/3g8RTWX+OJ/Z2V9I9w3nwGAyuy5meEkeqRdBRwsJl9Icx/HtjNzM7oWthdUzE5NJ34o8w2s8m9HU9XlHNsUN7xlXNsUN7xlXNs5aanEmzlds09nsI8HjCzg7uxe66Ojvo8t8T7oXHOOedcd8wHNsqaHwe83ddBeILGOeecc93xJLCFpE0k1QDHArf1dRAVU+TUCeU8yl45xwblHV85xwblHV85x9Zflds193gK83i6wczaJJ0B3E3UbPtqM3u+r+OomErBzjnnnHP5eJGTc8455yqeJ2icc845V/EqKkEjaSNJD0h6UdLzks4My9eWdK+kV8LrWln7nCfpVUkvSzqoRPFdJuklSc9KukXS6LB8gqSPJD0TpitLENuFkt7KiuGQrH365NoViO2mrLjmSXomLO+z6xbON1jSE5L+FeL7Xlhe8vuuQGwlv+cGEklHh+ufkZS3ua2kg8M98aqkc3sxnrz3Zoft5kmaG+6F2b0QR8HPq8gvwvpnJU3q6Rg6GU+DpCVZ348LejGWqyW9K+m5POv79Nr0C2ZWMROwPjApvB8B/BvYBvgxcG5Yfi7wo/B+G+BfQC2wCfAaUFWC+A4EBoXlP8qKbwLwXImv3YXA2Tm277Nrly+2Dtv8BLigr69bOJ+A4eF9NfA4sHs53HcFYiv5PTeQJmBrYEsgCUzOs01VuBc2BWrCPbJNL8WT897Msd08YEwvxVD08wKHAHeF+3h34PFe/BvFiacBuKOP7pl9gEn5vo99eW36y1RROTRmtsDMngrvlwIvAhsSdbE8M2w2EzgivJ8K/MnMWszsP8CrRF0092l8ZnaPmbWFzR4jaqPfpwpcu3z67NoVi02SgGOAG3vj/DHiMzNbFmarw2SUwX2XL7ZyuOcGEjN70cxeLrLZyu7hzSwFtHcP3xvy3Zt9Kc7nnQpcG+7jx4DRktYvYTx9xsweAt4vsElfXpt+oaISNNkkTQB2IvpFuq6ZLYDonyOwTthsQ+DNrN3mU/ifeG/Fl+1kolR3u00kPS3pQUl7lyi2M0KW5tVZWdMluXZ5rtvewDtm9krWsj69bpKqQpHXu8C9ZlY2912e2LKV/J5zQN/eF/nuzY4MuEfSHEnTeziGOJ+3L69J3HPVhyLcuyRt20uxxFGy/1+VqiL7oZE0HJgFnGVmH0Q/4HNvmmNZr7dT7xhf1vLzgTbg+rBoATDezJol7Qz8TdK22fv0dmySfgNcRHRdLiIq2jmZEly7fNcN+Ayr5870+XUzszSwY6iLcouk7Qps3qfXLldsZvYclMc9119Iug9YL8eq883s1jiHyLGsy/dFoXg6cZg9zextSesA90p6KeQc9IQ4n7cvvytxzvUU0fhHyxTVJ/wbsEUvxVNMWQwnUEkqLkEjqZron971ZvbXsPgdSeub2YKQJfduWN7n3THniQ9JJwCHAVPMogJSM2sBWsL7OZJeAz4G9HjlvHyxmdk7Wet/C9wRZvv02hW4boOATwE7ty/r6+uWzcwWS0oCB1NG912O2J4rh3uuPzGzA7p5iB69LwrFIynfvdnxGG+H13cl3UJULNNTCZo4n7cvvytFz5WdsDezOyX9WtIYMyvFwJVlMZxAJamoIqdQl+L3wItm9tOsVbcBJ4T3JwC3Zi0/VlKtpE2IUtpP9HV8kg4GvgUcbmbLs5aPlVQV3m8a4nu9j2PLLpM9Emivcd9n167A3xXgAOAlM5uftX2fXbes840O74e0x0QZ3Hf5YiuHe86toS+7h893b64kaZikEe3viSqS52xx00VxPu9twPGK7A4saS8q6wVF45G0XngeIWlXov+Rzb0UTzF9eW36h67UJC7VBOxFlOX2LPBMmA4B6oBG4JXwunbWPucT1Wx/GfhEieJ7lagstH3ZlWH7acDzRLXtnwI+WYLY/gjMDctvA9bv62uXL7aw7g/AaR2277PrFs63PfB0iO85VrW2Kvl9VyC2kt9zA2ki+jEwnyj36x3g7rB8A+DOrO0OIWrF9xpRUVVvxZPz3syOh6i1z7/C9HxvxJPr8wKntX+niYpVfhXWzyVPC7E+jOeMrO/HY8AevRjLjURFwK3h3jmllNemP0w+9IFzzjnnKl5FFTk555xzzuXiCRrnnHPOVTxP0DjnnHOu4nmCxjnnnHMVzxM0zjnnnKt4nqBxzjnnXMXzBE0fk7Ss+FbdOv6dkkaH6ctd2L9B0h3Ft1xt+yWS7syz/g+SjupsHJUoXIs9sua/Jum/kq4oZVzO9RZJF0o6W9L3JR0Qlu0t6XlJz0gaIumyMH9ZoWP0beSuP6q4oQ9cYWZ2CKwc5PHLwK/74LQPm9lhvXkCSYNs1ejR5aoBWAb8E8DMfiZpETC5lEE519vM7IKs2eOA/2dm1wBIOhUYa9GwG871Gs+hKQOSdpT0mKIRr29RGPFaUlLSjyQ9IenfCiMjSxoq6eaw/U2SHpc0OaybJ2kMcCmwWfiVdFnHnBdJV0g6Mbw/WNJLkh4hGjepfZthikbgflLR6MxTY3wWhWO/IOnvZI3yK2lnRSM8z5F0d/uwC5J2CZ+lKcTaPrDiiZL+LOl2ohGBc8ajaLTpy8LyZ8MDFEnrS3ooXIPnVGBkaUkHhvM/Fc45PCy/IBz3OUkzsrpF/2r4jM9K+lNIQJ4GfC2cz0exdv2SpPMlvaxocMwtw7I/SDpK0heAY4ALJF0v6TZgGPC4pE/HOHannoXOZfMETXm4FviWmW1P1MX1d7PWDTKzXYGzspZ/GVgUtr+IrIEbs5wLvGZmO5rZOflOLGkw8Fvgk8DerD567/nA/Wa2C7AfcJmiMV8KOZLoITcR+CKwRzhPNfBL4Cgz2xm4Grg47HMNUXff9UC6w/HqgRPMbP8C8ZxCNM7JLsAuwBcVjaH0WaIu6HcEdiAaAiDXNRgDfBs4wMwmEQ3U+PWw+goz28XMtgOGEA32CNH13Sn8DU4zs3nAlcDPwjV/uMh1cq7iKBqh/VhgJ6IfP7tkrzez3xENoXKOmR1nZocDH4XvxE0xTtHZZ6FzK3mRU4lJGgWMNrMHw6KZwJ+zNmkfeXoOMCG83wv4OYCZPSfp2W6EsBXwHzN7JcRzHTA9rDsQODyrfHswMB54scDx9gFuNLM08Lak+8PyLYHtgHtDJkcVsEDRwIojzOyfYbsbWJVoALjXzN4vEs+BwPZaVVdnFNGgi08CV4fE1N/M7Jk8Me8ObAM8GmKrAZrCuv0kfRMYCqxNNM7L7URjJ10v6W/A3wpcD+f6k72BWywMeBpyYHpEF5+Fzq3kCZry117unGbV30tdOE4bq+fIDc56n29ALwHTzOzlTp4r1/EEPB9yYVYtDFnKBXxYLJ5QDPQVM7t7jZNK+wCHAn+UdJmZXZsntnvN7DMd9h1MVAdpspm9KelCVl23Q4kSb4cD35G0bZHP4Vx/UaoBAHM9C51byYucSszMlgCLssqEPw88WGAXgEeIyqmRtA1R8U5HS4ERWfNvANtIqg2/hKaE5S8Bm0jaLMxn/1O/G/hKVr2RnWJ8pIeAY0O9lvWJioYgGnV6rKT6cKxqSdua2SJgqaTdw3bHFjh2vnjuBr4UcmKQ9LFQ32Zj4F0z+y3we2BSnuM+BuwpafOw/1BJH2NV4mVhqFNzVFifADYysweAbwKjgeGsec2d628eAo5U1HppBFFRdY/o4rPQuZU8ldv3hkqanzX/U+AE4EpJQ4HXgZOKHOPXwMxQ1PQ0UfHHkuwNzKxZ0qOhgu1dZnaOpJvDtq+E/TCzFZKmA3+XtJAosbRdOMxFwOXAsyERMY/Vi4NyuQXYn6j8+9+EB5KZpUKR0C9CgmpQOPbzRHVgfivpQyDZ8bNkyRfP74iyoJ8Ky98DjiBqdXSOpFai1kfH5zqomb2nqIL0jZJqw+Jvm9m/Jf02fJZ5REVYEBWXXRc+h4jqzSxWVHn5L4oqK3/F69G4/sbMnpJ0E1F9tDeAnr7HO/ssdG4lmZUq99B1laQqoDokRjYDGoGPmVmqBLE0AGd3p9m2pOFmtiy8PxdY38zO7JkISysklCab2RmljsU55/ozz6GpTEOBB0IRi4AvlSIxE6SA7STd2d4HThccKuk8ovvxDeDEngqulCR9jagp96xSx+Kcc/2d59C4AUXS40Bth8WfN7O5pYjHuYFC0vnA0R0W/9nMLs61vXOd5Qka55xzzlU8b+XknHPOuYrnCRrnnHPOVTxP0DjnnHOu4nmCxjnnnHMV7/8D7uqsMLTH2ecAAAAASUVORK5CYII=\n",
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {
+ "needs_background": "light"
+ },
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "fig, ax = plt.subplots(1, 2, figsize=(8, 3))\n",
+ "\n",
+ "ds_trajectory.plot.scatter(\n",
+ " x=\"longitude\", y=\"latitude\", c=\"k\",\n",
+ " s=9, ax=ax[0], label=\"traj\"\n",
+ ")\n",
+ "selected.plot.scatter(\n",
+ " x=\"lon2d\", y=\"lat2d\", c=\"r\",\n",
+ " s=9, ax=ax[0], label=\"grid\"\n",
+ ")\n",
+ "ax[0].legend()\n",
+ "ax[0].set_title(\"target and selected locations\")\n",
+ "\n",
+ "selected.plot.scatter(\n",
+ " x=\"diff_lon\", y=\"diff_lat\", hue=\"along_track\",\n",
+ " ax=ax[1]\n",
+ ")\n",
+ "ax[1].set_title(\"differences\")\n",
+ "\n",
+ "fig.tight_layout();"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3",
+ "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.8.8"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py
index ba127ed..8395007 100644
--- a/src/xoak/accessor.py
+++ b/src/xoak/accessor.py
@@ -224,16 +224,13 @@ def _get_pos_indexers(self, indices, indexers):
return pos_indexers
- def sel(
+ def query(
self, indexers: Mapping[Hashable, Any] = None, **indexers_kwargs: Any
) -> Union[xr.Dataset, xr.DataArray]:
- """Selection based on a ball tree index.
+ """Directly query the underlying tree index.
The index must have been already built using `xoak.set_index()`.
- It behaves mostly like :meth:`xarray.Dataset.sel` and
- :meth:`xarray.DataArray.sel` methods, with some limitations:
-
- Orthogonal indexing is not supported
- For vectorized (point-wise) indexing, you need to supply xarray
objects
@@ -253,11 +250,35 @@ def sel(
indices = self._query(indexers)
if not isinstance(indices, np.ndarray):
- # TODO: remove (see todo below)
+ # TODO: remove (see TODO in self.sel below)
indices = indices.compute()
pos_indexers = self._get_pos_indexers(indices, indexers)
+ return xr.Dataset(pos_indexers)
+
+ def sel(
+ self, indexers: Mapping[Hashable, Any] = None, **indexers_kwargs: Any
+ ) -> Union[xr.Dataset, xr.DataArray]:
+ """Selection based on an tree index.
+
+ The index must have been already built using `xoak.set_index()`.
+
+ It behaves mostly like :meth:`xarray.Dataset.sel` and
+ :meth:`xarray.DataArray.sel` methods, with some limitations:
+
+ - Orthogonal indexing is not supported
+ - For vectorized (point-wise) indexing, you need to supply xarray
+ objects
+ - Use it for nearest neighbor lookup only (it implicitly
+ assumes method="nearest")
+
+ This triggers :func:`dask.compute` if the given indexers and/or the index
+ coordinates are chunked.
+
+ """
+ pos_indexers = self.query(indexers, **indexers_kwargs)
+
# TODO: issue in xarray. 1-dimensional xarray.Variables are always considered
# as OuterIndexer, while we want here VectorizedIndexer
# This would also allow lazy selection
diff --git a/src/xoak/tests/conftest.py b/src/xoak/tests/conftest.py
index c131b68..e985904 100644
--- a/src/xoak/tests/conftest.py
+++ b/src/xoak/tests/conftest.py
@@ -23,7 +23,11 @@ def indexer_array_lib(request):
@pytest.fixture(
- params=[(('d1',), (200,)), (('d1', 'd2'), (20, 10)), (('d1', 'd2', 'd3'), (4, 10, 5))],
+ params=[
+ (('d1',), (200,)),
+ (('d1', 'd2'), (20, 10)),
+ (('d1', 'd2', 'd3'), (4, 10, 5)),
+ ],
scope='session',
)
def dataset_dims_shape(request):
@@ -31,7 +35,11 @@ def dataset_dims_shape(request):
@pytest.fixture(
- params=[(('i1',), (100,)), (('i1', 'i2'), (10, 10)), (('i1', 'i2', 'i3'), (2, 10, 5))],
+ params=[
+ (('i1',), (100,)),
+ (('i1', 'i2'), (10, 10)),
+ (('i1', 'i2', 'i3'), (2, 10, 5)),
+ ],
scope='session',
)
def indexer_dims_shape(request):
@@ -64,6 +72,34 @@ def query_brute_force(dataset, dataset_dims_shape, indexer, indexer_dims_shape,
return dataset.isel(indexers=pos_indexers)
+def query_brute_force_indexers(
+ dataset, dataset_dims_shape, indexer, indexer_dims_shape, metric='euclidean'
+):
+ """Find indexers for nearest neighbors using brute-force approach."""
+
+ # for lat/lon coordinate, assume they are ordered lat, lon!!
+ X = np.stack([np.ravel(c) for c in indexer.coords.values()]).T
+ Y = np.stack([np.ravel(c) for c in dataset.coords.values()]).T
+
+ if metric == 'haversine':
+ X = np.deg2rad(X)
+ Y = np.deg2rad(Y)
+
+ positions, _ = pairwise_distances_argmin_min(X, Y, metric=metric)
+
+ dataset_dims, dataset_shape = dataset_dims_shape
+ indexer_dims, indexer_shape = indexer_dims_shape
+
+ u_positions = list(np.unravel_index(positions.ravel(), dataset_shape))
+
+ pos_indexers = {
+ dim: xr.Variable(indexer_dims, ind.reshape(indexer_shape))
+ for dim, ind in zip(dataset_dims, u_positions)
+ }
+
+ return xr.Dataset(pos_indexers)
+
+
@pytest.fixture(scope='session')
def geo_dataset(dataset_dims_shape, dataset_array_lib):
"""Dataset with coords lon and lat on a grid of different shapes."""
@@ -93,7 +129,22 @@ def geo_indexer(indexer_dims_shape, indexer_array_lib):
@pytest.fixture(scope='session')
def geo_expected(geo_dataset, dataset_dims_shape, geo_indexer, indexer_dims_shape):
return query_brute_force(
- geo_dataset, dataset_dims_shape, geo_indexer, indexer_dims_shape, metric='haversine'
+ geo_dataset,
+ dataset_dims_shape,
+ geo_indexer,
+ indexer_dims_shape,
+ metric='haversine',
+ )
+
+
+@pytest.fixture(scope='session')
+def geo_expected_indexers(geo_dataset, dataset_dims_shape, geo_indexer, indexer_dims_shape):
+ return query_brute_force_indexers(
+ geo_dataset,
+ dataset_dims_shape,
+ geo_indexer,
+ indexer_dims_shape,
+ metric='haversine',
)
diff --git a/src/xoak/tests/test_s2_adapters.py b/src/xoak/tests/test_s2_adapters.py
index fd47a68..2c1ecb5 100644
--- a/src/xoak/tests/test_s2_adapters.py
+++ b/src/xoak/tests/test_s2_adapters.py
@@ -16,6 +16,14 @@ def test_s2point(geo_dataset, geo_indexer, geo_expected):
xr.testing.assert_equal(ds_sel.load(), geo_expected.load())
+def test_s2point_via_query(geo_dataset, geo_indexer, geo_expected):
+ geo_dataset.xoak.set_index(['lat', 'lon'], 's2point')
+ ds_indexer = geo_dataset.xoak.query(lat=geo_indexer.latitude, lon=geo_indexer.longitude)
+ ds_sel = geo_dataset.isel(ds_indexer)
+
+ xr.testing.assert_equal(ds_sel.load(), geo_expected.load())
+
+
def test_s2point_sizeof():
ds = xr.Dataset(coords={'lat': ('points', [0.0, 10.0]), 'lon': ('points', [-5.0, 5.0])})
points = np.array([[0.0, -5.0], [10.0, 5.0]])
diff --git a/src/xoak/tests/test_scipy_adapters.py b/src/xoak/tests/test_scipy_adapters.py
index a15f119..70776c1 100644
--- a/src/xoak/tests/test_scipy_adapters.py
+++ b/src/xoak/tests/test_scipy_adapters.py
@@ -13,6 +13,15 @@ def test_scipy_kdtree(xyz_dataset, xyz_indexer, xyz_expected):
xr.testing.assert_equal(ds_sel.load(), xyz_expected.load())
+def test_scipy_kdtree_via_indexer(xyz_dataset, xyz_indexer, xyz_expected):
+ xyz_dataset.xoak.set_index(['x', 'y', 'z'], 'scipy_kdtree')
+
+ indexers = xyz_dataset.xoak.query(x=xyz_indexer.xx, y=xyz_indexer.yy, z=xyz_indexer.zz)
+ ds_sel = xyz_dataset.isel(indexers)
+
+ xr.testing.assert_equal(ds_sel.load(), xyz_expected.load())
+
+
def test_scipy_kdtree_options():
ds = xr.Dataset(coords={'x': ('points', [1, 2]), 'y': ('points', [1, 2])})
diff --git a/src/xoak/tests/test_sklearn_adapters.py b/src/xoak/tests/test_sklearn_adapters.py
index 0a019c5..f4c5ec6 100644
--- a/src/xoak/tests/test_sklearn_adapters.py
+++ b/src/xoak/tests/test_sklearn_adapters.py
@@ -38,10 +38,13 @@ def test_sklearn_balltree_options():
assert ds.xoak._index._index_adapter._index_options == {'leaf_size': 10}
-def test_sklearn_geo_balltree(geo_dataset, geo_indexer, geo_expected):
+def test_sklearn_geo_balltree(geo_dataset, geo_indexer, geo_expected, geo_expected_indexers):
geo_dataset.xoak.set_index(['lat', 'lon'], 'sklearn_geo_balltree')
- ds_sel = geo_dataset.xoak.sel(lat=geo_indexer.latitude, lon=geo_indexer.longitude)
+ ds_indexers = geo_dataset.xoak.query(lat=geo_indexer.latitude, lon=geo_indexer.longitude)
+ xr.testing.assert_equal(ds_indexers.load(), geo_expected_indexers.load())
+
+ ds_sel = geo_dataset.xoak.sel(lat=geo_indexer.latitude, lon=geo_indexer.longitude)
xr.testing.assert_equal(ds_sel.load(), geo_expected.load())
@@ -52,4 +55,7 @@ def test_sklearn_geo_balltree_options():
# sklearn tree classes init options are not exposed as class properties
# user-defined metric should be ignored
- assert ds.xoak._index._index_adapter._index_options == {'leaf_size': 10, 'metric': 'haversine'}
+ assert ds.xoak._index._index_adapter._index_options == {
+ 'leaf_size': 10,
+ 'metric': 'haversine',
+ }