CliMAF tour

CliMAF is a framework for Climate Model Assessment.
It allows to combine diagnostic binaries, scripts and Python functions as easily as possible
It is documented on ReadTheDocs

1 - Init CliMAF, and resets its cache

In [1]:
from climaf.api import *
craz()
Climaf version = 0.12
Cache directory set to : ~/tmp/climaf_cache (use $CLIMAF_CACHE if set) 
Available macros read from ~/.climaf.macros are : []

Depending on the depth of the demo , one may tell CliMAF to be more or less verbose (both on stderr and on logging file). Levels are : critical < error < warning (stderr default) < info (log file default) < debug

In [2]:
#clog('info')
#clog_file('debug')

2 - Data location and datasets

A 'project' means a data organization and a number of specific facets.
It can have multiple data location roots.
Let us declare project 'example' with only one additional facet.
Note : All that is usually defined done in config files.

In [3]:
cproject("example" , ("frequency","monthly") )
print cprojects['example'].facets
[__init__   : classes.py, L. 98  ] : WARNING  : Redefining project example

['project', 'simulation', 'variable', 'period', 'domain', 'frequency']

For this project, we will use data in the distro.

In [4]:
dataloc(project="example",organization="generic",url=[\
        cpath+"/../examples/data/${simulation}/L/${simulation}SFXYYYY.nc", \
        cpath+"/../examples/data/${simulation}/A/${simulation}PLYYYY.nc"])
Out[4]:
<climaf.dataloc.dataloc instance at 0x7fc2f6738758>

Define a default value for some dataset facets

In [5]:
cdef("project","example")
cdef("period","1980-1981")

Define a dataset based on a given simulation, and one variable (this uses the facets default values)

In [6]:
tas_ds=ds(simulation="AMIPV6ALB2G", variable="tas")

You can also define a dataset restricted to a latlon box

In [7]:
tas_ds_r=ds(simulation="AMIPV6ALB2G", variable="tas", domain=[10,80,-50,40])

Have a look at dataset's symbolic expression (other patterns could apply , according to project's "DRS")

In [8]:
tas_ds_r
Out[8]:
ds('example.AMIPV6ALB2G.tas.1980-1981.[10, 80, -50, 40].monthly')

CliMAF can identify all the files relevant to the dataset (see CMIP5 DRS case further below, § 10)

In [9]:
tas_ds.baseFiles()
Out[9]:
'/cnrm/aster/data3/aster/vignonl/code/climaf/examples/data/AMIPV6ALB2G/L/AMIPV6ALB2GSFX1980.nc /cnrm/aster/data3/aster/vignonl/code/climaf/examples/data/AMIPV6ALB2G/L/AMIPV6ALB2GSFX1981.nc /cnrm/aster/data3/aster/vignonl/code/climaf/examples/data/AMIPV6ALB2G/L/AMIPV6ALB2GSFX1981.nc'

It can also hande other organizations : CMIP5 DRS (see below), EM, CAMIobs and OBS4MIPs at CNRM, OCMIP5 on IPSL's Ciclad

To TOC

3 - Using external scripts and binaries for implementing diagnostics

First diagnostic is basic climatology - Let us declare a 'CliMAF operator' for this diag, using an external binary and a few keywords conventions

In [10]:
cscript('time_mean' ,'cdo timmean ${in} ${out}' )
Out[10]:
CliMAF operator : time_mean

Declaring a 'CliMAF operator' actually creates a Python function

In [11]:
time_mean
Out[11]:
<function climaf.operators.time_mean>

A doc string is available for each 'operator function' as long a Rest syntax doc file exists with same name. Here, we just have the default doc.

In [12]:
help(time_mean)
Help on function time_mean in module climaf.operators:

time_mean(*args, **dic)
    CliMAF wrapper for command : cdo timmean ${in} ${out}


Symbolically apply the script to a dataset, using the 'CliMAF operator', and have a look at object result by displaying its expression in CliMAF Reference Syntax (CRS)

In [13]:
tas_avg=time_mean(tas_ds)
tas_avg
Out[13]:
time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'))

Define another operator, apply it on result above, and look at CRS

