Function: WCSPOS Purpose: WCSPOS prompts the user to define a position in a subset. Category: USER IO, PHYSICAL COORDINATES File: ioutils.py Author: M.G.R. Vogelaar Use: grids, world, pixels, subsetunits, errmes = WCSset.wcspos(subnum, # Integer key="POS=", # Character string mes="", # Character string userpos='') # String wcspos A method of class WCSset. subnum Index of current subset in subset arrray. Note that this is an index not a subset coordinate word. key Keyword for prompt. If you do not specify this parameter then by default it is "POS=" mes Message for prompt. If you do not specify this message, then an appropriate message will be created. userpos A string to replace the manual input. No prompt will appear. WCSPOS therefore is also suitable to use in a non interactive environment. If you want a prompt then this parameter should be omitted (or it must be an empty string). It returns: grids A 2 dimensional NumPy array with grids values for the entered position(s). Each row represents a position. Each element in a row is a coordinate. world A 2 dimensional NumPy array with world coordinates for the entered position(s) pixels A 2 dimensional NumPy array with pixel coordinates. Pixel coordinates following the FITS standard, e.g. they start at 1. These coordinates can be used as an interface to the Kapteyn Package. This is a Python based toolkit for plotting and conversion of world coordinates. subsetunits Units corresponding to axes in the current subset. The units of the first world coordinate in a position has the units of the first element in this list etc. errmes An error message if something was wrong. One should check this return value to decide whether the input was valid or not. Description: WCSPOS is a method that returns a user given position in grid-, world- and pixel coordinates. This can be a single position or a sequence of positions. The syntax of a position is described in the 'Syntax' section. WCSPOS can only be used in GIPSY tasks written in Python. That is because it uses a system of world coordinates based on the library WCSLIB for which we made a so called 'Python binding'. If one writes an application in Python than one should use the more versatile routines 'wcsinp', 'wcspos' and 'wcsbox' for the input of GIPSY sets (or FITS files), positions and axis limits. WCSPOS replaces the function GDSPOS. The input of sky systems is the main difference between these two routines. It has a more options to specify sky systems. It supports ICRS, elliptic terms of abberation for transformations to or from FK4, a rich set of projection systems etc. WCSPOS and GDSPOS support a compatible set of units, but WCSPOS can also deal with inverse units like 1/m. Notes: In GIPSY tasks we often use WCSPOS in combination with WCSINP and WCSBOX. In most situations the conversion between grids and world coordinates is independent of the subset. E.g. in a sequence of channel maps, each RA,Dec image represents the same coordinate system. However, when we create subsets with only one spatial axis, then the matching spatial coordinate depends on the selected subset (e.g. with XV maps). Therefore we need to enter the index of the subset we want to use. Example of a small GIPSY test task ++++++++++++++++++++++++++++++++++ Next example can be copied (cut & paste) to a file on disk. Give it the name 'postest' and make it executable with: > chmod u+x postest Then start GIPSY with: > gipsy and enter the name of your test application (postest). #!/usr/bin/env python import gipsy import ioutils gipsy.init() # A GIPSY set is created first. # Then a list with positions is processed as a closure test and as a test # for error messages. Finally, you are prompted in a loop (aborted by an # empty input i.e. 'enter' only) to enter your own test positions. # If you want to use your own GIPSY sets or FITS files, then # use method 'wcsinp()' instead of constructor 'WCSset()' setname = "wcspos_testset" command = "CREATE BUNIT= \ CDELT1= -1.200000000000E-03 CDELT2= 1.497160000000E-03 CDELT3= -7.812500000000E+04 \ CROTA2= \ CRPIX1= 5 CRPIX2= 5 CRPIX3= 5 \ CRVAL1= 1.787792000000E+02 CRVAL2= 5.365500000000E+01 CRVAL3= 1.415418199417E+09 \ CTYPE1= RA---TAN CTYPE2= DEC--TAN CTYPE3= FREQ-OHEL CTYPE4= \ CUNIT1= CUNIT2= CUNIT3= HZ \ DELETE= y \ DRVAL3= 1.050000000000E+03 DUNIT3= KM/S FREQ0= 1415418199.417 \ FUNCTION= rang(0.0,3.0) \ INSTRUME= WSRT \ NAXIS1= 10 NAXIS2= 10 NAXIS3= 10 \ OUTSET= %s"%(setname) gipsy.xeq(command) gipsyset = ioutils.WCSset(setname + " freq 0") gipsyset.printinfo() userpos = ["0 0", "eq 178.7792 eq 53.655", '{} header("crval1"), {} header("crval2")', "eq 178 ga 22"] for upos in userpos: print gipsyset.wcspos(0, userpos=upos) gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) if errm == '': gipsy.anyout("%s -> grids: %f %f -> world: %f %f (%s %s)" % (upos, gr[0][0], gr[0][1], wor[0][0], wor[0][1], units[0], units[1])) else: gipsy.anyout("Error message test: %s --> %s" % (upos, errm)) # Using a prompt to enter positions for the same GIPSY set: # Abort the loop if only enter is pressed. key = "POS=" while 1: grids, world, pix, units, errmes = gipsyset.wcspos(0, key, "") gipsy.cancel(key) if not len(grids): if errmes: gipsy.reject(key, errmes) else: break else: for g,w in zip(grids, world): s = str(g) + "=" + str(w) + ' ' + str(units) gipsy.anyout(s) gipsy.finis() # Run the module ioutils.py as a GIPSY task in test mode to see # the output of many input examples. Syntax ++++++ The syntax is described in the documentation of the 'positions' module in the Kapteyn Package. The web link is: http://www.astro.rug.nl/software/kapteyn/positions.html#position-syntax Real input examples for sky definitions ....................................... POS={eq} 178.7792 {} 53.655 # As a sky definition between curly brackets POS={} 178.7792 {} 53.655 # A world coordinate in the native sky system POS={eq,B1950,fk4} 178.12830409 {} 53.93322241 # With sky system, reference system and equinox POS={fk4} 178.12830409 {} 53.93322241 # With reference system only. POS={eq, B1950,fk4, J1983.5} 178.1283 {} 53.933 # With epoch of observation (FK4 only) POS={eq B1950 fk4 J1983.5} 178.1283 {} 53.933 # With space as separator POS=ga 140.52382927 ga 61.50745891 # Galactic coordinates POS=ga 140.52382927 {} 61.50745891 # Second definition copies from first POS=su 61.4767412, su 4.0520188 # Supergalactic POS=ec 150.73844942 ec 47.22071243 # Ecliptic POS={} 178.7792 6.0 # Mix world- and pixel coordinate POS=5.0, {} 53.655 # Mix with world coordinate in native system Notes: * Mixing sky definitions for one position is not allowed i.e. one cannot enter POS=ga 140.52382927 eq 53.655 * If you mix a pixel- and a world coordinate in a spatial system then this world coordinate must be defined in the native system, i.e. {} We can also specify positions in data structures with only one spatial axis and a non-spatial axis (e.g. position velocity diagrams). The grid coordinate of the missing spatial axis is derived from the subset (entered as an index with parameter 'subnum'). If one of the axes is a spectral axis, then one can enter world coordinates in a compatible spectral system: POS={} 53.655 1.415418199417E+09 hz # Spatial and spectral world coordinate POS={} 53.655 1.415418199417E+03 Mhz # Change Hz to MHz POS=53.655 deg 1.415418199417 Ghz # to GHz POS={} 53.655 vopt 1.05000000e+06 # Use spectral translation to enter optical velocity POS={} 53.655 , vopt 1050 km/s # Change units POS=10.0 , vopt 1050000 m/s # Combine with a pixel position POS={} 53.655 vrad 1.05000000e+06 # Radio velocity POS={} 53.655 vrad 1.05000000e+03 km/s # Radio velocity with different unit POS={} 53.655 FREQ 1.41541820e+09 # A Frequency POS={} 53.655 wave 21.2 cm # A wave length with alternative unit POS={} 53.655 vopt c_/285.51662 # Use speed of light constant to get number in m/s Note: For positions in a data structure with one spatial axis, the other (missing) spatial axis is identified by a grid coordinate derived from the current subset. This restricts the spatial world coordinates to their native WCS. We define a world coordinate in its native sky system with {} Notes: A sky definition needs not to be repeated for other coordinates than the first. The second definition therefore can be empty as in {} as in: POS={eq} 178.7792 {} 53.655 Only one definition is allowed in a position. If you mix definitions as in: POS={eq, B1950} 178.7792 ga 22.01234 Note: World coordinates followed by a unit, are supposed to be compatible with the Projection object (attribute of an object from class WCSset). So if you have a header with spectral type FREQ but with a spectral translation set to VOPT, then POS={} 53.655 1.415418199417E+09 hz is invalid, but POS=10.0 , vopt 1050000 m/s is ok and also POS={} 53.655 FREQ 1.415418199417e+09 is correct. Sexagesimal notation ..................... Read the documentation at function 'parsehmsdms()' in the positions module of the Kapteyn Package for the details. Here are some examples: POS=11h55m07.008s 53d39m18.0s POS={B1983.5} 11h55m07.008s {} -53d39m18.0s POS= -33d 0d Reading from file with function 'READCOL()' ............................................. We introduced a command 'readcol' to read a column from a file on disk. Function readcol() requires the filename enclosed by double quotes. At the GIPSY command line (Hermes) you need to use back-quotes to escape these double quotes. Here is an example: POS=`readcol("test123positions.txt", col=0) readcol("test123positions.txt", col=2)` !! Note the back quotes to prevent an unwanted interpretation of the input by the GIPSY user interface Hermes. Structure of output .................... In the first example of this document we presented a small program that uses method 'wcspos()' for processing positions. If we focus on the anyout statement that prints the result in the GIPSY log file then we that it is called in a loop (two strings that represent a position in the list 'userpos'). For each position 'wcspos()' returns a tuple with the following elements: gr: The entered position in grids. Returned is a two dimensional array. Each row represents a position. Each element in a row is one of the coordinates. The order of the coordinates is the same as the order of the axes in your data structure. wor: The entered position in world coordinates. See also remarks about grids. pix: The entered position in (FITS) pixels. See also remarks about grids. units: Tuple with units. One for each coordinate. errm: If empty then no error condition exists, else, an error occurred and one should take care of this. userpos = ["0 0", "eq 178.7792 eq 53.655"] for upos in userpos: gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) if errm == '': gipsy.anyout("%s -> grids: %f %f -> world: %f %f (%s %s)" % (upos, gr[0][0], gr[0][1], wor[0][0], wor[0][1], units[0], units[1])) The code of GIPSY's ioutils.py in $gip_tsk shows many examples. Task mapmaker.py shows how 'wcspos()' can be used in an interactive environment. Errors: ........ The position parser is flexible but there are some rules. If the input cannot be transformed into coordinates then an appropriate message will be returned. In some cases the error message seems not related to the problem but that seems inherent to parsers. If a number is wrong, the parser tries to parse it as a sky system or a unit. If it fails, it will complain about the sky system or the unit and not about the number. Updates: Apr 23, 2009: VOG, Document created. Jul 15, 2010: VOG, Made compatible with Kapteyn Package 2.0 (module positions.py). Extended the documentation. Feb 03, 2011: VOG Removed documentation that can be found in positions module of K.P. Improved example.