First, import the pyhemco.emissions Python module. The commands below allow to import the package without installing it.
In [1]:
import sys
sys.path.extend('/home/bovy/MyGitRepos/pyHEMCO') # replace with your path/to/pyHEMCO_git_repo
from pyhemco import emissions
In [2]:
input_content=!cat HEMCO_input_sample.txt
from IPython.display import HTML
HTML('<pre style="font-size:80%;">' + '<br>'.join(input_content) + '</pre>')
Out[2]:
Steps required for creating a new emission setup using pyhemco are detailed below.
GCField object. In this example, we first create new GCField objects for fields that will be included in the emission setup (a GEOS-Chem Python API currently in development will provide high-level functions for reading NetCDF files and importing the NetCDF variables as GCField objects).
In [3]:
noxscalar = emissions.GCField('NOXScalar',
var_name='NOXScalar',
filename='/path/to/data/emis/AnnualScalar/nc/NOx-AnnualScalar.geos.1x1.nc',
ndim=2,
unit='unitless')
geia_dow_nox = emissions.GCField('GEIA_DOW_NOX',
ndim=1, # does all scalar-based scale factors have 'xy' dims in HEMCO input file ?
unit='unitless',
data=[1.0706, 1.0706, 1.0706, 1.0706, 1.0706, 0.863, 0.784])
emep_mask = emissions.GCField('EMEP_MASK',
var_name='MASK',
filename='/path/to/data/emis/MASKS/nc/EMEP_mask.geos.1x1.nc',
ndim=2,
unit='xxx')
geia_no = emissions.GCField('GIEA_NO',
var_name='NO',
filename='/path/to/emis/GEIA/nc/GEIA_NOSO4.fullyear.geos.4x5.nc',
ndim=2,
unit='kg/m2/s')
emep_ship_no = emissions.GCField('EMEP_SHIP_NO',
var_name='NO',
filename='/path/to/data/emis/EMEP/nc/EMEP.updt.ship.geos.1x1.nc',
ndim=2,
unit='kg/m2/s')
Datafields that will be used by the HEMCO core and its extensions are created using the base_emission_field, scale_factor and mask functions, which basically add the metadata required by HEMCO to GCField objects.
Scale factors and masks IDs (fid) can be left undefined, although it must be set when saving the emission setup to a file. Note that mask fields are specific scale factor fields and therefore have some common metadata attributes.
Scale factors and masks can be assigned to a base emission field with the scale_factor argument of the base_emission_field function.
Emission metadata can be added to a GCField in place (copy=False, default) or to a copy of this GCField (copy=True). The latter option may be useful, for example, when we need to create multiple scale factors or masks from one GCField.
In [4]:
# copy GCField and add metadata
totfuel_thisyr = emissions.scale_factor(noxscalar, 'TOTFUEL_THISYR', '1985-2008/1/1/0', operator='*',
fid=1, copy=True)
totfuel_1985 = emissions.scale_factor(noxscalar, 'TOTFUEL_1985', '1985/1/1/0', operator='/',
fid=2, copy=True)
totfuel_2002 = emissions.scale_factor(noxscalar, 'TOTFUEL_2002', '2002/1/1/0', operator='/',
fid=3, copy=True)
# add metadata in place
emissions.scale_factor(geia_dow_nox, geia_dow_nox.name, '*/*/*/*', operator='*', fid=None) # fid will be set later
emissions.mask(emep_mask, emep_mask.name, '*/*/*/*', mask_window=[-30, 30, 45, 70], fid=1000)
emissions.base_emission_field(geia_no, geia_no.name, '1985/1-12/1/0', 'NO', 1, 1,
scale_factors=[totfuel_thisyr, totfuel_1985, geia_dow_nox])
emissions.base_emission_field(emep_ship_no, emep_ship_no.name, '1990-2007/1/1/0', 'NO', 1, 10,
scale_factors=[emep_mask])
EmissionExt object for each HEMCO extension to include in the setup (must be done at least for HEMCO Core!). Currently, the EmissionExt class only stores basic information about the extension (extension name, on/off switch, extension number). Base emissions fields can eventually be attached to an extension. Similar to scale factors, User can optionally give a value for eid (i.e., the extension number) or leave it undefined (None, default), though it is higly recommended to define a value. When settings are written to a file, a eid value must be defined for each extension (it can be done automatically, see below). Note that HEMCO Core must have eid set to 0 !
In [5]:
core_ext = emissions.EmissionExt('HEMCOCore', eid=0, base_emission_fields=[geia_no])
seaflux_ext = emissions.EmissionExt('SeaFlux', enabled=False, eid=101)
paranox_ext = emissions.EmissionExt('PARANOX', enabled=False, eid=102, base_emission_fields=[emep_ship_no])
lightnox_ext = emissions.EmissionExt('LightNOx', enabled=True, eid=103)
soilnox_ext = emissions.EmissionExt('SoilNOx', enabled=True, eid=None) # eid will be set later.
extra_settings keyword argument. It is also possible to define an empty setup and add extensions later (see section 'Editing an existing emission setup' below).
In [6]:
example_setup = emissions.Emissions([core_ext, seaflux_ext, paranox_ext, lightnox_ext, soilnox_ext],
description='An emission setup example (not complete)',
extra_settings={'Verbose': False})
empty_setup = emissions.Emissions([], description='An empty setup')
It is easy to access to any emission setting, either directly or through filtering. Iterating over settings is also straightforward. A few basic examples are given below.
Get all base emission fields (returns a collection of GCField objects, see class pyhemco.datatypes.ObjectCollection):
In [7]:
example_setup.base_emission_fields
Out[7]:
Get all scale factors and masks that are assigned to base emission fields (read-only collection of GCField objects):
In [8]:
example_setup.scale_factors
Out[8]:
Get all HEMCO extensions involved in the setup (read-only collection of EmissionExt objects):
In [9]:
example_setup.extensions
Out[9]:
Get the 'GIEA_NO' base emission field
In [10]:
giea_no = example_setup.base_emission_fields.get_object(name='GIEA_NO')
giea_no
Out[10]:
Base emission field metadata is stored as a dictionary in the 'emission_base' attribute of GCField:
In [11]:
giea_no.attributes['emission_base']
Out[11]:
Similarly, for scale factors and masks ('emission_scale_factor' attribute):
In [12]:
totfuel_thisyr = giea_no.emission_scale_factors.get_object(name='TOTFUEL_THISYR')
totfuel_thisyr.attributes['emission_scale_factor']
Out[12]:
Get all base emission fields that will be used by the 'PARANOX' extension:
In [13]:
paranox_ext = example_setup.extensions.get_object(name='PARANOX')
paranox_fields = paranox_ext.base_emission_fields
paranox_fields
Out[13]:
Get all masks:
In [14]:
filter_masks = lambda field: 'mask_window' in field.attributes['emission_scale_factor'].keys()
masks = example_setup.scale_factors.filter(filter_masks)
masks
Out[14]:
Which extension will use the 'GIEA_NO' base emission field?
In [15]:
example_setup.extensions.get_object(lambda ext: giea_no in ext.base_emission_fields)
Out[15]:
There is more filtering / sorting possibilities, see docstring of class pyhemco.datatypes.ObjectCollection.
In [16]:
acet_seawater = emissions.GCField('ACET_SEAWATER',
var_name='ACET',
filename='/path/to/data/emis/ACET/nc/ACET_seawater.generic.1x1.nc',
ndim=2,
unit='ngL')
emissions.base_emission_field(acet_seawater, acet_seawater.name, '2005/1/1/0', 'ACET', 1, 1,
extension=seaflux_ext)
seaflux_ext = example_setup.extensions.get_object(name='SeaFlux')
seaflux_ext.base_emission_fields.add(acet_seawater)
seaflux_ext.base_emission_fields
Out[16]:
In [17]:
# base emission fields are automatically updated:
example_setup.base_emission_fields
Out[17]:
Create and assign a new scale factor to the 'GEIA_NO' base emission field:
In [18]:
edgar_todnox = emissions.GCField('EDGAR_TODNOX',
var_name='NOXscale',
filename='/path/to/data/emis/EDGAR/nc/dailyScal.nc',
ndim=2,
unit='unitless')
emissions.scale_factor(edgar_todnox, edgar_todnox.name, '2000/1/1/1-24', operator='*', fid=25)
giea_no = example_setup.base_emission_fields.get_object(name='GIEA_NO')
giea_no.emission_scale_factors.add(edgar_todnox)
giea_no.emission_scale_factors
Out[18]:
In [19]:
# scale factors listed in setup are automatically updated:
example_setup.scale_factors
Out[19]:
Remove the mask assigned to the 'EMEP_SHIP_NO' base emission field:
In [20]:
emep_ship_no = example_setup.base_emission_fields.get_object(name='EMEP_SHIP_NO')
emep_ship_no.emission_scale_factors.get(name='EMEP_MASK').remove()
emep_ship_no.emission_scale_factors
Out[20]:
In [21]:
example_setup.scale_factors
Out[21]:
Remove all base emission fields attached to the 'PARANOX' HEMCO extension:
In [22]:
paranox_ext = example_setup.extensions.get_object(name='PARANOX')
paranox_ext.base_emission_fields.get_all().remove()
example_setup.base_emission_fields
Out[22]:
Change the ID (fid) of the 'EDGAR_TODNOX' scale factor:
In [23]:
edgar_todnox = example_setup.scale_factors.get_object(name='EDGAR_TODNOX')
edgar_todnox.attributes['emission_scale_factor']['fid'] = 26
edgar_todnox.attributes
Out[23]:
In [24]:
giea_no.emission_scale_factors.get_object(name='EDGAR_TODNOX').attributes
Out[24]:
In [25]:
# should raise an error because we left some IDs undefined above
example_setup.check_id()
resolve_id allows to set automatically new values for missing or duplicate IDs
In [26]:
geia_dow_nox = example_setup.scale_factors.get_object(name='GEIA_DOW_NOX')
geia_dow_nox.attributes['emission_scale_factor']['fid'] is None # there is a missing ID here (should return True)
Out[26]:
In [27]:
example_setup.resolve_id()
example_setup.check_id() # should not raise any error anymore
geia_dow_nox.attributes['emission_scale_factor']['fid'] # should return a valid ID
Out[27]:
In [ ]:
In [28]:
%pylab inline
pylab.rcParams['figure.figsize'] = (10.0, 8.0)
import networkx as nx
from networkx.readwrite import json_graph
import matplotlib.pyplot as plt
import json
Building the graph from emission setup
In [29]:
G = nx.DiGraph()
G.add_node("SETUP")
for ext in example_setup.extensions.sorted(lambda ext: ext.eid):
ext_node_label = "{} ({})".format(ext.name, ext.eid)
G.add_node(ext_node_label)
G.add_edge("SETUP", ext_node_label)
for bef in ext.base_emission_fields:
G.add_node(bef.name)
G.add_edge(ext_node_label, bef.name)
for sf in bef.emission_scale_factors:
# TODO: handle the case where a scale factor is assigned to
# several base emission fields (only one node?)
sf_fid = sf.attributes['emission_scale_factor']['fid']
sf_node_label = "{} ({})".format(sf.name, sf_fid)
G.add_node(sf_node_label)
G.add_edge(bef.name, sf_node_label)
Visualization using matplotlib
Requires:
In [30]:
plt.title(example_setup.description)
glayout = nx.pygraphviz_layout(G, prog="dot")
nx.draw(G, glayout, with_labels=True, arrows=False) # default parameters, output not very pretty
plt.show()
Visualization using d3.js
In [32]:
# serialize graph into d3 tree format and write to a json file
setup_dict = json_graph.tree_data(G, "SETUP")
json.dump(setup_dict, open('setup_example.json', 'w'))
In [33]:
%%html
<div id="tree-container"></div>
<style>
.node circle {
fill: #fff;
stroke: steelblue;
stroke-width: 1.5px;
}
.node {
font: 10px sans-serif;
}
.link {
fill: none;
stroke: #ccc;
stroke-width: 1.5px;
}
</style>
In [34]:
%%javascript
require.config({paths: {d3: "http://d3js.org/d3.v3.min",
jquery: "http://code.jquery.com/jquery-1.10.2.min.js"}
});
require(["d3", "jquery"], function(d3) {
var width = 500,
height = 400;
var tree = d3.layout.tree()
.size([height, width - 160]);
var diagonal = d3.svg.diagonal()
.projection(function(d) { return [d.y, d.x]; });
var svg = d3.select("#tree-container").append("svg")
.attr("width", width)
.attr("height", height)
.append("g")
.attr("transform", "translate(40,0)");
d3.json("files/setup_example.json", function(error, root) {
var nodes = tree.nodes(root),
links = tree.links(nodes);
var link = svg.selectAll(".link")
.data(links)
.enter().append("path")
.attr("class", "link")
.attr("d", diagonal);
var node = svg.selectAll(".node")
.data(nodes)
.enter().append("g")
.attr("class", "node")
.attr("transform", function(d) { return "translate(" + d.y + "," + d.x + ")"; })
node.append("circle")
.attr("r", 4.5);
node.append("text")
.attr("dx", function(d) { return d.children ? -8 : 8; })
.attr("dy", 3)
.style("text-anchor", function(d) { return d.children ? "end" : "start"; })
.text(function(d) { return d.id; });
});
d3.select(self.frameElement).style("height", height + "px");
});
Visualization with d3.js (advanced)
Pan, Zoom, Collapse
In [35]:
%%html
<div id="tree-container2"></div>
<style>
.node rect {
cursor: pointer;
fill: #fff;
fill-opacity: .5;
stroke: #3182bd;
stroke-width: 1.5px;
}
.node text {
font: 10px sans-serif;
pointer-events: none;
}
path.link {
fill: none;
stroke: #9ecae1;
stroke-width: 1.5px;
}
</style>
In [36]:
%%javascript
require.config({paths: {d3: "http://d3js.org/d3.v3.min",
jquery: "http://code.jquery.com/jquery-1.10.2.min.js"}
});
require(["d3", "jquery"], function(d3) {
var margin = {top: 30, right: 20, bottom: 30, left: 20},
width = 500 - margin.left - margin.right,
barHeight = 20,
barWidth = width * .8;
var i = 0,
duration = 400,
root;
var tree = d3.layout.tree()
.nodeSize([0, 20]);
var diagonal = d3.svg.diagonal()
.projection(function(d) { return [d.y, d.x]; });
var svg = d3.select("#tree-container2").append("svg")
.attr("width", width + margin.left + margin.right)
.append("g")
.attr("transform", "translate(" + margin.left + "," + margin.top + ")");
d3.json("files/setup_example.json", function(error, flare) {
flare.x0 = 0;
flare.y0 = 0;
update(root = flare);
});
function update(source) {
// Compute the flattened node list. TODO use d3.layout.hierarchy.
var nodes = tree.nodes(root);
var height = Math.max(500, nodes.length * barHeight + margin.top + margin.bottom);
d3.select("svg").transition()
.duration(duration)
.attr("height", height);
d3.select(self.frameElement).transition()
.duration(duration)
.style("height", height + "px");
// Compute the "layout".
nodes.forEach(function(n, i) {
n.x = i * barHeight;
});
// Update the nodes…
var node = svg.selectAll("g.node")
.data(nodes, function(d) { return d.id || (d.id = ++i); });
var nodeEnter = node.enter().append("g")
.attr("class", "node")
.attr("transform", function(d) { return "translate(" + source.y0 + "," + source.x0 + ")"; })
.style("opacity", 1e-6);
// Enter any new nodes at the parent's previous position.
nodeEnter.append("rect")
.attr("y", -barHeight / 2)
.attr("height", barHeight)
.attr("width", barWidth)
.style("fill", color)
.on("click", click);
nodeEnter.append("text")
.attr("dy", 3.5)
.attr("dx", 5.5)
.text(function(d) { return d.id; });
// Transition nodes to their new position.
nodeEnter.transition()
.duration(duration)
.attr("transform", function(d) { return "translate(" + d.y + "," + d.x + ")"; })
.style("opacity", 1);
node.transition()
.duration(duration)
.attr("transform", function(d) { return "translate(" + d.y + "," + d.x + ")"; })
.style("opacity", 1)
.select("rect")
.style("fill", color);
// Transition exiting nodes to the parent's new position.
node.exit().transition()
.duration(duration)
.attr("transform", function(d) { return "translate(" + source.y + "," + source.x + ")"; })
.style("opacity", 1e-6)
.remove();
// Update the links…
var link = svg.selectAll("path.link")
.data(tree.links(nodes), function(d) { return d.target.id; });
// Enter any new links at the parent's previous position.
link.enter().insert("path", "g")
.attr("class", "link")
.attr("d", function(d) {
var o = {x: source.x0, y: source.y0};
return diagonal({source: o, target: o});
})
.transition()
.duration(duration)
.attr("d", diagonal);
// Transition links to their new position.
link.transition()
.duration(duration)
.attr("d", diagonal);
// Transition exiting nodes to the parent's new position.
link.exit().transition()
.duration(duration)
.attr("d", function(d) {
var o = {x: source.x, y: source.y};
return diagonal({source: o, target: o});
})
.remove();
// Stash the old positions for transition.
nodes.forEach(function(d) {
d.x0 = d.x;
d.y0 = d.y;
});
}
// Toggle children on click.
function click(d) {
if (d.children) {
d._children = d.children;
d.children = null;
} else {
d.children = d._children;
d._children = null;
}
update(d);
}
function color(d) {
return d._children ? "#3182bd" : d.children ? "#c6dbef" : "#fd8d3c";
}
});
In [ ]: