{ "cells": [ { "cell_type": "markdown", "source": [ "# Grid Metrics\n", "Most modern circulation models discretize the partial differential equations needed to simulate the earth system on a logically rectangular grid. This means the grid for a single time step can be represented as a 3-dimensional array of cells. Even for more complex grid geometries like [here](https://xgcm.readthedocs.io/en/latest/grid_topology.html), subdomains are usually organized in this manner. A noteable exception are models with unstructured grids [example](https://en.wikipedia.org/wiki/Unstructured_grid), which currently cannot be processed with the datamodel of xarray and xgcm.\n", "\n", "Our grid operators work on the logically rectangular grid of an ocean model, meaning that e.g. differences are evaluated on the 'neighboring' cells in either direction, but even though these cells are adjacent, cells can have different size and geometry.\n", "\n", "In order to convert operators acting on the logically rectangular grid to physically meaningful output models need 'metrics' - information about the grid cell geometry in physical space.\n", "In the case of a perfectly [rectangular cuboid](https://en.wikipedia.org/wiki/Cuboid), the only metrics needed would be three of the edge distances. All other distances can be reconstructed exactly. Most ocean models have however slightly distorted cells, due to the curvature of the earth. To accurately represent the volume of the cell we require more metrics. \n", "\n", "Each grid point has three kinds of fundamental metrics associated with it which differ in the number of described axes:\n", "\n", "1. **Distances**: A distance is associated with a single axis (e.g. ('X',),('Y',) or ('Z',)). Each distance describes the distance from the point to either face of the cell associated with the grid point. \n", "2. **Areas**: An area is associated with a pair of axes (e.g. ('X', 'Y'), ('Y', 'Z') and ('X', 'Z')). Each grid point intersects three areas.\n", "3. **Volume**: The cell volume is unique for each cell and associated with all three axes (('X', 'Y', 'Z')).\n", "\n", "## Using metrics with xgcm\n", "Once the user assigns the metrics (given as coordinates in most model output) to the grid object, xgcm is able to automatically select and apply these to calculate e.g. derivatives and integrals from model data. \n", "\n", "
\n", "\n", "*Note*: xgcm does not currently check for alignment of missing values between data and metrics. The user needs to check and mask values appropriately\n", "\n", "
" ], "metadata": {} }, { "cell_type": "code", "execution_count": 1, "source": [ "import xarray as xr\n", "import numpy as np\n", "from xgcm import Grid\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ], "outputs": [ { "output_type": "stream", "name": "stderr", "text": [ "/Users/tim/anaconda3/lib/python3.7/site-packages/docrep/__init__.py:413: MatplotlibDeprecationWarning: \n", "The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.\n", " s = dedents(s)\n", "/Users/tim/anaconda3/lib/python3.7/site-packages/docrep/__init__.py:342: MatplotlibDeprecationWarning: \n", "The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.\n", " s = dedents('\\n' + '\\n'.join(lines[first:]))\n" ] } ], "metadata": {} }, { "cell_type": "code", "execution_count": 2, "source": [ "# hack to make file name work with nbsphinx and binder\n", "import os\n", "fname = '../datasets/mitgcm_example_dataset_v2.nc'\n", "if not os.path.exists(fname):\n", " fname = '../' + fname\n", " \n", "ds = xr.open_dataset(fname)\n", "ds" ], "outputs": [ { "output_type": "execute_result", "data": { "text/plain": [ "\n", "Dimensions: (XC: 90, XG: 90, YC: 40, YG: 40, Z: 15, Zl: 15, time: 1)\n", "Coordinates:\n", " * time (time) timedelta64[ns] 11:00:00\n", " maskC (Z, YC, XC) bool ...\n", " dyC (YG, XC) float32 ...\n", " hFacC (Z, YC, XC) float32 ...\n", " rA (YC, XC) float32 ...\n", " hFacS (Z, YG, XC) float32 ...\n", " Depth (YC, XC) float32 ...\n", " * YG (YG) float32 -80.0 -76.0 -72.0 -68.0 -64.0 ... 64.0 68.0 72.0 76.0\n", " * Z (Z) float32 -25.0 -85.0 -170.0 -290.0 ... -3575.0 -4190.0 -4855.0\n", " PHrefC (Z) float32 ...\n", " dyG (YC, XG) float32 ...\n", " rAw (YC, XG) float32 ...\n", " drF (Z) float32 ...\n", " * YC (YC) float32 -78.0 -74.0 -70.0 -66.0 -62.0 ... 66.0 70.0 74.0 78.0\n", " dxG (YG, XC) float32 ...\n", " * XG (XG) float32 0.0 4.0 8.0 12.0 16.0 ... 344.0 348.0 352.0 356.0\n", " iter (time) int64 ...\n", " maskW (Z, YC, XG) bool ...\n", " * Zl (Zl) float32 0.0 -50.0 -120.0 -220.0 ... -3280.0 -3870.0 -4510.0\n", " rAs (YG, XC) float32 ...\n", " rAz (YG, XG) float32 ...\n", " maskS (Z, YG, XC) bool ...\n", " dxC (YC, XG) float32 ...\n", " hFacW (Z, YC, XG) float32 ...\n", " * XC (XC) float32 2.0 6.0 10.0 14.0 18.0 ... 346.0 350.0 354.0 358.0\n", "Data variables:\n", " UVEL (time, Z, YC, XG) float32 ...\n", " VVEL (time, Z, YG, XC) float32 ...\n", " WVEL (time, Zl, YC, XC) float32 ...\n", " SALT (time, Z, YC, XC) float32 ...\n", " THETA (time, Z, YC, XC) float32 ...\n", " PH (time, Z, YC, XC) float32 ...\n", " Eta (time, YC, XC) float32 ...\n", "Attributes:\n", " Conventions: CF-1.6\n", " title: netCDF wrapper of MITgcm MDS binary data\n", " source: MITgcm\n", " history: Created by calling `open_mdsdataset(extra_metadata=None, ll..." ] }, "metadata": {}, "execution_count": 2 } ], "metadata": {} }, { "cell_type": "markdown", "source": [ "For mitgcm output we need to first incorporate partial cell thicknesses into new metric coordinates" ], "metadata": {} }, { "cell_type": "code", "execution_count": 3, "source": [ "ds['drW'] = ds.hFacW * ds.drF #vertical cell size at u point\n", "ds['drS'] = ds.hFacS * ds.drF #vertical cell size at v point\n", "ds['drC'] = ds.hFacC * ds.drF #vertical cell size at tracer point" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "To assign the metrics, the user has to provide a dictionary with keys and entries corresponding to the spatial orientation of the metrics and a list of the appropriate variable names in the dataset. " ], "metadata": {} }, { "cell_type": "code", "execution_count": 4, "source": [ "metrics = {\n", " ('X',): ['dxC', 'dxG'], # X distances \n", " ('Y',): ['dyC', 'dyG'], # Y distances \n", " ('Z',): ['drW', 'drS', 'drC'], # Z distances \n", " ('X', 'Y'): ['rA', 'rAz', 'rAs', 'rAw'] # Areas \n", "}\n", "grid = Grid(ds, metrics=metrics)" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "## Grid-aware (weighted) integration" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "It is now possible to integrate over any grid axis. For example, we can integrate over the Z axis to compute the \n", "discretized version of:\n", "\n", "$$\\int_{-H}^0 u dz$$\n", "\n", "in one line:" ], "metadata": {} }, { "cell_type": "code", "execution_count": 5, "source": [ "grid.integrate(ds.UVEL, 'Z').plot();" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "