Given a Numeric array of 2-D latitude-longitude
data, make a smooth, filled contour plot (with overlain
contour levels) with color bar and customized labeling.
Let's make some sample data (with some missing values)
and the longitude and latitude vectors.
(This is the same data and method
as in our simple example except
for the missing values.)
The grid is global and 15 x 15 degrees.
The data is in data and
consists of two "distorted bullseyes" centered around each pole
(positive over the North Pole and negative over the South Pole).
The latitude
lat vector goes from 90S to 90N, and the longitude
lon vector goes from -180E to 180E.
We set a few points in the mid-eastern
Pacific Ocean to "missing" by setting them equal
to a missing value value (1e+20):
import MA
import Numeric as N
import cdms, vcs
lat = N.arange(13) * 15.0 - 90.
lon = N.arange(25) * 15.0 - 180.0
data = N.outerproduct( N.sin(lat*N.pi/360.) \
, N.cos(lon*N.pi/360.) )
missing = 1e20
data[6:8,2:3] = missing
data[4:5,3:4] = missing
The data array is 13 x 24,
latitude x longitude. This means each
row corresponds to a constant latitude, and each
column is a constant longitude. We import MA since we'll
use that package to help us deal with missing values later
on.
We first open a VCS canvas and create our own graphics methods and templates. We'll be revising these two objects to customize our plot. Separate graphics methods are created for the isoline and isofill methods. Here we also create the template for the isofill plot; later on we'll create the template for the overlain contour lines based upon our customized isofill template:
v = vcs.init()
my_fill_gm = v.createisofill('my_fill_gm', 'default')
my_line_gm = v.createisoline('my_line_gm', 'default')
my_fill_tpl = v.createtemplate('my_fill_tpl', 'default')
Next we
create a custom color map
with a blue gradient for negative values and a red gradient for
positive values. The switch from blue to red occurs between
color index 127 (last blue) to 128 (first red);
this will be consistent with our use of
getcolors to set
color indices to contour levels later on in this example.
After making the color map, we use
setcolormap to make it the active color map:
my_cmap = v.createcolormap('my_cmap', 'default')
for icell in range(16,128):
rvalue = int( float(icell-16)/(127.-16.) * 75. )
gvalue = rvalue
bvalue = 100
my_cmap.index[icell] = [rvalue, gvalue, bvalue]
for icell in range(128,240):
rvalue = 100
gvalue = 75 - int( float(icell-128)/(239.-128.) \
* 75. )
bvalue = gvalue
my_cmap.index[icell] = [rvalue, gvalue, bvalue]
v.setcolormap('my_cmap')
A picture of what the color bar looks like can be found here.
In this example we will customize the size and position of the plot on the page. To make this easier, we make variables of the coordinates of a bounding box, and later on we will use those coordinates in calculating positioning of labels and titles. The bounding box will describe the corners of the plot data field (i.e. the origin and upper-right corner of the data field). Font heights are given in arbitrary units and are set for the overall title, the x-axis and y-axis titles, and the tick labels (for both axis). The tick length and the bounding box coordinates are in normalized coordinates. Remember, here we're just making temporary variables; the actual setting of plot elements occurs later on:
bbllx = 0.2 bblly = 0.3 bburx = 0.9 bbury = 0.9 title_fontht = 40 xytitle_fontht = 30 ticklabel_fontht = 30 ticklength = 0.01
Here we make variables that specify which font to use for the tick labels and all of the titles (overall and the axes). Later on we'll set the plot attributes controlling those fields to these values. Font code 1 is a monaco-looking font:
ticklabel_font = 1
titles_font = 1
Because font heights are in arbitrary units, to make our calculations easier we define a conversion factor that we multiply font height coordinates by in order to get the equivalent value in normalized coordinates. This "magic number" was obtained through trial-and-error:
fontht2norm = 0.0007
Now the complicated part begins. VCS knows where to put everything
based upon the template that's provided to it, but our current template
is bare-bones. What we have to do now is create secondary objects
that will define the attributes of our template to be what we want
it to be (e.g. font, font height, alignment).
First, let's create text-table and text-orientation
objects for the x- and y-axis tick labels. We use the standard default
objects for all except the text-orientation for the x-axis tick label,
which we want to be centered at the ticks, so we use the
'defcenter' object:
ttab_xticklabel = \
v.createtexttable( 'xticklabel_ttab', 'default' )
ttab_yticklabel = \
v.createtexttable( 'yticklabel_ttab', 'default' )
tori_xticklabel = \
v.createtextorientation( 'xticklabel_tori' \
, 'defcenter' )
tori_yticklabel = \
v.createtextorientation( 'yticklabel_tori' \
, 'default' )
ttab_xticklabel.font = ticklabel_font
ttab_yticklabel.font = ticklabel_font
ttab_xticklabel.color = 241
ttab_yticklabel.color = 241
tori_xticklabel.height = ticklabel_fontht
tori_yticklabel.height = ticklabel_fontht
tori_yticklabel.valign = 'half'
The default template settings for overall and x- and y-axis
titles are a bit plain. Here we'll create gussied-up text-table and
text-orientation objects to replace those fields in the template, and
then pass in the string values via the plot method's
xname etc. keywords; this has the benefit of keeping
everything in the template and accessed via one plot
command.
(See the complex
x-y line plot example for another way of placing titles.)
In the template, the x-axis title is element xname,
the y-axis title is yname, and the overall title
is element dataname
(although the overall title
value is passed in plot via the name keyword).
For both the overall title and x-axis title we want the text-orientation
to be horizontally centered. For the y-axis title we want the
text-orientation to be vertical and vertically centered.
We also set the font and font height.
ttab_title = \
v.createtexttable( 'title_ttab', 'default' )
tori_title = \
v.createtextorientation( 'title_tori' \
, 'defcenter' )
ttab_title.font = titles_font
ttab_title.color = 241
tori_title.height = title_fontht
ttab_xytitle = \
v.createtexttable( 'xytitle_ttab', 'default' )
ttab_xytitle.font = titles_font
ttab_xytitle.color = 241
tori_xtitle = \
v.createtextorientation( 'xtitle_tori' \
, 'defcenter' )
tori_ytitle = \
v.createtextorientation( 'ytitle_tori' \
, 'defcentup' )
tori_xtitle.height = xytitle_fontht
tori_ytitle.height = xytitle_fontht
my_fill_tpl.dataname.texttable = ttab_title
my_fill_tpl.dataname.textorientation = tori_title
my_fill_tpl.xname.texttable = ttab_xytitle
my_fill_tpl.xname.textorientation = tori_xtitle
my_fill_tpl.yname.texttable = ttab_xytitle
my_fill_tpl.yname.textorientation = tori_ytitle
We turn on the
titling elements using the priority keyword
and turn off some other text string fields that are on by default
in the default template:
my_fill_tpl.dataname.priority = 1
my_fill_tpl.xname.priority = 1
my_fill_tpl.yname.priority = 1
my_fill_tpl.mean.priority = 0
my_fill_tpl.max.priority = 0
my_fill_tpl.min.priority = 0
Here we define the coordinates of the data field, the box around
the data field,
the x-axis ticks (whose lengths are set to ticklength),
and the location of the x-axis tick labels.
We also replace the default text-table and text-orientation objects
for the x-axis tick labels with the secondary objects we've created
for our plot. Finally we turn off the upper x-axis and right y-axis
ticks, so that the only ticks there are are on the lower x-axis and
left y-axis:
my_fill_tpl.data.x1 = my_fill_tpl.box1.x1 = bbllx my_fill_tpl.data.y1 = my_fill_tpl.box1.y1 = bblly my_fill_tpl.data.x2 = my_fill_tpl.box1.x2 = bburx my_fill_tpl.data.y2 = my_fill_tpl.box1.y2 = bbury my_fill_tpl.xtic1.y1 = bblly - ticklength my_fill_tpl.xtic1.y2 = bblly my_fill_tpl.xlabel1.y = bblly - (3.0 * ticklength) my_fill_tpl.xlabel1.texttable = ttab_xticklabel my_fill_tpl.xlabel1.textorientation = tori_xticklabel my_fill_tpl.xtic2.priority = 0 my_fill_tpl.ytic2.priority = 0
Since latitude and longitude is usually given in degrees N/S and E/W, we reformat the x- and y-axis labels to look that way. The x-axis labels are done straightforwardly:
my_fill_gm.xticlabels1 = \
{ 0.0:"0"
, 45.0:"45E", -45.0:"45W"
, 90.0:"90E", -90.0:"90W"
, 135.0:"135E", -135.0:"135W"
, 180.0:"180E", -180.0:"180W" }
The y-axis labels are a little more complicated, since location of the y-axis title depends on how long the tick labels are. So we make the tick labels dictionary, calculate the length of the longest y-axis tick label, then set the y-axis ticks, location of their labels, and replace the text-table and text-orientation objects for the y-axis tick labels with our custom objects:
ylabels = { 0.0:"EQ"
, 30.0:"30N", -30.0:"30S"
, 60.0:"60N", -60.0:"60S"
, 90.0:"90N", -90.0:"90S" }
my_fill_gm.yticlabels1 = ylabels
longest_ylabels = max([ len(ylabels.values()[i]) \
for i in range(len(ylabels)) ])
my_fill_tpl.ytic1.x1 = bbllx - ticklength
my_fill_tpl.ytic1.x2 = bbllx
my_fill_tpl.ylabel1.x = bbllx - ticklength \
- ( longest_ylabels \
* ticklabel_fontht * fontht2norm )
my_fill_tpl.ylabel1.texttable = ttab_yticklabel
my_fill_tpl.ylabel1.textorientation = tori_yticklabel
Now we set the coordinates for the overall title, x-axis title, and y-axis title. Note that these template attributes are not set as lists, unlike in the case of a text-combined object:
my_fill_tpl.dataname.x = (bburx-bbllx)/2.0 + bbllx
my_fill_tpl.dataname.y = bbury + ( 1.7 * title_fontht \
* fontht2norm )
my_fill_tpl.xname.x = (bburx-bbllx)/2.0 + bbllx
my_fill_tpl.xname.y = my_fill_tpl.xlabel1.y \
- (1.7 * title_fontht * fontht2norm)
my_fill_tpl.yname.x = my_fill_tpl.ylabel1.x \
- ( 1.6 * title_fontht * fontht2norm )
my_fill_tpl.yname.y = (bbury-bblly)/2.0 + bblly
We calculate the maximum and minimum non-missing data values and calculate the contour levels. The colors that describe these contour levels are set so color indices 16 to 127 cover negative values and 128 to 239 cover positive values.
my_min = MA.minimum(MA.masked_values(data, missing))
my_max = MA.maximum(MA.masked_values(data, missing))
con_levels = vcs.mkscale(my_min, my_max, 16)
my_fill_gm.levels = con_levels
my_fill_gm.fillareacolors = vcs.getcolors(con_levels, split=1)
A little clean-up task: we make sure the color bar lettering is
black. To do so, we create a text-table object that has color
set to black and copy it to the template legend element:
ttab_black = v.createtexttable('ttab_black', 'default')
ttab_black.color = 241
my_fill_tpl.legend.texttable = ttab_black
Now that we have a "good" template to use when we make the isofill call, let's copy it to another template for use in the isoline call (so the contour lines will overlie the filled contours). We turn off text etc. on the isoline template to make sure only the contours are written out:
my_line_tpl = v.createtemplate('my_line_tpl', 'my_fill_tpl')
my_line_tpl.dataname.priority = 0
my_line_tpl.xname.priority = 0
my_line_tpl.yname.priority = 0
my_line_tpl.xlabel1.priority = 0
my_line_tpl.xlabel2.priority = 0
my_line_tpl.ylabel1.priority = 0
my_line_tpl.ylabel2.priority = 0
my_line_tpl.xtic1.priority = 0
my_line_tpl.xtic2.priority = 0
my_line_tpl.ytic1.priority = 0
my_line_tpl.ytic2.priority = 0
my_line_tpl.box1.priority = 0
We use the same contouring levels for the contour line plot as the filled contour plot and turn on contour labels in our isoline graphics method. The line style for the contour lines are set to dashed if the contour is negative and solid it is positive or 0; this is done through creating a line secondary object:
my_line_gm.level = con_levels
my_line_gm.label = 'y'
line_pos_con = v.createline('pos_con_line', 'default')
line_neg_con = v.createline('neg_con_line', 'default')
line_pos_con.type = 'solid'
line_neg_con.type = 'dash'
my_line_gm.line = \
N.where( N.array(con_levels) >= 0.0 \
, line_pos_con, line_neg_con ).tolist()
For some reason the x- and y-axis need to be specified as CDMS axis objects, rather than just vectors. This is done here:
lonAxis = cdms.createAxis(lon) latAxis = cdms.createAxis(lat)
We now render the plot, first plotting the filled contours (with all the labeling, the color bar, etc.) and then overplot the contour levels. Places where data is missing will be filled with white:
v.plot( data, my_fill_gm, my_fill_tpl \
, xaxis=lonAxis, yaxis=latAxis \
, continents=1 \
, xname='Longitude [deg]', yname='Latitude [deg]' \
, name='A Complex Contour Plot' )
v.plot( data, my_line_gm, my_line_tpl \
, xaxis=lonAxis, yaxis=latAxis )
We can write Postscript and GIF output of the plot with the following commands:
v.postscript('latlon_plot.ps')
v.gif('latlon_plot.gif', 'r')
And the plot we get is:
Notes: This discussion applies to CDAT 3.3.