Plotting ICON Topography

Here, we get started with the ICON data. We will do the following

• we will read the ICON grid file: It contains information on center and edge point of the triangular grid of ICON

• we will read the ICON extpar file: It contains e.g. surface height information (and additional all other possible external parameters

• we will learn to how plot an ICON 2D field onto a map

Load Python Libraries

%matplotlib inline

# system libs
import os, sys, glob

# array operators and netcdf datasets
import numpy as np
import xarray as xr

# plotting
import pylab as plt
import seaborn as sns
sns.set_context('talk')

# drawing onto a map
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader

Open Data with xarray

General Paths

exercise_path = '/work/bb1224/2023_MS-COURSE'

data_path = f'{exercise_path}/data'
icon_data_path = f'{data_path}/icon-lem'

This is the content of the data directory:

os.listdir(icon_data_path)

['grids-extpar', 'experiments', 'bc-init']

The grid & extpar data are stored under grids-extpar. There are subdirectories named after the grid used.

grid_name = 'lpz_r2'
grid_path = f'{icon_data_path}/grids-extpar/{grid_name}'

os.listdir( grid_path )

['external_parameter_icon_lpz_r2_dom02_DOM02_tiles.g2',
 'external_parameter_icon_lpz_r2_dom01_DOM01_tiles.extpar.log',
 'external_parameter_icon_lpz_r2_dom02_DOM02_tiles.extpar.log',
 'oflxd21.e23ad4422bde130aa73773d4104d9a6c000.zip',
 'lpz_r2_dom01_DOM01.nc',
 'gridgen.log',
 'gridgen.nml',
 'lpz_r2_dom01_DOM01.html',
 'external_parameter_icon_lpz_r2_dom01_DOM01_tiles.g2',
 'external_parameter_icon_lpz_r2_dom02_DOM02_tiles.nc',
 'bc-grid-lpz_r2_dom01_DOM01.grid.nc',
 'lpz_r2_dom02_DOM02.html',
 'external_parameter_icon_lpz_r2_dom01_DOM01_tiles.nc',
 'bc-grid-lpz_r2_dom02_DOM02.grid.nc',
 'README.TXT',
 'lpz_r2_dom02_DOM02.nc']

OK, this is rather croudy. The files with .ncare netcdf files which contain our data, i.e. grid and extpar.

Perhaps, you see DOM01 & DOM02 in the filenames? Yes?

• DOM01 is our simulation nest

• DOM02 is a second and inner nest which we don’t use for our exercises -> more advanced simulation -> curious?

Grid Files

grid1 = xr.open_dataset( f'{grid_path}/lpz_r2_dom01_DOM01.nc', chunks={} )

# you can see the content if you just comment the line below out
# grid_file1

Extpar Files

extpar1 = xr.open_dataset( f'{grid_path}/external_parameter_icon_lpz_r2_dom01_DOM01_tiles.nc', chunks={} )

Plotting the Topography

Shortcuts for Fieldnames

# read center lon/lat in radiant
clon_rad = grid1['clon']                # center longitude  / rad
clat_rad = grid1['clat']                # center latitutde  / rad

# convert to degrees
clon = np.rad2deg( clon_rad )
clat = np.rad2deg( clat_rad )

# read surface height
h    = extpar1['topography_c']      # topography height / m

Tricontour Plot

ICON is on a triangular grid. Hence, we need a special plotting function for this which is contained in matplotlib

plt.figure( figsize = (10,10))
plt.tricontourf( clon, clat, h, levels = 20, cmap = plt.cm.gist_earth, vmin = -800, vmax = 1200.)

plt.title( 'ICON DOM01', fontweight = 'bold');

OK, this looks already quite nice!

I guess you know * where Leipzig is located? * what the names of the mountain ranges are?

Plotting onto a Map

Wouldn’t it be much cooler to state or district borders in the topography map?

Yes! It would be much easier to locate some places …. In the following, we will use cartopy to draw the topography field onto a map. We will copy the contour plot from above and add some map details.

plt.figure( figsize = (10,10))

# this is our map projection
ax = plt.axes(projection = ccrs.PlateCarree())

plt.tricontourf( clon, clat, h, levels = 20, cmap = plt.cm.gist_earth, vmin = -800, vmax = 1200.)

# and here we draw country borders
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth = 2)


plt.title( 'ICON DOM01', fontweight = 'bold');

OK, next level: include also district borders.

This is a bit more complicated, because in Germany federal states borders are not known to cartopy. Hence, we need to first read this information from disk (no worries, data are provided).

# read shape file
shape_file = f'{data_path}/shapes/vg2500_bld.shp'
reader = shpreader.Reader( shape_file )

#  and create cartopy feature
states_data = list(reader.geometries())
states = cfeature.ShapelyFeature(states_data, ccrs.PlateCarree())

plt.figure( figsize = (10,10))
ax = plt.axes(projection = ccrs.PlateCarree())


ax.tricontourf( clon, clat, h, levels = 20,
                                cmap = plt.cm.gist_earth,
                                vmin = -800,
                                vmax = 1200.)

ax.add_feature(cfeature.BORDERS.with_scale('50m'),  linewidth = 2)

# and here we draw states borders
ax.add_feature(states, facecolor='none', edgecolor='black', linewidth = 0.4)

plt.title( 'ICON DOM01', fontweight = 'bold');

Tasks

Please explore the data with this notebook:

• What else is in the extpar file?

  • Please try to guess which data are there!

  • Take ICON Tutorial as reference!

• Please select one additional extpar field and plot this field with the methods described above!

Additional task:

• Can you reproduce the topography plots with DOM02?