notebooks/Interpolate_one_grid_into_another.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0332940b",
   "metadata": {},
   "source": [
    "# Interpolate values of one grid into another one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d5a7877",
   "metadata": {},
   "outputs": [],
   "source": [
    "# solve issue with autocomplete\n",
    "%config Completer.use_jedi = False\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e44900e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geotree import gtree\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "524f47c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiate gtree\n",
    "mytree = gtree()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc23ae73",
   "metadata": {},
   "source": [
    "## Define the first set of points or `base`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a011262",
   "metadata": {},
   "outputs": [],
   "source": [
    "npoints = 100\n",
    "lons = np.random.randint(-180, 180, npoints)\n",
    "lats = np.random.randint(-90, 90, npoints)\n",
    "depths = np.zeros(npoints)\n",
    "\n",
    "# assign some values (to each point)\n",
    "vals = np.zeros(npoints)\n",
    "vals[:(npoints//2)] = 0\n",
    "vals[(npoints//2):] = 1\n",
    "# random:\n",
    "# mu = 0\n",
    "# sigma = 1\n",
    "# vals = sigma*np.random.randn(npoints) + mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55161276",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add lons/lats/depths of the first set of points\n",
    "mytree.add_lonlatdep(lons=lons, \n",
    "                     lats=lats, \n",
    "                     depths=depths)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdba2d9b",
   "metadata": {},
   "source": [
    "## Define queries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47186edc",
   "metadata": {},
   "outputs": [],
   "source": [
    "q_npoints = 10000\n",
    "q_lons = np.random.randint(-180, 180, q_npoints)\n",
    "q_lats = np.random.randint(-90, 90, q_npoints)\n",
    "q_depths = np.zeros(q_npoints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cea3e198",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add lons/lats/depths of queries\n",
    "mytree.add_lonlatdep_query(lons=q_lons, \n",
    "                           lats=q_lats, \n",
    "                           depths=q_depths)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51469b69",
   "metadata": {},
   "source": [
    "## Assign values to the first set of points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fdc201a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# note: size of vals should be the same as the first set of points\n",
    "mytree.add_vals(vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55235242",
   "metadata": {},
   "source": [
    "## Interpolation \n",
    "\n",
    "compute the values of `queries` from the values of `base` points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e9f3ce0",
   "metadata": {},
   "outputs": [],
   "source": [
    "mytree.interpolate(num_neighs=2, method=\"kdtree\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2df814ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.figure(figsize=(15, 7))\n",
    "\n",
    "plt.scatter(q_lons, q_lats, \n",
    "            c=mytree.interp_vals, \n",
    "            marker=\"x\", \n",
    "            vmin=min(vals), vmax=max(vals),\n",
    "            label=\"queries\")\n",
    "\n",
    "plt.scatter(lons, lats,\n",
    "            c=vals, \n",
    "            marker=\"o\",\n",
    "            vmin=min(vals), vmax=max(vals), edgecolors=\"r\",\n",
    "            label=\"base\",\n",
    "            zorder=100)\n",
    "\n",
    "plt.legend(bbox_to_anchor=(0., 1.01, 1., .05), \n",
    "           loc=\"right\", ncol=2, \n",
    "           fontsize=16,\n",
    "           borderaxespad=0.)\n",
    "\n",
    "plt.title(f\"Method: {mytree.interp_method}\", size=16)\n",
    "plt.colorbar()\n",
    "plt.grid()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7725fada",
   "metadata": {},
   "source": [
    "## Project results using `cartopy`\n",
    "\n",
    "`cartopy` needs to be installed, see: https://scitools.org.uk/cartopy/docs/latest/index.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3954c7b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cartopy import crs\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "\n",
    "ax = fig.add_subplot(1, 1, 1, \n",
    "                     projection=crs.InterruptedGoodeHomolosine())\n",
    "ax.coastlines(color=\"black\")\n",
    "\n",
    "# plot the queries\n",
    "plt.scatter(\n",
    "    mytree.lons_q, \n",
    "    mytree.lats_q, \n",
    "    c=mytree.interp_vals,\n",
    "    transform=crs.PlateCarree(), \n",
    "    vmin=min(vals), vmax=max(vals), \n",
    "    marker=\"x\")\n",
    "\n",
    "# plot the base points\n",
    "plt.scatter(\n",
    "    mytree.lons,\n",
    "    mytree.lats, \n",
    "    c=mytree.vals,\n",
    "    transform=crs.PlateCarree(), \n",
    "    vmin=min(vals), vmax=max(vals), \n",
    "    edgecolors=\"r\")\n",
    "\n",
    "\n",
    "plt.title(f\"Method: {mytree.interp_method}\", size=16)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e598e97f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from cartopy import crs\n",
    "\n",
    "fig = plt.figure(figsize=(10, 7))\n",
    "\n",
    "ax = fig.add_subplot(1, 2, 1, \n",
    "                     projection=crs.Orthographic(central_longitude=0.0))\n",
    "ax.coastlines(color=\"black\")\n",
    "\n",
    "# plot the queries\n",
    "plt.scatter(\n",
    "    mytree.lons_q, \n",
    "    mytree.lats_q, \n",
    "    c=mytree.interp_vals,\n",
    "    transform=crs.PlateCarree(), \n",
    "    vmin=min(vals), vmax=max(vals), \n",
    "    marker=\"x\")\n",
    "\n",
    "# plot the base points\n",
    "plt.scatter(\n",
    "    mytree.lons, \n",
    "    mytree.lats, \n",
    "    c=mytree.vals,\n",
    "    transform=crs.PlateCarree(), \n",
    "    vmin=min(vals), vmax=max(vals), \n",
    "    edgecolors=\"r\")\n",
    "\n",
    "plt.title(f\"Method: {mytree.interp_method}\", size=16)\n",
    "plt.grid()\n",
    "\n",
    "# -------------------\n",
    "\n",
    "ax = fig.add_subplot(1, 2, 2, \n",
    "                     projection=crs.Orthographic(central_longitude=180.0))\n",
    "ax.coastlines(color=\"black\")\n",
    "\n",
    "# plot the queries\n",
    "plt.scatter(\n",
    "    mytree.lons_q, \n",
    "    mytree.lats_q, \n",
    "    c=mytree.interp_vals,\n",
    "    transform=crs.PlateCarree(), \n",
    "    vmin=min(vals), vmax=max(vals), \n",
    "    marker=\"x\")\n",
    "\n",
    "# plot the base points\n",
    "plt.scatter(\n",
    "    mytree.lons, \n",
    "    mytree.lats, \n",
    "    c=mytree.vals,\n",
    "    transform=crs.PlateCarree(), \n",
    "    vmin=min(vals), vmax=max(vals), \n",
    "    edgecolors=\"r\")\n",
    "\n",
    "plt.title(f\"Method: {mytree.interp_method}\", size=16)\n",
    "plt.grid()\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a6fe303",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}