Approximating metrics on a spherical grid.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spherical coordinates\n",
    "\n",
    "Latitude, $\\phi$, and longitude, $\\lambda$, determine position on the surface of a sphere of radius, $r$.\n",
    "\n",
    "The Cartesian coordinates $(x,y,z)$ of position $(\\lambda,\\phi)$ are\n",
    "\n",
    "\\begin{align}\n",
    "x &= r \\cos{\\lambda} \\cos{\\phi} \\\\\n",
    "y &= r \\sin{\\lambda} \\cos{\\phi} \\\\\n",
    "z &= r \\sin{\\phi}\n",
    "\\end{align}\n",
    "\n",
    "The spherical coordinate $(\\lambda,\\phi)$ in terms Cartesian coordinates $(x,y,z)$ are\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda &= \\tan^{-1}{ \\frac{y}{x} } \\\\\n",
    "\\phi &= \\sin^{-1}{ \\frac{z}{r} }\n",
    "\\end{align}\n",
    "\n",
    "The line element is given by\n",
    "\n",
    "\\begin{align}\n",
    "ds^2 &= h_\\lambda^2 d\\lambda^2 + h_\\phi^2 d\\phi^2\n",
    "\\end{align}\n",
    "\n",
    "where the scale factors are\n",
    "\n",
    "\\begin{align}\n",
    "h_\\lambda &= \\sqrt{\n",
    "\\left( \\frac{\\partial x}{\\partial \\lambda} \\right)^2\n",
    "+ \\left( \\frac{\\partial y}{\\partial \\lambda} \\right)^2\n",
    "+ \\left( \\frac{\\partial z}{\\partial \\lambda} \\right)^2 }\n",
    "\\\\\n",
    "h_\\phi &= \\sqrt{\n",
    "\\left( \\frac{\\partial x}{\\partial \\phi} \\right)^2\n",
    "+ \\left( \\frac{\\partial y}{\\partial \\phi} \\right)^2\n",
    "+ \\left( \\frac{\\partial z}{\\partial \\phi} \\right)^2 }\n",
    "\\end{align}\n",
    "\n",
    "The scale factors for the spherical projection simplify to\n",
    "\n",
    "\\begin{align}\n",
    "h_\\lambda &= r \\cos{\\phi} \\\\\n",
    "h_\\phi &= r\n",
    "\\end{align}\n",
    "so the line element is\n",
    "$$\n",
    "d{\\mathbf s} = \\hat{\\boldsymbol \\lambda} r \\cos{\\phi} d\\lambda + \\hat{\\boldsymbol \\phi} r d\\phi\n",
    "$$\n",
    "where the unit vectors are\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{\\boldsymbol \\lambda} &=\n",
    "\\left( \\begin{array}{c}\n",
    "- \\sin{\\lambda} \\\\ \\cos{\\lambda} \\\\ 0\n",
    "\\end{array} \\right)\n",
    "\\\\\n",
    "\\hat{\\boldsymbol \\phi} &=\n",
    "\\left( \\begin{array}{c}\n",
    "\\cos{\\lambda} \\sin{\\phi} \\\\ \\sin{\\lambda} \\sin{\\phi} \\\\ \\cos{\\phi}\n",
    "\\end{array} \\right)\n",
    "\\end{align}\n",
    "\n",
    "The area element is\n",
    "$$\n",
    "d{\\mathbf a} = r^2 \\cos{\\phi} d\\lambda d\\phi\n",
    "$$\n",
    "\n",
    "A line segment of constant longitude has length $\\Delta s = \\int h_\\phi d\\phi = r \\Delta \\phi$ where $\\Delta \\phi$ is the change in latitude for the segment.\n",
    "\n",
    "A line segment of constant latitude has length $\\Delta s = \\int h_\\lambda d\\lambda = ( r \\cos{\\phi} ) \\Delta \\lambda$ where $\\Delta \\lambda$ is the change in longitude for the segment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A spherical grid\n",
    "\n",
    "A spherical grid is a mesh made from lines of constant longitude uniformly spaced by $\\Delta \\lambda$ and lines of constant latitude spaced by $\\Delta \\phi$.\n",
    "\n",
    "The mesh parameters are the starting longitude, $\\lambda_0$, starting latitude, $\\phi_0$ and the spacings $\\Delta \\lambda$ and $\\Delta \\phi$.\n",
    "\n",
    "For $ni \\times nj$ cells there are $(ni+1)(nj+1)$ nodes.\n",
    "\n",
    "Positions of the nodes are\n",
    "\n",
    "\\begin{align}\n",
    "\\lambda_i = \\lambda_0 + i \\Delta \\lambda \\;\\;\\; i = 0, 1, \\ldots, ni \\\\\n",
    "\\phi_j = \\phi_0 + j \\Delta \\phi \\;\\;\\; j = 0, 1, \\ldots, nj\n",
    "\\end{align}\n",
    "\n",
    "The mesh element lengths are\n",
    "\n",
    "\\begin{align}\n",
    "\\Delta y_{i,j+1/2} = r \\left( \\phi_{j+1} - \\phi_j \\right) = r \\Delta \\phi\n",
    "& \\;\\; \\forall \\;\\; i = 0, 1, \\ldots, ni \\;\\; \\textrm{and} \\;\\; j = 0, 1, \\ldots, nj-1 \\\\\n",
    "\\Delta x_{i+1/2,j} = r \\cos{\\phi_j} \\left( \\lambda_{i+1} - \\lambda_i \\right) = r \\cos{\\phi_j} \\Delta \\lambda\n",
    "& \\;\\; \\forall \\;\\; i = 0, 1, \\ldots, ni-1 \\;\\; \\textrm{and} \\;\\; j = 0, 1, \\ldots, nj\n",
    "\\end{align}\n",
    "\n",
    "The cell areas are\n",
    "$$\n",
    "A_{i+1/2,j+1/2} = r^2 \\left( \\sin{\\phi_{j+1}} - \\sin{\\phi_j} \\right) \\Delta \\phi \\Delta \\lambda\n",
    "\\;\\; \\forall \\;\\; i = 0, 1, \\ldots, ni-1 \\;\\; \\textrm{and} \\;\\; j = 0, 1, \\ldots, nj-1\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical approximation of the scale factors\n",
    "\n",
    "Since $\\lambda$,$\\phi$ are orthogonal coordinates we can measure the scale factors\n",
    "\n",
    "\\begin{align}\n",
    "h_\\lambda &= \\left. \\frac{\\partial s}{\\partial \\lambda} \\right|_\\phi\n",
    "\\\\\n",
    "h_\\phi &= \\left. \\frac{\\partial s}{\\partial \\phi} \\right|_\\lambda\n",
    "\\end{align}\n",
    "\n",
    "The derivatives can be approximated numerically. For example, using second order centered differences\n",
    "\n",
    "\\begin{align}\n",
    "h_\\lambda &= \\left. \\frac{\\partial s}{\\partial \\lambda} \\right|_\\phi\n",
    "\\\\\n",
    "&= \\left. \\frac{\\partial s}{\\partial i} \\right|_\\phi \\left. \\frac{\\partial i}{\\partial \\lambda} \\right|_\\phi\n",
    "\\\\\n",
    "&\\approx \\frac{ \\left( \\frac{ s(i+\\epsilon,j) - s(i-\\epsilon,j) }{ 2\\epsilon } + O(\\epsilon^2) \\right)\n",
    "}{ \\left( \\frac{ \\lambda(i+\\epsilon,j) - \\lambda(i-\\epsilon,j) }{ 2\\epsilon } + O(\\epsilon^2) \\right) }\n",
    "\\\\\n",
    "&\\approx \\frac{ \\left( s(i+\\epsilon,j) - s(i-\\epsilon,j) \\right)\n",
    "}{ \\left( \\lambda(i+\\epsilon,j) - \\lambda(i-\\epsilon,j) \\right) } + O(\\epsilon^2)\n",
    "\\\\\n",
    "&\\approx \\frac{\n",
    "\\left( s(i+\\epsilon,j) - s(i-\\epsilon,j) \\right)\n",
    "}{ 2 \\epsilon \\Delta \\lambda } + O(\\epsilon^2)\n",
    "\\end{align}\n",
    "\n",
    "For small $\\epsilon$ the distance $s(i+\\epsilon,j) - s(i-\\epsilon,j)$ can be approximated by the great arc distance. Higher order differences can be substituted for accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SG:\n",
    "    def __init__(self, lon0, lat0, dlon, dlat, r=1):\n",
    "        \"\"\"A spherical grid starting at (lon0,lat0) with spacing (dlon,dlat)\n",
    "        given in degrees. r=1 by default.\"\"\"\n",
    "        self.r = r\n",
    "        self.lon0, self.lat0 = lon0, lat0\n",
    "        self.dlon, self.dlat = dlon, dlat\n",
    "        deg2rad = numpy.pi / 180.\n",
    "        self.dlam, self.dphi = deg2rad * dlon, deg2rad * dlat\n",
    "    def lon(self, j, i):\n",
    "        \"\"\"Returns longitude of nodes in degrees\"\"\"\n",
    "        return self.lon0 + self.dlon * numpy.array(i)\n",
    "    def lat(self, j, i):\n",
    "        \"\"\"Returns latitude of nodes in degrees\"\"\"\n",
    "        return self.lat0 + self.dlat * numpy.array(j)\n",
    "    def lam(self, j, i):\n",
    "        \"\"\"Returns longitude of nodes in radians\"\"\"\n",
    "        deg2rad = numpy.pi / 180.\n",
    "        return self.lon(j, i) * deg2rad\n",
    "    def phi(self, j, i):\n",
    "        \"\"\"Returns latitude of nodes in radians\"\"\"\n",
    "        deg2rad = numpy.pi / 180.\n",
    "        return self.lat(j, i) * deg2rad\n",
    "    def dy(self, j, i):\n",
    "        \"\"\"Returns length of latitude segments\"\"\"\n",
    "        phi = self.phi(j, i)\n",
    "        return self.r * ( phi[1:,:] - phi[:-1,:] )\n",
    "    def dx(self, j, i):\n",
    "        \"\"\"Returns length of longitude segments\"\"\"\n",
    "        phi, lam = self.phi(j, i), self.lam(j, i)\n",
    "        phi = 0.5 * ( phi[:,1:] + phi[:,:-1] )\n",
    "        dlam = lam[:,1:] - lam[:,:-1]\n",
    "        return self.r * numpy.cos(phi) * dlam\n",
    "    def area(self, j, i):\n",
    "        \"\"\"Returns areas of cells\"\"\"\n",
    "        phi, lam = self.phi(j, i), self.lam(j, i)\n",
    "        phi = 0.5 * ( phi[:,1:] + phi[:,:-1] )\n",
    "        dphi = phi[1:,:] - phi[:-1,:]\n",
    "        lam = 0.5 * ( lam[1:,:] + lam[:-1,:] )\n",
    "        dlam = lam[:,1:] - lam[:,:-1]\n",
    "        return self.r**2 * ( numpy.sin(phi[1:,:]) - numpy.sin(phi[:-1,:]) ) * dlam\n",
    "    def hi(self, j, i):\n",
    "        \"\"\"Returns scale factor h_lambda\"\"\"\n",
    "        phi = self.phi(j,i)\n",
    "        return self.r * numpy.cos( phi )\n",
    "    def hj(self, j, i):\n",
    "        \"\"\"Returns scale factor h_phi\"\"\"\n",
    "        phi = self.phi(j,i)\n",
    "        return self.r + 0 * phi\n",
    "    def great_arc_distance(self, j0, i0, j1, i1):\n",
    "        \"\"\"Returns great arc distance between nodes (j0,i0) and (j1,i1)\"\"\"\n",
    "        # https://en.wikipedia.org/wiki/Great-circle_distance\n",
    "        phi0, lam0 = self.phi(j0,i0), self.lam(j0,i0)\n",
    "        phi1, lam1 = self.phi(j1,i1), self.lam(j1,i1)\n",
    "        dphi, dlam = phi1 - phi0, lam1 - lam0\n",
    "        # Haversine formula\n",
    "        d = numpy.sin( 0.5 * dphi)**2 + numpy.sin( 0.5 * dlam)**2 * numpy.cos(phi0) * numpy.cos(phi1)\n",
    "        return self.r * 2. * numpy.arcsin( numpy.sqrt( d ) )\n",
    "    def numerical_hi(self, j, i, eps, order=6):\n",
    "        \"\"\"Returns a numerical approximation to h_lambda\"\"\"\n",
    "        reps = 1. / ( self.dlam * eps )\n",
    "        ds2 = self.great_arc_distance(j, i+eps, j, i-eps)\n",
    "        if order == 2: return 0.5 * ds2 * reps\n",
    "        ds4 = self.great_arc_distance(j, i+2.*eps, j, i-2.*eps)\n",
    "        if order == 4: return ( 8. * ds2 - ds4 ) * (1./12.) * reps\n",
    "        ds6 = self.great_arc_distance(j, i+3.*eps, j, i-3.*eps)\n",
    "        if order == 6: return ( 45. * ds2 - 9. * ds4 + ds6 ) * (1./60.) * reps\n",
    "        raise Exception('order not coded')\n",
    "    def numerical_hj(self, j, i, eps, order=6):\n",
    "        \"\"\"Returns a numerical approximation to h_phi\"\"\"\n",
    "        reps = 1. / ( self.dphi * eps )\n",
    "        ds2 = self.great_arc_distance(j+eps, i, j-eps, i)\n",
    "        if order == 2: return 0.5 * ds2 * reps\n",
    "        ds4 = self.great_arc_distance(j+2.*eps, i, j-2.*eps, i)\n",
    "        if order == 4: return ( 8. * ds2 - ds4 ) * (1./12.) * reps\n",
    "        ds6 = self.great_arc_distance(j+3.*eps, i, j-3.*eps, i)\n",
    "        if order == 6: return ( 45. * ds2 - 9. * ds4 + ds6 ) * (1./60.) * reps\n",
    "        raise Exception('order not coded')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lon = [[  0  60 120 180 240 300 360]\n",
      " [  0  60 120 180 240 300 360]\n",
      " [  0  60 120 180 240 300 360]\n",
      " [  0  60 120 180 240 300 360]\n",
      " [  0  60 120 180 240 300 360]]\n",
      "lat = [[  0   0   0   0   0   0   0]\n",
      " [ 45  45  45  45  45  45  45]\n",
      " [ 90  90  90  90  90  90  90]\n",
      " [135 135 135 135 135 135 135]\n",
      " [180 180 180 180 180 180 180]]\n",
      "dx / pi = [[ 3.33333333e-01  3.33333333e-01  3.33333333e-01  3.33333333e-01\n",
      "   3.33333333e-01  3.33333333e-01]\n",
      " [ 2.35702260e-01  2.35702260e-01  2.35702260e-01  2.35702260e-01\n",
      "   2.35702260e-01  2.35702260e-01]\n",
      " [ 2.04107800e-17  2.04107800e-17  2.04107800e-17  2.04107800e-17\n",
      "   2.04107800e-17  2.04107800e-17]\n",
      " [-2.35702260e-01 -2.35702260e-01 -2.35702260e-01 -2.35702260e-01\n",
      "  -2.35702260e-01 -2.35702260e-01]\n",
      " [-3.33333333e-01 -3.33333333e-01 -3.33333333e-01 -3.33333333e-01\n",
      "  -3.33333333e-01 -3.33333333e-01]]\n",
      "dy / pi = [[0.25 0.25 0.25 0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25 0.25 0.25 0.25]\n",
      " [0.25 0.25 0.25 0.25 0.25 0.25 0.25]]\n",
      "area / pi = [[ 0.23570226  0.23570226  0.23570226  0.23570226  0.23570226  0.23570226]\n",
      " [ 0.09763107  0.09763107  0.09763107  0.09763107  0.09763107  0.09763107]\n",
      " [-0.09763107 -0.09763107 -0.09763107 -0.09763107 -0.09763107 -0.09763107]\n",
      " [-0.23570226 -0.23570226 -0.23570226 -0.23570226 -0.23570226 -0.23570226]]\n",
      "hi = [[ 1.00000000e+00  1.00000000e+00  1.00000000e+00  1.00000000e+00\n",
      "   1.00000000e+00  1.00000000e+00  1.00000000e+00]\n",
      " [ 7.07106781e-01  7.07106781e-01  7.07106781e-01  7.07106781e-01\n",
      "   7.07106781e-01  7.07106781e-01  7.07106781e-01]\n",
      " [ 6.12323400e-17  6.12323400e-17  6.12323400e-17  6.12323400e-17\n",
      "   6.12323400e-17  6.12323400e-17  6.12323400e-17]\n",
      " [-7.07106781e-01 -7.07106781e-01 -7.07106781e-01 -7.07106781e-01\n",
      "  -7.07106781e-01 -7.07106781e-01 -7.07106781e-01]\n",
      " [-1.00000000e+00 -1.00000000e+00 -1.00000000e+00 -1.00000000e+00\n",
      "  -1.00000000e+00 -1.00000000e+00 -1.00000000e+00]]\n",
      "hj = [[1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]]\n"
     ]
    }
   ],
   "source": [
    "g = SG(0, 0, 60, 45)\n",
    "i,j = numpy.meshgrid( numpy.arange(7), numpy.arange(5) ) # 6 x 4 cells\n",
    "print( 'lon =', g.lon(j,i) )\n",
    "print( 'lat =', g.lat(j,i) )\n",
    "print( 'dx / pi =', g.dx(j,i) / numpy.pi)\n",
    "print( 'dy / pi =', g.dy(j,i) / numpy.pi )\n",
    "print( 'area / pi =', g.area(j,i) / numpy.pi )\n",
    "print( 'hi =', g.hi(j,i) )\n",
    "print( 'hj =', g.hj(j,i) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total area / 4 pi = 3.533949646070574e-17\n",
      "Lower hemisphere area / 2 pi = 0.9999999999999999\n",
      "Western hemisphere area / 2 pi = 1.766974823035287e-17\n",
      "Pole-to-pole distance / pi = 1.0\n",
      "Length of equator / 2 pi = 6.123233995736766e-17\n"
     ]
    }
   ],
   "source": [
    "print( 'Total area / 4 pi =', g.area(j,i).sum() / (4*numpy.pi) )\n",
    "print( 'Lower hemisphere area / 2 pi =',g.area(j,i)[:2,:].sum() / (2*numpy.pi) )\n",
    "print( 'Western hemisphere area / 2 pi =',g.area(j,i)[:,:3].sum() / (2*numpy.pi) )\n",
    "print( 'Pole-to-pole distance / pi =',g.dy(j,i)[:,0].sum() / numpy.pi )\n",
    "print( 'Length of equator / 2 pi =',g.dx(j,i)[2,:].sum() / (2*numpy.pi) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]\n",
      " [1. 1. 1. 1. 1. 1. 1.]]\n",
      "[[ 1.          1.          1.          1.          1.          1.        ]\n",
      " [ 0.97603415  0.97603415  0.97603415  0.97603415  0.97603415  0.97603415]\n",
      " [ 0.95492966  0.95492966  0.95492966  0.95492966  0.95492966  0.95492966]\n",
      " [-0.97603415 -0.97603415 -0.97603415 -0.97603415 -0.97603415 -0.97603415]\n",
      " [-1.         -1.         -1.         -1.         -1.         -1.        ]]\n"
     ]
    }
   ],
   "source": [
    "# Check great arc distance == dy\n",
    "print( g.great_arc_distance(j[:-1],i[:-1],j[1:],i[1:]) / g.dy(j,i) )\n",
    "# Check great arc distance <= dx\n",
    "print( g.great_arc_distance(j[:,:-1],i[:,:-1],j[:,1:],i[:,1:]) / g.dx(j,i) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  1.,  1.,  1.,  1.,  1.,  1.],\n",
       "       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.],\n",
       "       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.],\n",
       "       [-1., -1., -1., -1., -1., -1., -1.],\n",
       "       [-1., -1., -1., -1., -1., -1., -1.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compare numerical scale factor hi to true scale factor hi\n",
    "g.numerical_hi(j,i,1e-2,order=6) / g.hi(j,i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1., 1., 1., 1., 1., 1.],\n",
       "       [1., 1., 1., 1., 1., 1., 1.],\n",
       "       [1., 1., 1., 1., 1., 1., 1.],\n",
       "       [1., 1., 1., 1., 1., 1., 1.],\n",
       "       [1., 1., 1., 1., 1., 1., 1.]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compare numerical scale factor hj to true scale factor hj\n",
    "g.numerical_hj(j,i,1e-2,order=6) / g.hj(j,i)"
   ]
  }
 ],
 "metadata": {
  "_draft": {
   "nbviewer_url": "https://gist.github.com/56416874ef83801a11382071c75fb2b9"
  },
  "gist": {
   "data": {
    "description": "Documents/GitHub/NumericalMethodsForESMs/Untitled.ipynb",
    "public": false
   },
   "id": "56416874ef83801a11382071c75fb2b9"
  },
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}