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

Basic NCL examples

Some basic examples for using NCL that have proven helpful.

Passing arguments to NCL

An elegant and secure way to pass arguments to NCL is to set environment variables and then invoke NCL which reads the variables using the getenv function.

runncl.sh is a general shell script to pass arguments to NCL.

#!/bin/bash
#
# Universal wrapper script for ncl. 
# Pass arguments from the command line to environment variables
#
# version 0.1, Thierry Corti, C2SM ETH Zurich
# 

E_BADARGS=65

if [ ! -n "$1" ]
then
  echo "Usage: `basename $0` script.ncl argument1 argument2 etc."
  exit $E_BADARGS
fi  

# save number of arguments to environment variable NCL_N_ARG
export NCL_N_ARGS=$#

# save command line arguments to environment variable NCL_ARG_#
for ((index=1; index<=$#; index++))
do
  eval export NCL_ARG_$index=\$$index
done   

# run ncl
ncl $1

Usage:

   ./runncl.sh script.ncl argument1 "string argument"

Sample application:

#!/bin/bash
#
# Use "runncl" shell script to pass different arguments to a NCL script

for YEAR in {2000..2005}
do
  echo "---------- starting $YEAR"
  ./runncl.sh runwithargs.ncl $YEAR
done  

Same script, but now runncl.sh runs in parallel !

#!/bin/bash

for YEAR in {2000..2005}
do
  echo "---------- starting $YEAR"
  # start runncl.sh and send it to the background
  ./runncl.sh runwithargs.ncl $YEAR &
done  

In NCL, the arguments can be retrieved using "getenv", optionally accompanied by a check for the right number of arguments. Here is an example:

;================================================;
;  runwithargs.ncl
;================================================;
;
; Concepts illustrated:
;   - Demonstrate usage of runncl.sh shell script
;   - Determine number of arguments passed to the script
;   - Retrieve arguments
;

begin
  ; set number of arguments (optional)
  if(ismissing(getenv("NCL_N_ARGS"))) then 
     NCL_N_ARGS = 0  
  else
     NCL_N_ARGS = getenv("NCL_N_ARGS")
  end if

  ; test for correct number of arguments (optional)
  if (NCL_N_ARGS .lt. 2) then
    print("This script needs at least 2 command line arguments. Exit.")
    print("Usage: ./runncl.sh runwithargs.ncl $YEAR")
    exit()
  end if

  a1 = getenv("NCL_ARG_1")
  print("First Argument (script name) ="+a1)

  a2 = getenv("NCL_ARG_2")
  print("Second Argument (YEAR) ="+a2)


  ;
  ; Do what you want here
  ;


end

Map with country borders

By default, NCL country borders are pretty outdated. This can be solved by choosing the more recent, high resolution boundary data set, as now documented in the maps example 2 on the official NCL webpage

  res@mpDataBaseVersion = "Ncarg4_1"            ; choose more recent database
  res@mpDataSetName = "Earth..4"                ; high resolution 

More than country borders - Natural Earth

Natural Earth is a public domain map dataset with cultural and physical features. The script below demonstrates how to use this data in NCL.

naturalearth.png

;----------------------------------------------------------------------
; This example shows how to read geographic data 
; from Natural Earth shapefiles
; and plot them as polylines and polygons. 
;----------------------------------------------------------------------
; This particular example plots data for Switzerland.
;----------------------------------------------------------------------
; Download the shapefiles from http://www.naturalearthdata.com/
; Unzip to a directory 
;----------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
  NaturalEarthDir = "Data/"

;--- Open workstation and define colormap.
  wks_type = "newpdf"        ; use "newpdf" instead of "pdf" for smaller files
  wks = gsn_open_wks(wks_type,"naturalearth")
  gsn_define_colormap(wks,(/"white","black","Beige","LightBlue","tan4"/))

;---Set some resources for the map.
  res                 = True
  res@mpFillOn        = False
  res@gsnTickMarksOn  = False
  res@gsnDraw         = False
  res@gsnFrame        = False
  res@gsnMaximize     = True                ; blow up plot

;---Zoom in on Switzerland
  res@mpLimitMode     = "LatLon"
  res@mpMinLatF       = 44.5
  res@mpMaxLatF       = 51.0
  res@mpMinLonF       = 5.5
  res@mpMaxLonF       = 10.5

  lnres               = True       ; resources for polylines
  plres               = True       ; resources for polygons

  map = gsn_csm_map(wks,res)

;---Shapefiles 
  fnames = NaturalEarthDir \
          + (/"10m-admin-0-countries" \
             ,"10m_rivers_lake_centerlines","10m_rivers_europe" \
             ,"10m_lakes","10m_lakes_europe"   /) + ".shp"
  linecolors = (/"tan4","SkyBlue","SkyBlue","SkyBlue","SkyBlue"/)   
  fillcolors = (/"Beige","","","SkyBlue","SkyBlue"/)   
  thick = (/6,3,3,0,0/)

;---Loop through files that we want to read geographic information from.
  prims = True
  lines = True
  do n=0,dimsizes(fnames)-1
;---Open the shapefile.
    f = addfile(fnames(n),"r")

;---Read data off the shapefile
    segments = f->segments
    geometry = f->geometry
    segsDims = dimsizes(segments)
    geomDims = dimsizes(geometry)

;---Read global attributes  
    geom_segIndex = f@geom_segIndex
    geom_numSegs  = f@geom_numSegs
    segs_xyzIndex = f@segs_xyzIndex
    segs_numPnts  = f@segs_numPnts
    geometry_type = f@geometry_type

    numFeatures   = geomDims(0)

    lon = f->x
    lat = f->y
 
    ; put if statement outside the loop
    if (geometry_type.eq."polygon") then 
      plres@gsFillColor = fillcolors(n)
      plres@gsEdgesOn   = True       ; draw border around polygons      
      plres@gsEdgeColor = linecolors(n) 
      plres@gsEdgeThicknessF = thick(n)  

      ;---Section to draw polygons on map.
      do i=0, numFeatures-1  
        startSegment = geometry(i, geom_segIndex)
        numSegments  = geometry(i, geom_numSegs)
        do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT = startPT + segments(seg, segs_numPnts) - 1
          ;---This call adds the polygon.
          dumstr = unique_string("lines")
          lines@$dumstr$ = gsn_add_polygon(wks, map , lon(startPT:endPT), lat(startPT:endPT), plres)
        end do
      end do
    else
      lnres@gsLineThicknessF = thick(n)  
      lnres@gsLineColor = linecolors(n)
      ;---Section to draw polylines on map.
      do i=0, numFeatures-1  
        startSegment = geometry(i, geom_segIndex)
        numSegments  = geometry(i, geom_numSegs)
        do seg=startSegment, startSegment+numSegments-1
          startPT = segments(seg, segs_xyzIndex)
          endPT   = startPT + segments(seg, segs_numPnts) - 1
          ;---This call adds the line segment.
          dumstr = unique_string("primitive")
          prims@$dumstr$ = gsn_add_polyline(wks, map, lon(startPT:endPT), lat(startPT:endPT), lnres)
        end do
      end do
    end if

;---Clean up before we read in same variables again.
    delete(lat)
    delete(lon)
    delete(segments)
    delete(geometry)
    delete(segsDims)
    delete(geomDims)
  end do

;---We're done adding primitives, draw everything and advance frame.
  draw(map)
  frame(wks)
end

Reading Arc/Info ASCII Grid Files

NCL is rather weak when it comes to reading ASCII formats. Your best bet is to preformat your data into a very simple, regular ASCII file and take a look at some examples on the official website. Here is a function to read lat/lon Grid files using the asciiread command.

; ***********************************************************************
; Name: readAGR
; Description: Read an ESRI grid formated file
; For a data format description, see http://en.wikipedia.org/wiki/ESRI_grid
;
; NOTE: The function silently assumes a regular latitude/longitude grid and 
; sets the metadata accordingly. It is the user's responsibility to transform
; the grid or data if necessary.

undef("readAGR")

function readAGR(filename)
local raw, data, ncols, nrows, xllcorner, yllcorner, cellsize, NODATA_value, lon, lat 
begin

  ; read in all numbers 
  raw =  asciiread(filename, -1, "float")
 
  ; first 6 numbers are header
  ncols = floattoint(raw(0)) ; number of columns
  nrows = floattoint(raw(1)) ; number of rows
  xllcorner    = raw(2) ; x coord of lower left corner
  yllcorner    = raw(3) ; y coord of lower left corner
  cellsize     = raw(4) ; cellsize
  NODATA_value = raw(5) ; Fill value

  ; reshape the rest of the data
  data = onedtond(raw(6::),(/nrows,ncols/))
 
  ; create longitude
  lon = xllcorner + (ispan(1,ncols,1) - 0.5) * cellsize
  lon@units = "degrees_east"

  ; create latitude, reverse to fit data
  lat = yllcorner + (ispan(1,nrows,1) - 0.5) * cellsize
  lat = lat(::-1)
  lat@units = "degrees_north"

  ; add coorinates to data
  data!0 = "lat"
  data&lat = lat

  data!1 = "lon"
  data&lon = lon
  
  ; attach the remaining metadata
  data@_FillValue = NODATA_value
  data@cellsize = cellsize
  data@xllcorner = xllcorner
  data@yllcorner = yllcorner

  return(data)

end ; of function readAGR

If you have vector data instead of raster data, from Arc/Info, your best bet is probably to use the widely known shapefile format which can be read by NCL (see documention for more information).

Attach
png
maponly_2.png (40.24K)
version 1 uploaded by ThierryCorti on 20 Jan 2011 - 09:30
txt
version 1 uploaded by KatherineOsterried on 08 Aug 2018 - 11:10
png
version 1 uploaded by ThierryCorti on 01 Feb 2011 - 11:53
txt
version 1 uploaded by KatherineOsterried on 08 Aug 2018 - 11:10
txt
version 1 uploaded by KatherineOsterried on 08 Aug 2018 - 11:10
txt
version 1 uploaded by ThierryCorti on 15 Feb 2011 - 15:30
txt
version 1 uploaded by KatherineOsterried on 08 Aug 2018 - 11:10
spacer
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!