Getting started using NCL - example 4

Example 4 - streamline plots

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

This example reads in a GRIB file and creates three streamline plots, each using a different dataset. Resources are used to change aspects of each streamline plot.

To find out more about the GRIB file being opened, this example uses NCL inquiry functions like isatt, getfilevarnames, getfilevaratts, and getfilevardims to determine what variables and attributes are in the file. These functions are all documented in the "Built-in NCL functions and procedures" section of the NCL Reference Manual.

This example also shows how to write some of the GRIB data to a netCDF file.

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

gsun04n.ncl
gsn_code.ncl

and then type:

ncl < gsun04n.ncl

Output from example 4

Frame 1	Frame 2	Frame 3

(Click on any frame to see it enlarged.)

NCL code for example 4

  1. load "gsn_code.ncl"
  2. 
  3. begin
  4. 
  5.   data_dir = ncargpath("data")
  6.   grb_file = addfile(data_dir + "/grb/ced1.lf00.t00z.eta.grb","r")
  7. 
  8.   names = getfilevarnames(grb_file)  ; Get the variable names in the
  9.   print(names)                       ; GRIB file and print them out.
 10. 
 11.   atts = getfilevaratts(grb_file,names(0)) ; Get the variable attributes and
 12.   dims = getfilevardims(grb_file,names(0)) ; dimension names from the GRIB
 13.   print(atts)                              ; file and print them out.
 14.   print(dims)
 15. 
 16.   wks = gsn_open_wks("x11","gsun04n") ; Open an X11 workstation.
 17. 
 18. ;----------- Begin first plot -----------------------------------------
 19. 
 20.   resources             = True
 21. 
 22.   uvar = grb_file->U_GRD_6_ISBL(0,:,:) ; Store two GRIB variables to 
 23.   vvar = grb_file->V_GRD_6_ISBL(0,:,:) ; local NCL variables.
 24. 
 25.   if(isatt(uvar,"units").eq.True) ; Check if "uvar" has an attribute "units".
 26.     resources@tiMainString = "GRD_6_ISBL (u,v " + uvar@units + ")"
 27.   else
 28.     resources@tiMainString = "GRD_6_ISBL"
 29.   end if              ; The space is required between "end" and "if".
 30.   resources@tiMainFont    = "Times-Roman"
 31.   resources@tiXAxisString = "streamlines"
 32. 
 33.   plot = gsn_streamline(wks,uvar(::2,::2),vvar(::2,::2),resources) ; Create and
 34.                                                  ; draw a streamline plot.
 35. ;----------- Begin second plot -----------------------------------------
 36. 
 37.   delete(uvar) ; Delete uvar and vvar since you are going to use
 38.   delete(vvar) ; them again, but with different sized arrays.
 39. 
 40.   uvar = grb_file->U_GRD_6_TRO ; Store two different GRIB file variables
 41.   vvar = grb_file->V_GRD_6_TRO ; to uvar and vvar.
 42. 
 43.   if(isatt(uvar,"units"))
 44.     resources@tiMainString = "GRD_6_TRO (u,v " + uvar@units + ")"
 45.   else
 46.     resources@tiMainString = "GRD_6_TRO"
 47.   end if
 48. 
 49.   resources@tiXAxisFont   = "Times-Roman"  ; Change the default font used.
 50.   resources@tmXBLabelFont = "Times-Roman"
 51.   resources@tmYLLabelFont = "Times-Roman"
 52. 
 53.   resources@stLineColor = "green" ; Change streamlines to green.
 54. 
 55.   plot = gsn_streamline(wks,uvar(::2,::2),vvar(::2,::2),resources) ; Draw streamline plot.
 56. 
 57. ;----------- Begin third plot -----------------------------------------
 58. 
 59.   getvalues plot                        ; Retrieve some resource values.
 60.     "stArrowLengthF"    : arrowlength
 61.     "stMinLineSpacingF" : spacing
 62.   end getvalues
 63.     
 64.   resources@stMinLineSpacingF = spacing * 2.0     ; Set some resources based
 65.   resources@stArrowLengthF    = arrowlength * 2.0 ; on resources you retrieved.
 66.   resources@stLineColor       = "white"           ; Change line colors back to 
 67.   resources@stLineThicknessF  = 1.5               ; white.
 68. 
 69.   uvar = grb_file->U_GRD_6_GPML(0,:,:) ; Get new GRIB variables.
 70.   vvar = grb_file->V_GRD_6_GPML(0,:,:)
 71. 
 72.   if(isatt(uvar,"units"))
 73.     resources@tiMainString = "GRD_6_GPML (u,v " + uvar@units + ")"
 74.   else
 75.     resources@tiMainString = "GRD_6_GPML"
 76.   end if
 77. 
 78.   plot = gsn_streamline(wks,uvar(::2,::2),vvar(::2,::2),resources) ; Draw streamline plot.
 79. 
 80.   delete(uvar) ; Remove some variables we don't need any more.
 81.   delete(vvar)
 82.   delete(plot)
 83. 
 84. ;----------- Write GRIB data to netCDF file --------------------------
 85. 
 86.   cdf_filename = "ced1.lf00.t00z.eta.nc"
 87.   system("/bin/rm -f " + cdf_filename)
 88.   cdf_file = addfile(cdf_filename,"c") ; Create a new netCDF file.
 89. 
 90.   cdf_file@title = "data from a GRIB file" ; Add some global attributes to
 91.   cdf_file@date = systemfunc("date")       ; the netCDF file.
 92. 
 93.   do i = 0,dimsizes(names)-1
 94.     names_char = stringtochar(names(i))
 95.     if(names_char(0:3).eq."PRES") then   ; Only write variables that start
 96.       print(names(i))                    ; with "PRES".
 97.       cdf_file->$names(i)$ = grb_file->$names(i)$
 98.     end if
 99.     delete(names_char)
