{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Part 1: Deriving Derivatives on Triangular Grid"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"_To run this script you need output from [Selecting and Reindexing of Area of Interest](../cutting_grid_window.ipynb)._"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Theoretical Introduction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since we are working on a grid whitout any right angles, we have to come up with something new to calculate gradients or derivatives.
\n",
"Basically, the idea is this: We do a local two-dimensional Taylor evolution for a cell $i$, giving us a plane equation containing 3 unknowns. If we apply this plane equation to three points (here the three neighbors) in the neighborhood of the cell $i$, the equation can be solved. This procedure can be applied in matrix form to all cells simultaneously and thus calculate the derivative for a certain area in one go.\n",
"\n",
"To look at this from a mathematical point of view, we make use of linear algebra: if you have the coordinates of three points, you can use them to create a plane and thus a linear system of equations.
\n",
"This has the form:\n",
"\n",
"$$\\underline{\\underline{A}}~ \\underline{\\phi} = \\underline{p}.$$\n",
"\n",
"If the matrix $\\underline{\\underline{A}}$ is square and has full rank, then the system has a unique solution given by:\n",
"\n",
"$$\\underline{\\phi} = \\underline{\\underline{A}}^{-1}~\\underline{p},$$\n",
"where $\\underline{\\underline{A}}^{-1}$ is the inverse of $\\underline{\\underline{A}}$. (For matrices which are not square or of full rank the Moore-Penrose pseudoinverse can help here, and possibly be used to generalize this code.)\n",
"\n",
"We reach the respective centers of the neighboring cells via their coordinates and write down the equations in parametric form of a plane:\n",
"\n",
"$$\n",
"\\begin{matrix}\n",
"\\alpha(x_0, y_0) + \\beta (x_1-x_0) + \\gamma (y_1-y_0) = p_1(x_1, y_1)\\\\\n",
"\\alpha(x_0, y_0) + \\beta (x_2-x_0) + \\gamma (y_2-y_0) = p_2(x_2, y_2)\\\\\n",
"\\alpha(x_0, y_0) + \\beta (x_3-x_0) + \\gamma (y_3-y_0) = p_3(x_3, y_3)\n",
"\\end{matrix}\n",
"$$\n",
"\n",
"Perhaps you can already see the practical use of this approach; the derivative with respect to $x_i$ or $y_i$, with $i \\in \\{1,2,3\\}$ results in the coefficients $\\beta$ or $\\gamma$:\n",
"\n",
"$$\\frac{\\partial p_i}{\\partial x_i} = \\beta ~~\\text{and}~~ \\frac{\\partial p_j}{\\partial y_j} = \\gamma$$\n",
"This means that we have to calculate these coefficients to get the derivatives of the corresponding variable.\n",
"\n",
"We write down our equations in matrix notation:\n",
"\n",
"$$\n",
"\\underbrace{\n",
"\\begin{bmatrix}\n",
"1 & (x_1-x_0) & (y_1-y_0)\\\\\n",
"1 & (x_2-x_0) & (y_2-y_0)\\\\\n",
"1 & (x_3-x_0) & (y_3-y_0)\n",
"\\end{bmatrix}}_{\\underline{\\underline{A}}}\n",
"\\underbrace{\n",
"\\begin{bmatrix}\n",
"\\alpha\\\\\n",
"\\beta\\\\\n",
"\\gamma\n",
"\\end{bmatrix}}_{\\underline{\\phi}}\n",
"=\n",
"\\underbrace{\n",
"\\begin{bmatrix}\n",
"p_1(x_1, y_1)\\\\\n",
"p_2(x_2, y_2)\\\\\n",
"p_3(x_3, y_3)\n",
"\\end{bmatrix}}_{\\underline{p}}\n",
"$$\n",
"\n",
"As a final step, we calculate the inverse of $\\underline{\\underline{A}}$ and automatically obtain the coefficients $\\alpha, \\beta, \\gamma$ by calculating the matrix product $$\\underline{\\phi} = \\underline{\\underline{A}}^{-1}~\\underline{p}.$$\n",
"\n",
"Let's get our hands dirty now!"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import xarray as xr\n",
"import numpy as np"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Deriving Derivatives"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The individual output variables of the ICON simulations do not contain any coordinates apart from the height (only in 3D output) and the time. Therefore it is necessary to get the coordinates from the grid-file.
\n",
"We have created a `new_grid` via the [Selecting and Reindexing of Area of Interest](../cutting_grid_window.ipynb) script, which contains the same information as the global grid-file focused on our area of interest with the advantage of being smaller and easier to handle."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Range of area of interest:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# Max Plancks birthplace, Kiel, Schleswig-Holstein, Germany\n",
"left_bound = 9.88\n",
"right_bound = 10.38\n",
"top_bound = 54.57\n",
"bottom_bound = 54.07"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We import the `selected_indices`, hence the indices for `cells`, `vertices` and `edges`, just like the `new_grid` which contains the grid information of the area of interest."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
<xarray.Dataset>\n", "Dimensions: (cell: 70, nv: 3, vertex: 50, ne: 6,\n", " edge: 119, no: 4, nc: 2,\n", " max_stored_decompositions: 4, two_grf: 2,\n", " cell_grf: 14, max_chdom: 1, edge_grf: 24,\n", " vert_grf: 13)\n", "Coordinates:\n", " clon (cell) float64 0.1804 0.1805 ... 0.174\n", " clat (cell) float64 0.9453 0.946 ... 0.9488\n", " vlon (vertex) float64 0.1815 0.1807 ... 0.1717\n", " vlat (vertex) float64 0.9456 0.9466 ... 0.9482\n", " elon (edge) float64 0.1811 0.1814 ... 0.1724\n", " elat (edge) float64 0.9461 0.9471 ... 0.9487\n", " * cell (cell) int64 4282376 4282377 ... 4283130\n", " * vertex (vertex) int32 2144929 2144932 ... 2145286\n", " * edge (edge) int32 6427302 6427311 ... 6428420\n", "Dimensions without coordinates: nv, ne, no, nc, max_stored_decompositions,\n", " two_grf, cell_grf, max_chdom, edge_grf, vert_grf\n", "Data variables: (12/91)\n", " clon_vertices (cell, nv) float64 0.1794 0.1802 ... 0.173\n", " clat_vertices (cell, nv) float64 0.9457 0.9446 ... 0.9492\n", " vlon_vertices (vertex, ne) float64 9.969e+36 ... 9.969e+36\n", " vlat_vertices (vertex, ne) float64 9.969e+36 ... 9.969e+36\n", " elon_vertices (edge, no) float64 0.1807 0.1805 ... 0.1728\n", " elat_vertices (edge, no) float64 0.9466 0.946 ... 0.9485\n", " ... ...\n", " edge_dual_normal_cartesian_x (edge) float64 -0.6706 -0.7373 ... -0.7362\n", " edge_dual_normal_cartesian_y (edge) float64 -0.5071 0.4974 ... 0.4973\n", " edge_dual_normal_cartesian_z (edge) float64 0.5414 0.4572 ... 0.4589\n", " cell_circumcenter_cartesian_x (cell) float64 0.576 0.5755 ... 0.5739\n", " cell_circumcenter_cartesian_y (cell) float64 0.105 0.105 ... 0.1003 0.1009\n", " cell_circumcenter_cartesian_z (cell) float64 0.8107 0.8111 ... 0.8127\n", "Attributes: (12/43)\n", " title: ICON grid description\n", " institution: Max Planck Institute for Meteorology/Deutscher ...\n", " source: git@git.mpimet.mpg.de:GridGenerator.git\n", " revision: d00fcac1f61fa16c686bfe51d1d8eddd09296cb5\n", " date: 20180529 at 222250\n", " user_name: Rene Redler (m300083)\n", " ... ...\n", " topography: modified SRTM30\n", " subcentre: 1\n", " number_of_grid_used: 15\n", " history: Thu Aug 16 11:05:44 2018: ncatted -O -a ICON_gr...\n", " ICON_grid_file_uri: http://icon-downloads.mpimet.mpg.de/grids/publi...\n", " NCO: netCDF Operators version 4.7.5 (Homepage = http...
<xarray.DataArray 'pfull' (time: 1, height: 90, ncells: 20971520)>\n", "[1887436800 values with dtype=float32]\n", "Coordinates:\n", " * height (height) float64 1.0 2.0 3.0 4.0 5.0 ... 86.0 87.0 88.0 89.0 90.0\n", " * time (time) float64 2.02e+07\n", "Dimensions without coordinates: ncells\n", "Attributes:\n", " standard_name: air_pressure\n", " long_name: Pressure\n", " units: Pa\n", " param: 0.3.0\n", " CDI_grid_type: unstructured\n", " number_of_grid_in_reference: 1
<xarray.DataArray 'pfull' (height: 10, cell: 70)>\n", "array([[ 89579.93 , 89524.6 , 89488.8 , 89601.09 , 89637.8 ,\n", " 89505.664, 89588.016, 89515.9 , 89455.82 , 89628.78 ,\n", " 89769.734, 89697.72 , 89747.2 , 89601.88 , 89620.94 ,\n", " 89537.88 , 89568.84 , 89619.6 , 89473.53 , 89578.05 ,\n", " 89530.805, 89863.78 , 89874.59 , 89801.36 , 89865.234,\n", " 89873.83 , 89874.29 , 89878.98 , 89876.27 , 89785.62 ,\n", " 89815.086, 89723.125, 89753.766, 89685.4 , 89664.7 ,\n", " 89672.19 , 89590.06 , 89466.48 , 89633.836, 89566.055,\n", " 89569.05 , 89584.96 , 89607.36 , 89586.35 , 89625.45 ,\n", " 89485.46 , 89533.4 , 89723.61 , 89655.95 , 89720.44 ,\n", " 89754.84 , 89648.88 , 89516.305, 89732.47 , 89718.28 ,\n", " 89825.27 , 89695.81 , 89882.02 , 89862.805, 89879.836,\n", " 89692.67 , 89654.26 , 89836.68 , 89712.78 , 89724.66 ,\n", " 89881.46 , 89886.16 , 89770.2 , 89740.29 , 89750.305],\n", " [ 91697.42 , 91636.45 , 91596.43 , 91717.51 , 91757.08 ,\n", " 91612.414, 91704.62 , 91624.65 , 91559.664, 91747.234,\n", " 91900.195, 91821.46 , 91874.39 , 91716.66 , 91740.4 ,\n", " 91648.15 , 91683.46 , 91739.58 , 91580.266, 91695.84 ,\n", " 91642.55 , 92001.88 , 92012.78 , 91934.52 , 92003.15 ,\n", " 92011.06 , 92011.81 , 92015.88 , 92013.9 , 91915.88 ,\n", "...\n", " 101646.22 , 101715.98 , 102000.59 , 101901.94 , 101991.586,\n", " 102043.75 , 101875.74 , 101698.97 , 101997.73 , 101980.17 ,\n", " 102144.83 , 101936.79 , 102235.68 , 102205.15 , 102234.97 ,\n", " 101932.21 , 101877.91 , 102163.36 , 101973.68 , 101991.84 ,\n", " 102228.2 , 102237.72 , 102058.66 , 102013.766, 102033.27 ],\n", " [102309.24 , 102219.92 , 102155.516, 102326.03 , 102382.24 ,\n", " 102168.05 , 102313.945, 102194.7 , 102103.3 , 102368.69 ,\n", " 102579.664, 102468.74 , 102534.625, 102313.586, 102361.57 ,\n", " 102221.02 , 102278.79 , 102365.4 , 102133.664, 102308.28 ,\n", " 102222.21 , 102722.5 , 102731.53 , 102631.77 , 102719.91 ,\n", " 102719.28 , 102722.78 , 102720.305, 102726.414, 102594.2 ,\n", " 102638.445, 102489.87 , 102553.766, 102450.336, 102398.12 ,\n", " 102403.766, 102294.336, 102106.266, 102357.83 , 102246.2 ,\n", " 102257.36 , 102266.836, 102310.47 , 102292.81 , 102346.68 ,\n", " 102137.52 , 102208.1 , 102495.82 , 102396.266, 102486.766,\n", " 102539. , 102368.85 , 102190.805, 102492.61 , 102475.08 ,\n", " 102641.555, 102430.74 , 102734.87 , 102702.695, 102734.59 ,\n", " 102426.09 , 102371.79 , 102660.34 , 102468.66 , 102487.75 ,\n", " 102726.49 , 102736.54 , 102553.94 , 102508.5 , 102528.7 ]],\n", " dtype=float32)\n", "Coordinates:\n", " * height (height) float64 81.0 82.0 83.0 84.0 85.0 86.0 87.0 88.0 89.0 90.0\n", " * cell (cell) int64 4282376 4282377 4282378 ... 4283128 4283129 4283130