Python matplotlib example: unstructured ICON contour plot (python 3)#
Software requirements:
Python 3
Numpy
matplotlib
cartopy
xarray
python-cdo
ncells: 20480
run time: 0.932 s
Script DKRZ_example_ICON_tricontourf.py#
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
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cdo import *
cdo = Cdo()
def main():
t1 = time.time() #-- retrieve start time
#-- set title string
title = 'ICON tricontourf plot'
#-- define path, file and variable name
dir_name = '/Users/k204045/data/ICON/'
file_name = 'ta_ps_850.nc'
grid_name = 'grids/r2b4_amip.nc'
varName = 'ta'
fname = os.path.join(dir_name,file_name)
gname = os.path.join(dir_name,grid_name)
#-- extract subregion of the global grid with CDO
minlon, maxlon, minlat, maxlat = -179.5, 179.5, -89.5, 89.5
#-- make sure that nan, inf, and -inf are set to missing value
cdo_cmd = '-setmisstoc,-9999. -setgrid,'+ gname +' '+ fname +' '
Data = cdo.sellonlatbox('{0},{1},{2},{3}'.format(minlon,maxlon,minlat,maxlat),
input=cdo_cmd,
options='-r',
returnXDataset=True)
#-- get missing_value
missing_value = Data[varName].encoding['missing_value']
#-- get variable
var = Data[varName][0,0,:].values
var = var - 273.15
#-- get coordinates and convert radians to degrees
clon = np.rad2deg(Data.clon.values)
clat = np.rad2deg(Data.clat.values)
#-- 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()
#-- set colormap for nlevs
cmap = plt.get_cmap('Spectral_r', nlevs)
#-- 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.gridlines(draw_labels=True,
linewidth=0.5,
color='dimgray',
alpha=0.4,
zorder=2)
ax.coastlines(linewidth=0.5, zorder=2)
#-- plot the title string
plt.title(title)
#-- contour plot
cnf = ax.tricontourf(clon, clat, var,
vmin=varMin,
vmax=varMax,
levels=levels,
cmap=cmap,
extend='neither',
zorder=0)
#-- add a color bar
cbar_ax = fig.add_axes([0.2, 0.25, 0.6, 0.015], autoscalex_on=True) #-- x,y,w,h
cbar = fig.colorbar(cnf, cax=cbar_ax, orientation='horizontal')
plt.setp(cbar.ax.get_xticklabels()[::2], visible=False)
cbar.set_label('[deg C]')
#-- maximize and save the PNG file
plt.savefig('plot_ICON_tricontourf.png', bbox_inches='tight',dpi=150)
#-- get wallclock time
t2 = time.time()
print('Wallclock time: %0.3f seconds\n' % (t2-t1))
if __name__ == '__main__':
main()