In [14]:
cscript('space_mean' ,'cdo fldavg ${in} ${out}' )
tas_ga=space_mean(tas_avg)
tas_ga
Out[14]:
space_mean(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')))

To TOC

4 - From CliMAF objects to files or Masked Array

Computation actually occurs when wanting to 'export' objects in some way : here as a file using 'cfile'

In [15]:
cfile(tas_avg)
Done in 0.4 s with script computation for select(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')) 
Done in 0.1 s with script computation for time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')) 

Out[15]:
'/home/vignonl/tmp/climaf_cache/59/0.nc'

The file cache will be used whenever the same object is requested (or even for an included objet). Here, in addition, we request a copy of the cache file, at some location

In [16]:
cfile(tas_ds,"~/tmp/tasavg.nc")
Out[16]:
'/home/vignonl/tmp/tasavg.nc'

You can display the header information of a netCDF representation of a CliMAF object

In [17]:
ncdump(tas_ds)
netcdf c {
dimensions:
	lon = 256 ;
	lat = 128 ;
	time = UNLIMITED ; // (36 currently)
variables:
	float lon(lon) ;
		lon:standard_name = "longitude" ;
		lon:long_name = "longitude" ;
		lon:units = "degrees_east" ;
		lon:axis = "X" ;
	float lat(lat) ;
		lat:standard_name = "latitude" ;
		lat:long_name = "latitude" ;
		lat:units = "degrees_north" ;
		lat:axis = "Y" ;
	double time(time) ;
		time:standard_name = "time" ;
		time:units = "days since 1975-01-01 00:00:00" ;
		time:calendar = "standard" ;
	float tas(time, lat, lon) ;
		tas:long_name = "2m_Temperature" ;
		tas:units = "K" ;
		tas:grid_type = "gaussian" ;
		tas:_FillValue = 1.e+20f ;
		tas:missing_value = 1.e+20f ;
		tas:comment = "instantaneous values" ;

// global attributes:
		:CDI = "Climate Data Interface version 1.6.9 (http://mpimet.mpg.de/cdi)" ;
		:Conventions = "CF-1.0" ;
		:history = "Mon Jan 04 17:10:33 2016: cdo copy climaf_mcdo_LsxkP8/AMIPV6ALB2GSFX1980.nc climaf_mcdo_LsxkP8/AMIPV6ALB2GSFX1981.nc climaf_mcdo_LsxkP8/AMIPV6ALB2GSFX1981.nc climaf_mcdo_LsxkP8/mcdo.sh.nc\n",
			"Mon Jan 04 17:10:33 2016: cdo -f nc seldate,1980-01-01T00:00:00,1981-12-31T23:59:00 -selname,tas /cnrm/aster/data3/aster/vignonl/code/climaf/examples/data/AMIPV6ALB2G/L/AMIPV6ALB2GSFX1980.nc climaf_mcdo_LsxkP8/AMIPV6ALB2GSFX1980.nc\n",
			"Mon Feb 16 16:03:44 2015: cdo selvar,tas,tasmin,tasmax,hfls,hfss AMIPV6ALB2GSFX1980.nc AMIPV6ALB2GSFX1980e.nc\n",
			"Thu Apr 10 07:47:44 2014: ncrcat -n 12,2,1 AMIPV6ALB2GSFX198001.nc AMIPV6ALB2GSFX1980.nc\n",
			"created at time: Thu Apr 10 07:47:07 2014" ;
		:institution = "CNRM (Centre National de Recherches Meteorologiques, Meteo-France, Toulouse, France)" ;
		:title = "[EXPE]" ;
		:nco_openmp_thread_number = 1 ;
		:CDO = "Climate Data Operators version 1.6.9 (http://mpimet.mpg.de/cdo)" ;
		:CRS_def = "select(ds(\'example.AMIPV6ALB2G.tas.1980-1981.global.monthly\'))" ;
		:CliMAF = "CLImate Model Assessment Framework version 0.12 (http://climaf.rtfd.org)" ;
}

Out[17]:
ncdump(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'))

Another export format : MaskedArrays (others could apply, as e.g. CDAT's MV)

In [18]:
tas_avg_MA=cMA(tas_avg)

Looking at the result

In [19]:
type(tas_avg_MA)
Out[19]:
numpy.ma.core.MaskedArray
In [20]:
tas_avg_MA.shape
Out[20]:
(1, 128, 256)

To TOC

5 - Basic and advanced displays

Plot the data, first way : let us declare a 'CliMAF operator', for using binary 'ncview', which has no managed output

In [21]:
cscript('my_ncview' ,'ncview ${in}')
Out[21]:
CliMAF operator : my_ncview

Using this operator is straightforward

In [22]:
my_ncview(tas_avg)
Out[22]:
my_ncview(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')))

2nd plot : declare a simplified script plot with PNG output. Here in NCL with a number of args for graphics settings

In [23]:
cscript('plotmap',\
        "ncl "+ cpath +"/../scripts/plotmap.ncl infile=${in} plotname=${out} \
                    cmap=${color} vmin=${min} vmax=${max} vdelta=${delta} var=${var} \
                    title=${title} scale=${scale} offset=${offset} units=${units}",format="png")
Out[23]:
CliMAF operator : plotmap

Applying the operator (actually, let us use a richer one: 'plot', a general purpose plot operator for plotting 1D and 2D datasets - maps, cross-sections and profiles -)
The only effect is to define a figure object

In [24]:
map=plot(tas_avg, min=260, max=300, delta=4)
map
Out[24]:
plot(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')),delta=4,max=300,min=260)

Require the figure as a file; this leads to actually compute it if needed (i.e. launch the script)

In [25]:
cfile(map)
Done in 1.0 s with script computation for plot(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')),delta=4,max=300,min=260) 

Out[25]:
'/home/vignonl/tmp/climaf_cache/e8/8.png'

'cshow' will display the figure using the file cache [ Nothing will show on Web-only demo]

In [26]:
cshow(map)

To TOC

6 - Mixing CliMAF basics and Python (and/or IPython)

A more realistic way of configuring the NCL-based plots, using Python features for graphics settings on the basis of on plotted variable

First define an auxiliary function

In [27]:
def map_graph_attributes(var) :
    rep=dict()
    rep["offset"]=0.
    rep["scale"]=1.
    rep["color"]="BlueDarkRed18"
    if var=='tas' :
        rep["offset"]=-273.15 ;  rep["scale"]=1.0 ; rep["units"]="C"
        rep["min"]=-30  ; rep["max"]=30 ; rep["delta"]=2.
    return rep

This allows to supply graphics attributes which are better tuned to the variable to display

In [28]:
map2=plot(tas_avg,title="The Title",**map_graph_attributes(varOf(tas_avg)))
map2
Out[28]:
plot(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')),color='BlueDarkRed18',delta=2.0,max=30,min=-30,offset=-273.15,scale=1.0,title='The Title',units='C')

Let launch the computation of the figure

In [29]:
figfile=cfile(map2) ; print(figfile)
/home/vignonl/tmp/climaf_cache/ce/2.png

Done in 1.0 s with script computation for plot(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')),color='BlueDarkRed18',delta=2.0,max=30,min=-30,offset=-273.15,scale=1.0,title='The Title',units='C') 

Bonus : using Ipython Notebook features for inline display

In [30]:
from IPython.display import Image
Image(filename=figfile)
Out[30]:

To TOC

7 - Remapping data

For remapping a dataset to a named grid (known by cdo), let us use operator 'regridn'

In [31]:
dgr=regridn(tas_avg,cdogrid="r90x45") 
mapgr=plot(dgr, title="2° Grid", min=260, max=300, delta=4)
Image(filename=cfile(mapgr))
Done in 0.1 s with script computation for regridn(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')),cdogrid='r90x45') 
Done in 0.9 s with script computation for plot(regridn(time_mean(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly')),cdogrid='r90x45'),delta=4,max=300,min=260,title='2\xc2\xb0 Grid') 

Out[31]:

An other way is to regrid to the grid of another dataset, with 'regrid'; Let us use the dataset defined earlier on a latlon box

In [32]:
dgr2=regrid(tas_ds, tas_ds_r) 
ncview(dgr2)
Done in 0.3 s with script computation for select(ds('example.AMIPV6ALB2G.tas.1980-1981.[10, 80, -50, 40].monthly')) 
Done in 0.1 s with script computation for regrid(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),ds('example.AMIPV6ALB2G.tas.1980-1981.[10, 80, -50, 40].monthly')) 

Out[32]:
ncview(regrid(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),ds('example.AMIPV6ALB2G.tas.1980-1981.[10, 80, -50, 40].monthly')))

To TOC

8 - Computing an annual cycle

Let us compute an annual cycle, using swiss knife operator 'ccdo' (which applies cdo on a single dataset or object, with a cdo operator as argument), and look at it (nothing will show on Web only demo)

In [33]:
anncycle=ccdo(tas_ds,operator='ymonavg')
ncview(anncycle)
Done in 0.4 s with script computation for ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg') 

Out[33]:
ncview(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'))

Define the average annual cycle on a latlon box using operator 'llbox', plot it (January will show), and display the plot

In [34]:
extract=llbox(anncycle, latmin=30, latmax=60, lonmin=-30, lonmax=30)
map_extract=plot(extract, title='extract', min=260, max=300, delta=4)
Image(filename=cfile(map_extract))
Done in 0.1 s with script computation for llbox(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),latmax=60,latmin=30,lonmax=30,lonmin=-30) 
Done in 0.9 s with script computation for plot(llbox(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),latmax=60,latmin=30,lonmax=30,lonmin=-30),delta=4,max=300,min=260,title='extract') 

Out[34]:

Let's compute the space average of annual cycle

In [35]:
space_average=ccdo(extract,operator='fldavg')

Creating a figure with standard operator 'curves'

In [36]:
fig_avg=curves(space_average, title="Annual cycle")

Get the figure computed, and get its filename in CliMAF file cache

In [37]:
Image(filename=cfile(fig_avg))
Done in 0.1 s with script computation for ccdo(llbox(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),latmax=60,latmin=30,lonmax=30,lonmin=-30),operator='fldavg') 
Done in 0.8 s with script computation for curves(ccdo(llbox(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),latmax=60,latmin=30,lonmax=30,lonmin=-30),operator='fldavg'),title='Annual cycle') 

Out[37]:

Computation of seasonal average and plot bias

In [38]:
# obs
if atCNRM : obs=ds(project='cruts3',variable='tas', period='1980-1981',grid='T127')
else :
    if onCiclad : obs=ds(project='ref_ipsl', variable='tas', period='1980-1981', product='ERAI')
    else :
        print("I do not know how to find data on this machine")
        exit(0)
anncycle_obs=ccdo(obs,operator='ymonavg')
In [39]:
# seasonal average
JFM=ccdo(anncycle,operator='selmon,1,2,3')
JAS=ccdo(anncycle,operator='selmon,7,8,9')

JFM_obs=ccdo(anncycle_obs,operator='selmon,1,2,3')
JAS_obs=ccdo(anncycle_obs,operator='selmon,7,8,9')
In [40]:
# bias
JFM_regrid=regrid(JFM,JFM_obs)
JFM_bias=minus(JFM_regrid,JFM_obs)

JAS_regrid=regrid(JAS,JAS_obs)
JAS_bias=minus(JAS_regrid,JAS_obs)
In [41]:
# plot
JFM_map=plot(JFM,title='JFM',contours=1)
JAS_map=plot(JAS,title='JAS', contours=1)

#JFM_obs_map=plot(JFM_obs,title='JFM')
#JAS_obs_map=plot(JAS_obs,title='JAS')

JFM_bias_map=plot(JFM_bias,title='Bias JFM')
JAS_bias_map=plot(JAS_bias,title='Bias JAS')

multiplot=cpage([[JFM_map, JAS_map],[JFM_bias_map,JAS_bias_map]],heights=[0.2, 0.2])
Image(filename=cfile(multiplot))
Done in 0.2 s with script computation for ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,1,2,3') 
Done in 1.0 s with script computation for plot(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,1,2,3'),contours=1,title='JFM') 
Done in 0.2 s with script computation for ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,7,8,9') 
Done in 1.0 s with script computation for plot(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,7,8,9'),contours=1,title='JAS') 
						[selectGene : dataloc.py, L. 334 ] : WARNING  : Cannot yet filter files re. time using only file content. TBD
Done in 1.9 s with script computation for ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg') 
Done in 0.2 s with script computation for ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,1,2,3') 
Done in 0.1 s with script computation for regrid(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,1,2,3'),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,1,2,3')) 
Done in 0.1 s with script computation for minus(regrid(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,1,2,3'),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,1,2,3')),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,1,2,3')) 
Done in 1.0 s with script computation for plot(minus(regrid(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,1,2,3'),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,1,2,3')),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,1,2,3')),title='Bias JFM') 
Done in 0.2 s with script computation for ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,7,8,9') 
Done in 0.1 s with script computation for regrid(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,7,8,9'),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,7,8,9')) 
Done in 0.1 s with script computation for minus(regrid(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,7,8,9'),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,7,8,9')),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,7,8,9')) 
Done in 1.0 s with script computation for plot(minus(regrid(ccdo(ccdo(ds('example.AMIPV6ALB2G.tas.1980-1981.global.monthly'),operator='ymonavg'),operator='selmon,7,8,9'),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,7,8,9')),ccdo(ccdo(ds('cruts3.N/A.tas.1980-1981.global.T127'),operator='ymonavg'),operator='selmon,7,8,9')),title='Bias JAS') 

Out[41]: