Creating new coords from multiple experiments

In this tutorial, we are going to construct new coord objects and their axes from parameters that vary across a range of experiments. Numerical model runs often involve configuration (control) files that specify the value of parameters used in the runs. Here, we will look at two types of model runs, one where the Drake Passage is open, and one type where it is closed, equilibrated under a range of atmospheric CO2 values. The CO2 concentrations are specified in configuration files that we have saved alongside the Netcdf data in the experiment directories. First, import the modules:


In [1]:
import spacegrids as sg
import numpy as np
import matplotlib.pyplot as plt
import os
figsize(8,4)


Using netCDF4

Obtain a dictionary of available projects for easy reference:


In [2]:
D = sg.info()


Projects rooted in PROJECTS
----------
DP                my_project        tmp              

Initialize a project:


In [3]:
P = sg.project(D['DP'])


Creating exp object GM_C_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_1250ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_1000ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_500ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_1750ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_390ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_750ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_2000ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_1750ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_1500ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_1250ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_1000ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_2000ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_750ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object GM_C_500ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_390ppm_flat from /home/wim/PROJECTS/DP4/
Creating exp object DPclsd_C_1500ppm_flat from /home/wim/PROJECTS/DP4/
Trying to read mask using yt, xt for mask.
...no masks read in project init.

Bring axes names X, Y, Z, T into namespace:


In [4]:
for c in P['GM_C_flat'].axes:
    exec c.name + ' = c'

As can be seen from the names of created experiment objects, we have two types of experiments, the DP type and the GM type. The Drake Passage is closed in the DP type, whereas this ocean gateway is open in the GM type experiment. The experiment names also contain a CO2 concentration value (e.g. 750ppm), indicating the atmospheric CO2 concentration used in that run. In both experiment types, DP and GM, CO2 varies from 280ppm to 2000ppm. The Netcdf output data contains fields in equilibrium with each CO2 value. The CO2 values are not stored in the Netcdf files. Instead, they are specified in the original configuration files, named control.in, used in the experiments:


In [5]:
ls /home/wim/PROJECTS/DP4/GM_C_flat/


control.in*  mk.in*  time_mean.nc

The control.in file looks like this:

$ tail control.in &hlmix / &isopyc slmx=0.01, ahisop=4.e6, athkdf=4.e6, del_dm=0.4e-02, s_dm=0.1e-2 / &ppmix / &smagnl / &adv_q diffactor=50. / &co2 co2ccn=280. / &paleo pyear=1850. / &ice dampice=5., ice_yr=1850. / &veg crops_yr=1850. / &mtlm TIMESTEP=3600., INT_VEG=.true., VEG_EQUIL=.false. / $

The input parameter co2ccn is what we are looking for. Let's examine temperature:


In [6]:
P.load('temp')


OK. Fetched field temp for DPclsd_C_1750ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_1250ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_500ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_2000ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_2000ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_1000ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_1500ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_1750ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_750ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_750ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_1000ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_500ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_390ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_1500ppm_flat. 1 file found.
OK. Fetched field temp for DPclsd_C_flat. 1 file found.
OK. Fetched field temp for GM_C_390ppm_flat. 1 file found.
OK. Fetched field temp for GM_C_1250ppm_flat. 1 file found.

The project insert method takes a function that takes an exper object as input and outputs a list of (name,value) tuple pairs. These pairs are subsequently inserted into the params dictionary attribute of each experiment object by the project insert method. In our case here, we use a function defined in sg that uses the path attribute to find the configuration file for that experiment, parse it, and insert the name/ value pairs into the params dictionary of that experiment. sg.read_control_func is actually a function that yields a function that reads the control file specified in the argument to sg.read_control_func. sg.read_control_func('control.in') is then a function that is passed an exper object, and outputs a list of (name, value) pairs. Such functions can be passed to P.insert to be processed, where for each experiment in the project, the function is executed, and the resulting (name, value) pairs are inserted into the vars or params dictionaries of that experiment. We will encounter another such function, using regular expressions.


In [7]:
P.insert(sg.read_control_func('control.in')) # read, parse and insert parameters from control.in

The value that we saw in the control.in file now appears in the params attribute of P['GM_C_flat']:


In [8]:
P['GM_C_flat'].params['co2ccn'] # obtain the value of the 'co2ccn' param


Out[8]:
280.0

The param2gr project method creates a new field defined on a grid containing a new coord object constructed from a specified parameter:


In [9]:
bigfield_DP = P.param2gr('co2ccn' , lambda x:x['temp'], 'DP*')  # do this for the DP type

