Python plot data on a triangular grid with tripcolor#
This example shows how the ICON data can be displayed with Python on the triangular grid. Anyone who has tried to do this with Matplotlib’s matplotlib.pyplot.tripcolor() routine will have been surprised by the artfully nice visualizations. The so-called ‘wrapping of lines’ of the data once around is suppressed here by masking the triangles. The triangulation is generated using Matplotlib’s tri.Triangulation(), and for plotting the colored triangles we use the pyplot.tripcolor() function.
Software requirements:
Python 3
numpy
xarray
matplotlib
cartopy
Example script#
triangular_grid_with_tripcolor_ICON.py
#!/usr/bin/env python
# coding: utf-8
'''
-------------------------------------------------------------------------------
DKRZ example
Plotting ICON data sets on triangular grid
This example shows how the ICON data can be displayed with Python on the
triangular grid. Anyone who has tried to do this with Matplotlib's
`matplotlib.pyplot.tripcolor()` routine will have been surprised by the
artfully nice visualizations. The so-called 'wrapping of lines' of the data
once around is suppressed here by masking the triangles. The triangulation is
generated using Matplotlib's `tri.Triangulation()`, and for plotting
the colored triangles we use the `pyplot.tripcolor()` function.
Inspired by
- Matplotlib Tripcolor Demo
https://matplotlib.org/stable/gallery/images_contours_and_fields/tripcolor_demo.html
- easy.gems Triangular Meshes and Basic
https://easy.gems.dkrz.de/Processing/playing_with_triangles/tripcolor.html
ICON Model Tutorial
https://www.dwd.de/DE/leistungen/nwv_icon_tutorial/pdf_einzelbaende/icon_tutorial2024.pdf;jsessionid=D9D8D98D646668F16569C7F4136E8E32.live21072?__blob=publicationFile&v=3
Content
1. Open ICON data and grid files
2. Coordinates
3. Mask
4. Variable data
5. Grid area
6. Data range and colors
7. Triangulation
8. Plotting with coordinates in radians
9. Plotting on a map
10. Plotting using projections
10.1 Orthographic
10.2 Mollweide
10.3 Robinson
10.4 North polar stereographic
-------------------------------------------------------------------------------
2024 copyright DKRZ licensed under CC BY-NC-SA 4.0 <br>
(https://creativecommons.org/licenses/by-nc-sa/4.0/deed.en)
-------------------------------------------------------------------------------
'''
import time, os
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import matplotlib.colors as mcolors
import matplotlib.path as mpath
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def main():
plt.switch_backend('agg')
#-- Let's see how long it takes.
t1 = time.time() #-- retrieve start time
#-- Open ICON data and grid files
fname = '../../data/ta_ps_850.nc' #-- data file
gname = '../../data/r2b4_amip.nc' #-- grid file
varname = 'ta' #-- variable name
ds_data = xr.open_dataset(fname)
ds_grid = xr.open_dataset(gname)
#-- Have a look at the different shapes.
print(f'vlon/vlat: {ds_grid.vlon.shape}, {ds_grid.vlat.shape}')
print(f'vertex_of_cell: {ds_grid.vertex_of_cell.shape}')
print(f'variable: {ds_data[varname].shape}')
#-- Data and coordinates
#
# vlon, vla: longitude, latitude of the triangle vertices in radians
vlon = ds_grid.vlon
vlat = ds_grid.vlat
# Variable data
#
# Now, we can use the mask array to extract the data variable for only wanted
# cells. In the following the first time step will be used and the data will
# be converted from K to °C.
var = ds_data[varname].isel(time=0).squeeze()
var = var - 273.15
print(f'Variable shape: {var.size}')
# Grid area
#
# vertex_of_cell:
# The indices vertex of cell(:,i) denote the vertices that belong to the triangle i.
# The vertex of cell index array is ordered counter-clockwise for each cell.
voc = ds_grid.vertex_of_cell.T.values - 1
print(f'vertex_of_cell shape: {voc.shape}')
print(f'First indices vertex of cells: \n{voc[:3]}')
# Data range and colors
#
# Set the contour level range, increment, and labels.
var_min = -32.
var_max = 28.
var_int = 2
levels = np.arange(var_min, var_max+var_int, var_int)
nlevs = levels.size
labels = ['{:.2f}'.format(x) for x in levels]
print(f'levels: {levels}\nnumber of levels: {nlevs:3d}')
#-- Colors
#
# Define the color map for the contour levels. Assign the color indices of
# all cells in between the contour levels.
#
# The color assignment of the values time-consuming.
cmap = plt.cm.Spectral_r(range(nlevs))
colors = np.zeros([var.size, 4], np.float32) #-- assign color array for triangle colors
for m in range(nlevs-1):
vind = 0
for i in range(var.size-2):
if (var[i] >= levels[m] and var[i] < levels[m+1]):
colors[i,:] = cmap[m]
vind += 1
print(f'level {m:3d}: {levels[m]:>5.1f} < var > {levels[m+1]:>5.1f} -- {vind:5d} polygons considered')
# set values below and above value range
colors[np.where(var < var_min),:] = cmap[0]
colors[np.where(var >= var_max),:] = cmap[-1]
#-- Triangulation
#
# The triangulation for plotting can be generated with Matplotlib's
# `tri.Triangulation()` class. It takes the longitude and latitude vertices,
# and the vertex of cell indices to compute the triangles.
#
# To avoid wrapping of the lines, a mask is also created here and assigned to
# the triangulation object with the set_mask() method. See also
# https://stackoverflow.com/questions/52457964/how-to-deal-with-the-undesired-triangles-that-form-between-the-edges-of-my-geo
triang = tri.Triangulation(vlon.values.flatten(), vlat.values.flatten(), voc)
triangles = triang.triangles
# generate the mask
mlon = vlon.values.flatten()[triangles] - np.roll(vlon.values.flatten()[triangles], 1, axis=1)
mlat = vlat.values.flatten()[triangles] - np.roll(vlat.values.flatten()[triangles], 1, axis=1)
max_value = np.max(np.sqrt(mlon**2 + mlat**2), axis=1)
max_radius = 1
# assign mask to triang
triang.set_mask(max_value > max_radius)
#-- Plotting with coordinates in radians
#
# The vlon and vlat values are in radians and we use Matplotlib's
# `pyplot.tripcolor()` function to generate the plot showing the colored
# triangles.
fig, ax = plt.subplots(figsize=(10,6))
ax.grid()
plot = ax.tripcolor(triang, var.values,
cmap='Spectral_r',
color=colors,
edgecolor='k',
linewidth=0.05)
#-- add a color bar
cb = plt.cm.ScalarMappable(cmap='Spectral_r',
norm=plt.Normalize(vmin=var_min, vmax=var_max))
cb.set_array([])
cbar = plt.colorbar(cb,ax=ax,
orientation='horizontal',
ticks=levels,
boundaries=levels,
format='%0.0f',
shrink=0.6,
pad=0.06,
aspect=30)
plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
cbar.set_label('[deg C]')
ax.set_title('ICON: coordinates in radians', fontsize=16, weight='bold')
#-- save figure to PNG file
fig.savefig('plot_triangulation_example_ICON_in_radians.png',
bbox_inches='tight', facecolor='white', dpi=150);
#-- Plotting on a map
#
# Now we don't want to work with radians anyway, but with degrees and a map
# instead, which is why the vlon and vlat values have to be converted to
# degrees. Of course, the triangles have to be recalculated but we can reuse
# the boolean mask from above.
triang2 = tri.Triangulation(np.rad2deg(vlon).values.flatten(), np.rad2deg(vlat).values.flatten(), voc)
triangles2 = triang2.triangles
# assign boolean mask from above
triang2.set_mask(max_value > max_radius)
# There are two keywords used by plotting routines for the different coordinate
# reference systems (CRS), `projection` and `transform`, that are important to
# understand. One is the **data CRS** which is set with the `transform` keyword,
# and the **map projection CRS** which is set with the `projection` keyword.
#
# Read more about 'Understanding the transform and projection keywords' in the
# Cartopy documentation:
# https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html
#
# Here, the ICON data CRS can be set to PlateCarree (cylindrical equidistant projection).
data_crs = ccrs.PlateCarree()
# The map projection in Cartopy for a cylindrical equidistant projection is
# called PlateCarree.
#
# Set the map projection CRS
projection = ccrs.PlateCarree()
# Generate the plot
fig, ax = plt.subplots(figsize=(14,7), subplot_kw=dict(projection=projection))
ax.coastlines()
ax.gridlines(draw_labels=True)
plot = ax.tripcolor(triang2, var.values,
cmap='Spectral_r',
color=colors,
edgecolor='k',
linewidth=0.05,
transform=data_crs)
# add a color bar
cb = plt.cm.ScalarMappable(cmap='Spectral_r',
norm=plt.Normalize(vmin=var_min, vmax=var_max))
cb.set_array([])
cbar = plt.colorbar(cb,ax=ax,
orientation='horizontal',
ticks=levels,
boundaries=levels,
format='%0.0f',
shrink=0.6,
pad=0.06,
aspect=30)
plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
cbar.set_label('[deg C]')
ax.set_title('ICON: coordinates in degrees', fontsize=16, weight='bold')
# save figure to PNG file
fig.savefig('plot_triangulation_example_ICON_in_degrees.png',
bbox_inches='tight', facecolor='white', dpi=150);
#-- Plotting using projections
#
# If a projection other than `ccrs.PlateCarree()` is used, the `ax.set_global()`
# or `ax.set_extent()` must be set and the `transform=data_crs` has to be added
# to the `ax.tripcolor()` call.
#
# Read more about 'Understanding the transform and projection keywords' in the
# Cartopy documentation:
# https://scitools.org.uk/cartopy/docs/latest/tutorials/understanding_transform.html
#
# Set the data CRS
data_crs = ccrs.PlateCarree()
#-- Orthographic
projection = ccrs.Orthographic()
# Create plot
fig, ax = plt.subplots(figsize=(8,8), subplot_kw=dict(projection=projection))
ax.set_global() #-- has to be set
ax.coastlines()
ax.gridlines(draw_labels=True)
plot = ax.tripcolor(triang2, var.values,
cmap='Spectral_r',
color=colors,
edgecolor='k',
linewidth=0.05,
transform=data_crs)
# add a color bar
cb = plt.cm.ScalarMappable(cmap='Spectral_r',
norm=plt.Normalize(vmin=var_min, vmax=var_max))
cb.set_array([])
cbar = plt.colorbar(cb,ax=ax,
orientation='horizontal',
ticks=levels,
boundaries=levels,
format='%0.0f',
shrink=0.7,
pad=0.06,
aspect=30)
plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
cbar.set_label('[deg C]')
ax.set_title('ICON: with Orthographic projection', fontsize=16, weight='bold')
# save figure to PNG file
fig.savefig('plot_triangulation_example_ICON_ortho_proj.png',
bbox_inches='tight', facecolor='white', dpi=150);
#-- Mollweide
#
# Set the map projection CRS
projection = ccrs.Mollweide()
# Create the plot
fig, ax = plt.subplots(figsize=(8,8), subplot_kw=dict(projection=projection))
ax.set_global() #-- has to be set
ax.coastlines()
ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)
plot = ax.tripcolor(triang2, var.values,
cmap='Spectral_r',
color=colors,
edgecolor='k',
linewidth=0.05,
transform=data_crs)
# add a color bar
cb = plt.cm.ScalarMappable(cmap='Spectral_r',
norm=plt.Normalize(vmin=var_min, vmax=var_max))
cb.set_array([])
cbar = plt.colorbar(cb,ax=ax,
orientation='horizontal',
ticks=levels,
boundaries=levels,
format='%0.0f',
shrink=0.7,
pad=0.04,
aspect=30)
plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
cbar.set_label('[deg C]')
ax.set_title('ICON: with Mollweide projection', fontsize=16, weight='bold')
# save figure to PNG file
fig.savefig('plot_triangulation_example_ICON_mollweide_proj.png',
bbox_inches='tight', facecolor='white', dpi=150);
#-- Robinson
#
# Set the map projection CRS
projection = ccrs.Robinson()
# Create the plot
fig, ax = plt.subplots(figsize=(8,8), subplot_kw=dict(projection=projection))
ax.set_global() #-- has to be set
ax.coastlines()
ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)
plot = ax.tripcolor(triang2, var.values,
cmap='Spectral_r',
color=colors,
edgecolor='k',
linewidth=0.05,
transform=data_crs)
# add a color bar
cb = plt.cm.ScalarMappable(cmap='Spectral_r',
norm=plt.Normalize(vmin=var_min, vmax=var_max))
cb.set_array([])
cbar = plt.colorbar(cb,ax=ax,
orientation='horizontal',
ticks=levels,
boundaries=levels,
format='%0.0f',
shrink=0.7,
pad=0.04,
aspect=30)
plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
cbar.set_label('[deg C]')
ax.set_title('ICON: with Robinson projection', fontsize=16, weight='bold')
# save figure to PNG file
fig.savefig('plot_triangulation_example_ICON_robinson_proj.png',
bbox_inches='tight', facecolor='white', dpi=150);
#-- Run time
t2 = time.time() #-- retrieve end time
print(f'Run time: {t2 - t1:3.2f}s')
if __name__ == '__main__':
main()
Plot result:
