gipsypython recipes

Python recipes for GIPSY

Introduction

This page is dedicated to the programmer who likes to examine practical examples to extract whatever is useful to use in his/her own application. The Python binding for GIPSY allows a programmer to write both simple and complicated GIPSY programs with the ease of writing a COLA script. COLA, COntrol LAnguage for GIPSY has been used as a rapid development environment to write utilities and small test programs for starting a sequence of GIPSY applications. This Python binding is much more powerful than COLA.

On the other hand, you can write real GIPSY programs in Python as you could with Fortran, Sheltran or C with the advantage of all the available modules and packages that you can use in your own Python environment.

First step

Assuming you installed the Python binding for GIPSY, copy the next example to a file. Save it to disk using an arbitrary file name and make sure that the script is executable (chmod u+x myprogram.py).
Start the program using one of these methods:

Start the GIPSY user interface (type: gipsy) and start your new task in the usual way by typing its name then press enter.
For jobs that don't need keyword interaction (e.g. batch jobs) you can also start your task via ngipsy at the Unix prompt. Output will be appended to the file GIPSY.LOG in the directory where you started your application. Example:
ngipsy hello1.py

Note that Python code after the finis() call will not be executed because this routine includes an implicit exit. In the next piece of code we created a familiar program...

hello1.py

#!/usr/bin/env python
from gipsy import *

init()
anyout("GIPSY")
finis()

Note that you can avoid problems with name spaces if you import GIPSY functionality in a different way. To illustrate that we present the alternative:

hello2.py

#!/usr/bin/env python
import gipsy

gipsy.init()
gipsy.anyout("GIPSY")
gipsy.finis()

How to get help

If you need more information about GIPSY's Python binding, you can consult the documentation at http://www.astro.rug.nl/~gipsy/pguide/pygipsy.html.
In this documentation you will encounter references to dc2 documents. These are documents which describe GIPSY functions for programmers. They can be found in your GIPSY tree ($gip_sub).
Questions can also be mailed to: gipsy@astro.rug.nl or helpdesk@astro.rug.nl

Starting and controlling GIPSY tasks

Starting existing GIPSY tasks

The primary goal of introducing Python for GIPSY was to offer a modern alternative for the batch language COLA. The most important functionality of COLA is the possibility to start other GIPSY tasks with predefined keywords. Therefore, The first function we made available in our Python environment was 'xeq' which does the same thing. The next example is an application that starts another task (DISK) which is called with a pre-defined keyword (INSET=;). If you are familiar with GIPSY tasks then you will know that a semi colon in a keyword stands for 'enter' and is used to terminate a keyword loop. So with this, rather useless, program you start GIPSY task DISK which will not prompt you with a keyword INSET=

xeq.py

#!/usr/bin/env python
from gipsy import *

init()
xeq("DISK INSET=;")

finis()

Two more elaborate examples of batch processing are given below. Both run a simple task which does nothing more than 'thinking' about a given set during some random time. It is not a perfect task: it has a failure rate of about 10%. The code is given below.

think

#!/usr/bin/env python

from gipsy import *
import time, random

init()
inset = usertext('inset=')
delay = 2.0*random.random()+1.0
anyout('Start thinking about set "%s".' % inset)
time.sleep(delay)                                # simulate processing
if random.random()<0.1:                          # simulate 10% failure
   error(4, 'couldn\'t think about set "%s" anymore' % inset)
anyout('%s thought about set "%s" during %g seconds.'
        % (myname(), inset, delay))
finis()

The first example obtains an input set name and a list of subset numbers from the user. Then it runs the task think for all given subsets of the set. When think fails, an exception is raised which is caught by the script, which notes the subset number and eventually reports for which subsets the problem occurred.

batch1

#!/usr/bin/env python
from gipsy import *

init()

inset  = usertext("INSET=", "Input set name")
work   = userint("SUBSETS=", "Subset numbers", nmax=1000)

failed = []                           # list of subsets which failed

for subset in work:
   try:
      xeq('think inset=%s %d' % (inset, subset)) # start and wait for completion
   except:                                       # in case of failure ...
      failed.append(subset)                      # ... register subset number

if failed:
   anyout("Subset(s) failed: %s" % repr(failed))

finis()

Now suppose that you have a multi-core or multiple CPU machine and that the operations on subsets are independent of each other. Then you may want to run multiple parallel instances of think. GIPSY's event mechanism makes this possible. The following example shows how.
Again an input set name and a list of subset numbers are obtained from the user. Also how many parallel tasks should be run. The first tasks are started in a loop by calling the function runjob(). For every task started, runjob() schedules a keyword event which will be handled by the callback function jobdone() which will note possible failures and call runjob() to start the next instance of think if there is still work to be done. When the first tasks have been started, the script calls mainloop(). Mainloop() then will call jobdone() whenever a task has exited. It will return when all subsets have been processed.

batch2

#!/usr/bin/env python
from gipsy import *

def runjob():
   if work:
      subset  = work.pop(0)             # get next subset to be processed
      job     = inactive.pop(0)         # obtain free job identifier ...
      active.append(job)                # ... and put it in the active list
      jobname = 'job%d' % job           # job name (and job keyword)
      KeyCallback(jobdone, jobname, subset=subset, job=job)
                                        # schedule task completion callback
      xeq('%s(think) inset=%s %d' % (jobname, inset, subset), jobname)
                                        # start task and continue

def jobdone(cb):
    active.remove(cb.job)               # move job identifier from active ...
    inactive.append(cb.job)             # ... to inactive list
    fate = userint(cb.key)              # read task's exit status ...
    if fate<0:                          # ... and in case of failure ...
       failed.append(cb.subset)         # ... register subset number
    runjob()                            # start next job
    cb.deschedule()                     # deschedule this job's callback

init()

inset  = usertext("INSET=", "Input set name")
work   = userint("SUBSETS=", "Subset numbers", nmax=1000)
npar   = userint("NPAR=", "Number of parallel jobs")

failed   = []                           # list of subsets which failed
active   = []                           # job identifiers in use
inactive = range(npar)                  # available job identifiers

while inactive:                         # start the first 'npar' jobs
   runjob()

mainloop()                              # event loop: returns when there ...
                                        # ... are no scheduled callbacks left
if failed:
   anyout("Subset(s) failed: %s" % repr(failed))

finis()

User interaction

Your application can interact with the user with GIPSY user input routines (see User input documentation). The input syntax is similar to the input routines in C. The values for 'default' are

0: no default
1: default
2: hidden
4: exact number

The general form is:

myvalues = usertext(keyword="INPUT=", message="Enter text string", default=0, defval=None, nmax=None)

Read a text string from the user. 'nmax' is the maximum number of characters to be returned (default 80 or the length of 'defval')

myvalues = userint(keyword="INPUT=", message=None, default=0, defval=0, nmax=1)
myvalues = userreal(keyword="INPUT=", message=None, default=0, defval=0.0, nmax=1)
myvalues = userdble(keyword="INPUT=", message=None, default=0, defval=0.0, nmax=1)

Read one or more integer/real/double precision numbers from the user. No more that 'nmax' numbers can be given. If no input is given, the value of argument 'defval' will be returned. If nmax=1, a single value will be returned and if nmax>1 a list containing the value(s), even if only one value is given.

myvalues = userlog(keyword="INPUT=", message=None, default=0, defval=0, nmax=1)