100.   end do                       ; The space is required b/w "end" and "do".
101. end

Explanation of example 4

Line 1:

load "gsn_code.ncl"

Load the NCL script that contains the functions and procedures (the ones that start with "gsn_") that are used in this example. load in NCL works much like include works in C and Fortran 90 programs.

Line 3:

begin

Begin the NCL script.

Lines 5-6:

  data_dir = ncargpath("data")
  grb_file = addfile(data_dir + "/grb/ced1.lf00.t00z.eta.grb","r")

Open a GRIB file with addfile.

Note that you cannot write GRIB files, but you can only read them.

Lines 8-9:

  names = getfilevarnames(grb_file)
  print(names)

Retrieve the list of variable names in the GRIB file with the getfilevarnames function and print them out. names is an array of strings, each string being a variable name. This function is useful if you don't know ahead of time what variables are in the data file you opened with addfile.

Lines 11-14:

  atts = getfilevaratts(grb_file,names(0))
  dims = getfilevardims(grb_file,names(0))
  print(atts)
  print(dims)

Using the getfilevaratts and getfilevardims functions, retrieve the list of attributes and dimension names for the first variable in the GRIB file and print them out (just to show how this is done). These functions are useful if you need to know about a variable's attributes and dimension sizes.

Line 16:

  wks = gsn_open_wks("x11","gsun04n")

Open an X11 workstation for which to draw the streamline plots.

Line 18:

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

Draw a basic streamline plot with a title.

Line 20:

  resources             = True

Get ready to set up a list of resources for changing the look of the plot. See example 1 for an explanation on how resources are set. Streamline plot resources are part of the "StreamlinePlot" group and all start with the letters "st". They are documented in the "StreamlinePlot resource descriptions" section of the NCAR Graphics Reference Manual.

Streamline plots use the same data field as vector plots, so the data associated with a streamline plot is called a "VectorField."

Lines 22-23:

  uvar = grb_file->U_GRD_6_ISBL(0,:,:)
  vvar = grb_file->V_GRD_6_ISBL(0,:,:)

Retrieve the variables U_GRD_6_ISBL and V_GRD_6_ISBL from the GRIB file. The first dimension of these variables is a level, so use "(0,:,:)" to select the first level.

Lines 25-29:

  if(isatt(uvar,"units").eq.True)
    resources@tiMainString = "GRD_6_ISBL (u,v " + uvar@units + ")"
  else
    resources@tiMainString = "GRD_6_ISBL"
  end if

NCL contains "if-else-end if" statements just like C and Fortran with the exceptions that if statements in NCL do not need a then, there is no elseif equivalent, and you need a space between end and if since NCL won't recognize endif. To compare two expressions with each other, you can use ".eq.", ".lt.", ".ge.", etc., which is similar to Fortran syntax. For more information about if statements, see the "NCL Statements" section of the NCL Reference Manual.

In the if statement above, the function isatt returns True if "units" is an attribute of the variable uvar and False otherwise. If isatt is True in the if statement above, then the "units" attribute is used to create part of the main title. Otherwise, just the name of the variable is used as the title. Since isatt returns a logical value, you could have written the if statement without the .eq.True part:

  if(isatt(uvar,"units"))

Lines 30-31:

  resources@tiMainFont    = "Times-Roman"
  resources@tiXAxisString = "streamlines"

