CosmoState¶
A high-level python module for reading GRIB1 and GRIB2 COSMO output files based on the grib_api/eccodes python module.¶
Grib_api
: Standard library for reading/writing grib files developed at ECMWF,
replaces DWDs grib library
eccodes : successor of grip_api
Why not use CosmoUtils? No support for GRIB2 (based on wgrib)...
Why not use grib_api python module directly? It is very low-level, records (=2D-slices of >=3D variables) are accessed sequentially. CosmoState allows the user to access model variables without having to deal with low-level technical details...
%matplotlib inline
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
Getting started¶
from cosmo_state import CosmoState
Get cosmo_state.py from http://www.meteo.physik.uni-muenchen.de/~leo.scheck/cosmo_state.py and put it somewhere in your $PYTHONPATH
Requirements: grib_api or eccodes python module.
Listing file contents¶
g2=CosmoState('/home/l/Leonhard.Scheck/userdata/grib_examples/grib2_fcst/de_2016060112.grib2') # operational COSMO-DE forecast in GRIB2 format
g2.list_variables(narrow=True) # narrow -> do not print long variable description
To generate this overview, cosmo_state.py can also be run directly from the command line with a grib file name as argument (like a grib_ls replacement). You may need to do "chmod +x cosmo_state.py" before...
Accessing variables and meta information¶
plt.imshow( g2['RELHUM_2M'], origin='lower' ) # variable is identified by 'short name' and read from file when you access it
g2['u'].shape # access the first variable with short name 'u' -- this happens to be u on pressure levels
type(g2['u']), g2['u'].max() # variables = numpy array, remain stored in the object (will not be read every time you access it)
g2['u:M'].shape # append ':M' to get the u on model levels
g2['u:Z'].shape, g2['u:P'].shape # ':Z' means geometric height, ':P' pressure levels (see source code for all types)
Excerpt from cosmo_state.py source:
# level types used in Grib1 and Grib2 files + level type classes:
# M = model levels
# P = pressure l.
# Z = geometrical height
# S = surface
# 0 = mean sea level
# L = special level, e.g. cloud top
# X = satellite or other vertically integrating obs.
level_types = { 'M':['hybrid','hybridLayer','generalVertical','generalVerticalLayer'],
'P':['isobaricInhPa'],
'Z':['depthBelowLand','depthBelowLandLayer','heightAboveSea','heightAboveGround'],
'S':['surface'],
'0':['meanSea'],
'L':['isothermZero','cloudTop','cloudBase','nominalTop'],
'U':['unknown'],
'X':['TOA','meanLayer','entireAtmosphere','isobaricLayer'] }
Different "type of level" names in GRIB1 and GRIB2, different names for levels and layers (half-levels) -> use classes M, P, Z...
g2.meta['RELHUM_2M'] # meta information is stored in .meta, can be accessed once the variable has been accessed
g2.meta['u:P']['levels'] # u:P is available on which pressure levels?
s,pres,n = g2.distribution( 'pres:M', 't:M' )
plt.plot( s['mean'], pres/1e2, 'k' )
plt.plot( s['mean']+0.5*s['std'], pres/1e2, '--r')
plt.plot( s['mean']-0.5*s['std'], pres/1e2, '--r')
plt.xlabel( g2.meta['t:M']['name'] )
plt.ylabel( g2.meta['pres:M']['name'] )
plt.ylim((1000,0))
A more realistic example¶
Common situation: Not all grib files contain all variables, e.g. constants are stored only in first output file -> specify second grib file that will be searched for variables not found in the first file
If you know which variables you will need, specify them using the "preload" argument. This will be faster...
Please note: cosmo_state is not yet optimized for speed (work in progress), but at least with "preload" it is quite fast.
cs = CosmoState('/home/userdata/leo.scheck/bacy_ahutt/cosmo_letkf/data/1000.05b/20140516060000/laf20140516080000.det',
constfile='/home/userdata/leo.scheck/bacy_ahutt/cosmo_letkf/data/1000.05b/20140516090000/lff20140516080000.det',
preload=['QC','RLAT','RLON'],verbose=True)
plt.figure(2,figsize=(8,8))
plt.gray()
plt.imshow( np.log10(np.maximum(cs['QC'].max(axis=2),1e-4)), origin='lower' )
plt.contour( cs['RLAT'], levels=np.arange(40,60,2.5), colors='#009900')
plt.contour( cs['RLON'], levels=np.arange(0,30,2.5), colors='#009900')
Time and ensemble dimensions¶
ens=CosmoState('/home/l/Leonhard.Scheck/userdata/grib_examples/grib2_ensemble/2016052800.grib2') # a grib2 file containing 2 variables for 20 members and 25 times
ens['tp'].shape # access all times and members
ens['tp@2'].shape # access time step 2 only (but all members)
ens['tp@2#10'].shape # access time step 2 of member 10
ens['tp#10'].shape # access all time steps of member 10
plt.plot( ens['tp#10'].mean(axis=(0,1)), color='#cc0033', label='member #10' )
plt.plot( ens['tp'].mean(axis=(0,1,3)), color='#0099ff', label='ensemble mean' )
plt.ylabel( '%s [%s]'%(ens.meta['tp']['name'],ens.meta['tp']['units']))
plt.legend(loc="upper left")