Read one or more boolean values (logicals) from the user. No more that 'nmax' values can be given. If no input is given, the value of argument 'defval' will be returned. If nmax=1, a single value will be returned and if nmax>1 a list containing the value(s), even if only one value is given.

myvalues = userchar(keyword="INPUT=", message=None, default=0, defval=None, width=1, nmax=1)
myvalues = usercharu(keyword="INPUT=", message=None, default=0, defval=None, width=1, nmax=1)
myvalues = usercharl(keyword="INPUT=", message=None, default=0, defval=None, width=1, nmax=1)

Read a list of character strings all-case/uppercase/lowercase from the user. No more that 'nmax' strings can be given and each string can be 'width' characters wide. If no input is given, the value of argument 'defval' will be returned. If nmax=1, a single value will be returned and if nmax>1 a list containing the value(s), even if only one value is given.

myvalues = userangle(keyword="INPUT=", message=None, default=0, defval=0.0, nmax=1)

The userangle routine is documented in userangle.dc2

Note that all arguments have a default.

userinput.py

#!/usr/bin/env python
from gipsy import *

init()

# Syntax: x = userreal(keyword="INPUT=", message=None, default=0, defval=0.0, nmax=1)
# example 1: All arguments are default

x = userreal()
anyout( str(x) )


# example 2: Specify everything

x = userreal(keyword="NUMBER2=", message="Give one value: [10.3]", 
             default=1, defval=10.3, nmax=1 )
anyout( str(x) )


# example 3: Specify everything, but in fixed order so that you can omit
# the argument names

x = userreal("NUMBER3=", "Give one value: [10.3]", 1, 10.3, 1 )
anyout( str(x) )


# example 4: Get a list of n numbers. Now 'x' is a list so you can
# index the individual numbers


x = userreal("NUMBER4=", "Give max. 10 values", 1, (1.3), 10 )
anyout( str(x) )
anyout( "Number = %f" % x[0] )


# example 5: Get a list of n numbers. with defaults

x = userreal("NUMBER5=", "Give max. 10 values", 1, (1,2,3,4), 10 )
anyout( str(x) )
anyout( "First Number = %f" % x[0] )

finis()

Note that if you allow for an input of more than one number, and you want to supply defaults, the list of defaults needs parentheses even if you supply only one default!

Next also demonstrate the use of keyword cancelling and rejecting.
Note that the last loop is endless and needs to be aborted with ctrl-C.

cancelreject.py

#!/usr/bin/env python
from gipsy import *

init()

# example 1: Repeat over user input until program agrees with input
# The reject shows a message and cancels the keyword input so it can
# be entered again

ok = False
keyword = "NUMBER1="
while not ok:
   x = userreal(keyword, "Enter a number:", 0, 0, 1)   
   anyout( str(x) )
   if x >= 0 :
      ok = True
   else:
      reject( keyword, "Negative value not allowed" )


# example 2: Repeat over user input. Cancel the keyword.
# Note that we are using the default value for the keyword because
# we did not specify the parameter keyword=

attempt = 1
while True:
   v = userreal(nmax=1, default=1, message='Enter number (attempt %d):'% attempt)
   anyout(str(v))
   cancel()
   attempt += 1

finis()

User input and applying exceptions

Usually you will extend your program by accepting user input with GIPSY keywords. In many cases you will perform some sort of check on this input. Here is an example that uses Python exceptions to control the flow of a program. If a set exists then a set object is created. Else, an exception is raised, an error message is retrieved and printed and the program will be aborted because we supplied a level 4 to the GIPSY error function. If you supply an empty name or non existing set, then the output of this program will be "GDSError on line 10."

except.py

#!/usr/bin/env python
from gipsy import *

init()

name = usertext(message="Set name")

try:
   s = Set(name)
except GipsyError, msg:
   error(4,str(msg))

del s

finis()

Calling some GIPSY tasks with keyword values

A more useful program is a Python script that will convert a fits file to a GIPSY set and then display that file. It asks for the name of your fits file and the name of the GIPSY set. Then it calls the GIPSY task RFITS with keywords to convert the fits file to GDS and it calls the task VIEW to display the result.

fitsviewxeq.py

#!/usr/bin/env python
from gipsy import *

init()

fitsfile = "file000001.mt"             # Set a default name for an fits file
fitsfile = usertext("FITSFILE=",
                    "Enter name of fits file:  [%s]" % fitsfile, defval=fitsfile)

# If the fits file you entered does not exist, then ens this program 
# in the GIPSY way e.g. call the 'finis' function
if not os.path.exists(fitsfile):
   anyout(">>>>>>> WARNING: Cannot find %s <<<<<<<<" % fitsfile)
   finis()


# Get a name for a GIPSY set to store the fits data in gds format
outputset = "gipsyset"
outputset = usertext("GIPSYSET=",
                     "Enter name of output set: [%s]" % outputset, 1, outputset)
try:
   s = Set(outputset)
except:
   pass
else:
   s.delete()
   del s


xeq("RFITS AUTO=Y FITSFILE=%s; INFILES=0; OUTSET=%s" % (fitsfile, outputset) )
xeq("VIEW INSET=%s CLIP=" % outputset)

finis()

GIPSY sets and subsets

How about input of set and subsets?

GIPSY sets and subsets are described in detail in the document Data Structure. but here we will summarize some definitions.

Data in GIPSY is stored in a so called SET. A set consists of two components: an image component and a descriptor component. These components are stored in separate files, but in GIPSY they are always used together.
Within a set subsets can be defined. A SUBSET is a part of a set with a dimensionality not greater than the dimensionality of the set.
A FRAME is a part of a subset.
Associated with the image data there exists a hierarchical system of DESCRIPTOR (or "header") items. A descriptor item is identified by the combination of its name and its subset identification.
AXIS PERMUTATION array. Axes in a set have a fixed order. Input of set and subsets returns a list of axis numbers ordered by subset axes followed by the numbers of the 'repeat' axes'. Repeat axes are axes along which the subsets are repeatedly extracted.

Next we show a picture of a Set AURORA with axes x, y and z named by RA, DEC and FREQ.

In the next example we get the set, subset(s) input with routine usertext and use this input to create the set object. First some input examples:

AURORA FREQ
Return planes in RA and DEC for each grid on the FREQ axis
AURORA FREQ -10:10
Return planes in RA and DEC for grids -10 to 10 on the FREQ axis
AURORA FREQ 2 4 5
Return planes in RA and DEC for grids 2 4 and 5 on the FREQ axis
AURORA RA 10
Return one plane with axis DEC and FREQ at RA grid 10
AURORA DEC 5 6 FREQ 10
Return two lines in RA at grids DEC 5, FREQ 10 and DEC 6, FREQ 10
AURORA RA 5 DEC 5 FREQ 10
Return one pixel at grid RA 5, DEC 5 and FREQ 10

Test this with an arbitrary set with this little example program. If you don't have a set, create one using the code in one of the recipes below. The program shows the list of coordinate words associated with subsets and the axis permutation array. Think of coordinate words as being id's of subsets. The anyout parmater should be checked before using them because if a header does not have an entry for CRVAL, then NULL is returned which cannot be formatted as a floating point number.

setproperties.py

#!/usr/bin/env python
from gipsy import *

init()

