{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "570a34d8-7958-4de7-b4b1-76d369ce4146",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "### GRIB: nearest gridpoint"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddb547b2-1d39-4c50-9777-19591bd0f6fa",
   "metadata": {},
   "source": [
    "First we ensure the example file is available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "87e4f240-b145-4c4d-9906-38db9a7b3a99",
   "metadata": {},
   "outputs": [],
   "source": [
    "import earthkit.data as ekd\n",
    "from earthkit.geo import nearest_point_haversine, nearest_point_kdtree, GeoKDTree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8ff42074-f9af-4f1b-a7c3-3d19ae54cf43",
   "metadata": {},
   "outputs": [],
   "source": [
    "ekd.download_example_file(\"test.grib\")\n",
    "ds = ekd.from_source(\"file\", \"test.grib\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9a78d17-7f63-4944-ab3e-22349954e601",
   "metadata": {},
   "source": [
    "#### Using the KDTree method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "159cd5bd-6d6f-4d59-aeb6-c5faaf985cfe",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "We define a reference point and get the index of the nearest gridpoint from the first field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8270dd2a-7ed6-4c7d-95d7-848e31465a9b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([102]), array([218417.94491761]))"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latlon = ds[0].to_latlon()\n",
    "lat = latlon[\"lat\"]\n",
    "lon = latlon[\"lon\"]\n",
    "\n",
    "p_ref = (51.45, -0.97)\n",
    "# distance is in m\n",
    "idx, distance = nearest_point_kdtree(p_ref, (lat, lon))\n",
    "idx, distance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c5ee441-7439-42d2-849b-6825d69d36e2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "With the resulting index we can get the value at the nearest gridpoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d2388b9e-9895-4ab0-be80-69e14baa76c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([280.68066406])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = ds[0].values[idx]\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6db1df33-2b17-4a0b-bd68-bfabe16c544f",
   "metadata": {},
   "source": [
    "The same method works with multiple reference points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e4b12a90-1429-41a1-bf5d-a17b89861b00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([102, 144, 116])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_ref = [\n",
    "    [51.45, 44.49, 50.73], \n",
    "    [-0.97, 18.34, -17.1]\n",
    "]\n",
    "idx, distance = nearest_point_kdtree(p_ref, (lat, lon))\n",
    "idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "f6f28de5-3556-4d34-8c9c-b03c4017c71e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([280.68066406, 296.4765625 , 285.05761719])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = ds[0].values[idx]\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bde851d-f9b3-4cdf-b483-b92ed5a87222",
   "metadata": {},
   "source": [
    "#### Using a KDTree object"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52c393e5-19d2-4505-a577-b6601a32624c",
   "metadata": {},
   "source": [
    "The problem with *nearest_point_kdtree()* is that the KDTree is rebuilt at each call (since it is not cached at the moment). This can be a costly operation for large grids. To overcome this difficulty we can create a *GeoKDTree* object and call the nearest point computations directly on it multiple times. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8cc6bd2a-3cd4-44a1-8049-1d8f715e6d26",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([102, 144, 116]),\n",
       " array([218417.94491761, 120066.38079566, 235683.22276184]))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tree = GeoKDTree(lat, lon)\n",
    "idx, distance = tree.nearest_point(p_ref)\n",
    "idx, distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "080d37b3-1b1b-4c99-b2c6-0ce26fbd4915",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([280.68066406, 296.4765625 , 285.05761719])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = ds[0].values[idx]\n",
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "3c5d0b56-4695-4c7a-a63c-390a7794347f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([144]), array([120066.38079566]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_ref = (44.49,18.34)\n",
    "idx, distance = tree.nearest_point(p_ref)\n",
    "idx, distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "16f872ce-5924-4462-a504-9739a11b64c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([296.4765625])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = ds[0].values[idx]\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e35f59b4-586c-42b5-a24f-5a8c50db98aa",
   "metadata": {},
   "source": [
    "#### Using the Haversine method"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cc90875-8ab5-4a27-b760-db0099fab25d",
   "metadata": {},
   "source": [
    "The nearest gridpoint can also be determined with the Haversine method, which is a \"brute-force\" approach since it computes all the distances between the reference point and the other points using the Haversine distance formula, then finds the points with the shortest distance."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24fd9b22-becd-404d-a240-0bde0d952d16",
   "metadata": {},
   "source": [
    "The *nearest_point_haversine()* method can be used exactly in the same way as *nearest_point_kdtree()*.\n",
    "\n",
    "We define a **reference point** and get the index of the nearest gridpoint from the first field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2578089a-c895-44a0-be5d-e6f5d6f30998",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([102]), array([218417.94491761]))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latlon = ds[0].to_latlon()\n",
    "lat = latlon[\"lat\"]\n",
    "lon = latlon[\"lon\"]\n",
    "\n",
    "p_ref = (51.45, -0.97)\n",
    "# distance is in m\n",
    "idx, distance = nearest_point_haversine(p_ref, (lat, lon))\n",
    "idx, distance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5f1d0e8-41da-4a8a-b3a6-ea8e0c50b420",
   "metadata": {},
   "source": [
    "With the resulting index we can get the value at the nearest gridpoint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e8b8cc37-efab-4870-bc54-20ebb6c8faa7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([280.68066406])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = ds[0].values[idx]\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68f4738a-48f7-488d-8534-631c6422ae27",
   "metadata": {},
   "source": [
    "This time we use 3 reference points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4809ecd5-837e-4d36-b411-cced61c306f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([102, 144, 116])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_ref = [\n",
    "    [51.45, 44.49, 50.73], \n",
    "    [-0.97, 18.34, -17.1]\n",
    "]\n",
    "idx, dist = nearest_point_haversine(p_ref, (lat, lon))\n",
    "idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "eb68fa7b-dc43-4b71-a843-8454ced5b4cb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([280.68066406, 296.4765625 , 285.05761719])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v = ds[0].values[idx]\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb92e03b-3ab7-4470-a423-ad3fca9262cb",
   "metadata": {},
   "source": [
    "### "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev_ecc",
   "language": "python",
   "name": "dev_ecc"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}