Log In | Users | Register
Edit | Attach | New | Raw | Diff | Print | | Tools
You are here: Data » DocTools » ToolsVis » VisNCL » VisNCLFAQ


This page contains a collection of the most frequently asked questions about the COSMO library and NCL in general. If you don't find your answer here, maybe you will find it on the official FAQ of the NCL software. Please feel free to add your own FAQ entries!

How to put a copyright note in the lower left corner of a plot? NEW

This code snippet will add the plot production date and "(C) MeteoSwiss?" in the lower left corner of the plot

resc                  = True
copyright             = systemfunc("date -u") 
resc@txFontHeightF    = 0.01
resc@txFontColor      = "blue"
gsn_text_ndc(wks,copyright+" / ~F35~c~F21~ COSMO.MeteoSwiss.MeteoSwiss",.15,.015,resc)

How to use accents and umlauts? NEW

Supposed your used font is F21(helvetica) (see font table), the following code does the trick:

 = o~H-14F35~H~H3F21~
 = a~H-14F35~H~H3F21~
 = u~H-14F35~H~H3F21~

This code uses the horizontal colon from font F35 and shifts it over the preceding character. The same holds for french accents, except that ~H~ has to be replaced by ~A~, ~B~ or ~C~, depending on the accent you want to place. (Make sure you set ~ as TextFuncCode in your .hluresfile)

How can I read an L2E file? NEW

The COSMO library facilitates the reading of L2E files. An example script would look like this...

load "$NCL_COSMOLIB/src/l2e.ncl"
  fname = "OBS-3-ext.pol"      ; name of L2E file to read
  nprof = 1                    ; number of profile in file to extract
  prof = l2e_read(fname,nprof) ; read data from L2E file
  printVarSummary(prof)        ; print data

How can I include a function from a Fortran77/90 code?

Including a function which you have already programmed in Fortran77/90 in NCL is very easy. More details on how to go about this are given on the page describing the WRAPIT command. Let's first start with a Fortran 77 example.

Fortran 77

If you have a subroutine to compute quadratic polynomial values in a file ex01.f, you will have to add statements (NCLFORSTART, NCLEND) around the subroutine declaration. Here's an example...

      subroutine cquad (a, b, c, nq, x, quad)
      real x(nq), quad(nq)
      do 10 i=1,nq
        quad(i) = a*x(i)**2 + b*x(i) + c
   10 continue


Now you will have to compile this code and pack it into a shared object library that can be used by NCL. In order to do that on buin, you will have to issue the following commands...

> module switch PrgEnv-pgi PrgEnv-gnu
> WRAPIT -gf ex01.f

If everything goes well, you will have a ex01.so file in the same directory. Now the function can be accessed very easily from NCL. An example script ex01.ncl could be of the following form...

external EX01 "./ex01.so"
  nump = 3
  x = (/ -1., 0.0, 1.0 /)
  qval = new(nump,float)              
  EX01::cquad(-1., 2., 3., nump, x, qval)
  print("Polynomial value = " + qval)

Fortran 90

The procedure for a Fortran 90 subroutine is somewhat different. Say you have a subroutine to compute quadratic polynomial values in a file ex01.f90. Here's an example...

subroutine cquad(a,b,c,nq,x,quad)
  implicit none
  integer, intent(in)  :: nq
  real,    intent(in)  :: a,b,c,x(nq)
  real,    intent(out) :: quad(nq)
  integer              :: i

  quad = a*x**2+b*x+c
end subroutine cquad

Now, you will have to create an separate file ex01.stub containing the following declarations...

      subroutine cquad(a,b,c,nq,x,quad)
      real a,b,c
      integer nq
      dimension x(nq),quad(nq)

Now you will have to compile this code and pack it into a shared object library that can be used by NCL. In order to do that on buin, you will have to issue the following commands...

> module switch PrgEnv-pgi PrgEnv-gnu
> WRAPIT -gf ex01.stub ex01.f90

If everything goes well, you will have a ex01.so file in the same directory. The usage of the function is now exactly the same as in the Fortran 77 case above.

How can I make a panel plot?

The basic rules for making panel plots are the following:

  1. Simply make a series of plots and panel them using the gsn_panel procedure (see documentation)
  2. Be sure to add the resource gsnMaximize = False to all your jmb_map and jmb_contour calls so that all plots have the same size.
  3. Add the jmbPanel = True resource to all calls to jmb_overlays so that the plots are not drawn and have the same size.
  4. Make sure that the base plots (the first plot in the call to jmb_overlays) correspond to different plots. For example, if your first plot is a map from jmb_map, you will have to draw a separate map for each panel.