while True:
   input = usertext("SET=", "Set, (subsets)")
   cancel("SET=")
   s = Set(input)
   anyout("set opened with: %s" % s.spec)
   anyout("set name       : %s" % s.name)
   anyout("subsets        : %s" % str(s.subsets))
   anyout("axperm         : %s" % str(s.axperm()))
   anyout("Axis in fixed order of set:" )

   for ax in s.axes():
      anyout( "Name = %-10s from %5d to %5d  grid 0 is %10g %-10s" % (
        ax.name, ax.slo, ax.shi, ax.crval, ax.cunit ) )

   anyout("Axis in order of user input:" )
   for a in s.axes( True, 'inside' ):
      anyout( "%s (subset axis)" % a.name )

   for a in s.axes( True, 'outside' ):
      anyout( "%s (repeat axis)" % a.name )

   del s

finis()

As a bonus we add a piece of code illustrating the use of the axis permutation array as an alternative to the previous example. It fits in the program body of the example just above the line with "del s". The code illustrates the use of the axes permutation method in two modes. First we use it to loop over the axis numbers in the subset to get the axes names. Second we loop over the the so called repeat axes to get the name of the axes along which the subsets are repeated. Note that we showed in the previous example that we do not really need the axis permutation array. The code becomes more clear if you use the axes() method.

setpropbonus.py

   n = s.ndims( subset=True )
   anyout( "dimension of the subsets is: %s" % n )

   # use the axis permutation 'k' to get the axis name
   for k in s.axperm('inside'):
      anyout("subset axis name         : %s" % str(s.axname(k)))
   for k in s.axperm('outside'):
      anyout("repeat axis name         : %s" % str(s.axname(k)))

Calculate the mean of subsets using GIPSY input syntax

In the previous example we showed how to obtain a list of coordinate words that represents subsets. Now we want to calculate the mean of the data in these subsets. Looping over subsets can be done with two methods. The method subimage() needs a coordinate word to create a slice of the set data. Its shape is set by the number of axes in the subset and the length of these axes. The other method is called subimages() and is an example of a generator that gets subset data without a coordinate word as a parameter. It is very suitable for looping over subsets.

This way we can loop over all subsets and calculate the mean of the data in these subsets using a method mean() that is available to arrays. The recipe shows that we can write a real analysis application with only a few lines of code.

meansubimages1.py

#!/usr/bin/env python
from gipsy import *

init()

s = Set(usertext("INSET=", "Set,subset(s): "))

for subset in s.subsets:
   anyout( "mean of subset" + " = " + str(s.subimage(subset).mean()) )

finis()

However if you want to access individual subsets, then you can get slices using the subimage() method, which has a coordinate word as input parameter. Now you don't need to loop over you subsets in a fixed order. Here is the example:

meansubimages2.py

# (due to changes in the gipsy module, this code is now
#  equivalent to the previous example's)
#
for cword in s.subsets:
   subsetslice = s.subimage(cword)
   anyout( "mean of subset" + " = " + str(subsetslice.mean()) )

Accessing the properties of a set (II)

In many GIPSY programs one needs to know some properties of the set. In the next example we demonstrate how to access set properties like axis lengths etc. Note that the argument of the range method is '0'. It is a coordinate word that represents the entire set (set level). For accessing header items on set level (=0) you can omit the zero. Header items can be removed from the set using the Python value None.
This script also writes some info in the header of the input set. You have to tell this explicitely to the method that opens the set with the parameter write=True

headerinfo.py

#!/usr/bin/env python
from gipsy import *

init()

s = Set(usertext(message="Setname"), write=True)

dim = s.naxis
anyout( "Dimension of set: %d" % dim )

# Get coordinate words of first data point in set and of the last one
lo, hi = s.range(0)

# Print the axis name and the range in grids for all the axes in the set
for i in range(dim):
   anyout("%s: from %d to %d" % (s.axname(i), s.grid(i, lo), s.grid(i, hi)))

# Print the value of a header item (at set level)
t = s[0, 'INSTRUME']
anyout( str(t) )


s[0, 'OBSERVER'] = 'Me'        # Change the value of a header item.
s['TEMPHEAD'] = None           # Delete item 'TEMPHEAD' on set level (=0)

del s                          # Close the set

finis()

Creating a set from scratch

Never has it been so easy to create a set from scratch than with this Python interface. Create the Set object with the parameter create=True and create the axes with the extend method. Then fill the image that this new set represents with some numbers and as a bonus, extract a plane.

createset.py

#!/usr/bin/env python
from gipsy import *

init()

name = usertext(message="Set name")

s = Set(name, create=True)      # Mode 'create' creates a new set

s.extend("RA",   10, 5)         # Create an axis. Parameters: Name, length, origin
s.extend("DEC",  15, 4)
s.extend("FREQ",  0,  3)

i = s.image                     # Now you can create its data
anyout( str(s.slo) )
anyout( str(s.shi) )

i[:] = range(5)                 # Fill with numbers

x = i[1,:,]                     # Get a slice

anyout(str(x))
anyout(str(x.shape))

del s

finis()

Calculate the mean of subsets in a set, an alternative method

This example is written for a 3d cube and we want a list of the mean values of the data in the subsets without using the GIPSY syntax for the specification of subsets i.e. we extract the data from the cube directly. If we create the set object 's' we can access all its data in the 'image' member which is a multidimensional array with the right shape. Note that in Python the first array index corresponds to the last axis in the GIPSY subset. Array slicing is used to get separate subsets (for a GIPSY RA-DEC-FREQ cube, these are the subsets with RA and DEC axis).

meanwithoutsubsets.py

#!/usr/bin/env python
from gipsy import *

init()                  # Contact Hermes

# Get the name of the GIPSY set by prompting with a default keyword INPUT=
name = usertext(message="Set name")

s = Set(name)           # Create set object s

i = s.image             # Extract the image part of the set object

k = i.shape[0]          # Length of first element = last axis

# Loop over all subsets
for j in range(0,k):
   m = i[j]             # Extract a subset along the last axis. Note that
                        # i[j] is an abbreviation of the array slice i[j,:,:]
   me = m.mean()        # Get the mean of this subset

   anyout( "mean of channel "+str(j)+" = "+str(me) )


del s                   # Release the object

finis()                 # Disconnect Hermes

Writing data to an existing set

Writing data to a set directly is easy. The Set object has an image part with the right dimensions. If we have a three dimensional data cube with dimensions 1024x512x64 then a matrix M has indices M[0..63, 0..511, 0..1023] In Python you can use some shortcuts to address rows, planes or the entire cube. In the next example we use the M[:] notation to set the entire input set to a value entered by the user. Do not forget to tell the set object that it is allowed to change its image data.

writedatatoset.py

#!/usr/bin/env python
from gipsy import *

init()

name = usertext(message="Set name")

s = Set(name, write=True)

i = s.image

val = userreal("value=", "Enter value to fill set")
anyout( str(val) )

# Entire set is set to this value
i[:] = val

del s
finis()

Defining a frame (box input)

We defined a frame as a part of a subset. If you refer to a subset as a slice through your data, then it easy to imagine that it can be useful to have a method that defines a part of this subset. GIPSY users recognize the input of a box. Before you can enter a box, a set object needs to be created, because a box is related to a set. Only then you can enter a box in physical coordinates. Read more about box input in Specifying a frame in the programmers guide. The setbox method of the set object has a string as input parameter. This string contains the definition of the box, in this case entered by the user with usertext.

Before the box method is used, the set has two attributes blo and bhi which contain the coordinates of the entire subset. After the setbox method these values are changed to what was specified in its argument, in this case a string, entered with usertext. As explained before, the subimage method extracts a data slice from the set. The method has an extra argument box=True or box=False. The first is a default. Then any slice (=frame) will have the limits as set by the box input. Otherwise, (with box=False) you get the entire slice. Below the example source, we will describe the code in detail. The program is an endless loop which can be aborted with ctrl-C.

