Triangular Meshes and Basic Plotting#

[1]:

%matplotlib inline

[2]:

import matplotlib.pylab as plt

anatomy of a triangular mesh#

Triangular meshes are often represented as a collection of corner points (vertices) and a collection of triangles. Triangles are typically represented as three indices from within the vertex arrays.

Let’s create four vertices at the corner of the unit square and fill the area with two triangles. Each of the two triangles will cary a value (think of surface temperature, some fluxes etc…).

[3]:

xs = [0, 1, 0, 1]
ys = [0, 0, 1, 1]
tris = [[0, 1, 2], [1, 2, 3]]
values = [1, 2]

Matplotlib’s ``tripcolor` <https://matplotlib.org/stable/gallery/images_contours_and_fields/tripcolor_demo.html>`__ takes all of these pieces and creates a plot out of it.

[4]:

plt.tripcolor(xs, ys, tris, values)
plt.colorbar()

[4]:

<matplotlib.colorbar.Colorbar at 0x2b8079edc790>
../../_images/Processing_playing_with_triangles_tripcolor_6_1.png

application to model data#

Let’s fetch some data from the intake catalog, attach the gridfile to it and select a lat-lon region for plotting.

[5]:

import intake
import numpy as np
import xarray as xr
import yaml
from gridlocator import (
    merge_grid,
)  # gridlocator.py can be found further below at 'playing with triangles'

[6]:

col_url = "/pf/k/k202134/DYAMOND/Processing/full_catalog.json"
col = intake.open_esm_datastore(col_url)

[7]:

res = col.search(
    variable="prw",
    project="NextGEMS",
    model="ICON-SAP-5km",
    ensemble_member="dpp0052",
    operation="inst",
)
# res

[8]:

d = res.to_dataset_dict(cdf_kwargs={"chunks": {"time": 1}})


--> The keys in the returned dictionary of datasets are constructed as follows:
        'project.institute.model.experiment.domain.frequency.grid.level_type.ensemble_member.operation'
100.00% [1/1 00:00<00:00] 
[9]:

ds = merge_grid(next(iter(d.values())))
ds

[9]:

<xarray.Dataset>
Dimensions:                         (time: 1969, cell: 20971520, nv: 3, vertex: 10485762, ne: 6, edge: 31457280, no: 4, nc: 2, max_stored_decompositions: 4, two_grf: 2, cell_grf: 14, max_chdom: 1, edge_grf: 24, vert_grf: 13)
Coordinates:
  * time                            (time) float64 2.02e+07 ... 2.02e+07
    clon                            (cell) float64 ...
    clat                            (cell) float64 ...
    vlon                            (vertex) float64 ...
    vlat                            (vertex) float64 ...
    elon                            (edge) float64 ...
    elat                            (edge) float64 ...
Dimensions without coordinates: cell, nv, vertex, ne, edge, no, nc, max_stored_decompositions, two_grf, cell_grf, max_chdom, edge_grf, vert_grf
Data variables: (12/92)
    prw                             (time, cell) float32 dask.array<chunksize=(1, 20971520), meta=np.ndarray>
    clon_vertices                   (cell, nv) float64 ...
    clat_vertices                   (cell, nv) float64 ...
    vlon_vertices                   (vertex, ne) float64 ...
    vlat_vertices                   (vertex, ne) float64 ...
    elon_vertices                   (edge, no) float64 ...
    ...                              ...
    edge_dual_normal_cartesian_x    (edge) float64 ...
    edge_dual_normal_cartesian_y    (edge) float64 ...
    edge_dual_normal_cartesian_z    (edge) float64 ...
    cell_circumcenter_cartesian_x   (cell) float64 ...
    cell_circumcenter_cartesian_y   (cell) float64 ...
    cell_circumcenter_cartesian_z   (cell) float64 ...
Attributes:
    uuidOfHGrid:             0f1e7d66-637e-11e8-913b-51232bb4d8f9
    references:              see MPIM/DWD publications
    history:                 /work/mh0287/k203123/GIT/icon-aes-dyw3/bin/icon ...
    intake_esm_varname:      ['prw']
    source:                  git@gitlab.dkrz.de:icon/icon-aes.git@83f3dcef81e...
    title:                   ICON simulation
    institution:             Max Planck Institute for Meteorology/Deutscher W...
    grid_file_uri:           http://icon-downloads.mpimet.mpg.de/grids/public...
    number_of_grid_used:     15
    CDI:                     Climate Data Interface version 1.8.3rc (http://m...
    Conventions:             CF-1.6
    intake_esm_dataset_key:  NextGEMS.MPIM-DWD-DKRZ.ICON-SAP-5km.Cycle1.atm.3...

Here, we select a 1° x 1° region such that we’ll have a chance to actually see some triangles. In ICON terminology, a triangle is a “cell”. So we can select our desired triangles based on the cell-latitude and cell-logitude.

[10]:

mask = (
    (ds.clat.values > np.deg2rad(0))
    & (ds.clat.values < np.deg2rad(1))
    & (ds.clon.values > np.deg2rad(-60))
    & (ds.clon.values < np.deg2rad(-59))
)

For each of the selected cell, we need the indices of the three corners of the triangles (the tris array from above). In ICON grid files, this array is called vertex_of_cell. It also is transposed and uses Fortran-indexing (which starts counting by one). Let’s load this array for our selected region and correct for the index offset:

[11]:

voc = ds.vertex_of_cell.T[mask].values - 1

We also want to know the bounding box of our triangles to properly set the extent of our plot. We could use the range which we used above to select our triangles, but then we couldn’t apply the example to other selections. Furthermore, here we are interested in the bounding box of the triangle corners not the triangle centers.

So what we’ll do is to find out which vertices actually are used by any of the selected triangles:

[12]:

used_vertices = np.unique(voc)

Now we can compute the extent of our plot:

[13]:

lat_min = ds.vlat[used_vertices].min().values
lat_max = ds.vlat[used_vertices].max().values
lon_min = ds.vlon[used_vertices].min().values
lon_max = ds.vlon[used_vertices].max().values

Now, we are all set and can do our little plot.

[14]:

plt.tripcolor(ds.vlon, ds.vlat, voc, ds.prw.isel(time=100, cell=mask).values)
plt.xlim(lon_min, lon_max)
plt.ylim(lat_min, lat_max)
plt.colorbar()

[14]:

<matplotlib.colorbar.Colorbar at 0x2b80afa4b8e0>
../../_images/Processing_playing_with_triangles_tripcolor_22_1.png