Constructing Field Plots with Python and Diana

This document aims to show how to use Diana to create plots from Python scripts and from within the ipython notebook.

Installation

To get access to the Python modules that expose Diana's plotting features, you need to install the python-diana package from the local package repository. Either use the package manager to do this, or enter the following in a terminal:

apt-get install python-diana

Some examples can also be found in the python-diana-examples package.

Importing Modules

You can use Diana's functionality from Python by importing the bdiana and plotcommands modules from the metno package. If you are using an old version of ipython (earlier than version 1.0) then it is also useful to import the embed function from the metno.ipython_extensions module.



In [1]:

    
# Import the modules needed to use Diana.
from metno import bdiana, plotcommands
# Import an extension that lets us embed an image into a worksheet.
from metno.ipython_extensions import embed

Initialisation and Setup

The bdiana module contains the BDiana class. This provides a convenient interface to Diana's functions and also wraps several initialisation steps into one object. To start using Diana, we simply create an instance of this class.



In [2]:

    
# Create a BDiana object - this provides a simpler interface than calling Diana functions directly.
b = bdiana.BDiana()

By default, Diana doesn't know anything about data files or other resources. The locations of data files, definitions of plot types, palettes, and other pieces of information are stored in one or more setup files. We choose one of these and use the setup method to let Diana know where to get data from.



In [3]:

    
# Load a setup file - this describes where Diana finds data files and other resources.
b.setup("/etc/diana/setup/diana.setup-COMMON")

The method will return True if the setup file was loaded correctly.

Getting Model and Field Information

Now we have a working interface to Diana, we can ask for information about the models available to us by using the getModels method.



In [4]:

    
models = b.getModels()
# Uncomment the following lines to see the list of models known to Diana.
models = list(models)
models.sort()
models

We will use the WAM 4km model.



In [5]:

    
model = "WAM.4km"

To find the fields available for this model, we call the getFields method with the model.



In [6]:

    
fields = b.getFields(model)
# Uncomment the following lines to see the list of fields available.
fields = list(fields)
fields.sort()
fields

We will use the mean wave direction field.



In [7]:

    
field = "mean_wave_direction"

Before we can plot anything, we need to know which times we have field data for. To do this, we call the getFieldTimes method with the model and field as arguments.



In [8]:

    
times = b.getFieldTimes(model, field)
# Show the last five times.
times[-5:]

Having obtained the model, field and time information, we can now create a plot. To do this, we need to create some plot commands that describe how the plot will appear.

Creating a Plot

Diana creates plots by reading and executing a series of plot commands. To create a plot, we need to create these commands. First, we write a plot command to create a map by creating a Map object and setting some options. These are explained below.



In [9]:

    
m = plotcommands.Map()
m.setOption("map", "Gshhs-Auto")
m.setOption("backcolour", "white")
m.setOption("land", "on")
m.setOption("land.colour", "flesh")

The map option describes the map data to use. "Gshhs-Auto" is a good general map to use.
The backcolour option describes the colour to use for the plot background. Named colours or colour components can be used to specify the colour. Use "0:0:0:0" for a transparent background.
The land option determines whether the land is filled.
The land.colour describes the fill colour used.

We write a plot command for the field by creating a Field object and setting the model, plot and plottype options. In the example below, we pass these to the class initialisation method instead of setting them separately. Either way can be used.



In [10]:

    
f = plotcommands.Field(model = model, plot = field, plottype = "direction")
#f.setOption("model", model)
#f.setOption("plot", field)
#f.setOption("plottype", "direction")

Normally, plot commands are contained in an input file that describes what kind of output to produce, as well as other syntax to ensure that Diana understands exactly what is to be plotted. For convenience, the BDiana class provides the setPlotCommands method that lets us ignore those details. We pass the plot commands to this method to set up the plot.



In [11]:

    
b.setPlotCommands([m, f])
# Print the actual plot commands for clarity.
print m.text()
print f.text()

The last thing we need to do before asking Diana to plot the figure is to set the plot time. We choose the last one available from the collection we obtained earlier.



In [12]:

    
b.setPlotTime(times[-1])

A return value of True indicates that the specified time was acceptable. If False was returned then the time is not available. This shouldn't happen in this example.

Now we are able to create an image for the plot we specified. This one will be 500 pixels wide and 500 pixels high.



In [13]:

    
im = b.plotImage(500, 500)

When using an ipython notebook, it is useful to be able to preview the image we obtained. The embed function does this for us.



In [14]:

    
embed(im)

Plotting a Specific Area

By default, the area plotted will be the area covered by the fields used in the plot. To plot a specific area, we need to specify an area in the list of plot commands. This is done using the Area plot command class.



In [15]:

    
a = plotcommands.Area(name = "Norge")

Now we can create the plot again using this area.



In [16]:

    
b.setPlotCommands([a, m, f])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im)

Plotting Multiple Fields

To plot more than one field, simply create plot commands for each field and combine them in the list of plot commands. You will first need to make sure that there is data available for the time you wish to create a plot for.



In [17]:

    
field2 = "maximum_wave_height"
times2 = b.getFieldTimes(model, field2)
if times[-1] in times2:
    print times[-1], "is available in both fields."
else:
    print times[-1], "is not available for", field2

If the time is available in both fields then we can create a plot command and include it in the list we pass to the setPlotCommands method. In this case, we will include it before the command to show the mean wave direction.



In [18]:

    
f2 = plotcommands.Field(model = model, plot = field2, plottype = "contour")
f2.setOption("palettecolours", "standard")
f2.setOption("line.interval", 1)
b.setPlotCommands([a, m, f2, f])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im)