Using these simple rules you can create nice looking panel plots. Here's an example...



  ; open files
  f = addfile("lm_coarse/lfff00000000.grb", "r")
  c = addfile("lm_coarse/lfff00000000c.grb", "r")

  ; read data
  u10 = jmb_getvar(f, "U_10M")
  v10 = jmb_getvar(f, "V_10M")
  jmb_getgrid(f, c, u10)
  jmb_getgrid(f, c, v10)

  ; convert units
  u10 = u10 * 3.6
  u10@units = "km/h"
  v10 = v10 * 3.6
  v10@units = "km/h"

  ; open graphic port
  wks = gsn_open_wks("pdf", "test")
  rct = jmb_set_ct(wks, "uv_17lev", False)

  ; plot map
  r = True
  r@gsnMaximize = False
  r@mpProjection = "Stereographic"
  mp1 = jmb_map(wks, u10, r)
  r@tmYLOn = False
  mp2 = jmb_map(wks, v10, r)

  ; contour winds
  r = rct
  r@gsnMaximize = False
  r@lbLabelBarOn = False
  r@jmbAnnotationsOn = False
  cn1 = jmb_contour(wks, u10, r)
  cn2 = jmb_contour(wks, v10, r)

  ; make plots
  r = True
  r@jmbPanel = True
  pl1 = jmb_overlays(wks, (/mp1,cn1/), r)
  pl2 = jmb_overlays(wks, (/mp2,cn2/), r)

  ; make paneling
  r = True
  r@gsnPanelLabelBar = True
  r@txString = "Wind Components at 10m ("+u10@units+")"
  r@gsnPanelFigureStrings = (/"U","V"/)
  r@gsnPanelFigureSTringsFontHeightF = 0.1
  r@amJust = "BottomLeft"
  r@amOrthogonalPosF = -0.43
  r@amParallelPosF = -0.48
  gsn_panel(wks, (/pl1,pl2/), (/1,2/), r)

How can I get syntax highlighting for NCL scripts?

Editor enhancements including syntax highlighting is available for the most popular editors. Check out the corresponding webpage for Emacs, JED, nedit and VIM.

Comment: If you use nedit and syntax highlighting does not highlight "comments" perform the following steps: open nedit, select Preferences > Default Settings > Syntax Highlighting > Recognition Patterns, then select Language Mode "NCL", select Pattern "Comment", tick Matching "Highlight between starting and ending REs", add "$" to Ending Regular Expression, press Apply.

How can I loop over a complete COSMO run?

Here's a simple example of how you can loop over hourly output from a COSMO model run.

dir = "lm_coarse"
files = systemfunc("cd "+dir+"; ls lfff[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]")
nfiles = dimsizes(files)
do i = 0, nfiles-1
  print("Reading file "+dir+"/"+files(i))
  ffile = addfile(dir+"/"+files(i)+".grb","r")

  ; ... do some more processing here

end do

Note that variables in NCL are, once they have been assigned or declared, not overwritten by a consecutive assignement. Thus, in a loop you typicall have to delete variables at the end of the loop.

How can I set the paper size to A4?

If you are using the gsn routines, simply add the following lines to your resources...

a4_height    = 29.70
a4_width     = 16.00
cm_per_inch  =  2.54
resm@gsnMaximize           = False ; maximize size of plot
resm@gsnPaperWidth         = a4_width/cm_per_inch
resm@gsnPaperHeight        = a4_height/cm_per_inch
resm@gsnPaperMargin        = 0.5/cm_per_inch
;resm@gsnPaperOrientation   = "Portrait"

Note that if using the COSMO Library you have to attach this resource to the lowermost plot. This corresponds to the first plot specified in the call to jmb_overlays.

How can I reorder/permute the dimensions of a field?

Say you have a temperature field named T with the named dimesions

T!0 = "time"
T!1 = "lat"
T!2 = "lon"
you can easily reorder the dimensions of the field using coordinate subscripting
reordered_T = T(lon|:,lat|:,time|:)

How can I draw high-resolution country outlines over a COSMO-2 contour plot?

The COSMO library has a command for plotting geographical information on top of a plot. See the jmb_geography command above.

How can I draw a contour plot using the official MO probability colormap?

