Selecting and Reindexing of Area of Interest#

Cells#

[5]:

window_cell = (
    (grid.clat >= np.deg2rad(bottom_bound))
    & (grid.clat <= np.deg2rad(top_bound))
    & (grid.clon >= np.deg2rad(left_bound))
    & (grid.clon <= np.deg2rad(right_bound))
).values

(window_cell_indices,) = np.where(window_cell)
window_cell_indices

[5]:

array([4282376, 4282377, 4282378, 4282400, 4282401, 4282402, 4282403,
       4282404, 4282405, 4282407, 4282408, 4282409, 4282410, 4282411,
       4282412, 4282413, 4282414, 4282415, 4282420, 4282421, 4282422,
       4282480, 4282481, 4282482, 4282483, 4282484, 4282485, 4282486,
       4282487, 4282488, 4282489, 4282490, 4282491, 4282493, 4282497,
       4282501, 4282508, 4282510, 4282511, 4282512, 4282513, 4282514,
       4282515, 4282516, 4282517, 4282518, 4282519, 4282520, 4282521,
       4282522, 4282523, 4282527, 4282558, 4282936, 4282938, 4282939,
       4282941, 4282960, 4282961, 4282962, 4282967, 4282968, 4282969,
       4282970, 4282971, 4282972, 4282973, 4283128, 4283129, 4283130])

What we have shown above are all the indices of the cells within the green window of the native grid. Since we selected them via np.where() we obtained the indices in 0-based python thinking. We now do the same for the vertices and edges that appear in our selected window:

Vertices#

We select with the function .isel() and the vertex_of_cell information from the grid all the vertices of the cells we have cut out one step before. We also sort the indices of the vertices and delete all duplicates vianp.unique(). Since the ICON code is written in Fortran, i.e. the integer-values are 1-based, we subtract 1 to get back to python thinking, hence 0-based.

[6]:

window_vertex_indices = (
    np.unique(grid.vertex_of_cell.isel(cell=window_cell_indices).values) - 1
)
window_vertex_indices

[6]:

array([2144929, 2144932, 2144933, 2144935, 2144937, 2144938, 2144939,
       2144953, 2144954, 2144955, 2144956, 2144957, 2144958, 2144959,
       2144960, 2144961, 2144962, 2144963, 2144968, 2144978, 2144983,
       2144984, 2145002, 2145003, 2145004, 2145005, 2145006, 2145007,
       2145008, 2145010, 2145011, 2145014, 2145020, 2145021, 2145022,
       2145023, 2145024, 2145025, 2145026, 2145217, 2145218, 2145220,
       2145222, 2145223, 2145224, 2145231, 2145237, 2145238, 2145239,
       2145286], dtype=int32)

Edges#

Same as for the vertices, we select with the function .isel() and the edge_of_cell information from the grid all the edges of the cells we have cut out two steps before. We also sort the indices of the edges and delete all duplicates via np.unique(). We subtract 1 to get back to python thinking.

[7]:

window_edge_indices = (
    np.unique(grid.edge_of_cell.isel(cell=window_cell_indices).values) - 1
)
window_edge_indices

[7]:

array([6427302, 6427311, 6427312, 6427313, 6427314, 6427315, 6427316,
       6427317, 6427318, 6427352, 6427353, 6427354, 6427355, 6427356,
       6427357, 6427358, 6427359, 6427360, 6427361, 6427362, 6427363,
       6427365, 6427366, 6427367, 6427368, 6427369, 6427370, 6427371,
       6427372, 6427373, 6427374, 6427375, 6427376, 6427377, 6427385,
       6427387, 6427388, 6427389, 6427390, 6427425, 6427426, 6427482,
       6427483, 6427484, 6427485, 6427486, 6427487, 6427488, 6427489,
       6427490, 6427491, 6427492, 6427493, 6427494, 6427495, 6427496,
       6427497, 6427498, 6427499, 6427500, 6427501, 6427504, 6427507,
       6427508, 6427513, 6427516, 6427527, 6427528, 6427529, 6427531,
       6427532, 6427533, 6427534, 6427535, 6427536, 6427537, 6427538,
       6427539, 6427540, 6427541, 6427542, 6427543, 6427544, 6427545,
       6427546, 6427547, 6427548, 6427549, 6427550, 6427553, 6427602,
       6428151, 6428152, 6428160, 6428161, 6428162, 6428164, 6428165,
       6428166, 6428167, 6428170, 6428202, 6428203, 6428204, 6428205,
       6428206, 6428207, 6428208, 6428213, 6428214, 6428215, 6428216,
       6428217, 6428218, 6428219, 6428410, 6428418, 6428419, 6428420],
      dtype=int32)

