{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unstructured grid\n",
    "=================\n",
    "\n",
    "The interpolation of this object is based on a\n",
    "[R*Tree](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.html#pyinterp.RTree)\n",
    "structure. To begin with, we start by building this object. By default, this\n",
    "object considers the WGS-84 geodetic coordinate system. But you can define\n",
    "another one using the class\n",
    "[System](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.geodetic.System.html#pyinterp.geodetic.System)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot\n",
    "import numpy\n",
    "import pyinterp\n",
    "\n",
    "mesh = pyinterp.RTree()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we will insert points into the tree. The class allows you to add points\n",
    "using two algorithms. The first one, called\n",
    "[packing](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.packing.html#pyinterp.RTree.packing),\n",
    "will enable you to enter the values in the tree at once. This mechanism is the\n",
    "recommended solution to create an optimized in-memory structure, both in terms\n",
    "of construction time and queries. When this is not possible, you can insert new\n",
    "information into the tree as you go along using the\n",
    "[insert](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.packing.html#pyinterp.RTree.insert)\n",
    "method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIZE = 2000\n",
    "X0, X1 = 80, 170\n",
    "Y0, Y1 = -45, 30\n",
    "lons = numpy.random.uniform(low=X0, high=X1, size=(SIZE, ))\n",
    "lats = numpy.random.uniform(low=Y0, high=Y1, size=(SIZE, ))\n",
    "data = numpy.random.random(size=(SIZE, ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Populates the search tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh.packing(numpy.vstack((lons, lats)).T, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the tree is created, you can interpolate data with three algorithms:\n",
    "\n",
    "- [Inverse Distance Weighting](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.inverse_distance_weighting.html#pyinterp.RTree.inverse_distance_weighting)\n",
    "    or IDW\n",
    "- [Radial Basis Function](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.radial_basis_function.html#pyinterp.RTree.radial_basis_function)\n",
    "    or RBF\n",
    "- [Window Function](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.radial_basis_function.html#pyinterp.RTree.window_function)\n",
    "\n",
    "Yon can also search the\n",
    "[nearest neighbors](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.RTree.query.html#pyinterp.RTree.query)\n",
    "on the tree.\n",
    "\n",
    "---\n",
    "**Note**\n",
    "\n",
    "When comparing an RBF to IDW, IDW will never predict values higher than\n",
    "the maximum measured value or lower than the minimum measured value.\n",
    "However, RBFs can predict values higher than the maximum values and\n",
    "lower than the minimum measured values.\n",
    "\n",
    "The window function restricts the analyzed data set to a range near the\n",
    "point of interest. The weighting factor decreases the effect of points\n",
    "further away from the interpolated section of the point.\n",
    "\n",
    "---\n",
    "\n",
    "In this example, we will under-sample the source grid at 1/32 degree\n",
    "over an area of the globe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "STEP = 1 / 32\n",
    "mx, my = numpy.meshgrid(numpy.arange(X0, X1 + STEP, STEP),\n",
    "                        numpy.arange(Y0, Y1 + STEP, STEP),\n",
    "                        indexing=\"ij\")\n",
    "\n",
    "idw, neighbors = mesh.inverse_distance_weighting(\n",
    "    numpy.vstack((mx.ravel(), my.ravel())).T,\n",
    "    within=False,  # Extrapolation is forbidden\n",
    "    k=11,  # We are looking for at most 11 neighbors\n",
    "    radius=600000,\n",
    "    num_threads=0)\n",
    "idw = idw.reshape(mx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolation with RBF method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbf, neighbors = mesh.radial_basis_function(\n",
    "    numpy.vstack((mx.ravel(), my.ravel())).T,\n",
    "    within=False,  # Extrapolation is forbidden\n",
    "    k=11,  # We are looking for at most 11 neighbors\n",
    "    radius=600000,\n",
    "    rbf=\"thin_plate\",\n",
    "    num_threads=0)\n",
    "rbf = rbf.reshape(mx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolation with a Window Function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wf, neighbors = mesh.window_function(\n",
    "    numpy.vstack((mx.ravel(), my.ravel())).T,\n",
    "    within=False,  # Extrapolation is forbidden\n",
    "    k=11,\n",
    "    radius=600000,\n",
    "    wf=\"parzen\",\n",
    "    num_threads=0)\n",
    "wf = wf.reshape(mx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize our interpolated data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = matplotlib.pyplot.figure(figsize=(10, 20))\n",
    "ax1 = fig.add_subplot(311)\n",
    "pcm = ax1.pcolormesh(mx, my, idw, cmap='jet', shading='auto', vmin=0, vmax=1)\n",
    "ax1.set_title(\"IDW interpolation\")\n",
    "ax2 = fig.add_subplot(312)\n",
    "pcm = ax2.pcolormesh(mx, my, rbf, cmap='jet', shading='auto', vmin=0, vmax=1)\n",
    "ax2.set_title(\"RBF interpolation\")\n",
    "ax3 = fig.add_subplot(313)\n",
    "pcm = ax3.pcolormesh(mx, my, wf, cmap='jet', shading='auto', vmin=0, vmax=1)\n",
    "ax3.set_title(\"Window function interpolation\")\n",
    "fig.colorbar(pcm, ax=[ax1, ax2, ax3], shrink=0.8)\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "9dbeda838dae5ee2fb6c16b710bbeb61ff7c35dcfc6a6c798a91ef3c0857de81"
  },
  "kernelspec": {
   "display_name": "Python 3.8.12 64-bit ('conda-forge': conda)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}