Set some resources to change the default font of the title and to label the X axis. Remember: predefined strings, like those listed in the font table, are case-insensitive, and you could have specified the above font with "times-roman" or "TIMES-ROMAN" or any another combination of uppercase and lowercase characters.

Line 33:

  plot = gsn_streamline(wks,uvar(::2,::2),vvar(::2,::2),resources)

Create and draw a streamline plot of the 2-dimensional arrays uvar and vvar. The first argument of the gsn_streamline function is the workstation variable returned from the previous call to gsn_open_wks. The next two arguments represent the 2-dimensional vector field to be plotted (they must have the same dimensions and can be of type float, double, or integer). The last argument is a logical value indicating whether you have set any resources. The resources you set here are just for titling the plot and for labeling the X axis.

The syntax "(::2)" is being used to "thin" the streamlines; that is, only every other data point is passed on to the plotting function.

The default streamline plot drawn contains labeled tick marks. Since no ranges have been defined for the X and Y axes in this plot, the range values default to 0 to n-1, where n is the number of points in that dimension.

Line 35:

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

Create a streamline plot using a different dataset, drawing the streamlines in color, and changing the font of the titles and labels in the plot.

Lines 37-38:

  delete(uvar)
  delete(vvar)

Delete the variables uvar and vvar since you are going to use these variables again. This is only necessary if you plan to assign data to uvar and vvar that is of a different type or dimensionality than when you originally initialized these two variables.

Lines 40-41:

  uvar = grb_file->U_GRD_6_TRO
  vvar = grb_file->V_GRD_6_TRO

Retrieve the variables U_GRD_6_TRO and V_GRD_6_TRO from the GRIB file. These two variables are 2-dimensional arrays, so you don't need to select a subset of either one.

Lines 43-47:

  if(isatt(uvar,"units"))
    resources@tiMainString = "GRD_6_TRO (u,v " + uvar@units + ")"
  else
    resources@tiMainString = "GRD_6_TRO"
  end if

Again, check if the variable in question has an attribute called "units", and if so, use it to create part of the main title.

Lines 49-51:

  resources@tiXAxisFont   = "Times-Roman"
  resources@tmXBLabelFont = "Times-Roman"
  resources@tmYLLabelFont = "Times-Roman"

Set some resources for changing the font of the X axis and the labels for the tick marks. Tick mark resources are part of the "TickMark" group and start with the letters "tm". They are documented in the "TickMark resource descriptions" section of the NCAR Graphics Reference Manual. The "XB" in tmXBLabelFont stands for "X bottom" and the "YL" in tmYLLabelFont stands for "Y left."

There are ways for setting fonts globally rather than having to set each font individually, and these are covered in example 8.

Line 53:

  resources@stLineColor = "green"

Change the color of the streamlines to green using the named colors capability of NCL. There are 650 valid color names you can use, but your color map must contain the named color you want to use or else NCL will use the closest match. For example, if your color map contains the following RGB values:

(/(/0.,0.,0./),(/1.,1.,1./),(/1.,0.,0./),(/0.,1.,0./),(/0.,0.,1./)/)

and you reference the color "CornflowerBlue" which is represented by the RGB value (/0.39215,0.58431,0.92941/), then the closest match in your color map is (/0.,0.,1./), and you will get a regular blue color instead of cornflower blue which is lighter. To assure that your color map contains the named colors you want to reference, add them specifically to your color map. GSUN example 7 shows how to set a color map using named colors.

For more information on named colors, see the section "Creating your own color map using named colors and RGB values" in the "Basics" chapter.

Line 55:

  plot = gsn_streamline(wks,uvar(::2,::2),vvar(::2,::2),resources)

Draw the new streamline plot using the new resources you just set.

Line 57:

;----------- Begin third plot -----------------------------------------

Create a streamline plot using a different dataset, changing the thickness and color of the streamlines, changing the length of the streamline arrows, and changing the spacing between streamlines.

Lines 59-62:

  getvalues plot
    "stArrowLengthF"    : arrowlength
    "stMinLineSpacingF" : spacing
  end getvalues

As noted in the first example, resources have default values. Some resources are set dynamically according to any number of things, like the size of the plot, the minimum and maximum values of the data, etc. There may be cases where you need to retrieve the value of one or more of these dynamic resources so you can change them or use their values to change other resources. The resources stArrowLengthF and stMinLineSpacingF are two examples of dynamic resources whose values you can retrieve so you can change them.

