6 Writing csh scripts

 6.1 How do I get pixel positions using the cursor?
 6.2 How do I get real world co-ordinate positions using the cursor?
 6.3 How do I overplot contours from one image on to another?
 6.4 How do I use scientific notation in bc?
 6.5 My file has been converted to NDF. How do I access FITS header keywords?
 6.6 How do I create an NDF file from an ASCII file?
 6.7 How to make a simple GUI

An excellent introduction to writing csh scripts can be found in the C-shell Cookbook (SC/4).

6.1 How do I get pixel positions using the cursor?

Seemingly a trivial problem with a simple solution this turns out to be slightly more complicated than you would initially expect. Presuming we need an integer pixel position to give to another application we are using in our script, we can use the following code block to do so.

  # Grab x,y position.
  cursor showpixel=true style="Colour(marker)=3" plot=mark \
        maxpos=1 marker=2 device=xwin frame="PIXEL" > /tmp/cursor_lock
  
  # Wait for cursor input.
  while ( ! -e /tmp/cursor_lock )
    sleep 1
  end
  rm -f /tmp/cursor_lock
  
  # Retrieve the position from the parameter file.
  set pos = ‘parget lastpos cursor‘
  
  # Get the pixel co-ordinates and convert to grid indices.  The
  # exterior NINT replaces the bug/feature -0 result with the desired 0.
  set xpix = ‘calc exp="nint(nint($pos[1]+0.5))" prec=_REAL‘
  set ypix = ‘calc exp="nint(nint($pos[2]+0.5))" prec=_REAL‘

Here we run the cursor application, requesting it to return co-ordinates in the PIXEL Frame, creating a lock file to stop further execution of the script until the user has clicked on the graphics display window. We then delete the lock file and retrieve the cursor position using parget. cursor returns each co-ordinates as a floating-point value, whereas we require an integer pixel index. We therefore use calc to convert from co-ordinates to indices by adding 0.5 and finding the nearest integer.

This functionality has been encapsulated within the script $DATACUBE_DIR/getcurpos.csh.

6.2 How do I get real world co-ordinate positions using the cursor?

A similar, but slightly easier, problem is to grab real world co-ordinates from a cursor click on a graphics display window.

  # Grab one position.
  cursor showpixel=true style="Colour(marker)=3" plot=mark \
        maxpos=1 marker=2 device=xwin > /tmp/cursor_lock
  
  # Wait for cursor input.
  while ( ! -e /tmp/cursor_lock )
    sleep 1
  end
  rm -f /tmp/cursor_lock
  
  # Retrieve WCS co-ordinates.
  set pos = ‘parget lastpos cursor‘
  set x = $pos[1]
  set y = $pos[2]

Here we run the cursor application, again creating a lock file to stop further execution of the script. We then delete the lock file and retrieve the cursor position using parget.

This functionality has been encapsulated within the script $DATACUBE_DIR/getcurpos.csh.

6.3 How do I overplot contours from one image on to another?

The following code will overplot the contours of one dataset on to the greyscale image of another of the same spatial size using the Kappa display and contour tasks.

  # Display greyscale image.
  display in=file1 device=xwin mode=per percentiles=[15,98] \
         axes=yes lut=$KAPPA_DIR/grey_lut margin=!
  
  # plot contours
  contour ndf=file2 device=xwin clear=no mode=equi axes=no \
         ncont=${numcont} pens=’colour=2’ margin=!

This is useful in many circumstances, for instance, plotting the contours of a white-light image over a passband image, or velocity map.

6.4 How do I use scientific notation in bc?

If you want to make use of the bc utility to carry out floating-point calculations in csh scripts you may come across a problem with scientific notation. In many cases applications return scientific notation in the form 3.0E+08, or 3.0E-0.8, while bc requires these numbers to be of the form 3.0*10^08, or 3.0*10^-08. The following example segment of code takes two numbers in the first format and passes them to bc in the correct manner so that it can do some arithmetic with them. bc will return a floating-point number (not in scientific notation).

  # first number
  set num1 = ‘echo ${num1} | sed ’s/E/\\*10\\^/’ | sed ’s/+//’‘
  
  # second number
  set num2 = ‘echo ${num2} | sed ’s/E/\\*10\\^/’ | sed ’s/+//’‘
  
  # answer
  set num3 = ‘echo "scale = 15; ${num1}-${num2}" | bc‘

The scale specifies the number of decimal places.

An alternative to using bc inside your scripts is the Kappa calc command which can evaluate many arbitrary mathematical expressions, as in this extract.

  # calculate velocity
  set delta_lambda = \
      ‘calc exp="’${centre_fit} - ${line_centre}’" prec=_double‘
  
  set velocity = \
      ‘calc exp="’${delta_lambda}/${line_centre}’" prec=_double‘
  
  set velocity = \
      ‘calc exp="’${velocity}*3.0E+08’" prec=_double‘

here we evaluate the Doppler velocity of an emission line using three calls to calc, although the velocity could have been derived in a single expression.

6.5 My file has been converted to NDF. How do I access FITS header keywords?

You may need to access information that was contained in the FITS header to carry out your data analysis, having converted your data to NDF format to make use of the Starlink software the FITS header keywords are still accessible as part of an NDF extension. Kappa has a number of tools specifically written to handle FITS keywords (see SUN/95 for details).

The recommended way to find the value of a FITS header keyword is by using the fitsval application. For instance, you could obtain the value of the FITS keyword CRVAL3 (see Section 4.3) as follows.

  % fitsval ifu_file CRVAL3
  7302.2918645
  %

