Python unstructured ICON triangles plot (python 3)

Software requirements:

• Python 3

• Numpy

• matplotlib

• cartopy

• xarray

  • ncells: 20480

Example script

DKRZ_ICON_triangles.py

#!/usr/bin/env python
# coding: utf-8
'''
DKRZ example 

Draw ICON data on original grid (triangles).

Content

- read ICON model data
- use variable 'ta'
- use cell vertices
- draw the triangles with PolyCollection
- add a colorbar
- save to PNG 

-------------------------------------------------------------------------------
2021 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 numpy.ma as ma

import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.tri as tri
from   matplotlib.collections import PolyCollection

import cartopy.crs as ccrs
import cartopy.feature as cfeature

def main():
    plt.switch_backend('agg')

    t1 = time.time()    #-- retrieve start time

    #-- set title string
    title = 'ICON triangles plot'

    #--  define data and grid path, and variable name
    varName   = 'ta'

    fname = '../../data/ta_ps_850.nc'
    gname = '../../data/r2b4_amip.nc'
    
    ds     = xr.open_dataset(fname)
    dsgrid = xr.open_dataset(gname)

    #-- get variable
    var = ds[varName][0,0,:].values
    var = var - 273.15

    #-- get coordinates and convert radians to degrees
    clon = np.rad2deg(dsgrid.clon.values)
    clat = np.rad2deg(dsgrid.clat.values)
    clon_vertices = np.rad2deg(dsgrid.clon_vertices.values)
    clat_vertices = np.rad2deg(dsgrid.clat_vertices.values)

    ncells, nv = clon_vertices.shape[0], clon_vertices.shape[1]

    #-- set contour levels, labels
    varMin, varMax, varInt = -32, 28, 2
    levels = np.arange(varMin, varMax+varInt, varInt)
    nlevs  = levels.size
    labels = ['{:.2f}'.format(x) for x in levels]

    #-- print information to stdout
    print('')
    print('Cells:            %6d ' % clon.size)
    print('Variable min/max: %6.2f ' % np.nanmin(var)+'/'+' %.2f' % np.nanmax(var))
    print('Contour  min/max: %6.2f ' % varMin+'/'+' %.2f' % varMax)
    print('')

    #-- set projection
    projection = ccrs.PlateCarree()

    #-- create figure and axes instances; we need subplots for plot and colorbar
    fig, ax = plt.subplots(figsize=(10,10), subplot_kw=dict(projection=projection))

    ax.set_global()

    #-- plot land areas at last to get rid of the contour lines at land
    ax.coastlines(linewidth=0.5, zorder=2)
    ax.gridlines(draw_labels=True, linewidth=0.5, color='dimgray', alpha=0.4,
                 zorder=2)

    #-- plot the title string
    plt.title(title)

    #-- define color map
    cmap     = plt.get_cmap('Spectral_r', nlevs)        #-- read the color map
    cmaplist = [i for i in range(cmap.N)]               #-- color bar indices
    ncol     = len(cmaplist)                            #-- number of colors
    colors   = np.ones([ncells,4], np.float32)       #-- assign color array for triangles

    print('levels:      ',levels)
    print('nlevs:       %3d' %nlevs)
    print('ncol:        %3d' %ncol)
    print('')

    #-- set color index of all cells in between levels
    for m in range(0,ncol-1):
        vind = []
        for i in range(0,ncells-2, 1):
            if (var[i] >= levels[m] and var[i] < levels[m+1]):
               colors[i,:] = cmap(cmaplist[m])
               vind.append(i)
        print('set colors: finished level %3d' % m ,
              ' -- %5d ' % len(vind) ,
              ' polygons considered')
        del vind

    colors[np.where(var < varMin),:]  = cmap(cmaplist[0])
    colors[np.where(var >= varMax),:] = cmap(cmaplist[ncol-1])    
    
    #-- create the triangles
    clon_vertices = np.where(clon_vertices < -180., clon_vertices + 360., clon_vertices)
    clon_vertices = np.where(clon_vertices >  180., clon_vertices - 360., clon_vertices)

    triangles = np.zeros((ncells, nv, 2), np.float32)

    for i in range(0, ncells, 1):
        triangles[i,:,0] = np.array(clon_vertices[i,:])
        triangles[i,:,1] = np.array(clat_vertices[i,:])

    print('')
    print('--> triangles done')

    #-- create polygon/triangle collection
    coll = PolyCollection(triangles, array=None, fc=colors, edgecolors='black',
                          linewidth=0.05, transform=ccrs.Geodetic(), zorder=0)
    ax.add_collection(coll)

    print('--> polygon collection done')

    #-- add a color bar
    cb = plt.cm.ScalarMappable(cmap=cmap,
                               norm=plt.Normalize(vmin=varMin, vmax=varMax))
    cb.set_array([])
    cbar = plt.colorbar(cb,
                        orientation='horizontal',
                        ticks=levels,
                        boundaries=levels,
                        format='%0.0f',
                        shrink=0.8,
                        pad=0.04,
                        aspect=30,
                        )
    plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
    cbar.set_label('[deg C]')

    #-- maximize and save the PNG file
    plt.savefig('plot_ICON_triangles.png', bbox_inches='tight',dpi=150)

    #-- get wallclock time
    t2 = time.time()
    print('Wallclock time:  %0.3f seconds' % (t2-t1))
    print('')


if __name__ == '__main__':
    main()