Here, the parameter is co2ccn and the second argument is a function acting on an exper object (experiment), simply yielding the temp field contained by that exper: in general, this function must take an experiment argument and yield a field object. A simple wildcard filter is applied in the last argument to obtain only results for the DP type runs. Similar for the GM type runs:


In [10]:
bigfield_GM = P.param2gr('co2ccn', lambda x:x['temp'], 'GM*') # do this for the GM type

This new bigger field bigfield_DP, constructed from temp fields across a range of experiments has a 4 dimensional grid that includes the newly constructed co2ccn_crd coordinate:


In [11]:
bigfield_DP.gr


Out[11]:
(co2ccn_crd, zt, yt, xt)

This new coord contains the CO2 values that we saw in the experiments:


In [12]:
bigfield_DP.gr[0][:]


Out[12]:
array([  280.,   390.,   500.,   750.,  1000.,  1250.,  1500.,  1750.,
        2000.])

Both bigfield fields contain temperature values over a range of CO2 concentrations. We can examine how the average bottom ocean temperatures (grid cells 17 to 19) depend on CO2 for the DP open and closed type experiments:


In [13]:
h1, = sg.plot(bigfield_DP[Z,17:]/(X*Y*Z) + sg.zero_kelvin)
h2, = sg.plot(bigfield_GM[Z,17:]/(X*Y*Z) + sg.zero_kelvin)
lbl = plt.ylabel('C')
lbl = plt.xlabel('co2 ccn')
tle = plt.title('Bottom ocean temperature for DP open and closed')
lgnd = plt.legend([h1,h2],["DP closed","DP open"],loc=2)
plt.grid()


 spacegrids/plotting.py:400: UserWarning:Label units not added: not a string. 

Deep ocean temperature in the DP closed type is more sensitive!

The file names (and therefore experiment names) used contain sufficient information to reconstruct the co2 values used in those runs. In some situations, we may not have access the configuration files used, and file names is all we can go by. Spacegrids contains several regular expression based methods (of the project class) that allow us to do this. First, the following regular expression extracts the co2 value from the experiment (file) names:


In [14]:
regexp = '[\w_]+_(\d+)ppm[\w_]+'

The function sg.parse_fname_func takes this regular expression, and the intended parameter name, to use it to yield another function that extracts the parameter from the experiment (file) name using the regexp. This is similar to the file reading function discussed above. The resulting function is passed an exper object, and outputs a list of (name, value) pairs. Such functions can be passed to P.insert to be processed, where for each experiment in the project, the function is executed, and the resulting (name, value) pairs are inserted into the vars or params dictionaries of that experiment:


In [15]:
P.insert(sg.parse_fname_func(regexp,parname='co2'))

This yields the following (unordered) list of co2 values, this time extracted from the experiment names (their filenames), without reading the configuration files:


In [16]:
[E.params['co2'] for E in P.expers.values() if 'DP' in E.name ]


Out[16]:
[1750.0, 1250.0, 2000.0, 750.0, 1000.0, 500.0, 390.0, 1500.0, None]

The 280ppm value is absent (in fact, it corresponds to the None entry). This is because this value is not encoded in the 280 ppm DP experiment name DPclsd_C_flat. This can be fixed by specifying the nomatch_fill argument:


In [17]:
P.insert(sg.parse_fname_func(regexp,parname='co2', nomatch_fill = '280ppm'))

In [18]:
[E.params['co2'] for E in P.expers.values() if 'DP' in E.name ]


Out[18]:
[1750.0, 1250.0, 2000.0, 750.0, 1000.0, 500.0, 390.0, 1500.0, '280ppm']

These co2 parameters can again be used to construct big_field_DP:


In [19]:
big_field_DP = P.param2gr('co2' , lambda x:x['temp'], name_filter = 'DP*')

We have created the same big field again as before, but now we have constructed it by using the information contained in the file naming:


In [20]:
h1, = sg.plot(bigfield_DP[Z,17:]/(X*Y*Z) + sg.zero_kelvin)
lbl = plt.ylabel('C')
lbl = plt.xlabel('co2 ccn')


A shorthand for this regexp-based procedure is contained in the pattern2gr method:


In [21]:
bigfield_DP=P.pattern2gr(fld_name = 'temp',pattern = regexp,parname='co2', name_filter = 'DP*', nomatch_fill = 280.)

In [22]:
bigfield_DP.gr[0][:]


Out[22]:
array([  280.,   390.,   500.,   750.,  1000.,  1250.,  1500.,  1750.,
        2000.])

In [ ]: