# Part 2: Verification of Derivative Function#

To run this script you need output from Selecting and Reindexing of Area of Interest.

In Part 1: Deriving Derivatives on Triangular Grid we introduced a derivative() function to calculate derivatives on a triangular grid. Since it is difficult to tell with the naked eye whether our function did what we wanted it to do, we create an artificially generated field that has a constant gradient along one axis and remains constant along the other axis.

[1]:

import numpy as np
import xarray as xr


Range of area of interest:

[2]:

# Max Plancks birthplace, Kiel, Schleswig-Holstein, Germany
left_bound = 9.88
right_bound = 10.38
top_bound = 54.57
bottom_bound = 54.07


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.

[3]:

selected_indices = xr.open_dataset(
f"../selected_indices_region_{bottom_bound}-{top_bound}_{left_bound}-{right_bound}.nc"
)
new_grid = xr.open_dataset(
f"../new_grid_region_{bottom_bound}-{top_bound}_{left_bound}-{right_bound}.nc"
)

new_grid

[3]:

<xarray.Dataset>
Dimensions:                         (cell: 70, nv: 3, vertex: 50, ne: 6,
edge: 119, no: 4, nc: 2,
max_stored_decompositions: 4, two_grf: 2,
cell_grf: 14, max_chdom: 1, edge_grf: 24,
vert_grf: 13)
Coordinates:
clon                            (cell) float64 0.1804 0.1805 ... 0.174
clat                            (cell) float64 0.9453 0.946 ... 0.9488
vlon                            (vertex) float64 0.1815 0.1807 ... 0.1717
vlat                            (vertex) float64 0.9456 0.9466 ... 0.9482
elon                            (edge) float64 0.1811 0.1814 ... 0.1724
elat                            (edge) float64 0.9461 0.9471 ... 0.9487
* cell                            (cell) int64 4282376 4282377 ... 4283130
* vertex                          (vertex) int32 2144929 2144932 ... 2145286
* edge                            (edge) int32 6427302 6427311 ... 6428420
Dimensions without coordinates: nv, ne, no, nc, max_stored_decompositions,
two_grf, cell_grf, max_chdom, edge_grf, vert_grf
Data variables: (12/91)
clon_vertices                   (cell, nv) float64 0.1794 0.1802 ... 0.173
clat_vertices                   (cell, nv) float64 0.9457 0.9446 ... 0.9492
vlon_vertices                   (vertex, ne) float64 9.969e+36 ... 9.969e+36
vlat_vertices                   (vertex, ne) float64 9.969e+36 ... 9.969e+36
elon_vertices                   (edge, no) float64 0.1807 0.1805 ... 0.1728
elat_vertices                   (edge, no) float64 0.9466 0.946 ... 0.9485
...                              ...
edge_dual_normal_cartesian_x    (edge) float64 -0.6706 -0.7373 ... -0.7362
edge_dual_normal_cartesian_y    (edge) float64 -0.5071 0.4974 ... 0.4973
edge_dual_normal_cartesian_z    (edge) float64 0.5414 0.4572 ... 0.4589
cell_circumcenter_cartesian_x   (cell) float64 0.576 0.5755 ... 0.5739
cell_circumcenter_cartesian_y   (cell) float64 0.105 0.105 ... 0.1003 0.1009
cell_circumcenter_cartesian_z   (cell) float64 0.8107 0.8111 ... 0.8127
Attributes: (12/43)
title:                    ICON grid description
institution:              Max Planck Institute for Meteorology/Deutscher ...
source:                   git@git.mpimet.mpg.de:GridGenerator.git
revision:                 d00fcac1f61fa16c686bfe51d1d8eddd09296cb5
date:                     20180529 at 222250
user_name:                Rene Redler (m300083)
...                       ...
topography:               modified SRTM30
subcentre:                1
number_of_grid_used:      15
history:                  Thu Aug 16 11:05:44 2018: ncatted -O -a ICON_gr...
NCO:                      netCDF Operators version 4.7.5 (Homepage = http...
[4]:

%matplotlib inline
import matplotlib.pylab as plt


We create the artificially generated field:

[5]:

vlon = np.rad2deg(new_grid.vlon)
voc = new_grid.vertex_of_cell.T.values - 1

fig, ax = plt.subplots(figsize=(10, 6))
im = plt.tripcolor(vlon, vlat, voc, test_field, cmap=plt.cm.get_cmap("seismic"))
cbar = fig.colorbar(im)

cbar.set_label("arb. unit")
ax.set_title("Testfield")
ax.set_xlabel("Longitude / deg")
ax.set_ylabel("Latitude / deg")

plt.show()


We just copied the function derivative() from Part 2:

[6]:

def derivative(grid, data):
neighbors = (grid.neighbor_cell_index.values - 1).T
valid = np.all(neighbors >= 0, axis=-1)

p = data.values[..., neighbors[valid]]
ones = np.ones_like(neighbors_lon)
A = np.stack(
(
ones,
neighbors_lon - cell_lon[:, np.newaxis],
neighbors_lat - cell_lat[:, np.newaxis],
),
axis=2,
)
A_inv = np.linalg.inv(A)

alpha, beta, gamma = np.einsum("...ij,...j->i...", A_inv, p)

return alpha, beta, gamma, valid


We plug in our test_field and calculate the derivatives in both directions:

[7]:

test_0, dtestdlon, dtestdlat, valid = derivative(new_grid, test_field)


Let’s plot it and we see that we do not see much … But wait, eureka ! Take a look at the range of the colorbar. We are at $$1\mathrm{e}{-13}$$

[8]:

voc_valid = voc[valid]

fig, ax = plt.subplots(figsize=(10, 6))
im = plt.tripcolor(vlon, vlat, voc_valid, dtestdlon, cmap=plt.cm.get_cmap("seismic"))
cbar = fig.colorbar(im)

cbar.set_label("arb. unit")
ax.set_title("Longitudinal Derivative of Testfield")
ax.set_xlabel("Longitude / deg")
ax.set_ylabel("Latitude / deg")

plt.show()


…so we simply round to the 12th digit after the decimal point and get a white area with the corresponding color scheme. What we get is a derivative approximately equal to 1 for the whole area of interest.

[9]:

dtestdlon_rounded = np.round(dtestdlon, 12)

fig, ax = plt.subplots(figsize=(10, 6))
im = plt.tripcolor(
vlon, vlat, voc_valid, dtestdlon_rounded, cmap=plt.cm.get_cmap("seismic")
)
cbar = fig.colorbar(im)

cbar.set_label("arb. unit")
ax.set_title("Longitudinal Derivative of Testfield")
ax.set_xlabel("Longitude / deg")
ax.set_ylabel("Latitude / deg")

plt.show()


This is exactly what we expected. Because, if you have been paying attention, you will have seen that we have chosen as test_field the longitudinal coordinate of the corresponding cell, namely test_field = np.rad2deg(new_grid.clon). This means that for every 1 deg change in position along the west-east axis, the value of our test field has changed by exactly 1. The corresponding derivation along the chosen axis then inevitably results in 1 - voilà !

The same can of course still be done and checked along the latitudinal axis. For this purpose, a corresponding test_field must be generated, but we leave this to you as a small task. Enjoy !