Resource values in NCL are retrieved using the getvalues statement. The getvalues statement is like a block statement, since it has a begin statement (getvalues plotid) and an end statement (end getvalues). Inside the block is a list of all the resources whose values you want to retrieve (each resource should be in double quotes), followed by a colon, followed by a variable to put the value in. The variable plotid is the one returned from a call to one of the gsn_* plotting routines, and the one that the resources apply to (in other words, if you had also created a vector plot with a call to gsn_vector, you couldn't retrieve Streamline plot resources from this plot). In this case, the plotting routine is gsn_streamline and the variable name is plot. More information on getvalues can be found in the "NCL Statements" section of the NCL Reference Manual.

Lines 64-67:

  resources@stMinLineSpacingF = spacing * 4.0
  resources@stArrowLengthF    = arrowlength * 2.0
  resources@stLineColor       = "white"
  resources@stLineThicknessF  = 1.5

Now that the values for stMinLineSpacingF and stArrowLengthF have been retrieved and stored to the variables spacing and arrowlength respectively, you can use these values to change these resources. The spacing between streamlines is increased by a factor of four, and the arrow length of the streamlines is doubled. The last two resources change the line color back to white and increase the line thickness by 50 percent (the default is 1).

Lines 69-70:

  uvar = grb_file->U_GRD_6_GPML(0,:,:)
  vvar = grb_file->V_GRD_6_GPML(0,:,:)

Retrieve the first level of variables U_GRD_6_GPML and V_GRD_6_GPML from the GRIB file.

Lines 72-76:

  if(isatt(uvar,"units"))
    resources@tiMainString = "GRD_6_GPML (u,v " + uvar@units + ")"
  else
    resources@tiMainString = "GRD_6_GPML"
  end if

Change the main title to reflect the new data that you are using to create a streamline plot.

Line 78:

  plot = gsn_streamline(wks,uvar(::2,::2),vvar(::2,::2),resources)

Draw the new streamline plot, using the new resources you just set.

Lines 80-82:

  delete(uvar)
  delete(vvar)
  delete(plot)

Clean up by removing some variables that you don't need anymore.

Line 84:

;----------- Write GRIB data to netCDF file --------------------------

Write GRIB variables that start with "PRES" to a netCDF file. Writing files in the way that is demonstrated here is easy to understand, but it can be very inefficient for large files containing lots of variables and attributes. For information on how to efficiently write large GRIB files to netCDF, see the "Efficiency issues" section in the "Basics" chapter.

Lines 86-88:

  cdf_filename = "ced1.lf00.t00z.eta.nc"
  system("/bin/rm -f " + cdf_filename)
  cdf_file = addfile(cdf_filename,"c")

Open a netCDF file to write the GRIB data to. First remove the file if it already exists, and then use addfile to create it. NCL knows to create a netCDF file by the ".nc" suffix.

Lines 90-91:

  cdf_file@title = "data from a GRIB file"  
  cdf_file@date = systemfunc("date")

Add global attributes to the netCDF file called "title" and "date". The function systemfunc executes a system command and returns the output in a string. It is used here to put a date stamp in the file.

Line 93:

  do i = 0,dimsizes(names)-1

Do loops in NCL are similar to those in Fortran. They have an integer start value, an integer end value, an optional integer stride value, and must be terminated with end do (enddo does not work). For more information on do loops, see the "NCL Statements" section of NCL Reference Manual.

In the do loop above, you are looping through the names of the variables that you retrieved from the GRIB file. Remember, dimsizes returns the dimension sizes of an array, and since names is a 1-dimensional array, dimsizes returns its length as a scalar integer.

Line 94:

    names_char = stringtochar(names(i))

To do any kind of parsing on an NCL string variable, you must first convert it to a character array (strings in NCL are considered a single entity). You can use the stringtochar function to accomplish this.

Lines 95-98:

    if(names_char(0:3).eq."PRES") then
      print(names(i))
      cdf_file->$names(i)$ = grb_file->$names(i)$
    end if

Now that you have a character array, you can access each character as a separate array element. In the four lines above, then, the first four characters of names_char are compared with the string "PRES". If they are equal, then this variable gets printed out and written to the netCDF file that you just opened.

If you don't know the names of the variables in a file that has been opened with addfile, the syntax cdf_file->$names(i)$ allows you to reference file variables using another variable that contains the file variable names as strings.

For example, if names(0) is equal to "PRES_101_SFC", and "PRES_101_SFC" is a variable in the GRIB file, then the following two lines are equivalent:

  cdf_file->$names(0)$ = grb_file->$names(0)$
  cdf_file->PRES_101_SFC = grb_file->PRES_101_SFC

Line 99:

    delete(names_char)

Since you're going to use names_char again, and it most likely will be a different size the next time around in the do loop, delete it first.

Line 100:

  end do

End the do loop. Just like with end if, there must be a space between end and do since NCL won't recognize enddo.

Line 101:

end

End the NCL script.