This document aims to show how to work with projections in Diana from Python scripts and from within the ipython notebook.
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.
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. The pyproj module is used to describe projections in a
way that Diana understands. We also import the datetime module so that we can pass a time
to Diana without needing to obtain one from a data file.
In [1]:
from metno import bdiana, metlibs, plotcommands
from metno.ipython_extensions import embed
import datetime, pyproj
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
and load a setup file that describes where various resources are located.
In [2]:
b = bdiana.BDiana()
b.setup("/etc/diana/diana.setup-COMMON")
To start with, we create a simple geographic projection.
In [3]:
p = pyproj.Proj(proj="latlong", datum="WGS84")
x1, y1 = p(-180, -90)
x2, y2 = p(180, 90)
We create a plot command to describe the map.
In [4]:
m = plotcommands.Map(map="Gshhs-Auto", backcolour="white", land="on")
m.setOption("land.colour", "flesh")
m.setOption("lat", "on")
m.setOption("lon", "on")
We define an area using the above projection by setting the proj4string option and
defining a view with the rectangle option.
In [5]:
a = plotcommands.Area(proj4string = p.srs, rectangle = [x1, x2, y1, y2])
Then we set the plot commands, use the current time and plot the image.
In [6]:
b.setPlotCommands([a.text(), m.text()])
b.setPlotTime(datetime.datetime.now())
embed(b.plotImage(500, 250))
Let's try another projection. The list of projections supported by the pyproj module are
defined by the underlying PROJ.4 library. A list of projections can be found at
http://www.remotesensing.org/geotiff/proj_list/. We will try the orthographic projection.
In [7]:
p = pyproj.Proj(proj="ortho", datum="WGS84")
x1, y1 = p(-180, -90)
x2, y2 = p(180, 90)
We create an Area plot command using the projection by setting the proj4string option,
and we define a view rectangle with the rectangle option. The new plot shows the new
projection.
In [8]:
a = plotcommands.Area(proj4string = p.srs, rectangle = [x1, x2, y1, y2])
b.setPlotCommands([a.text(), m.text()])
b.setPlotTime(datetime.datetime.now())
embed(b.plotImage(250, 250))
By default, this projection is centered on the zero latitude and zero longitude. We can
modify this by adding PROJ.4 options to the Proj object when we create it. Since we are
specifying the visible area using latitude and longitude coordinates, the values we use for
this must be reasonable for the projection in use. Typically, these will be equidistant from
the coordinates passed to Proj.
In [9]:
p = pyproj.Proj(proj="ortho", datum="WGS84", lat_0=60, lon_0=15)
x1, y1 = p(-10, 35)
x2, y2 = p(40, 85)
Regenerating the plot from the plot commands shows the new projection.
In [10]:
a = plotcommands.Area(proj4string = p.srs, rectangle = [x1, x2, y1, y2])
b.setPlotCommands([a.text(), m.text()])
b.setPlotTime(datetime.datetime.now())
embed(b.plotImage(250, 250))
In [11]:
cx = (x1 + x2)/2
cy = (y1 + y2)/2
w = x2 - x1
h = y2 - y1
zoomed_rectangle = [cx - w/4, cx + w/4, cy - h/4, cy + h/4]
In [12]:
a = plotcommands.Area(proj4string = p.srs, rectangle = zoomed_rectangle)
b.setPlotCommands([a.text(), m.text()])
b.setPlotTime(datetime.datetime.now())
embed(b.plotImage(250, 250))
Note that the size of the image we create tends to reflect the size of the area we want to display. Diana does not stretch plots to fit the output, but instead scales the output to ensure that the smallest edge of the plot is fully visible.
In [13]:
embed(b.plotImage(500, 250))
This can be confusing if you assumed that the edges of the image coincide with particular locations shown in the plot. If the rectangle's aspect ratio does not match the aspect ratio of the output, any assumptions like this do not apply.