Getting started using NCL - example 8 - contour and multi-curve XY plot

Example 8 - contour and multi-curve XY plot

Example 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11

This example reads an ASCII file of random seismic data, interpolates it to a 2-dimensional grid, and then creates a contour plot of this grid and an XY plot of slices of the grid. The interpolation is done using Natgrid, a natural neighbor gridding package (based on Dave Watson's nngridr package).

This example also makes use of a resource file to set all fonts in the plot to "Times-Roman" (unless otherwise specified). It also shows how to write the graphical output to an NCGM file and a PostScript file, in addition to an X11 window.

To run this example, you must download the following files:

gsun08n.ncl
gsun08n.res
gsn_code.ncl

and then type:

ncl < gsun08n.ncl

Output from example 8

Frame 1	Frame 2

(Click on any frame to see it enlarged.)

NCL code for example 8

  1. load "gsn_code.ncl"
  2. 
  3. begin
  4. 
  5.   ascii_filename = "$NCARG_ROOT/lib/ncarg/data/asc/seismic.asc"
  6. 
  7.   seismic = asciiread(ascii_filename,(/52,3/),"float")
  8. 
  9.   x = seismic(:,0)  ; Column 1 of file contains X values.
 10.   y = seismic(:,1)  ; Column 2 of file contains Y values.
 11.   z = seismic(:,2)  ; Column 3 of file contains Z values.
 12. 
 13.   numxout = 20     ; Define output grid for call to "natgrids".
 14.   numyout = 20
 15.   xmin    = min(x)
 16.   ymin    = min(y)
 17.   xmax    = max(x)
 18.   ymax    = max(y)
 19.   xc      = (xmax-xmin)/(numxout-1)
 20.   yc      = (ymax-ymin)/(numyout-1)
 21.   xo      = xmin + ispan(0,numxout-1,1)*xc
 22.   yo      = ymin + ispan(0,numyout-1,1)*yc
 23. 
 24.   zo = natgrids(x, y, z, xo, yo)  ; Interpolate.
 25. 
 26.   xo@long_name = "x values"   ; Define some attributes.
 27.   yo@long_name = "y values"
 28.   zo@long_name = "Depth of a subsurface stratum"
 29. 
 30.   xwks   = gsn_open_wks( "x11","gsun08n") ; Open an X11 workstation.
 31.   cgmwks = gsn_open_wks("ncgm","gsun08n") ; Open an NCGM workstation.
 32.   pswks  = gsn_open_wks(  "ps","gsun08n") ; Open a PS workstation.
 33. 
 34.   cmap = (/(/1., 1., 1./), (/0., 0., 0./), (/1., 0., 0./), (/1., 0., .4/), \
 35.            (/1., 0., .8/), (/1., .2, 1./), (/1., .6, 1./), (/.6, .8, 1./), \
 36.            (/.2, .8, 1./), (/.2, .8, .6/), (/.2, .8, 0./), (/.2, .4, .0/), \
 37.            (/.2, .4, .4/), (/.2, .4, .8/), (/.6, .4, .8/), (/.6, .8, .8/), \
 38.            (/.6, .8, .4/), (/1., .6, .8/)/)
 39. 
 40.   gsn_define_colormap(  xwks,cmap) ; Define a color map for each of
 41.   gsn_define_colormap(cgmwks,cmap) ; the three workstations you
 42.   gsn_define_colormap( pswks,cmap) ; just opened.
 43. 
 44. ;----------- Begin first plot -----------------------------------------
 45. 
 46.   resources                       = True
 47.   resources@sfXArray              = xo            ; X axes data points
 48.   resources@sfYArray              = yo            ; Y axes data points
 49. 
 50.   resources@tiMainString          = zo@long_name  ; Main title
 51.   resources@tiMainFont            = "Times-Bold"
 52.   resources@tiXAxisString         = "x values"    ; X axis label.
 53.   resources@tiYAxisString         = "y values"    ; Y axis label.
 54. 
 55.   resources@cnFillOn              = True     ; Turn on contour fill.
 56.   resources@cnInfoLabelOn         = False    ; Turn off info label.
 57.   resources@cnLineLabelsOn        = False    ; Turn off line labels.
 58. 
 59.   resources@lbOrientation         = "Horizontal" ; Draw it horizontally.
 60.   resources@lbPerimOn             = False        ; Turn off perimeter.
 61.   resources@pmLabelBarDisplayMode = "Always"     ; Turn on a label bar.
 62.   resources@pmLabelBarSide        = "Bottom"     ; Change location of 
 63.                                                  ; label bar.
 64.   resources@vpYF = 0.9   ; Change Y location of plot.
 65. 
 66.   zo!0 = "i"  ; Name the dimensions of "zo".
 67.   zo!1 = "j"
 68. 
 69.   contour = gsn_contour(xwks,zo(j|:,i|:),resources) ; Draw a contour plot.
 70. 
 71. ;----------- Begin second plot -----------------------------------------
 72. 
 73.   delete(resources)   ; Start with a new list of resources.
 74. 
 75.   resources                     = True
 76.   resources@tiMainString        = ":F26:slices"   ; Define a title.
 77. 
 78.   resources@xyLineColors        = (/2,8,10,14/)   ; Define line colors.
 79.   resources@xyLineThicknessF    = 3.0             ; Define line thickness.
 80. 
 81.   resources@pmLegendDisplayMode    = "Always"     ; Turn on the drawing
 82.                                                   ; of a legend.
 83.   resources@pmLegendZone           = 0            ; Change the location
 84.   resources@pmLegendOrthogonalPosF = 0.31         ; of the legend
 85.   resources@lgJustification        = "BottomRight"
 86. 
 87.   resources@pmLegendWidthF         = 0.4          ; Change width and
 88.   resources@pmLegendHeightF        = 0.12         ; height of legend.
 89. 
 90.   resources@pmLegendSide           = "Top"        ; Change location of 
 91.   resources@lgPerimOn              = False        ; legend and turn off
 92.                                                   ; the perimeter.
 93. 
 94.   resources@xyExplicitLegendLabels = (/"zo(i,2)","zo(i,4)","zo(i,6)","zo(i,8)"/)
 95. 
 96.   resources@vpYF      = 0.90  ; Change size and location of plot.
 97.   resources@vpXF      = 0.18
 98.   resources@vpWidthF  = 0.74
 99.   resources@vpHeightF = 0.74
100. 
101.   resources@trYMaxF = 980  ; Set the maximum value for the Y axes.
102. 
103.   xy = gsn_xy(xwks,xo,zo(j|2:8:2,i|:),resources) ; Draw an XY plot.
104. 
105. ;----------- Draw to other workstations  ------------------------------
106. 
107.   NhlChangeWorkstation(contour,cgmwks)  ; Change the workstation that the
108.   NhlChangeWorkstation(xy,cgmwks)       ; contour and XY plot is drawn to.
109.   draw(contour)                         ; Draw the contour plot to the new
110.   frame(cgmwks)                         ; workstation and advance the frame.
111.   draw(xy)                              ; Draw the XY plot to the new
112.   frame(cgmwks)                         ; workstation and advance the frame.
113. 
114.   NhlChangeWorkstation(xy,pswks)        ; Do the same for the PostScript
115.   NhlChangeWorkstation(contour,pswks)   ; workstation.
116.   draw(contour)
117.   frame(pswks)
118.   draw(xy)
119.   frame(pswks)
120. 
121.   delete(xy)       ; Clean up.
122.   delete(contour)
123. end
124.

Explanation of example 8

Lines 5-11:

  ascii_filename = "$NCARG_ROOT/lib/ncarg/data/asc/seismic.asc"

  seismic = asciiread(ascii_filename,(/52,3/),"float")

  x = seismic(:,0)
  y = seismic(:,1)
  z = seismic(:,2)

The ASCII file "seismic.asc" contains seismic data taken from Dave Watson's book, "nngridr." Dave notes that this data set is a well-known data set from the geological literature. See Davis, J.C., 1973, "Statistics and data analysis in geology": John Wiley & Sons, New York, 550p.

The data are in three columns where the first column contains X values, the second column contains Y values, and the third column contains Z values. Each column has 52 values.

First read in these values with the function asciiread into a 2-dimensional float array called seismic that is dimensioned 52 x 3, and then extract each "column" of data to 1-dimensional float arrays x, y, and z.

Lines 13-22:

  numxout = 20
  numyout = 20
  xmin    = min(x)
  ymin    = min(y)
  xmax    = max(x)
  ymax    = max(y)
  xc      = (xmax-xmin)/(numxout-1)
  yc      = (ymax-ymin)/(numyout-1)
  xo      = xmin + ispan(0,numxout-1,1)*xc
  yo      = ymin + ispan(0,numyout-1,1)*yc

Define an output grid for the interpolation function you call in the next step. Let's assume you want the output grid to be over the same range of values as the input grid, only at equally spaced intervals, and that you want it to be dimensioned 20 x 20.

First use the NCL built-in functions min and max to get the minimum and maximum values of the x and y arrays (min and max take an array of any dimensionality and return the minimum and maximum values of that array).

Then use the function ispan to get a range of evenly spaced integer values. The first argument of ispan is an integer starting value, the second argument is an integer ending value, and the third argument is the stride.

Using the minimum/maximum values and the ispan function, create the output grid arrays xo and yo.

Line 24:

  zo = natgrids(x, y, z, xo, yo)

Interpolate the random data to a 2-dimensional grid using natgrids. The first two arguments of natgrids are 1-dimensional float arrays of the same size, containing the X and Y coordinates of the input data points. The third argument is another 1-dimensional float array (of the same size as the X and Y arrays) containing the values of the input function. The last two arguments are 1-dimensional float arrays (of any size, numxout and numyout in this case) containing the X and Y coordinates of the output data grid. The function natgrids returns a 2-dimensional float array dimensioned numxout by numyout, so in this case, zo is dimensioned 20 x 20.

Lines 26-28:

  xo@long_name = "x values"
  yo@long_name = "y values"
  zo@long_name = "Depth of a subsurface stratum"

Define an attribute called "long_name" for xo, yo, and zo. By default, if the X and/or Y arrays passed to gsn_xy have "long_name" as an attribute, then the value of this attribute is used to put a label on the corresponding X and/or Y axis (this takes effect in the second plot of this example).

Lines 30-32:

  xwks   = gsn_open_wks( "x11","gsun08n")
  cgmwks = gsn_open_wks("ncgm","gsun08n")
  pswks  = gsn_open_wks(  "ps","gsun08n")

To have the graphical output go to multiple workstations, you need to call gsn_open_wks for each one. Here you are creating one of each kind of workstation available: an X11 window, an NCGM file, and a PostScript file.

In this example, the second argument of gsn_open_wks comes into play in two ways. First, it is used to determine the name of the NCGM and PostScript files that are created (with ".ncgm" and ".ps" appended, respectively). Second, it is used to specify the name of a resource file (with ".res" appended) to read in (if the file exists). In this case, the NCGM, PostScript, and resource files are called "gsun08n.ncgm," "gsun08n.ps," and "gsun08n.res." The resource file is discussed below.

Lines 34-42:

  cmap = (/(/1., 1., 1./), (/0., 0., 0./), (/1., 0., 0./), (/1., 0., .4/), \
           (/1., 0., .8/), (/1., .2, 1./), (/1., .6, 1./), (/.6, .8, 1./), \
           (/.2, .8, 1./), (/.2, .8, .6/), (/.2, .8, 0./), (/.2, .4, .0/), \
           (/.2, .4, .4/), (/.2, .4, .8/), (/.6, .4, .8/), (/.6, .8, .8/), \
           (/.6, .8, .4/), (/1., .6, .8/)/)

  gsn_define_colormap(  xwks,cmap)
  gsn_define_colormap(cgmwks,cmap)
  gsn_define_colormap( pswks,cmap)

Define a color table with 18 entries. For each workstation that you want the color table to take effect for, you have to call the gsn_define_colormap procedure.

Details on creating your own color map are covered in example 2.

Line 44:

;----------- Begin first plot -----------------------------------------

Draw a contour plot of zo, the 2-dimensional array returned from natgrids. Label the axes of the contour plot, fill the contour levels, and put a label bar at the bottom.

Lines 47-48:

  resources@sfXArray              = xo
  resources@sfYArray              = yo

Set some ScalarField resources so that the X and Y axis values are the values of the output grid that you defined, rather than index values (which is the default).

Lines 50-53:

  resources@tiMainString          = zo@long_name
  resources@tiMainFont            = "Times-Bold"
  resources@tiXAxisString         = "x values"
  resources@tiYAxisString         = "y values"

Set some resources for putting a title at the top of the contour plot and for labeling the X and Y axes. The font for the main title is changed to index 26 ("Times-Bold").

Lines 59-62:

  resources@lbOrientation         = "Horizontal"
  resources@lbPerimOn             = False
  resources@pmLabelBarDisplayMode = "Always"
  resources@pmLabelBarSide        = "Bottom"

Turn on the drawing of a label bar. By default, the label bar is vertical and centered at the right side of the plot. Let's assume you want to draw it at the bottom and horizontally. Since the drawing of a label bar associated with a plot is controlled by the PlotManager, the PlotManager also controls the location of it. Thus, use PlotManager resources to change the location of the label bar.

The PlotManager is described in more detail in example 2.

Line 64:

  resources@vpYF = 0.9

Drawing a label bar at the bottom causes the plot to go out of the viewport window in this case, so adjust the Y position of the plot on the viewport.

Lines 66-67:

  zo!0 = "i"
  zo!1 = "j"

Name the dimensions of zo. This is useful later when you need to swap the indices of zo.

Line 69:

  contour = gsn_contour(xwks,zo(j|:,i|:),resources)

Create and draw a contour plot of zo. Note that natgrids returns zo(i,j) as the interpolated value at the coordinate (x(i),y(j)) for 0 < i < numxout and 0 < j < numyout. For the purpose of drawing a contour plot, however, NCL identifies the fastest varying dimension (the second dimension of a 2-dimensional array in NCL) with the Cartesian X-axis. So, if it is desired to have the x values along the Cartesian X-axis, you need to swap the x and y indices before contouring. You named the dimensions of zo in lines 66-67, so use the named subscripting capabilities of NCL to swap the indices of zo with the syntax (j|:,i|:). See the section "Named subscripting" in the "Basics" chapter for more information on this syntax.

Note on setting fonts: All the fonts in this contour plot (except the title font) are "Times-Roman." By using a resource file called "gsun08n.res" containing the single entry:

*Font : Times-Roman

you are able to set the font to "Times-Roman" for all characters drawn in the plot. Any resources that you set within the NCL script override what is set in a resource file. For example, in line 51, you set tiMainFont to "Times-Bold," so this font is used for the main title instead of "Times-Roman". The name of the resource file must have the same name as the second string passed to gsn_open_wks, with a ".res" appended.

By default, NCL looks for the resource file in the same directory as where NCL is being executed from. If you want to change this default directory, then read the section "Resource files" in the "Going beyond the basics" chapter. Also, as noted in the "Resource files" section, when you set resources in a resource file, don't put double quotes around predefined strings. Use double quotes only when you actually want the quotes to be part of the string itself.

A more detailed resource file is used in example 9.

Line 71:

;----------- Begin second plot -----------------------------------------

Take some slices of the 2-dimensional zo data and generate an XY plot with four of these slices.

Line 73:

  delete(resources)

Since you previously used the resources variable to define some ContourPlot resources, you need to start with a new list of resources for the XyPlot resources, and thus should delete the variable before redefining it.

Lines 78-79:

  resources@xyLineColors        = (/2,8,10,14/)
  resources@xyLineThicknessF    = 3.0

Set some XY plot resources to define the line colors and thickness for each line.

Lines 81-91:

  resources@pmLegendDisplayMode    = "Always"

  resources@pmLegendZone           = 0
  resources@pmLegendOrthogonalPosF = 0.31
  resources@lgJustification        = "BottomRight"

  resources@pmLegendWidthF         = 0.4
  resources@pmLegendHeightF        = 0.12

  resources@pmLegendSide           = "Top"
  resources@lgPerimOn              = False

Turn on the drawing of a legend. By default, the legend is centered below the plot. To draw the legend in the upper right corner and inside the XY plot, you need to use PlotManager resources to control where the legend is drawn.

By default, the legend appears outside the plot. If you set pmLegendZone to 0, then this allows you to put the legend inside the plot. Setting the resources pmLegendOrthogonalPosF, pmLegendSide, and lgJustification to 0.31, "Top", and "BottomRight", respectively moves the legend to the inside upper right corner of the XY plot.

Legend resources start with "lg" and are defined in the "Legend resource descriptions" section of the NCAR Graphics 4.1 Reference Manual.

Line 94:

  resources@xyExplicitLegendLabels = (/"zo(i,2)","zo(i,4)","zo(i,6)","zo(i,8)"/)

Set this resource to define labels for the four curves in the legend.

Lines 96-99:

  resources@vpYF      = 0.90
  resources@vpXF      = 0.18
  resources@vpWidthF  = 0.74
  resources@vpHeightF = 0.74

Adjust the size and location of the plot so it doesn't run out of the viewport.

Line 101:

  resources@trYMaxF = 980

Change the maximum value of the upper bound of the Y axis to 980.0. The actual maximum of the data is 942.131, which you can verify with print(zo). This is done to allow room for the legend you created above to fit inside the top of the plot.

Line 103:

  xy = gsn_xy(xwks,xo,zo(j|2:8:2,i|:),resources)

To pass multiple curves to gsn_xy, you must pass a 2-dimensional array where the first dimension is the number of curves and the second dimension is the number of points. Assuming that each point of zo is represented by zo(i,j) and that you want to take slices of the data at j = 2, 4, 6, and 8 for all i, then use the syntax (j|2:8:2,i|:).

Note: You may see a warning message that looks like the following:

warning:XyPlotSetValues:Setting xyComputeYMax to False because trYMaxF was specified

You can safely ignore this warning. To get rid of it, add the following line before the call to gsn_xy:

  resources@xyComputeYMax = False

Line 105:

;----------- Draw to other workstations  ------------------------------

So far you've drawn the graphical output to the X workstation only. Now draw the contour and XY plot to the other two workstations that you opened in lines 31 and 32.

Lines 107-112:

  NhlChangeWorkstation(contour,cgmwks)
  NhlChangeWorkstation(xy,cgmwks)
  draw(contour)
  frame(cgmwks)
  draw(xy)
  frame(cgmwks)

To draw the contour and XY plot to the NCGM workstation that you opened, you need to change the current workstation to the NCGM workstation, and then redraw the plots and advance the frame. The gsn_* plotting functions call draw and frame for you (unless otherwise specified), so that's why you didn't need to call draw and frame when the plots were drawn to the X11 window (via the calls to gsn_contour and gsn_xy).

Call the procedure NhlChangeWorkstation to change the workstation being drawn to. The first argument is the plot that you want to draw (returned from a previous call to one of the gsn_* plotting functions), and the second argument is the new workstation you want to draw the plot to (returned from a previous call to gsn_open_wks).

After you've changed the workstation being drawn to, redraw the plots and advance the frame with calls to draw and frame. If you don't call frame after you call draw, then the two plots are drawn in the same frame, one on top of the other.

The name of the NCGM file is the second string that you passed to the call to gsn_open_wks, with the suffix ".ncgm" appended. In this case, the NCGM file name is "gsun08n.ncgm". If you pass an empty string for the second argument to gsn_open_wks, then the NCGM file is named "gsnapp.ncgm".

For information on how to view and edit the NCGM file, see the section "Viewing and editing your CGMs and raster images" in Chapter 4 of the NCAR Graphics Fundamentals document. If you are running under X Windows, then you can quickly view your NCGM file with the following command:

ctrans -d X11 name.ncgm

After the first frame is drawn, click on it with the left mouse button to advance the frame or to exit ctrans.

Lines 114-119:

  NhlChangeWorkstation(xy,pswks)
  NhlChangeWorkstation(contour,pswks)
  draw(contour)
  frame(pswks)
  draw(xy)
  frame(pswks)

To draw the plots to the PostScript workstation that you opened, do the same thing that you did for the NCGM workstation. The name of the PostScript file is the second string that you passed to the call to gsn_open_wks, with the suffix ".ps" appended. In this case, the PostScript file name is "gsun08n.ps". If you pass an empty string for the second argument to gsn_open_wks, then the PostScript file is named "gsnapp.ps".