Open a workstation (wks) and add the following lines to your plot resources...

resources@gsnSpreadColorStart   = 2
resources@gsnSpreadColorEnd     = 11
resources@gsnSpreadColors       = True              ; spread color map
resources@cnMinLevelValF        = 10                ; set min contour level
resources@cnMaxLevelValF        = 90                ; set max contour level
resources@cnLevelSpacingF       = 10                ; set contour spacing

Define the colormap like...

cmap = (/"(/1,1,1/)", \
         "(/0,0,0/)", \
         "(/.839216,0.890196,0.929412/)", \
         "(/0.709804,0.788235,1/)", \
         "(/0.556863,0.698039,1/)", \
         "(/0.498039,0.588235,1/)", \
         "(/0.666667,0.807843,0.388235/)", \
         "(/0.909804,0.960784,0.619608/)", \
         "(/1,0.980392,0.0784314/)", \
         "(/1,0.819608,0.129412/)", \
         "(/1,0.639216,0.0392157/)", \

Note that the first two colors in a colormap always correspond to the foreground and background colors and are ignored when making a plot. Colors you use have to be part of the defined colormap, if not, NCL takes the color of the colormap that is closest to your choice. The resulting plot could look like that...


How can I mask/replace values in a multidimensional array?

NCL has a powerful function call where (see NCL Documentation for details).

For example, to limit values of a arbitrary n-dimensional array to 90 you can do the following

x = where(x.gt.90.0, 90.0, x)
You can use more conditions and combine them using logical expressions, for example
x = where(x .ge. 90.0 .and. ismissing(x), 90.0, x)

A typical application would be to mask out sea points from a COSMO field. This could be accomplished using the following example script

f = addfile("lm_fine/lfff00000000c.grb", "r")
h = jmb_getvar(f, "HSURF")
frl = jmb_getvar(f, "FR_LAND")
h = where(frl .lt. 0.5, h@_FillValue, h)

How can I flag certain values as missing values?

Use the mask function of NCL. For example if you want to flag all rh values at 850hPa which lie "below" topography, use the following:

rh850 = mask(rh850,fi850.gt.hsurf,1)

where fi850 is the geopotential divided by g and hsurf is the topography height. Alternatively you can use

rh850 = mask(rh850,ps.gt.850,1)

if you have the surface pressure available.

Which NCL map transformation corresponds to the COSMO grid transformation?

Suppose you want to draw COSMO model output using one of the NCL _map functions but you would like your data to appear "untransformed", i.e. in a rectangular, Cartesian domain. In order to achieve this, you will have to choose an NCL map transformation which exactly matches the one of COSMO. This can be easily done using the CylindricalEquidistant projection (which is actually the default).

lat = ...         ; 2d-array of geographical latitudes
lon = ...         ; 2d-array of geographical longitudes
pollat = 43.0     ; latitude of rotated pole
pollon = -170.0   ; longitude of rotated pole
r = True
r@cnFillOn = True
r@cnLinesOn = False
r@gsnSpreadColors = True
r@mpDataBaseVersion = "MediumRes"
r@mpProjection = "CylindricalEquidistant"
r@mpCenterLonF = 180.0+pollon
r@mpCenterLatF = 90.0-pollat
r@mpLimitMode = "Corners"
r@mpLeftCornerLatF = lat(0,0)
r@mpLeftCornerLonF = lon(0,0)
r@mpRightCornerLatF = lat(s(0)-1,s(1)-1)
r@mpRightCornerLonF = lon(s(0)-1,s(1)-1)
plt = gsn_csm_contour_map(wks,f,r)

Note that the COSMO library automatically chooses this transformation if you supply geo-referenced data to the jmb_map or jmb_contour plot commands.

How can I switch off paging for the NCL print command?

Simply add the following line to your .cshrc file.

setenv PAGER "more"

Make sure that these variables are set nowhere later on in your .cshrc script.

-- OliverFuhrer - 28 Jan 2009

ex_panel.png (122.56K)
version 2 uploaded by OliverFuhrer on 05 Feb 2009 - 17:15 ...
version 1 uploaded by OliverFuhrer on 05 Feb 2009 - 17:10
version 1 uploaded by OliverFuhrer on 27 Jan 2009 - 10:04
This site is managed by the Center for Climate Systems Modeling (C2SM).
This site is powered by FoswikiCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors. Ideas, requests, problems? Send feedback!