notebooks/interpolate_sphere.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test spherical spatial interpolators using Franke 3D evaluation function  \n",
    "### Test functions\n",
    "- http://www.sfu.ca/~ssurjano/franke2d.html  \n",
    "- http://www.sciencedirect.com/science/article/pii/S037704270100485X  \n",
    "\n",
    "### QR Factorization (Fornberg)\n",
    "- http://epubs.siam.org/doi/abs/10.1137/060671991  \n",
    "- http://epubs.siam.org/doi/abs/10.1137/09076756X  \n",
    "\n",
    "### Initial nodes\n",
    "- http://web.maths.unsw.edu.au/~rsw/Sphere/  \n",
    "- http://math.boisestate.edu/~wright/montestigliano/nodes.zip  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.io\n",
    "import numpy as np\n",
    "import spatial_interpolators as spi\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Franke's 3D evaluation function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def franke(x,y,z):\n",
    "\tF1 = 0.75*np.exp(-((9.*x-2.)**2 + (9.*y-2.)**2 + (9.0*z-2.)**2)/4.)\n",
    "\tF2 = 0.75*np.exp(-((9.*x+1.)**2/49. + (9.*y+1.)/10. + (9.0*z+1.)/10.))\n",
    "\tF3 = 0.5*np.exp(-((9.*x-7.)**2 + (9.*y-3.)**2 + (9.*z-5)**2)/4.)\n",
    "\tF4 = 0.2*np.exp(-((9.*x-4.)**2 + (9.*y-7.)**2 + (9.*z-5.)**2))\n",
    "\tF = F1 + F2 + F3 - F4\n",
    "\treturn F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Franke's evaluation function at nodal points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using max_determinant nodal points from\n",
    "# http://math.boisestate.edu/~wright/montestigliano/nodes.zip\n",
    "N = 324\n",
    "xd = scipy.io.loadmat('md{0:05d}.mat'.format(N))\n",
    "x,y,z = xd['x'][:,0],xd['x'][:,1],xd['x'][:,2]\n",
    "# compute functional values at nodes\n",
    "f = franke(x,y,z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Franke's evaluation function at grid points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert node coordinates to lat/lon\n",
    "lon,lat,_ = spi.spatial.to_sphere(x,y,z)\n",
    "\n",
    "# calculate output points (standard lat/lon grid)\n",
    "dlon = 5.0\n",
    "dlat = 5.0\n",
    "gridlon = np.arange(0,360+dlon,dlon)\n",
    "gridlat = np.arange(90,-90-dlat,-dlat)\n",
    "LON,LAT = np.meshgrid(gridlon,gridlat,indexing='ij')\n",
    "x,y,z = spi.spatial.to_cartesian(LON,LAT,a_axis=1.0,flat=0.0)\n",
    "# calculate functional values at output points\n",
    "FI = franke(x, y, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interpolate to spherical grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interpolate with radial basis functions\n",
    "m = 'gaussian'\n",
    "RBF = spi.sph_radial_basis(lon,lat,f,LON,LAT,method=m,smooth=0.0001,\n",
    "\tepsilon=1.0)\n",
    "QR = spi.sph_radial_basis(lon,lat,f,LON,LAT,method=m,epsilon=0.9,QR=True)\n",
    "# interpolate with spherical splines\n",
    "t = 2.\n",
    "SPL = spi.sph_spline(lon,lat,f,LON,LAT,tension=t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create output plot showing interpolation and points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot interpolation and real-interpolated\n",
    "fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(num=1, ncols=2, nrows=2,\n",
    "\tsharex=True, sharey=True, figsize=(9,4.9))\n",
    "cmap = plt.cm.get_cmap('Spectral_r').copy()\n",
    "cmap.set_bad('w',0.)\n",
    "\n",
    "ax1.scatter(lon, lat, c=f, vmin=f.min(), vmax=f.max(),\n",
    "\tcmap=cmap, edgecolors='w')\n",
    "ax1.imshow(FI.transpose(),extent=(0,360,-90,90),\n",
    "\tvmin=f.min(), vmax=f.max(), cmap=cmap)\n",
    "ax2.imshow(RBF.transpose(),extent=(0,360,-90,90),\n",
    "\tvmin=f.min(), vmax=f.max(), cmap=cmap)\n",
    "ax2.scatter(lon, lat, c=f, vmin=f.min(), vmax=f.max(),\n",
    "\tcmap=cmap, edgecolors='none')\n",
    "ax3.imshow(QR.transpose(),extent=(0,360,-90,90),\n",
    "\tvmin=f.min(), vmax=f.max(), cmap=cmap)\n",
    "ax3.scatter(lon, lat, c=f, vmin=f.min(), vmax=f.max(),\n",
    "\tcmap=cmap, edgecolors='none')\n",
    "ax4.imshow(SPL.transpose(),extent=(0,360,-90,90),\n",
    "\tvmin=f.min(), vmax=f.max(), cmap=cmap)\n",
    "ax4.scatter(lon, lat, c=f, vmin=f.min(), vmax=f.max(),\n",
    "\tcmap=cmap, edgecolors='none')\n",
    "\n",
    "# for each axis\n",
    "for ax in [ax1,ax2,ax3,ax4]:\n",
    "    # no ticks on the x and y axes\n",
    "\tax.get_xaxis().set_ticks([])\n",
    "\tax.get_yaxis().set_ticks([])\n",
    "\t# set x and y limits (global)\n",
    "\tax1.set_xlim(0,360)\n",
    "\tax1.set_ylim(-90,90)\n",
    "\n",
    "# set titles\n",
    "ax1.set_title('Franke Evaluation Function')\n",
    "ax2.set_title('RBF {0}'.format(m.capitalize()))\n",
    "ax3.set_title('RBF {0} with QR Factorization'.format(m.capitalize()))\n",
    "ax4.set_title('Spline with Tension {0:0.0f}'.format(t*100))\n",
    "# subplot adjustments\n",
    "fig.subplots_adjust(left=0.02,right=0.98,bottom=0.01,top=0.95,\n",
    "\twspace=0.02,hspace=0.11)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.6 64-bit",
   "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}