A Complex Latitude-Longitude Contour Plot


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.

Sample Data

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]

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.gif('latlon_plot.gif', 'r')

And the plot we get is:

Full Source Code For Example

Notes: This discussion applies to CDAT 3.3.

Return to the Tips and Examples index page.

Updated: February 16, 2004 by Johnny Lin <email address>. License.