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...
    ICON_grid_file_uri:       http://icon-downloads.mpimet.mpg.de/grids/publi...
    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)
vlat = np.rad2deg(new_grid.vlat)
voc = new_grid.vertex_of_cell.T.values - 1
test_field = np.rad2deg(new_grid.clon)

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()
../../../_images/Processing_playing_with_triangles_Gradient_on_Triangular_Grid_test_derivative_10_0.png

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)

    cell_lon = np.rad2deg(new_grid.clon.values[valid])
    cell_lat = np.rad2deg(new_grid.clat.values[valid])
    neighbors_lon = np.rad2deg(new_grid.clon.values[neighbors[valid]])
    neighbors_lat = np.rad2deg(new_grid.clat.values[neighbors[valid]])

    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()
../../../_images/Processing_playing_with_triangles_Gradient_on_Triangular_Grid_test_derivative_16_0.png

…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()
../../../_images/Processing_playing_with_triangles_Gradient_on_Triangular_Grid_test_derivative_18_0.png

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 !