<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Handling CTM datasets
This tutorial shows how to load and save datasets, and also how it can be manipulated.
The pygchem.datasets module provides functions to load and save CTM datasets
Step1:  NOTE
Step2:  Loading datasets
The function pygchem.datasets.load allows to load the content of a BPCH file or a netCDF file into a list of iris.cube.Cube objects. This function is just an alias of the iris.load function. Note that it doesn't load all the data into memory but enough information to select, manipulate and retreive the data.
See the Iris documentation
Step3:  Simple (unconstrained) loading
The line below loads the BPCH file corresponding to the first month of the simulation
Step4:  The line below print the list of the 20 lasts data fields of the list (name, units, dimensions and coordinates). Note that the field name may be reformatted during the loading so that it is consistent with the name set by the GAMAP's netCDF writing routines. Note also that the units may also be reformatted as required by the udunits2 library (used by Iris).
Step5:  The same function can be used for loading the content of the netCDF files. The netCDF files available in this simulation were created by the 'BPCH2COARDS' GAMAP routine. As there are differences between the COARDS Conventions and the CF Conventions, we have to use a callback function 'gamap_bpch2coards' to properly load the content of the netCDF file into Iris cubes
Step6:  PyGChem also provides a callback for loading netCDF files created by the GAMAP routine 'BPCH2NC'
Step7:  Constrained loading
For large datasets, it is possible to constrain the load to match specific metadata or to load only a subset of the data. It is recommended as the loading is faster with given constraints.
The following example loads only the "IJ_AVG_S__O3" variable
Step8:  We can specify multiple variables
Step9:  It is also possible to define more advanced constraints. For example, to load all "IJ-AVG-$" diagnostics
Step11:  A more advanced example, combining constraints and extracting data subsets
Step12:  Loading multiple files
It is possible to load multiple files at once. The pygchem.datafields.load function will try to merge the fields when it is possible, i.e., when the multiple fields describing the same phenomenom overlap in space and/or time, so that it reduces at the minimum the number of loaded fields.
The following example load the "IJ_AVG_S__O3" diagnostic for the entire 1-year simulation
Step13:  Merging fields may take a long time. If speed matters, it is still possible to load the fields without any merging
Step14:  Apply contraints to the field list after loading
We can extract data subsets and/or specific fields from a loaded list of fields
Step15:  To select only one field (cube), the extract_strict method can be used
Step16:  Saving datasets
Datasets can be written to a netCDF file (with the CF-conventions) using the function save of the datafields module (an alias to the iris.save function).
Step17:  A text representation of (the header information of) the written netCDF file using the ncdump utility (provided with the netCDF4 package)
Step18:  Loading the written file using the load function
Step19:  It is also possible to write the datasets to the BPCH format, using the low-level function write_bpch in the module pygchem.io.bpch (not yet documented).
Accessing data and metadata
String representation
Print a string representation of the data field (cube)
Step20:  Field name(s)
Get or set the name(s) of the field with the following function and/or attributes
Step21:  PyGChem considers that the GEOS-Chem variable name (category + tracer) is a standard name, although it is not CF-compliant (i.e., not listed in the standard name table of the udunits package).
Step22:  long_name is the full name of the diagnostic
Step23:  var_name is the (netCDF) variable name
Step24:  Attributes
Any attribute of the data field is stored in the following dictionary
Step25:  Field units
Get the field units
Step26:  It is easy to change the units of the field (data values are re-computed accordingly)
Step27:  Coordinates
Retreiving a coordinate by name
Step28:  Coordinate data and metadata
Step29:  Data
The field data can be accessed with the data attribute. It returns a Numpy array.
Manipulating datasets
It is very easy to manipulate data fields (Iris cubes)
Indexing
Indexing a cube return a new cube.
Step30:  Note that another way to extract a subset is by applying one or more constraints on the cube (see above).
Iterating
The slice method allow iterating over layers or time slices (return a Python generator). The slices can be 1-dimensional or n-dimensional.
The example below generate time slices (returns a 3-d cube for each time slice).
Step31:  Collapsing data dimensions (statistics)
The example below calculate the sum over the vertical levels
Step32:  Basic cube mathematics
The example below calculates the total columns of the tracer for all grid cells and all time slices
Step33:  Plotting datasets
Iris provides some modules for basic dataset plotting. It is built on top of matplotlib, and it uses the cartopy package for map projections.
Step34:  Plot the NO2 total columns for the first time slice
Step35:  A Hovmoller diagram example
Step36:  The CTM diagnostics
For loading datasets, Pygchem reads the information stored in the diaginfo.dat and tracerinfo.dat files. The module pygchem.diagnostics provides a basic API for reading / writing these files and handling their contents.
Step37:  To load a couple of files
Step38:  A CTMDiagnosticInfo object contains all information stored in those files. The attributes categories and diagnostics contains each record (line) in diaginfo.dat and tracerinfo.dat, respectively
Step39:  These attributes behave like a Python list, with added key reference and database lookup-like capabilities. Each item of the list coorespond to a record.
Step40:  It is aslo possible to filter the data (queries)
Step41:  We can add or remove entries
Step42:  Exporting to diaginfo and tracerinfo files
<ASSISTANT_TASK:>
Python Code:
from pygchem import datasets
Explanation: Handling CTM datasets
This tutorial shows how to load and save datasets, and also how it can be manipulated.
The pygchem.datasets module provides functions to load and save CTM datasets:
End of explanation
bmk_root = '/home/bovy/geoschem'
Explanation: NOTE: In the PyGChem development version, this module was previously named datafields and has been recently renamed to datasets (4-12-2014).
if
from pygchem import datasets
doesn't work, try
import pygchem.datafields as datasets
to use this notebook.
The Iris backend
PyGChem uses Iris (http://scitools.org.uk/iris/) for handling CTM datasets. Iris provides a data abstraction layer that is based on the Climate and Forecast (CF) data model (http://cfconventions.org/) and which allows to isolate analysis and visualisation code from data format specifics. It natively supports read/write access to the (CF)netCDF format. PyGChem currently adds support to the BPCH format and the netCDF files created with the GAMAP routines BPCH2NC and BPCH2COARDS (read access only).
With Iris, a data field (data + metadata) relative to a particular phenomenon is represented by a iris.cube.Cube object (more info on Iris data strutures here). Iris also supports basic data manipulation operations, such as arithmetic, interpolation, and statistics. This tutorial shows only a few examples. For a more detailled description of these features, see the Iris user guide.
Datasets used in this tutorial
This tutorial uses the datasets created by the 1-year GEOS-Chem benchmark simulations. These are available for download via anonymous FTP from ftp://ftp.as.harvard.edu/gcgrid/geos-chem/.
For running the code cells below, you must specify the root folder of your local copy of these benchmarks:
End of explanation
%cd {bmk_root}/1yr_benchmarks/v10-01/v10-01c/Run1
Explanation: Loading datasets
The function pygchem.datasets.load allows to load the content of a BPCH file or a netCDF file into a list of iris.cube.Cube objects. This function is just an alias of the iris.load function. Note that it doesn't load all the data into memory but enough information to select, manipulate and retreive the data.
See the Iris documentation: loading iris cubes for more details.
cd into the 1-year benchmark (v10-01c Run1)
End of explanation
filename = 'bpch/ctm.bpch.v10-01c-geosfp-Run1.20120801'
dataset = datasets.load(filename)
Explanation: Simple (unconstrained) loading
The line below loads the BPCH file corresponding to the first month of the simulation:
End of explanation
print dataset[-20:]
Explanation: The line below print the list of the 20 lasts data fields of the list (name, units, dimensions and coordinates). Note that the field name may be reformatted during the loading so that it is consistent with the name set by the GAMAP's netCDF writing routines. Note also that the units may also be reformatted as required by the udunits2 library (used by Iris).
End of explanation
filename = 'netcdf/v10-01c-geosfp-Run1.20120801.nc'
clb = datasets.load_callbacks['gamap_bpch2coards']
dataset = datasets.load(filename, callback=clb)
print dataset[-20:]
Explanation: The same function can be used for loading the content of the netCDF files. The netCDF files available in this simulation were created by the 'BPCH2COARDS' GAMAP routine. As there are differences between the COARDS Conventions and the CF Conventions, we have to use a callback function 'gamap_bpch2coards' to properly load the content of the netCDF file into Iris cubes:
End of explanation
% cd {bmk_root}/1yr_benchmarks/v9-02/v9-02r/geos5/Run0
filename = 'netcdf/v9-02r-geos5-Run0.20050101.nc'
clb = datasets.load_callbacks['gamap_bpch2nc']
dataset = datasets.load(filename, callback=clb)
print dataset[-20:]
%cd {bmk_root}/1yr_benchmarks/v10-01/v10-01c/Run1
Explanation: PyGChem also provides a callback for loading netCDF files created by the GAMAP routine 'BPCH2NC':
End of explanation
filename = 'netcdf/v10-01c-geosfp-Run1.20120801.nc'
clb = datasets.load_callbacks['gamap_bpch2coards']
dataset = datasets.load(filename, "IJ_AVG_S__O3",
                        callback=clb)
print dataset
Explanation: Constrained loading
For large datasets, it is possible to constrain the load to match specific metadata or to load only a subset of the data. It is recommended as the loading is faster with given constraints.
The following example loads only the "IJ_AVG_S__O3" variable:
End of explanation
dataset = datasets.load(filename, ["IJ_AVG_S__O3", "IJ_AVG_S__NO"],
                        callback=clb)
print dataset
Explanation: We can specify multiple variables:
End of explanation
import iris
check_ij_avg = lambda cube: cube.name().startswith("IJ_AVG_S")
ij_avg = iris.Constraint(cube_func=check_ij_avg)
dataset = datasets.load(filename, ij_avg,
                        callback=clb)
print dataset
Explanation: It is also possible to define more advanced constraints. For example, to load all "IJ-AVG-$" diagnostics:
End of explanation
def lon_subset(cell):
   
   return True or False as to whether the cell
   center in question should be kept
   
   return cell > 0. and cell < 20.
lon_cst = iris.Constraint(longitude=lon_subset)
dataset = datasets.load(filename,
                        "IJ_AVG_S__O3" & lon_cst,
                        callback=clb)
print dataset
# note the reduced grid-size for the longitude
Explanation: A more advanced example, combining constraints and extracting data subsets:
End of explanation
# note the wildcard character in the filename
# (UNIX expressions are supported)
filename = 'bpch/ctm.bpch.v10-01c-geosfp-Run1.*'
diagnostics = ["BXHGHT_S__BXHEIGHT",
               "BXHGHT_S__N(AIR)",
               "IJ_AVG_S__NO2"]
dataset = datasets.load(filename, diagnostics)
print dataset
# note the additional time dimension
Explanation: Loading multiple files
It is possible to load multiple files at once. The pygchem.datafields.load function will try to merge the fields when it is possible, i.e., when the multiple fields describing the same phenomenom overlap in space and/or time, so that it reduces at the minimum the number of loaded fields.
The following example load the "IJ_AVG_S__O3" diagnostic for the entire 1-year simulation:
End of explanation
dataset_nomerge = datasets.load_raw(filename, diagnostics)
print dataset_nomerge
Explanation: Merging fields may take a long time. If speed matters, it is still possible to load the fields without any merging:
End of explanation
dataset_lon_subset = dataset.extract(lon_cst)
print dataset_lon_subset
Explanation: Apply contraints to the field list after loading
We can extract data subsets and/or specific fields from a loaded list of fields:
End of explanation
no2_avg = dataset.extract_strict("IJ_AVG_S__NO2")
print no2_avg
Explanation: To select only one field (cube), the extract_strict method can be used:
End of explanation
outfile = 'netcdf/test.nc'
datasets.save(dataset, outfile)
Explanation: Saving datasets
Datasets can be written to a netCDF file (with the CF-conventions) using the function save of the datafields module (an alias to the iris.save function).
End of explanation
!ncdump -h netcdf/test.nc
Explanation: A text representation of (the header information of) the written netCDF file using the ncdump utility (provided with the netCDF4 package):
End of explanation
print datasets.load('netcdf/test.nc')
Explanation: Loading the written file using the load function:
End of explanation
print no2_avg
Explanation: It is also possible to write the datasets to the BPCH format, using the low-level function write_bpch in the module pygchem.io.bpch (not yet documented).
Accessing data and metadata
String representation
Print a string representation of the data field (cube):
End of explanation
no2_avg.name()
Explanation: Field name(s)
Get or set the name(s) of the field with the following function and/or attributes:
End of explanation
no2_avg.standard_name
Explanation: PyGChem considers that the GEOS-Chem variable name (category + tracer) is a standard name, although it is not CF-compliant (i.e., not listed in the standard name table of the udunits package).
End of explanation
no2_avg.long_name
Explanation: long_name is the full name of the diagnostic
End of explanation
no2_avg.var_name
Explanation: var_name is the (netCDF) variable name
End of explanation
no2_avg.attributes
Explanation: Attributes
Any attribute of the data field is stored in the following dictionary
End of explanation
no2_avg.units
Explanation: Field units
Get the field units:
End of explanation
no2_avg.convert_units('ppmv')
print no2_avg
Explanation: It is easy to change the units of the field (data values are re-computed accordingly):
End of explanation
lat_coord =  no2_avg.coord('latitude')
Explanation: Coordinates
Retreiving a coordinate by name:
End of explanation
lat_coord.points
lat_coord.bounds
lat_coord.units
Explanation: Coordinate data and metadata:
End of explanation
print no2_avg
# Get the first element of the 1st and last dimensions (time and model level number)
no2_avg_t0_l1 = no2_avg[0, :, :, 0]
print no2_avg_t0_l1
Explanation: Data
The field data can be accessed with the data attribute. It returns a Numpy array.
Manipulating datasets
It is very easy to manipulate data fields (Iris cubes)
Indexing
Indexing a cube return a new cube.
End of explanation
no2_avg_time_slices = no2_avg.slices(['longitude', 'latitude', 'model_level_number'])
for s in no2_avg_time_slices:
    print s
Explanation: Note that another way to extract a subset is by applying one or more constraints on the cube (see above).
Iterating
The slice method allow iterating over layers or time slices (return a Python generator). The slices can be 1-dimensional or n-dimensional.
The example below generate time slices (returns a 3-d cube for each time slice).
End of explanation
import iris.analysis
no2_avg_sum_levels = no2_avg.collapsed('model_level_number', iris.analysis.SUM)
print no2_avg_sum_levels
Explanation: Collapsing data dimensions (statistics)
The example below calculate the sum over the vertical levels
End of explanation
# extract the data fields (cubes) needed to compute the tracer columns
box_heights = dataset.extract_strict("BXHGHT_S__BXHEIGHT")
n_air = dataset.extract_strict("BXHGHT_S__N(AIR)")
# convert units back to ppbv for the NO2 tracer
no2_avg.convert_units('ppbv')
# calculate the columns
no2_avg_columns = (box_heights * n_air * no2_avg).collapsed('model_level_number',
                                                            iris.analysis.SUM)
# set name convert units to count/cm2 (count is used for #molecules)
no2_avg_columns.rename("NO2 columns")
no2_avg_columns.convert_units('count/cm2')
# string repr
print no2_avg_columns
Explanation: Basic cube mathematics
The example below calculates the total columns of the tracer for all grid cells and all time slices:
End of explanation
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import iris.quickplot as qplt
%matplotlib inline
Explanation: Plotting datasets
Iris provides some modules for basic dataset plotting. It is built on top of matplotlib, and it uses the cartopy package for map projections.
End of explanation
fig = plt.figure(figsize=(10, 8))
qplt.contourf(no2_avg_columns[0], 15)
plt.gca().coastlines()
Explanation: Plot the NO2 total columns for the first time slice:
End of explanation
import matplotlib.dates as mdates
no2_hovmoller = no2_avg_columns.collapsed('latitude',
                                          iris.analysis.MEAN)
fig = plt.figure(figsize=(10, 8))
qplt.contourf(no2_hovmoller, 20)
# fine control over time axis ticks and labels
plt.gca().yaxis.set_major_locator(mdates.MonthLocator())
plt.gca().yaxis.set_major_formatter(mdates.DateFormatter('%m-%Y'))
plt.gca().set_ylabel("Time")
Explanation: A Hovmoller diagram example
End of explanation
from pygchem import diagnostics
Explanation: The CTM diagnostics
For loading datasets, Pygchem reads the information stored in the diaginfo.dat and tracerinfo.dat files. The module pygchem.diagnostics provides a basic API for reading / writing these files and handling their contents.
End of explanation
dinfo = diagnostics.CTMDiagnosticInfo(diaginfo_file='diaginfo.dat',
                                      tracerinfo_file='tracerinfo.dat')
Explanation: To load a couple of files:
End of explanation
dinfo.categories
Explanation: A CTMDiagnosticInfo object contains all information stored in those files. The attributes categories and diagnostics contains each record (line) in diaginfo.dat and tracerinfo.dat, respectively
End of explanation
# get the 1st category (a Record like object)
cat_ij_avg = dinfo.categories[0]
cat_ij_avg
cat_ij_avg.name
cat_ij_avg.offset
# convert the record object to a dict
cat_ij_avg.to_dict()
Explanation: These attributes behave like a Python list, with added key reference and database lookup-like capabilities. Each item of the list coorespond to a record.
End of explanation
# select a category based on its name (key)
dinfo.categories.select_item("NS-FLX-$")
# select a diagnostic (tracer) based on its number (key)
dinfo.diagnostics.select_item(11)
# select categories based on other attributes
dinfo.categories.select(offset=3000)
# advanced selection
dinfo.diagnostics.select(lambda d: d.unit == 'ppbC' and d.number > 10)
Explanation: It is aslo possible to filter the data (queries):
End of explanation
new_tracer = diagnostics.CTMDiagnostic(9999, 'NEW', full_name='a new tracer')
dinfo.diagnostics.append(new_tracer)
dinfo.diagnostics[-1]
# select the new tracer added to the list 
s = dinfo.diagnostics.select(9999)
# remove the selected entry
s.selection_remove()
dinfo.diagnostics[-1]
Explanation: We can add or remove entries:
End of explanation
dinfo.save_diaginfo('diaginfo_test.dat')
dinfo.save_tracerinfo('tracerinfo_test.dat')
Explanation: Exporting to diaginfo and tracerinfo files:
End of explanation
<END_TASK>