ICON/Triangular grid

This example shows how to draw the edges of an ICON triangular mesh with easyclimate.plot.icon.plot_triangular_grid.

The helper reads ICON cell-center and vertex-boundary coordinates from the dataset, then draws the triangular edges on either plain Matplotlib axes or Cartopy map projections.

import xarray as xr
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import easyclimate as ecl

from easyclimate.plot.icon import plot_triangular_grid

Open a sample ICON pressure-level file. The file contains the native ICON cell-center coordinates and triangular cell bounds needed by the mesh helper.

data = ecl.open_tutorial_dataset("icon_native_pl_20080909T000000Z")
data

icon_native_pl_20080909T000000Z.nc ━━━━━━━ 100.0% • 3.5/3.5 • 26.8     • 0:00:00
                                                    MB        MB/s

<xarray.Dataset> Size: 6MB
Dimensions:    (ncells: 81920, vertices: 3, time: 1, plev: 1)
Coordinates:
    clon       (ncells) float64 655kB ...
    clat       (ncells) float64 655kB ...
  * time       (time) float64 8B 2.008e+07
  * plev       (plev) float64 8B 8.5e+04
Dimensions without coordinates: ncells, vertices
Data variables:
    clon_bnds  (ncells, vertices) float64 2MB ...
    clat_bnds  (ncells, vertices) float64 2MB ...
    temp       (time, plev, ncells) float32 328kB ...

xarray.Dataset

• Dimensions:
  • ncells: 81920
  • vertices: 3
  • time: 1
  • plev: 1
• Coordinates: (4)
  • clon
    (ncells)
    float64
    ...

    standard_name :

    longitude

    long_name :

    center longitude

    units :

    radian

    bounds :

    clon_bnds

    [81920 values with dtype=float64]

  • clat
    (ncells)
    float64
    ...

    standard_name :

    latitude

    long_name :

    center latitude

    units :

    radian

    bounds :

    clat_bnds

    [81920 values with dtype=float64]

  • time
    (time)
    float64
    2.008e+07

    standard_name :

    time

    units :

    day as %Y%m%d.%f

    calendar :

    proleptic_gregorian

    axis :

    T

    array([20080909.])

  • plev
    (plev)
    float64
    8.5e+04

    standard_name :

    air_pressure

    long_name :

    pressure

    units :

    Pa

    positive :

    down

    axis :

    Z

    array([85000.])

• Data variables: (3)
  • clon_bnds
    (ncells, vertices)
    float64
    ...

    [245760 values with dtype=float64]

  • clat_bnds
    (ncells, vertices)
    float64
    ...

    [245760 values with dtype=float64]

  • temp
    (time, plev, ncells)
    float32
    ...

    standard_name :

    air_temperature

    long_name :

    Temperature

    units :

    K

    param :

    0.0.0

    CDI_grid_type :

    unstructured

    number_of_grid_in_reference :

    1

    [81920 values with dtype=float32]

Draw the triangular mesh directly on the current Matplotlib axes.

plot_triangular_grid(data)

<matplotlib.collections.LineCollection object at 0x7a8799db7080>

Plot the global ICON mesh on a shifted PlateCarree projection. Passing a Cartopy GeoAxes lets the helper add the mesh with geographic coordinates.

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.PlateCarree(180))
ax.coastlines(resolution="110m", linewidth=0.6)

plot_triangular_grid(data, ax = ax, transform=ccrs.PlateCarree());
ax.set_global()
ax.gridlines(draw_labels=True, alpha = 0)

ax.set_title("Global ICON mesh 1")

Text(0.5, 1.0453342943342954, 'Global ICON mesh 1')

The same native triangular mesh can be displayed with a Robinson projection.

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.Robinson(0))
ax.coastlines(resolution="110m", linewidth=0.6)

plot_triangular_grid(data, ax = ax, transform=ccrs.PlateCarree());
ax.set_global()
ax.gridlines(draw_labels=True, alpha = 0)

ax.set_title("Global ICON mesh 2")

Text(0.5, 1.0453342943342951, 'Global ICON mesh 2')

Orthographic projection is useful for inspecting the global grid from a regional viewing angle.

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.Orthographic(central_longitude=110.0, central_latitude=70.0))
ax.coastlines(resolution="110m", linewidth=0.6)

plot_triangular_grid(data, ax = ax, transform=ccrs.PlateCarree())
ax.set_global()
ax.gridlines(draw_labels=True, alpha = 0)

ax.set_global()
ax.set_title("Global ICON mesh 3")

Text(0.5, 1.0381854681723834, 'Global ICON mesh 3')

Restrict the plotted cells to a regional longitude-latitude window. The helper selects cells by their centers and draws only the visible triangular edges.

fig = plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.PlateCarree(180))
ax.coastlines(resolution="50m", linewidth=0.6, color="b")

plot_triangular_grid(
    data,
    ax = ax,
    transform=ccrs.PlateCarree(),
    lon_min = 90,
    lon_max = 150,
    lat_min = 10,
    lat_max = 50,
);
ax.set_extent([90, 150, 10, 50], crs = ccrs.PlateCarree())
ax.gridlines(draw_labels=True, alpha = 0)

ax.set_title("East Asia ICON mesh 1")

Text(0.5, 1.0453342943342956, 'East Asia ICON mesh 1')