This could be used in a script to take appropriate action along with other Kappa FITS manipulation tools.

  # Get the input file name.
  echo -n "NDF input file: "
  set infile = $<
  
  # Check to see if the CTYPE3 keyword exists.
  set status = ‘fitsexist ${infile} CTYPE3‘
  
  # If the keyword exists...
  if ( ${status} == "TRUE" ) then
  
    # get the value of CTYPE3.
    set ctype3 = ‘fitsval ${infile} CTYPE3‘
  
    # Warn the user if its value is not LAMBDA.
    if ( ${ctype3} != "LAMBDA" ) then
       echo "Warning: Axis 3 not type LAMBDA"
    endif
  
  # The keyword does not exist.
  else
  
    # Warn the user that the keyword is missing.
    echo "Warning: Axis 3 type keyword missing"
  
  endif

Here we use Kappa fitsexist comtmnd to check that the the keyword we are interested in exists, then the fitsval command to query its value and act on the information.

6.6 How do I create an NDF file from an ASCII file?

If it possible to build an NDF file from a flat ASCII text file using the Convert application ascii2ndf and then use Kappa and Datacube applications to modify the structure and contents of the NDF extensions until it has the correct specifications.

  # Create a basic NDF from an ASCII file.
  ascii2ndf in=${datfile} out=${tmpfile} shape=[${dims[1]},${dims[2]}] \
           maxlen=1024 type=’_double’
  
  # Set bad pixels to the magic value VAL__BADD.
  setmagic in=${tmpfile} out=${outfile} repval=-9999.99
  rm -f "${tmpfile}.sdf"
  
  # Set the origin of the output file.
  setorigin ndf=${outfile} origin=[${lbnd[1]},${lbnd[2]}]
  
  # Attach the variance array.
  ascii2ndf in=${varfile} comp="Variance" out=${outfile} \
              shape=[${dims[1]},${dims[2]}] maxlen=1024 type=’_double’
  
  # We have a similar shaped NDF from which we want to clone the WCS
  # and AXIS information and attach to our newly created ${outfile}.
  
  # Clone the AXIS information from an similar shaped NDF.
  setaxis ndf=${outfile} like=${likefile}
  
  # Clone the WCS information.  This will be done incorrectly if the
  # AXIS structures does not exits before the WCS extension is cloned
  wcscopy ndf=${outfile} like=${likefile}

Here we take a flat ASCII data file $datfile and create an NDF of dimensions $dims[1] ×$dims[2], with data type _DOUBLE, using ascii2ndf. We then cal upon setmagic to flag all pixels that have the value − 9999.99 in the NDF with the standard bad ‘magic’ value, in this case VAL__BADD. setorigin sets the pixel origin to ($lbnd[1],$lbnd[2]). We then use ascii2ndf again to attach a variance array, from the ASCII flat file $varfile, to our newly created NDF.

For both invocations of ascii2ndf we make use of the MAXLEN parameter. This is the maximum record length (in bytes) of the records within the input text file, the default value being 512. If you attempt to generate an NDF from a file where many of the entries are double-precision numbers, it might be necessary to set this value higher than the default value otherwise some records (lines) may become truncated leading to ‘stepping’ effects within your output NDF.

After including our variance array, we attach AXIS information using the setaxis application, and incorporate WCS information with wcscopy application, both from Kappa. We clone this information from an NDF whose world co-ordinate information is the same as our newly created NDF. If we wanted to avoid copying the AXIS information (or if the NDF from where we were cloning had no WCS information) we could make use of the wcscopy command’s TR parameter to provide a transformation matrix.

Alternatively, we could make use of the Kappa setaxis command to create an AXIS structure within the NDF being cloned obtained from its existing WCS extension via the $likefile parameter. setaxis should be invoked for each axis. Then we copying the AXIS and WCS structures to our new NDF.

  # Create an AXIS structure from a WCS extension
  setaxis ndf=${likefile} mode=wcs comp=Centre dim=1
  setaxis ndf=${likefile} mode=wcs comp=Centre dim=2
  
  # Copy the AXIS structure to our new NDF.
  setaxis ndf=${outfile} like=${likefile}
  
  # Copy the WCS extension to our new NDF
  wcscopy ndf=${outfile} like=${likefile}

The above usage of setaxis is sometimes needed when including non-WCS compliant legacy applications, such as fitgauss, in scripts, as these legacy tasks do recognise the AXIS structure.

6.7 How to make a simple GUI

The Xdialog program is designed as a drop-in replacement for the dialog and cdialog programs. It can convert a simple terminal based script into a program with an X-Windows (GUI) interface. The requirements to install XDialog is that you must have the X11 libraries and GTK+ libraries (version 1.2.×) installed on your machine. GTK+ comes installed by default with most recent Linux distributions (being necessary to run the GNOME desktop), however, it can be compiled under both Solaris and Tru64 UNIX, and with the adoption of GNOME as the standard SUN desktop will increasing come installed by default on platforms other than Linux.

pict

Figure 20: An XDialog script based on velmap; here it asks for an input file.

pict

Figure 21: The same XDialog script later in the run; here we are prompted for the central wavelength of the line.

Most shell scripts are easily converted to use Xdialog, for example we quickly modified the velmap script to use it, see Figures 20 and 21. More information about XDialog can be found at http://xdialog.dyns.net/.