{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Create interpolator objects\n",
    "=====================\n",
    "\n",
    "In this example, we are going to build the basic objects allowing to\n",
    "carry out interpolations.\n",
    "\n",
    "Before starting, we will examine the properties of a Cartesian grid and\n",
    "the different classes associated with these objects.\n",
    "\n",
    "The first step is to open the NetCDF file and load the data. We use here\n",
    "the NetCDF4 library to detail the different steps, but we will see that\n",
    "we can automate the steps described below using the xarray objects\n",
    "library.\n",
    "\n",
    "Step-by-step creation of grids\n",
    "-----------------------------------"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "import timeit\n",
    "import netCDF4\n",
    "import pandas\n",
    "import pyinterp\n",
    "import pyinterp.backends.xarray\n",
    "import pyinterp.tests\n",
    "import numpy\n",
    "import xarray\n",
    "\n",
    "\n",
    "with netCDF4.Dataset(pyinterp.tests.grid3d_path()) as ds:\n",
    "    lon, lat, time, time_units, tcw = ds.variables[\n",
    "        \"longitude\"][:], ds.variables[\"latitude\"][:], ds.variables[\n",
    "            \"time\"][:], ds.variables[\"time\"].units, ds.variables[\"tcw\"][:]\n",
    "    time = numpy.array(netCDF4.num2date(time, time_units),\n",
    "                       dtype=\"datetime64[us]\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This regular 3-dimensional grid is associated with three axes:\n",
    "\n",
    "-   longitudes,\n",
    "-   latitudes and\n",
    "-   time.\n",
    "\n",
    "To perform the calculations quickly, we will build three objects that\n",
    "will be used by the interpolator to search for the data to be used.\n",
    "Let's start with the y-axis representing the latitude axis."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "y_axis = pyinterp.Axis(lat)\n",
    "y_axis"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For example, you can search for the closest point to 0.12 degrees north\n",
    "latitude."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "y_axis.find_index([0.12])"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, the x-axis representing the longitudinal axis. In this case, the\n",
    "axis is an axis representing a 360 degree circle."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "x_axis = pyinterp.Axis(lon, is_circle=True)\n",
    "x_axis"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The values -180 and 180 degrees represent the same point on the axis."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "x_axis.find_index([-180]) == x_axis.find_index([180])"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Finally, we create the time axis"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "t_axis = pyinterp.TemporalAxis(time)\n",
    "t_axis"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As these objects must communicate in C++ memory space, we use objects\n",
    "specific to the library much faster than other data models and manage\n",
    "the axes representing a circle. For example if we compare these objects\n",
    "to Pandas indexes:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "values = lon[10:20] + 1 / 3\n",
    "index = pandas.Index(lon)\n",
    "print(\"pandas.Index: %f\" % timeit.timeit(\n",
    "    \"index.searchsorted(values)\", globals=dict(index=index, values=values)))\n",
    "print(\"pyinterp.Axis %f\" % timeit.timeit(\n",
    "    \"x_axis.find_index(values)\", globals=dict(x_axis=x_axis, values=values)))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This time axis is also very efficient compared to the pandas index."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "index = pandas.Index(time)\n",
    "values = time + numpy.timedelta64(1, \"ns\")\n",
    "print(\"pandas.Index: %f\" % timeit.timeit(\n",
    "    \"index.searchsorted(values)\", globals=dict(index=index, values=values)))\n",
    "print(\"pyinterp.Axis %f\" % timeit.timeit(\n",
    "    \"t_axis.find_index(values)\", globals=dict(t_axis=t_axis, values=values)))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Before constructing the tensor for pyinterp, we must begin to organize\n",
    "the tensor data so that it is properly stored in memory for pyinterp.\n",
    "-   The shape of the tensor must be (len(x_axis), len(y_axis),\n",
    "    len(t_axis))"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "tcw = tcw.T"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "-   The undefined values must be set to nan."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "tcw[tcw.mask] = float(\"nan\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we can build the object handling the regular 3-dimensional grid.\n",
    "\n",
    "---\n",
    "**Note**\n",
    "\n",
    "Grid data are not copied, the Grid3D class just keeps a reference on the\n",
    "handled array. Axis data are copied for non-uniform axes, and only\n",
    "examined for regular axes.\n",
    "\n",
    "---"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "grid_3d = pyinterp.Grid3D(x_axis, y_axis, t_axis, tcw)\n",
    "grid_3d"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "xarray\n",
    "-------\n",
    "\n",
    "The construction of these objects manipulating the\n",
    "[regular grids](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.backends.xarray.RegularGridInterpolator.html#pyinterp.backends.xarray.RegularGridInterpolator)\n",
    "can be done more easily using the [xarray](http://xarray.pydata.org/) library\n",
    "and [CF](https://cfconventions.org/) convention usually found in NetCDF files."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "interpolator = pyinterp.backends.xarray.RegularGridInterpolator(\n",
    "    xarray.open_dataset(pyinterp.tests.grid3d_path()).tcw)\n",
    "interpolator.grid"
   ],
   "outputs": [],
   "metadata": {}
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}