boxtest.py

#!/usr/bin/env python
from gipsy import *

init()

s = Set(usertext("SET=", "Set, subset"), write=True)
anyout("blo= %s, bhi= %s" % (s.blo, s.bhi))

while True:
   b = usertext("BOX=", "Give box", defval='', default=1)
   cancel("BOX=")
   s.setbox(b)
   i = s.subimage(s.subsets[0])
   anyout(str(i.shape))
   anyout("blo= %s, bhi= %s" % (s.blo, s.bhi))
   v = userreal("VALUE=", "Give value")
   i[:] = v

finis()

The program prompts the user for a set and subset. It creates the set object s (allowing to change its data with write=True and prints the tupels that correspond to the grid limits of the entire set. Then it starts a loop. In this loop, the user is prompted for a box. The user enters a box and an array with image data corresponding to the subset is created with method subimage(). Then the user is prompted for a data value and then all pixels in the box are set to this value. You can check the result by starting task SLICEVIEW to display your subset. If you press its reload button after you entered box and value it will display the new result. Be aware of the fact that your set is changed so make a copy of it first if you want to keep the original.

Here are some examples related to box input:

PC projection centre, effectively grid position (0,0,...,0).
AC axis centre. If the length of axis i is NAXISi and the reference pixel is CRPIXi, then the i-th coordinate is given by the expression NAXISi / 2 - CRPIXi.

BOX= -3 -3 4 4 specifies a 2-dimensional frame containing 64 pixels from lower left to upper right

BOX= PC D 8 8 box with size 8x8 pixels with PC as center

BOX= * 3 0 0 * 45 0 0 D 6 arcmin 6 arcmin box centered at 3h0m0s and 45d0m0s with sizes in RA and Dec of 6 minutes of arc. The units of the sizes must be compatible with the axis units in the header of the set.

Writing history to a set header

It is useful to store some text in a set header to describe how a set was changed or used. There are two methods to write such information. In the next examples the whistory() method is once called without parameters. Then it writes the current date and time and the current task parameters to history in the descriptor of the set at top level. If you specify a string, then only that string is stored. The wcomment() method is used if you want to store comments.

historytest.py

s.whistory()
s.whistory('GIPSY (C) 2003')
s.wcomment('Processed on Albirumi')

A small program to copy data to a new set

Lets write a program that reads data from an input GIPSY set, processes this data and writes the result to a new set. As an bonus we want to manipulate the data as function of the data value and the position in the box.

Then we need to solve the a couple of problems:

Create a useful default string for the BOX keyword
Create an output set with the same properties as the input
Create a loop over all the subsets in the input set and write the manipulated data to the corresponding subsets in the output set
Convert the grid positions to positions in the array

An output set in this case is a copy of the input set. The copy() method copies the input set and takes the subset definitions and the box into account. Therefore, in the default of the copy() method it is possible that your output set can be smaller than the input set in any dimension.

It is possible to access image data in a loop using the shape of the array in which the data is stored. For a N x M array the first element is accessed by array[0,0] and the last element by array[M-1,N-1] (i.e. index of 'fastest' axis is last). Now it is possible to loop over all positions in x and y and retreive the pixel value at that point. Using an expression entered by the user, we use Python's eval function to calculate a new pixel value as function of its position and previous pixel value. Test the progam with e.g. EXPR=x+y+data.
Note that the dimension of the input must be a 2-dimensional (sub)set. The program does not check that your input has the right dimensionality.

copyset.py

#!/usr/bin/env python
from gipsy import *

init()
setin = Set(usertext("INSET=", "Set, subset:"))
anyout("flo= %s, fhi= %s" % (setin.slo, setin.shi))

# Prepare to enter a box

boxmes = 'Enter box in '
for k in setin.axperm('inside'):
    boxmes += str(setin.axname(k)) + ' '
boxmes += '   [entire subset]'
b = usertext("BOX=", boxmes, defval='', default=1)
setin.setbox(b)
anyout( "blo=%s bhi=%s" % (setin.blo,setin.bhi) )

outname = usertext("OUTSET=", "Name of output set: ")
setout=setin.copy( outname )

userexpr = usertext("EXPR=", "Enter expression f(x,y,data): ")

# Start to loop over all subsets. You need both the coordinate words 
# of the input set and the output set.

for s in range( len(setin.subsets) ):
   Min  = setin.subimage(setin.subsets[s])
   Mout = setout.subimage(setout.subsets[s])

   for y in range(setin.blo[1],setin.bhi[1]+1):
      iy = y - setin.blo[1]
      for x in range(setin.blo[0],setin.bhi[0]+1):
         ix = x - setin.blo[0]
         data = Min[iy,ix]
         Min[iy,ix] = eval(userexpr)

   Mout[:] = Min                     # Copy subset data to output set
   setout.wminmax(setout.subsets[s]) # Update header items DATAMIN,DATAMAX,NBLANK

setout.whistory()                    # Save history to header
del setin
del setout

finis()

Coordinates

Definitions

In GIPSY we have different coordinates. Here are the definitions.

File coordinates, which are directly related to the position in the image file. In GDS file coordinates are often represented as coordinate words. A COORDINATE WORD is an integer value into which the position of any subset (pixels, lines, planes etc.) can be coded.
GRID COORDINATES, related to the position with respect to the symbolic axes
PHYSICAL or astronomical COORDINATES

Converting grid coordinates to physical coordinates

Many times a GIPSY application needs to know the physical value that correspond to a location in grid coordinates. For example you want to know what the RA and DEC are in a RA-DEC subset at coordinate 0,0. Or you want to know the physical coordinate on the repeat axis/axes for the current subset in a loop. The next code is a recipe for converting a grid position in a subset to a physical value and to convert the grid position of the subset on the repeat axis (or axes) to a physical value.

What are these physical values? In the header of the GIPSY set you find header items that define a physical coordinate system. Item CRPIX indicates which pixel along an axis is grid 0. CRVAL is the physical value that corresponds to this grid 0 in units given by CUNIT. Then you need a step size CDELT which is a size in physical coordinates of one pixel. Finally you need a sky system and a projection system or something similar to do proper transformations. The projection is read from header item CTYPE (e.g. RA-NCP).

Lets introduce two new methods for the transformation between grid coordinates and physical coordinates. First is tophys which converts grids to physical coordinates annd second togrid which transforms physical coordinates to grids. The methods accept a number of parameters: a coordinate tupel, a coordinate word and a mode which can have values 'inside' and 'outside'. This mode determines what is transformed, either the coordinates inside the subset or coordinates outside the subset, i.e. the coordinates along the repeat axis/axes. For the last option you need to specify the grid coordinate also because the information about the repeat axes is stored not completely stored into the coordinate word if the repeat axis needs a matching axis. For example, given our RA-DEC-FREQ cube aurora we can specify subsets as: INSET=aurora dec -10:10 Then the grids along the repeat axis (dec) are extracted from the coordinate word, but it need a matching RA (otherwise a default is taken) to be able to transform to physical coordinates. This all is reflected into the two lines:

   p = s.tophys( x, cword, 'inside' )
   po = s.tophys( x, cword, 'outside' )

The routines are based on the C-function COTRANS but are much more comfortable to use because you don't need to bother about axis permutations as you should have done when you used COTRANS. Basically our next recipe has the same functionality as the GIPSY task COORDS.

In the recipe we present two methods to get the grid on a repeat axis that represents the position of the subset. In C we used to find this position using the axis permutation array and a grid function that converts a coordinate word to a grid. The second method uses the axes() method to get the properties of an axis. One of its attributes is 'blo' which represents the lower limit of a box in a subset or the grid position of a subset for a repeat axis (then also blo == bhi). But the method needs a coordinate to do the right conversion.

cotrans.py

#!/usr/bin/env python
from gipsy import *

init()

input = usertext("SET=", "Set, (subsets)")
s = Set(input)

ndim = s.ndims( subset=True )
x = userdble( "GRIDS=", "Enter position in subset (%d numbers): " % ndim, nmax=ndim )
anyout( "grid = %s" % str(x) )

axnames = ''
for a in s.axes( True, 'inside' ):
   axnames += a.name + '(s) '

for a in s.axes( True, 'outside' ):
   axnames += a.name + '(r) '

anyout( axnames )

for cword in s.subsets:
   p = s.tophys( x, cword, 'inside' )
   po = s.tophys( x, cword, 'outside' )

   #  Get the grids of the position of the subset on the repeat axes
   #  Method 1, use axis permutation array and grid() method
   #  st = ''
   #  for k in s.axperm('outside'):
   #     st += str(s.grid(k,cword))+ ' '

   # Method 2, use the axes() method for this subset. The repeat axes
   # need a coordinate word.
   st = ''
   for a in s.axes( cword, 'outside' ):
       st += str(a.blo) + ' '


   anyout( "%s repeat: (%s)=%s" % (str(p), st, str(po)) )

del s


finis()

Tables, columns and rows

Introduction

A GIPSY table resides in a GDS set. A set can contain as many tables as necessary and every table can contain as many columns as necessary. Columns can be created and deleted at any moment. For programmers details about tables see: GIPSY Table documentation.

There are 4 examples to illustrate the use of the table methods. The first example describes how to create a table.

createtable.py

#!/usr/bin/env python
from gipsy import *

init()

s = Set(usertext("set=",message="Setname"), write=True)
t = usertext("tab=",message="Table name")
c = usertext("col=",message="Column name")

# Create the column and prepare to store numbers of type DBLE.
s.crecol(t, c, "DBLE")

finis()

Writing column data

In the second example we demonstrate how you store some data in a column.

writecolumn.py

#!/usr/bin/env python
from gipsy import *

init()
# Get set name, table and column name
s = Set(usertext("set=",message="Setname"), write=True)
t = usertext("tab=",message="Table name")
c = usertext("col=",message="Column name")

# Get some doubles from user
data = userdble("data=", message="data", nmax=20)

# Write the data to that column at set level (0) and
# do not append but start writing at the begin.

s[0, t, c, :] = data

finis()

List the columns of a given table in a set

Assuming you know the table name (and level), this example shows how to obtain the name of the columns in that table.

listcolumns.py

#!/usr/bin/env python
from gipsy import *

init()

set   = Set(usertext("set=",message="Setname"))
table = usertext("tab=",message="Table name")

cols = set.tabinq(table)

for col in cols:
   ctype, ccomm, cunts, rows = set.colinq(table,col)
   anyout("%s: Type= %s, Comment= %s, Units= %s Rows = %d"
        % (col,      ctype,       ccomm,     cunts,    rows) )

finis()

Reading the data in a column

The next step is to display all the data in a column, given the set and th name of a table. In this example we get the data of a table on set level (0).

readingcolumndata.py

#!/usr/bin/env python
from gipsy import *

init()

s = Set(usertext("set=",message="Setname"))
t = usertext("tab=",message="Table name")
c = usertext("col=",message="Column name")


for l in s[0,t,c,:]:
   anyout(str(l))

finis()

Finally we present a small program that demonstrates how to find tables, read columns and delete tables.

table.py

#!/usr/bin/env python
import gipsy
import numpy

gipsy.init()


setname = gipsy.usertext("INSET=", message="Enter set name")
myset = gipsy.Set(setname, write=True)

# Let's write a new table and column first to append to the header
tabname = "AMIGA"
colname = "ID"
coltype = 'DBLE'
subset = 0
myset.crecol(tabname, colname, coltype, subset=subset) #, colcomm="Just testing", colunits="kHz")
data = numpy.linspace(0,10)
myset[subset, tabname, colname, :] = data


# Generate a list of tables. This routine will find ALL tables, i.e.
# it does not depend on the given subset level
# So with 'tablis', you will find different subset levels. They are not
# necessarily the same as the subsets you entered with INSET=
# This is kind of a pitfall. If you need to access tables on a certain
# subset level, you have to filter (see next section)

# 'tablis' needs a maximum. The default is 100 (it is a C routine which uses
# dynamic memory allocation. If there are more
# than 'maxtables', we try again with an increased number
maxtables = 200
attempts = 0
while attempts >= 0:
   try:   
      tablist = myset.tablis(nitems=maxtables)
      attempts = -1
   except:
      maxtables += maxtables/2
      attempts += 1

# Now filter the table list on subset levels
filteredtables = []
for tabl in tablist:
   tablename, subset = tabl
   s = "Found %s at level %d"%(tablename,subset)
   gipsy.anyout(s)
   if subset in myset.subsets:
      filteredtables.append(tabl)

# Sort on subset first and then on name
from operator import itemgetter
sortedtables = sorted(filteredtables, key=itemgetter(1,0))

gipsy.anyout("Tables on the required subset level(s):")
for tnam, tlev in sortedtables:
   gipsy.anyout("%s   (%d)"%(tnam,tlev))


# For all the tables, get the information about its columns
for tabl in sortedtables:
   tablename = tabl[0]
   subset = tabl[1]
   gipsy.anyout("\nTABLE: %s"%tablename)
   cols = myset.tabinq(tablename, subset=subset)      # Get the columns for this table at this subset level
   s = "%10s %10s %10s %10s %-20s" % ("Columns", "Type", "Units", "Rows",  "Comment")
   gipsy.anyout(s)
   gipsy.anyout("="*len(s))
   for col in cols:
      ctype, ccomm, cunits, rows = myset.colinq(tablename, col, level=subset)
      gipsy.anyout("%10s %10s %10s %10d %-20s" % (col, ctype, cunits, rows, ccomm))
      

# Let's assume we want the contents of all the columns in table with index 0
if len(sortedtables) >= 1:
   table = sortedtables[0]
   tablename, subset = table
   gipsy.anyout("\n\nContents of table %s:" % tablename)
   cols = myset.tabinq(tablename, subset=subset)
   for col in cols:
      rows = myset[subset, tablename, col, :]
      gipsy.anyout("%s = %s" % (col, str(rows)))


key = "DELTAB="
mes = "Enter name of table to delete ... [skip deleting]:"
tablename = gipsy.userchar(key, mes, default=1, defval="")

if tablename:
   # Find this name in the list of table tuples
   deltabs = [x for x in sortedtables if x[0]==tablename.upper()]
   if deltabs:
      tname, subset = deltabs[0]
      myset.deltab(tname, subset)
      gipsy.anyout("I deleted table [%s] at subset level [%d]"%(tname,subset))
   else:
      gipsy.anyout("Did not find this table to delete!")

key = "DELTABLIST="
mes = "Delete all tables in given subset(s) ... Y/[N]:"
deltablist = gipsy.userlog(key, mes, default=1, defval=False)
if deltablist:
   for tablename, subset in sortedtables:
      myset.deltab(tablename, subset)

key = "DELALLTAB="
mes = "Delete all tables on all levels ... Y/[N]:"
delalltab = gipsy.userlog(key, mes, default=1, defval=False)
if delalltab:
   for tablename, subset in tablist:
      myset.deltab(tablename, subset)
      
   
gipsy.finis()

Plotting

We developed KPLOT (the Kapteyn Plot module) for creating graphical representations of your data (images, contour plots, gray scales, line plots, scatter diagrams etc.). Information is found in the Kplot tutorial pages.
Note: unfortunately, the current version of KPLOT cannot be used at the moment, so the following example will not work.

Plot an image of a 2d GIPSY (sub)set

In the following code we extract a GIPSY image and plot it with Kplot. An important problem to address is that the box of an image is always an integer box. Grid positions usually point to the centre of a pixel, so to be able to plot the entire box, we have to increase the box with half a pixel on either side of the plot boundaries to get a correct world coordinate system. Also the origin of a pixel can be moved by a fraction of a pixel. GIPSY and FITS images have a header item 'CRPIX' that defines this shift (i.e. the non integer part). Panel method setWorldImage() does it all in one call. It sets the world coordinate system, using the supplied arguments.

plotimage.py

#!/usr/bin/env python
import gipsy
import math
import Kplot
from math import *

gipsy.init()

# Get GIPSY set and enter box

sname = gipsy.usertext("SET=", "Set, 1 subset:")
s = gipsy.Set(sname)
gipsy.anyout("flo= %s, fhi= %s" % (s.slo, s.shi))

# Use GIPSY way to end program if dimensionality is wrong
if s.ndims( subset=True ) <> 2:
   gipsy.error(4,str("Subset must be 2-dimensional"))


boxmes = 'Enter box in '
for a in s.axes( True, 'inside' ):
   boxmes += a.name + ' '
boxmes += '   [entire subset]'

b = gipsy.usertext("BOX=", boxmes, defval='', default=1)
s.setbox(b)
gipsy.anyout( "blo=%s bhi=%s" % (s.blo,s.bhi) )

cw = s.subsets[0]
levels = gipsy.userreal( "LEVELS=", "Enter <10 levels between [%g,%g]" %
                         (s[cw,'datamin'], s[cw,'datamax']), default=0, nmax=10 )


M = s.subimage( cw )                # Get the image data


# Next is the Kplot part. Open a device, set a panel,
# prepare th GIPSY image for use in Kplot. Tell the 
# image object that it's origin is at CRPIX (as stored
# in the GIPSY axis object)

dev = Kplot.createDevice()
panel = dev.createPanel()
image = Kplot.Image( M )
#image.origin = ( -s.blo[0]+0.5, -s.blo[1]+0.5 )
panel.setWorldImage( image )


# Get values for the offset due to non integer crpix
# This is not a final solution to the problem. A new
# method will be created in the near future.

crp1 = s.axes( True, 'inside')[0].crpix
crp2 = s.axes( True, 'inside')[1].crpix
fr1 = crp1 - floor(crp1)

if (fr1 > 0.5):
   fr1 -= 1.0
fr2= crp2 - floor(crp2)
if (fr2 > 0.5):
   fr2 -= 1.0

d1 = fr1 - s.blo[0]
d2 = fr2 - s.blo[1]

box = Kplot.Box(drawgrid=False, xinterval=1, yinterval=1)
box.matrix = Kplot.Matrix.translate((d1,d2))


gipsy.anyout( str(levels) )
contourmap = Kplot.Contours( image, levels=levels, color="red")

panel.add( image, box, contourmap )

dev.output()

endkey = 'PRESSKEY='
gipsy.userlog( endkey, message=" ", default=1 )
gipsy.cancel( endkey )
dev.close(prompt=False)                      # We do not want a Kplot prompt


del s

gipsy.finis()

Use existing task GPLOT to create plots

As an alternative to Kplot you can use the PGPLOT based facilities provided by the GIPSY task GPLOT. GPLOT is still intensively used by GIPSY users and especially this group will benefit by using Python to write GPLOT input files. The next recipe is an example how you can build and use a GPLOT input file. In the recipe pages of GPLOT you will find numerous examples of GPLOT input files. GPLOT has only one important keyword (COMMAND=) which accepts GPLOT commands. These commands can be stored in a text file and executed all together with

GPLOT COMMAND=input myfile.txt

or with some other file name. You can close GPLOT with

COMMAND=quit

If we combine these to entries in a call to the xeq function then we get in Python:

gipsy.xeq( "GPLOT COMMAND=input %s;quit" % name )

and name is the name of your GPLOT input file on disk. We choose to create an unique file name with

name = tempfile.mktemp(dir='.')

which creates a file name that is unique in the current directory. A class 'Gplot' hides the complexity of this by defining some methods so you can concentrate on the second par of the code. You add a GPLOT command to a list with method addcommand(). If you want to plot the result then call method plotondevice(). This method has two arguments. The first is 'device'. This is a string that represents the GPLOT device (X11 = default, PCPSFILE, LCPSFILE etc.). The second is a filename. If omitted a temporary file is created in which the GPLOT commands are stored. After executing GPLOT, the file is removed from disk. If you entered a file name then a file with that name is created and it is NOT removed afterwards so you are able to inspect it or use it directly as input file for GPLOT. First we show you the code of the Gplot class.

Gplotclass.py

#!/usr/bin/env python
import gipsy


import os
import tempfile
class Gplot:
   """
------------------------------------------------------------------
Example of a class to facilitate the use of an existing
GIPSY task. In this example we made methods for
writing and executing a GPLOT input file (see GPLOT documentation
in GIPSY)
Methods:
Gplot.addcommand( str ) : str is a string which represents a Gplot 
                          command.
Gplot.plotondevice( device, gplotinputfile ): device is a string 
                          representing a valid GPLOT device (X11, 
                          PCPSFILE etc.). gplotinputfile is a 
                          file name. If omitted, a tempory file
                          is created else a file with the given 
                          name is kept on disk.
------------------------------------------------------------------
"""
   def __init__( self ):
      self.commandlist = []
   def addcommand( self, command ):
      self.commandlist.append( command )
   def plotondevice( self, device="X11", gplotinputfile = '' ):
      if gplotinputfile == '':
         name = tempfile.mktemp(dir='.')       # Get a unique name for a file
         remfile = True
      else:
         name = gplotinputfile
         remfile = False
      f = open( name, 'w' )
      f.write( "device %s\n" % device )
      for com in self.commandlist:
         f.write( com + '\n' )

      f.close()

      # Start GPLOT and load the input file. Provide the COMMAND=
      # keyword with value 'quit' to finish GPLOT.

      gipsy.xeq( "GPLOT COMMAND=input %s;quit" % name )

      if remfile:
         os.remove( name )                     # Remove the file from disk

A set, subset is asked and the program checks whether the dimension of the set is 2. If not then the program is aborted with GIPSY function 'error'. The minimum and maximum data values are read from the header of the (sub)set and used in a prompt to get the values for the wanted contour levels. Note that we store the input in a string and not in an array of reals using userreal as input function.

gplot.py

#!/usr/bin/env python
import gipsy
from Gplot import *


gipsy.init()

# For a PostScript file, write something like PCPSFILE/myfile.ps
defdev = "X11"
dev = gipsy.usertext( "DEVICE=", "Enter device:  [%s]" % defdev, 1, defdev )

setname = gipsy.usertext( "INSET=", "Enter set, subset", 0 )
s = gipsy.Set( setname )

# Use GIPSY way to end program if dimensionality is wrong
if s.ndims( subset=True ) <> 2:
   gipsy.error(4,str("Subset must be 2-dimensional"))

cw = s.subsets[0]
levels = gipsy.usertext( "LEVELS=", "Enter <10 levels between [%g,%g]" %
                         (s[cw,'datamin'], s[cw,'datamax']), 0 )

# Write the GPLOT input file

ob = Gplot()
ob.addcommand( "inset %s" % setname )
ob.addcommand( "box" )
ob.addcommand( "xsize 120" )
ob.addcommand( "ysize 120" )
ob.addcommand( "colour cyan" )
ob.addcommand( "levels %s" % levels )
ob.addcommand( "contour" )
ob.addcommand( "colour yellow" )
ob.addcommand( "font roman" )
ob.addcommand( "frame" )
ob.plotondevice( dev )

# Alternatives:
#ob.plotondevice( dev, gplotinput="mygplotinput.txt" )
#obplotondevice()                                      # X11 is default device
#ob.plotondevice( device='X11', gplotinput="mygplotinput.txt" )

gipsy.finis()

Modules

Using standard modules to list GIPSY sets

In GIPSY there is a task DISK that lists all the GIPSY sets in your working directory. Let's use Python to create a similar task, but now with the option to enter a directory from which all sets are listed including those in subdirectories. First you try to find a method that gets all the files in an directory and its subdirectories. This method is os.walk() so we need to import the os module. A GIPSY set must have a .descr extension so we need to find a method to match. This can be done with fnmatch() from the fnmatch module. Finally we need to know if the image part (ends on .image) exists. For this we use the glob() method. Further, os.stat() gives us file info like file size. The next program uses all these methods to create a handy utility. Notes:

We did not include a routine that checks whether a file is really a GIPSY descriptor.
The os.walk() method does not follow symbolic links.

findsets.py

#!/usr/bin/env python
import os
import fnmatch
import gipsy
import glob


# Find all GIPSY sets in all sub directories starting at a given path
# Note: Symbolic links to directories are not treated as subdirectories


gipsy.init()
walkpath = gipsy.usertext( "PATH=", 
                           "Enter path to start searching:   [current dir]", 
                           default=1, defval='.', nmax=512 ) 
filetupel = os.walk( walkpath )
for dir in filetupel:
   for filenam in dir[2]:
      if fnmatch.fnmatch( filenam, "*.descr" ):
         base = filenam.replace( ".descr", "" )
         descrpart = dir[0]+'/'+ base + '.descr'
         impart = dir[0]+'/'+ base + '.image'
         sized = os.stat(descrpart)[6] / 1024.0
         if glob.glob( impart ): 
            sizei = os.stat(impart)[6] / 1024.0
         else:
            sizei = 0.0
         gipsy.anyout( "%s decr=%d im=%d (kb)" %(dir[0]+'/'+ base ,sized, sizei) )


gipsy.finis()

Using PYFITS to create a fits image viewer

PYFITS is a Python module that enables a programmer to read, write and process FITS files. Next recipe is a small demonstration program that converts a FITS file to GIPSY format and displays the data with GIPSY task SLICEVIEW. Note that we are using import gipsy and import pyfits to avoid problems with name spaces. The program asks for the name of a FITS file on disk. Then it wants to know what the name of the GIPSY set will be. A minimal set of header items is read from the FITS file and is used to construct the GIPSY set. The interface for the data is very simple. One can copy the image data with very few lines of code. The minimum and maximum of the image data are calculated and stored in the header of the new set. Then the GIPSY viewer SLICEVIEW is called to view the data for 2-dimensional data only.

fits2gds.py

#!/usr/bin/env python
import gipsy
import pyfits
import tempfile


gipsy.init()

fitsfile = gipsy.usertext("FITSFILE=", "Enter name of fits file: " )
fimg = pyfits.open( fitsfile )           # Open a fits file

# name = gipsy.usertext("SETNAME=", message="Set name: ")
name = tempfile.mktemp(dir='.')
gipsy.anyout( name )
s = gipsy.Set(name, create=True)         # Create a set with this name


# Next we want to read the minimum set of axes properties
# so that we can build a GIPSY set. We use the get() method
# to access a header item because a default value can be 
# returned if an item doed not exist. For example DSS fits 
# files does not have a ctype, so we substitute axis names 
# 'X', 'Y', 'Z', 'A', 'B' etc.

prihdr = fimg[0].header                  # Assuming the data is in first hdu 
naxis = prihdr['naxis']
for i in range(naxis):
   # Generate 'X', 'Y', 'Z', 'A', 'B' etc.
   defnam = chr(divmod(ord('X')-ord('A')+i,26)[1] + ord('A'))
   istr = str(i+1)                       # Start counting at 1, NAXIS1 etc.
   len = prihdr.get('naxis'+istr, 1)
   ori = prihdr.get('crpix'+istr, 1)
   nam = prihdr.get('ctype'+istr, defnam)
   gipsy.anyout( 'GIPSY: Creating axis %s %s %s' % (nam,ori,len ) )
   s.extend( nam, ori, len )
   s[0,'cunit'+istr] = prihdr.get('cunit'+istr, '??')   
   s[0,'crval'+istr] = prihdr.get('crval'+istr, 1)   
   s[0,'cdelt'+istr] = prihdr.get('cdelt'+istr, 1)   

gdsimage = s.image

# Copy the contents 
# Pitfall: Do not use the statement: gdsimage = fitsdata
# which would fail to copy the data because you lost your reference
# to s.image

gdsimage[:] = fimg[0].data

# Update the header with the minimum and maximum values of the data
s.wminmax()


if naxis == 2: 
   gipsy.xeq( "SLICEVIEW INSET=%s ggiopt=fixcol" % name )

s.delete()

gipsy.finis()

This program shows how fast you can develop an application but we are not ready yet. For a robust task you need to check whether your FITS file exists and that your GIPSY set does not already exists. Futher, we assumed that the image data was stored in the first hdu of the FITS file.

GUI Programming

Writing a GIPSY task in Python with a graphical user interface is relatively easy. Such a task usually consists of three components:

a collection of functions which will react to events caused by the user's interaction with the GUI elements;
code which creates and positions the GUI elements;
code which connects the events with the event-handling functions.

Most GUI elements are associated with a user input keyword which will be changed according to the user's actions on the element. Consequently it is usually sufficient if the task handles the keyword events caused by these keyword changes. An exception is the handling of cursor information from a plot window.
More information about GIPSY's graphical user interface and event handling can be found here.

Interactive graph example

This annotated program draws a graph of the function y = A.sin(B.x). The user can modify the values of A and B with either an input field or a slider ('gauge'). Any change causes the graph to be re-drawn.

gui-graph.py

#!/bin/env python

from gipsy import *
import math

# ========================= Initial Parameters ===============================
NPTS = 100
a  = 0.7; b = 1.0
xmin =  0.0; xmax = 6.2831853;
ymin = -1.1; ymax = 1.1
incr = (xmax-xmin)/NPTS
xvals = map(lambda x: xmin+x*incr, range(NPTS))   
yvals = None;

# ========================== Functions =======================================

# ---------------------- keyword setter --------------------------------------
def triggerkey(key, value=None):
   wkey("%s%s" % (key, usertext(key, default=2, defval=value)))

# ---------------------- graph drawer ----------------------------------------
def drawplot():
   global yvals
   if yvals: pgline(xvals, yvals)                    # erase any previous graph
   yvals = map(lambda x: a*math.sin(b*x), xvals)     # compute new y-values
   pgline(xvals, yvals)                              # draw graph

# ---------------------- event handler for QUIT= keyword ---------------------
def quit(cb):
   if userlog(cb.key):
      wkey(cb.key)                                   # set back to False
      finis()                                        # terminate program

# ---------------------- event handler for both A= and B= keywords -----------
def number(cb):
   globals()[cb.var] = userreal(cb.key)              # assign to global variable
   drawplot()                                        # (re-) draw graph

# ====================== Main Program ========================================

init()                                               # program initialization

base = GgiBase()                                     # Ggi initialization

button = base.Button("QUIT=", "Terminate program")   # create button
text1  = base.TextField("A=", None, 10)              #   ,,   text input field
text2  = base.TextField("B=", None, 10)              #   ,,   text input field
gauge1 = base.Gauge("A=", None, 200, -1.0, 1.0)      #   ,,   slider
gauge2 = base.Gauge("B=", None, 200, 0.5, 5.0)       #   ,,   slider
plot   = base.PlotField("MYPLOT", 400, 300)          #   ,,   PGPLOT device
label  = base.Label("y = A*sin(B*x)")                #   ,,   text label

gauge1.setLabel(" ", 1)                              # remove standard ...
gauge2.setLabel(" ", 1)                              # ... label from sliders

button.setPosition(0, None, 0, None)                 # top left in window
plot.setPosition(0, None, 10, button)                # below button
label.setPosition(-(plot.width+label.width)/2, plot, # centered above plot
                  -plot.height-label.height-5, plot)
text1.setPosition(0, None, 10, plot)                 # below plot
text2.setPosition(0, None, 0, text1)                 # below text1
gauge1.setPosition(0, text1, 10, plot)               # right of text1
gauge2.setPosition(0, text2,  0, text1)              # right of text2

base.realize()                                       # activate Ggi
plot.open()                                          # open PGPLOT device
pgenv(xmin, xmax, ymin, ymax)                        # draw frame
plot.xor(True)                                       # XOR drawing mode
pgsci(3)                                             # green graph colour

KeyCallback(quit,   "QUIT=")                         # register QUIT= callback
KeyCallback(number, "A=", var='a')                   #    ,,    A=       ,,
KeyCallback(number, "B=", var='b')                   #    ,,    B=       ,,

triggerkey("A=", a)                                  # initialize keywords ...
triggerkey("B=", b)                                  # ... to draw first graph

mainloop()                                           # action!

Simple image viewer

The following example shows how a simple image viewer for GIPSY sets can be constructed. It consists of a plot window into which the image will be drawn, a composite element which enables the user to manipulate the colour map, and some menus and buttons for specifying the data set etc.

gui-image.py

#!/usr/bin/env python
from gipsy import *

NCI=200                                   # number of color indices

# ========================== Functions =======================================

# ---------------- event handler for FILE= keyword ---------------------------
def file_cb(cb):
   choice = userint(cb.key, defval=-1)
   if choice==0:
      GgiBase.base.Inset('SETNAME=', 'SET=', 'BOX=').ndims(2)
   elif choice==1:
      if GgiBase.base.verify('Really quit?', 'Yes, quit.', 'No, continue!'):
         finis()

# ---------------- event handler for SET= keyword ----------------------------
def inset_cb(cb):
   global set
   try:
      set.close()                         # close any previously opened set
   except:
      pass
   set = Set(usertext(cb.key))            # open set

# ---------------- event handler for BOX= keyword ----------------------------
def box_cb(cb):
   set.setbox(usertext(cb.key, defval=''))
   wkey('SHOW=YES')                       # cause SHOW= keyword event

# ---------------- event handler for SHOW= keyword ---------------------------
def show_cb(cb):
   if userlog(cb.key):
      image = set.subimage(set.subsets[0])        # (numpy-) image array
      if len(image.shape)!=2:                     # check dimensionality:
         reject('SET=', 'Must be 2-dimensional')  # wrong: reject input
         return
      idim = image.shape[1]                       # horizontal size
      jdim = image.shape[0]                       # vertical size
      tr = [-0.5/idim, 1.0/idim, 0.0, -0.5/jdim, 0.0, 1.0/jdim] # pgimag matrix
      imgmin = image.min(); imgmax = image.max()  # image min- and max-values
      pgsimi(userlog('INTERPOL='))                # image interpolation?
      pgimag(image, idim, jdim, 1, idim, 1, jdim, imgmin, imgmax, tr) # draw
      coled.limits(imgmin, imgmax)                # adjust color wedge's ...
      coled.units(str(set['BUNIT']))              # ... limits and units
      show.active(True)                           # enable 'SHOW' button
   wkey(cb.key)                                   # reset keyword (and button)

# ====================== Main Program ========================================

init()                                    # Open Hermes communication

base  = GgiBase()                         # Ggi initialization

file  = base.Menu('FILE=', 'Options menu', ['Inset', 'Quit'])          # menu
intp  = base.Button('INTERPOL=', 'Interpolation ON/OFF').setLabel('INTP') # btn
show  = base.Button('SHOW=', 'Show image').active(False)               # button
plotf = base.Form('PLOTFORM', 4)          # Frame for PGPLOT window
coled = base.ColorEditor(16, NCI)         # Color manipulation element        
logo  = base.Logo()                       # The GIPSY logo

file.setPosition(0, None, 0, None)        # top left
intp.setPosition(0, file, 0, None)        # right to 'file'
plotf.setPosition(0, None, 10, file)      # below 'file'
coled.setPosition(0, None, 5, plotf)      # below 'plotf'
logo.setPosition(-logo.width, plotf, -logo.height, coled) # bottom right corner
show.setPosition(-show.width, plotf, 0, None)             # top right corner

plot = plotf.PlotField('PLOTTER', 400, 400)  # PGPLOT window
plot.setPosition(0, None, 0, None)        # top left in 'plotf'
plot.colors(16+NCI+1, False)              # allocate colors

base.realize()                            # initialize GUI
plot.open()                               # open PGPLOT device
pgscir(16, 16+NCI-1)                      # set color index range

KeyCallback(inset_cb, 'SET=')             # register keyword handler
KeyCallback(box_cb,   'BOX=')             #   ,,        ,,     ,,
KeyCallback(file_cb, 'FILE=')             #   ,,        ,,     ,,
KeyCallback(show_cb, 'SHOW=')             #   ,,        ,,     ,,

mainloop()                                # Action!

PC	projection centre, effectively grid position (0,0,...,0).
AC	axis centre. If the length of axis i is NAXISi and the reference pixel is CRPIXi, then the i-th coordinate is given by the expression NAXISi / 2 - CRPIXi.
BOX= -3 -3 4 4	specifies a 2-dimensional frame containing 64 pixels from lower left to upper right
BOX= PC D 8 8	box with size 8x8 pixels with PC as center
BOX= * 3 0 0 * 45 0 0 D 6 arcmin 6 arcmin	box centered at 3h0m0s and 45d0m0s with sizes in RA and Dec of 6 minutes of arc. The units of the sizes must be compatible with the axis units in the header of the set.