Constructing New Grid with Selected Cells, Vertices and Edges#

Wow, that’s already great ! We have received a lot of information in the form of indices in individual arrays about our green window. We merge them into one dataset so that everything is compact:

[8]:

selected_indices = xr.Dataset(
    {
        "cell": ("cell", window_cell_indices),
        "vertex": ("vertex", window_vertex_indices),
        "edge": ("edge", window_edge_indices),
    }
)

selected_indices

It could be that we need more variables for future calculations. Therefore we create a dictionary with further interesting variables, which we reindex as a precaution.

[9]:

vars_to_renumber = {
    "cell": [
        "adjacent_cell_of_edge",
        "cells_of_vertex",
        "neighbor_cell_index",
        "cell_index",
    ],
    "vertex": ["vertex_of_cell", "edge_vertices", "vertices_of_vertex"],
    "edge": ["edge_of_cell", "edges_of_vertex"],
}

We now come to the heart of this script: the reindexing. Several things happen here, which is why it is best to define a function. This function reindex_grid() needs 3 inputs and returns 1 output. The inputs are the original, complete grid and the parts of the grid that should be reindexed, hence indices and vars_to_renumber. The indices define the cells, vertices and edges. The vars_to_renumber are all variables that we are still interested in and can be composed of cells, vertices and edges. Output of our function will be a new_grid containing all indices and variables for our green window around the birthplace of Max Planck in such a way that everything starts counting at 0.

Let’s go through it step by step:

Line 1: We define a function wiht 3 input variables.

Line 2: We define as new_grid the area in the old grid that contains the indices we selected at the beginning of this script for the cells, vertices and edges. For this we use the .load() function, which loads the 17GB file into memory and processes it there: this is a little faster.

Line 3: We open a for-loop that accesses the coordinates and the entries of the array selected_indices.

Line 4: We open an array, which is only filled with -2 (exceptional value like nan but as an integer) in the original, old grid length and call it renumbering.

Line 5: We start counting at 0 at the index positions of the long renumbering array, which belong to the indices of the selected dark green area, until we have reached the length of the short, previously selected array. So what we get is an array with the length of the original grid dimension (20971520 cells, 10485762 vertices, 31457280 edges), which contains a value other than -2 only at that position within the array which is inside the dark green selected window.

Line 6: we open another for loop over the remaining variables vars_to_renumber to be reindexed.

Line 7: For the variables stored in the dictionary of a particular dimension (cell, vertex, edge), we take one item and access it in new_grid (line 2) and subtract 1 to work in python 0-based system; this is done in the square brackets on the right side of the equal sign. We use this to select the valid position in the renumbering array but in total we add 1 to output the new_grid in the same 1-based thinking as the original grid.

Line 8: We output the new_grid.

[10]:

def reindex_grid(grid, indices, vars_to_renumber):
    new_grid = grid.load().isel(
        cell=indices.cell, vertex=indices.vertex, edge=indices.edge
    )
    for dim, idx in indices.coords.items():
        renumbering = np.full(grid.dims[dim], -2, dtype="int")
        renumbering[idx] = np.arange(len(idx))
        for name in vars_to_renumber[dim]:
            new_grid[name].data = renumbering[new_grid[name].data - 1] + 1
    return new_grid

After long theory we want to use our function and create the actual new_grid:

[11]:

new_grid = reindex_grid(grid, selected_indices, vars_to_renumber)
new_grid

Let’s see if everything worked as we wanted it and choose reindexed variable like .vertex_of_cell and sort it:

[12]:

np.unique(new_grid.vertex_of_cell)

[12]:

array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16, 17,
       18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,
       35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50])

Voilà ! That worked and we are done and have now built a new grid-file tailored to the area of our interest, which was provided with a new indexing.

For further processing, the two datasets selected_indices and new_grid can be saved to have them quickly accessible for further calculations.

[13]:

selected_indices.to_netcdf(
    f"selected_indices_region_{bottom_bound}-{top_bound}_{left_bound}-{right_bound}.nc",
    mode="w",
)
new_grid.to_netcdf(
    f"new_grid_region_{bottom_bound}-{top_bound}_{left_bound}-{right_bound}.nc",
    mode="w",
)

Selecting and Reindexing of Area of Interest#

The elegant way of dealing with large output#

Importing the Grid-File#

Cells#

Vertices#

Edges#

Constructing New Grid with Selected Cells, Vertices and Edges#