Fine-Tuning Plot Options

Some plot types look better with certain plot options enabled. The direction plot type, for example, looks a lot better in image plots when anti-aliasing is enabled.

Here's a close-up of part of the previous plot.



In [19]:

    
embed(im.copy(200, 300, 50, 50).scaled(200, 200))

If we set the antialias plot option in the first field, the arrows appear smoother but we don't affect the way the rest of the plot is displayed.



In [20]:

    
f.setOption("antialiasing", True)
b.setPlotCommands([a, m, f2, f])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im.copy(200, 300, 50, 50).scaled(200, 200))

Note that the frame around the field area is also anti-aliased. The frame itself can be disabled in plots. Note that there are two frames being plotted here, so we need to disable them both.



In [21]:

    
f.setOption("frame", 0)
f2.setOption("frame", 0)
b.setPlotCommands([a, m, f2, f])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im.copy(200, 300, 50, 50).scaled(200, 200))

Creating PDF Files

While images are fine to display on the screen and can be easily embedded in Web pages, for publications it is better to create files in PDF format. The BDiana class provides a method for this.



In [22]:

    
b.plotPDF(500, 500, "/tmp/plot.pdf")

You can discard the return value from this call if you want. The output file can be found in the location passed to the plotPDF method call.

Adding Labels

It's useful to include some information about the data we are plotting. To do this, we add Label plot commands to the list we pass to the BDiana object. For example, we can create a label that shows the name of the model and field.



In [23]:

    
l1 = plotcommands.Label()
l1.setOption("data")
b.setPlotCommands([a, m, f2, f, l1])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im)

The Label plot command accept a number of different options. We can see some of these by calling its getOptions method.



In [24]:

    
l1.getOptions()

Let's try moving the label to the top-left of the plot and changing the font size.



In [25]:

    
l1.setOption("valign", "top")
l1.setOption("fontsize", 8)
b.setPlotCommands([a, m, f2, f, l1])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im)

We can also show the associated table for the maximum wave height plot, but we need to ensure that this plot generates the information needed for this by setting its table option to 1 and setting the line.values option.



In [26]:

    
f2 = plotcommands.Field(model = model, plot = field2, plottype = "fill_cell")
f2.setOption("frame", 0)
f2.setOption("antialiasing", True)
f2.setOption("palettecolours", "standard")
f2.setOption("table", 1)
f2.setOption("line.values", [0,1,2,3,4,5,6,7,8,9,10])
f2.text()

The command to show the table is a label that uses the anno option with a value which contains some special syntax. We place the label in the bottom-left corner of the plot. It is necessary to set the polystyle option to either none or border.



In [27]:

    
l2 = plotcommands.Label()
l2.setOption("anno", "<table,fcolour=white>")
l2.setOption("halign", "left")
l2.setOption("valign", "bottom")
l2.setOption("polystyle", "none")
l2.setOption("fontsize", 8)

b.setPlotCommands([a, m, f2, f, l1, l2])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im)

Using Projections

In addition to specifying an area to plot, we can also specify the projection used to display the data. To do this, we need to import the pyproj module.



In [28]:

    
import pyproj

A projection object is created by defining which projection and datum to use. We use the projection object to obtain a rectangle in the projected coordinate system whose bottom-left and top-right corners are specified as latitude and longitude values.



In [29]:

    
p = pyproj.Proj(proj="utm", datum="WGS84")
x1, y1 = p(-15, 25)
x2, y2 = p(35, 75)

With the projection and rectangle defined, we can create a plot command



In [30]:

    
a2 = plotcommands.Area(proj4string = p.srs, rectangle = (x1, x2, y1, y2))

b.setPlotCommands([a2, m, f2, f, l1, l2])
b.setPlotTime(times[-1])
im = b.plotImage(500, 500)
embed(im)

Adding Fields

Usually, all the fields you want to use are listed in the setup file that is passed to the BDiana object's setup method. However, you may want to load another field that is not listed there. To achieve this, the BDiana object provides an addField method that performs the work of loading the field into Diana.

For example, we can obtain a data file from a THREDDS server and examine its contents. We must provide a model name so that we can reference the data.



In [31]:

    
url = "http://thredds.met.no/thredds/dodsC/myocean/siw-tac/siw-fmi-baltic/2013/03/ice_thic_baltic_201303311400.nc"
b.addModel("baltic_sea_ice_thickness", url)

The new model can now be found in the list of available models.



In [32]:

    
models = list(b.getModels())
"baltic_sea_ice_thickness" in models

The fields provided by the model can be found as usual.



In [33]:

    
fields = list(b.getFields("baltic_sea_ice_thickness"))
fields.sort()
fields

We plot one of these fields in the usual way.



In [34]:

    
a = plotcommands.Area(name = "Europa")
m = plotcommands.Map(map = "Gshhs-Auto", backcolour = "white", land = "on")
m.setOption("land.colour", "flesh")
f = plotcommands.Field(model = "baltic_sea_ice_thickness",
                       plot = "ice_thickness", plottype = "fill_cell")
times = b.getFieldTimes("baltic_sea_ice_thickness", "ice_thickness")

b.setPlotCommands([a, m, f])
b.setPlotTime(times[-1])
embed(b.plotImage(500, 500))

Sometimes, we only want to use the data from a particular file, and so we want to remove all other information about other models when we load it. To do this, we call the addField with the clear argument set to True.



In [35]:

    
b.addModel("baltic_sea_ice_thickness", url, clear = True)

Now, the list of available models only contains the model we loaded.



In [36]:

    
b.getModels()