Python matplotlib ICON R02B09 datashader

Software requirements:

• Python 3

• matplotlib

• cartopy

• datashader

Example script

ICON_R02B09_fast_datashader_example.py

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

The example script shows a fast way to plot large ICON files e.g. R02B09.

Number of cells     20971520
Number of vertices  3
Number of levels    77

Thanks to Florian Ziemen/DKRZ and Lukas Kluft/MPI-M.

Content

- read ICON data
- use datashader to create the plot
- save to PNG

-------------------------------------------------------------------------------
2022 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, re
import xarray as xr
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from   cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

import datashader as dsh
from datashader.mpl_ext import dsshow

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

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

    #-- input data and grid files
    icon_file = '../../data/dpp0029_atm_3d_t_ml_20200301T000000Z.nc'
    grid_file = '../../data/icon_grid_0015_R02B09_G.nc'

    #-- open ICON and grid dataset
    ds_icon = xr.open_dataset(icon_file)
    ds_grid = xr.open_dataset(grid_file)

    #-- read variable ta
    var = ds_icon['ta'][:,0,:]

    #-- generate the date string used for the output plot file name of the
    #-- selected timestep
    if (re.search('.%f', ds_icon.time.units.split(' ')[2])):
        date_format = ds_icon.time.units.split(' ')[2].split('.')[0]
    date = pd.to_datetime(ds_icon.time.values, format=date_format)
    date = str(date.strftime('%Y-%m-%d')[0])
    print(date)

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

    #-- transform radians to geodetic coordinates
    coords = projection.transform_points(
                            ccrs.Geodetic(),
                            np.rad2deg(ds_grid.clon.values),
                            np.rad2deg(ds_grid.clat.values))

    #-- create data frame of the variable values and the geodetic coordinates.
    df = pd.DataFrame({'val': var.isel(time=0).values, 'x': coords[:,0], 'y': coords[:,1]})

    #-- choose colormap
    colormap = 'RdYlBu_r'

    #-- create the plot
    fig, ax = plt.subplots(figsize=(30,30), facecolor='white', subplot_kw={"projection": projection})
    #fig.canvas.draw_idle()

    ax.add_feature(cfeature.COASTLINE, linewidth=0.7)

    gl = ax.gridlines(crs=ccrs.PlateCarree(),
                      draw_labels=True,
                      color='black',
                      linewidth=0.5,
                      alpha=0.7,
                      x_inline=False)
    gl.xlocator   = mticker.FixedLocator(range(-180, 180+1, 60))
    gl.xformatter = LONGITUDE_FORMATTER

    artist = dsshow(df,
                    dsh.Point('x', 'y'),
                    dsh.mean('val'),
                    vmin=245,
                    vmax=300,
                    cmap=colormap,
                    plot_width=1600,
                    plot_height=1400,
                    ax=ax)

    fig.colorbar(artist, label='ta [K]', shrink=0.3, pad=0.02)

    plt.title('ICON R02B09 - ta', fontsize=18)
    plt.figtext(0.75, 0.33, "© 2022 DKRZ", ha="right", fontsize=5)
    plt.savefig('plot_ICON_ta_'+date+'_r2b9_dsshow.png', bbox_inches='tight', dpi=150)

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


if __name__ == '__main__':
    main()

Plot result: