{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2D interpolation\n",
    "=============\n",
    "\n",
    "Interpolation of a two-dimensional regular grid.\n",
    "\n",
    "Bivariate\n",
    "----------\n",
    "\n",
    "Perform a\n",
    "[bivariate](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.bivariate.html#pyinterp.bivariate)\n",
    "interpolation of gridded data points.\n",
    "\n",
    "The distribution contains a 2D field `mss.nc` that will be used in this\n",
    "help. This file is located in the `src/pyinterp/tests/dataset` directory at the\n",
    "root of the project.\n",
    "\n",
    "---\n",
    "**Warning**\n",
    "\n",
    "This file is an old version of the sub-sampled quarter step MSS CNES/CLS. Please do not use it for scientific purposes, download the latest updated high-resolution version instead\n",
    "[here](https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/mss.html).\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy.crs\n",
    "import matplotlib\n",
    "import matplotlib.pyplot\n",
    "import numpy\n",
    "import pyinterp\n",
    "import pyinterp.backends.xarray\n",
    "import pyinterp.tests\n",
    "import xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to load the data into memory and create the\n",
    "interpolator object:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xarray.open_dataset(pyinterp.tests.grid2d_path())\n",
    "interpolator = pyinterp.backends.xarray.Grid2D(ds.mss)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will then build the coordinates on which we want to interpolate our\n",
    "grid:\n",
    "\n",
    "---\n",
    "**Note**\n",
    "\n",
    "The coordinates used for interpolation are shifted to avoid using the\n",
    "points of the bivariate function.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mx, my = numpy.meshgrid(numpy.arange(-180, 180, 1) + 1 / 3.0,\n",
    "                        numpy.arange(-89, 89, 1) + 1 / 3.0,\n",
    "                        indexing='ij')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grid is\n",
    "[interpolated](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.backends.xarray.Grid2D.bivariate.html#pyinterp.backends.xarray.Grid2D.bivariate)\n",
    "to the desired coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mss = interpolator.bivariate(\n",
    "    coords=dict(lon=mx.ravel(), lat=my.ravel())).reshape(mx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the original grid and the result of the interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = matplotlib.pyplot.figure(figsize=(10, 8))\n",
    "ax1 = fig.add_subplot(\n",
    "    211, projection=cartopy.crs.PlateCarree(central_longitude=180))\n",
    "lons, lats = numpy.meshgrid(ds.lon, ds.lat, indexing='ij')\n",
    "pcm = ax1.pcolormesh(lons,\n",
    "                     lats,\n",
    "                     ds.mss.T,\n",
    "                     cmap='jet',\n",
    "                     shading='auto',\n",
    "                     transform=cartopy.crs.PlateCarree(),\n",
    "                     vmin=-0.1,\n",
    "                     vmax=0.1)\n",
    "ax1.coastlines()\n",
    "ax1.set_title(\"Original MSS\")\n",
    "ax2 = fig.add_subplot(212, projection=cartopy.crs.PlateCarree())\n",
    "pcm = ax2.pcolormesh(mx,\n",
    "                     my,\n",
    "                     mss,\n",
    "                     cmap='jet',\n",
    "                     shading='auto',\n",
    "                     transform=cartopy.crs.PlateCarree(),\n",
    "                     vmin=-0.1,\n",
    "                     vmax=0.1)\n",
    "ax2.coastlines()\n",
    "ax2.set_title(\"Bilinear Interpolated MSS\")\n",
    "fig.colorbar(pcm, ax=[ax1, ax2], shrink=0.8)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Values can be interpolated with several methods: *bilinear*, *nearest*,\n",
    "and *inverse distance weighting*. Distance calculations, if necessary,\n",
    "are calculated using the [Haversine\n",
    "formula](https://en.wikipedia.org/wiki/Haversine_formula).\n",
    "\n",
    "Bicubic\n",
    "---------\n",
    "\n",
    "To interpolate data points on a regular two-dimensional grid. The\n",
    "interpolated surface is smoother than the corresponding surfaces\n",
    "obtained by bilinear interpolation. Spline functions provided by\n",
    "[GSL](https://www.gnu.org/software/gsl/) achieve bicubic interpolation.\n",
    "\n",
    "---\n",
    "**Warning**\n",
    "\n",
    "\n",
    "When using this interpolator, pay attention to the undefined values.\n",
    "Because as long as the calculation window uses an indefinite point, the\n",
    "interpolator will compute indeterminate values. In other words, this\n",
    "interpolator increases the area covered by the masked values. To avoid\n",
    "this behavior, it is necessary to\n",
    "[pre-process](https://pangeo-pyinterp.readthedocs.io/en/latest/auto_examples/ex_fill_undef.html) the grid to\n",
    "delete undefined values.\n",
    "\n",
    "---\n",
    "\n",
    "The interpolation\n",
    "[bicubic](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.backends.xarray.Grid2D.bicubic.html#pyinterp.backends.xarray.Grid2D.bicubic)\n",
    "function has more parameters to define the data frame used\n",
    "by the spline functions and how to process the edges of the regional\n",
    "grids:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mss = interpolator.bicubic(coords=dict(lon=mx.ravel(), lat=my.ravel()),\n",
    "                           nx=3,\n",
    "                           ny=3).reshape(mx.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "**Warning**\n",
    "\n",
    "\n",
    "The grid provided must have strictly increasing axes to meet the\n",
    "specifications of the GSL library. When building the grid, specify the\n",
    "[increasing_axes](https://pangeo-pyinterp.readthedocs.io/en/latest/generated/pyinterp.grid.Grid2D.html#pyinterp.grid.Grid2D)\n",
    "option to flip the decreasing axes and the grid automatically. For example:\n",
    "\n",
    "```python\n",
    "interpolator = pyinterp.backends.xarray.Grid2D(\n",
    "    ds.mss, increasing_axes=True)\n",
    "```\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = matplotlib.pyplot.figure(figsize=(10, 8))\n",
    "ax1 = fig.add_subplot(\n",
    "    211, projection=cartopy.crs.PlateCarree(central_longitude=180))\n",
    "pcm = ax1.pcolormesh(lons,\n",
    "                     lats,\n",
    "                     ds.mss.T,\n",
    "                     cmap='jet',\n",
    "                     shading='auto',\n",
    "                     transform=cartopy.crs.PlateCarree(),\n",
    "                     vmin=-0.1,\n",
    "                     vmax=0.1)\n",
    "ax1.coastlines()\n",
    "ax1.set_title(\"Original MSS\")\n",
    "ax2 = fig.add_subplot(212, projection=cartopy.crs.PlateCarree())\n",
    "pcm = ax2.pcolormesh(mx,\n",
    "                     my,\n",
    "                     mss,\n",
    "                     cmap='jet',\n",
    "                     shading='auto',\n",
    "                     transform=cartopy.crs.PlateCarree(),\n",
    "                     vmin=-0.1,\n",
    "                     vmax=0.1)\n",
    "ax2.coastlines()\n",
    "ax2.set_title(\"Bicubic Interpolated MSS\")\n",
    "fig.colorbar(pcm, ax=[ax1, ax2], shrink=0.8)\n",
    "fig.show()"
   ]
  }
 ],
 "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
}