{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Example calculation of finite difference operators on a uniformly spaced A-grid using xgcm\n", "* Primary example is our Reanalysis Data\n", "* We handle the convergence of the longitude lines analytically in \n", "order to simplify the metric terms\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import xarray as xr\n", "from matplotlib import pyplot as plt\n", "%matplotlib inline\n", "plt.rcParams['figure.figsize'] = (8,4)\n", "import xgcm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### First we will import a sample vector field (u10,v10)\n", "Note, this will only work on our machines which see the kage drives" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Dimensions: (X: 360, Y: 181)\n", "Coordinates:\n", " * X (X) float32 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 ...\n", " * Y (Y) float32 90.0 89.0 88.0 87.0 86.0 85.0 84.0 83.0 82.0 81.0 ...\n", "Data variables:\n", " si10 (Y, X) float32 dask.array\n", " u10 (Y, X) float32 dask.array\n", " v10 (Y, X) float32 dask.array" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "files = '/net/kage/d5/datasets/ERAInterim/monthly/Surface/*10.nc'\n", "ds1 = xr.open_mfdataset(files,decode_times=False)\n", "ds = ds1.mean('T')\n", "#ds.to_netcdf('ERAInt_surf_clim.nc')\n", "#ds = xr.open_dataset('/home/pangeo/data/ERAInt_surf_clim.nc')\n", "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The reanalysis data is on an 'A-grid' (all data on the same grid)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We need to define the 'half point' grids\n", "N.B. The Y grid is going from North to South" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "Dimensions: (X: 360, Xh: 360, Y: 181, Yh: 180)\n", "Coordinates:\n", " * X (X) float32 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 ...\n", " * Y (Y) float32 90.0 89.0 88.0 87.0 86.0 85.0 84.0 83.0 82.0 81.0 ...\n", " * Xh (Xh) float64 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 ...\n", " * Yh (Yh) float64 89.5 88.5 87.5 86.5 85.5 84.5 83.5 82.5 81.5 80.5 ...\n", "Data variables:\n", " si10 (Y, X) float32 dask.array\n", " u10 (Y, X) float32 dask.array\n", " v10 (Y, X) float32 dask.array" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "delx = ds.X[1]-ds.X[0]; dely = ds.Y[1]-ds.Y[0]\n", "xh = ds.X + delx/2.0; ds['Xh'] = ('Xh',xh)\n", "yh = ds.Y[0:-1] + dely/2.0; ds['Yh'] = ('Yh',yh)\n", "\n", "ds" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "-" } }, "source": [ "## Next we create our xgcm Grid:\n", "Yes, I know it is not periodic in Y, but it is easier (don't need to give boundary=...) and we just won't look at the derivatives near the poles\n", "\n", "### Recall the possible positions of our original A-grid coordinates and our new half grid coordinates:\n", "\n", "![positions](http://byrd.ldeo.columbia.edu/images/positions.png)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "X Axis (periodic):\n", " * center X (360) --> right\n", " * right Xh (360) --> center\n", "Y Axis (periodic):\n", " * center Y (181) --> inner\n", " * inner Yh (180) --> center" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "grid = xgcm.Grid(ds, coords={'X': {'center': 'X', 'right': 'Xh'},\n", " 'Y': {'center': 'Y', 'inner': 'Yh'}}, periodic=['X','Y'])\n", "grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### and our metric terms:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(\n", " array(111000.), \n", " array(-111000.))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coslat = np.cos(np.deg2rad(ds.Y))\n", "coslath = np.cos(np.deg2rad(ds.Yh))\n", "meterPdegree = 111000.0\n", "\n", "dxm = delx*meterPdegree\n", "dym = dely*meterPdegree\n", "dxm,dym" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now that we have our grids and metrics defined, the hard part is over. \n", "\n", "### Let's use xgcm to calculate our standard differential operators: \n", "* divergence\n", "* gradient\n", "* curl " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Divergence\n", "\n", "\n", "\n", "The horizontal divergence of $\\vec u = (u,v)$ in spherical coordinates is\n", "\n", "$$ \\nabla \\cdot \\vec u = \\frac{1}{\\cos y} \\left( \\frac{\\partial u}{\\partial x} + \\\n", " \\frac{\\partial (v \\cos y )}{\\partial y} \\right)$$\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# calculate the divergence of the 10m winds, (u10,v10)\n", "# ORDER OF THE diff and interp MATTER:\n", "# Method I (discouraged) - calculate the derivatives and then put them where needed\n", "dudx = grid.diff(ds.u10,axis='X')/dxm \n", "dvdy = grid.diff(coslat*ds.v10,axis='Y')/dym\n", "divergence = (grid.interp(dudx,axis='X')+ grid.interp(dvdy,axis='Y'))/coslat\n", "ds['divergence1'] = divergence" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# same calculation, \n", "# Method II (preferred) - interpolate the quantities to edge of trapezoid and then compute derivatives\n", "uh = grid.interp(ds.u10,axis='X')\n", "vh = grid.interp(ds.v10,axis='Y')\n", "dudx = grid.diff(uh,axis='X')/dxm \n", "dvdy = grid.diff(coslath*vh,axis='Y')/dym\n", "divergence = (dudx + dvdy)/coslat\n", "ds['divergence'] = divergence" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "div1 = ds.sel(X=240,Y=slice(30,-30)).divergence1\n", "div2 = ds.sel(X=240,Y=slice(30,-30)).divergence\n", "div_diff = (div2 - div1)*1000\n", "div1.plot(marker='o', label=\"Method I\")\n", "div2.plot(label=\"Method II\")\n", "div_diff.plot(label=\"1000*difference\")\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "del ds['divergence1']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Gradient\n", "\n", "The gradient of a scalar field, $F$, in spherical coordinates is\n", "\n", "$$ \\nabla F = (\\frac{1}{cos y}\\frac{\\partial F}{\\partial x} , \\\n", " \\frac{\\partial F}{\\partial y} )$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# calculate the gradient of the mean 10m wind speed, si10\n", "ds.si10.plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Fhx = grid.interp(ds.si10,axis='X')\n", "Fhy = grid.interp(ds.si10,axis='Y')\n", "ds['wnspgradx'] = grid.diff(Fhx,axis='X')/(dxm*coslat)\n", "ds['wnspgrady'] = grid.diff(Fhy,axis='Y')/dym" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#plt.quiver(ds.X, ds.Y[5:-5], gradx[0,5:-5], grady[0,5:-5])\n", "ds.sel(X=240,Y=slice(30,-30)).wnspgradx.plot()\n", "ds.sel(X=240,Y=slice(30,-30)).wnspgrady.plot()\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Curl (Vorticity)\n", "\n", "The curl of $\\vec u = (u,v)$ in spherical coordinates is\n", "\n", "$$\\zeta = \\nabla \\times \\vec u = \\frac{1}{\\cos y} \\left ( \\frac{\\partial v}{\\partial x} - \\frac{\\partial (u \\cos y)}{\\partial y} \\right\n", " )$$\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "uh = grid.interp(ds.u10,axis='Y')\n", "vh = grid.interp(ds.v10,axis='X')\n", "dudy = grid.diff(uh*coslath,axis='Y')/dym\n", "dvdx = grid.diff(vh,axis='X')/dxm\n", "ds['curl'] = (dvdx - dudy)/coslat\n", "ds" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds.curl[5:-5].plot(vmin = -1e-5, vmax = 1e-5)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "gist_info": { "gist_id": null, "gist_url": null }, "kernelspec": { "display_name": "Py3 pangeo", "language": "python", "name": "pangeo" }, "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.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }