#!/usr/bin/env python #---------------------------------------------------------------------- # FILE: ioutils.py # PURPOSE: Provides three functions for user I/O # 1) wcsinp() for the input of GIPSY sets or FITS files # 2) wcsbox() for the input of axis limits in a selected data slice # 3) wcspos() for the input of positions in grids or world coordinates # AUTHOR: M.G.R. Vogelaar, University of Groningen, The Netherlands # DATE: April 20, 2009 # UPDATE: Jun 05, 2010 # Jan 13, 2011 # VERSION: 0.8 # # (C) University of Groningen # Kapteyn Astronomical Institute # Groningen, The Netherlands # E: gipsy@astro.rug.nl #---------------------------------------------------------------------- # Imports from math import * # For Pythons 'eval()' function from re import split as re_split from string import upper as string_upper from sys import stdout # For flushing a buffer from numpy import nan as unknown # new name for a NaN from numpy import asarray, zeros, floor from numpy import array2string # Used to format output import gipsy # Provides classes and bindings to GIPSY routines from kapteyn import wcs # The WCSLIB binding from kapteyn import positions # The coordinate parser (expressions, units) from os import listdir, stat, path # Used to get names of files on disk from fnmatch import fnmatch # Match file name against pattern def nint(x): """-------------------------------------------------------------- Purpose: Calculate a nearest integer compatible with then definition used in GIPSY's coordinate routines. Inputs: x- A floating point number to be rounded to the nearest integer Returns: The nearest integer for 'x'. Notes: This definition adds a rule for half-integers. This rule is implemented with function floor() which implies that the left side of a pixel, in a sequence of horizontal pixels, belongs to the pixel while the right side belongs to the next pixel. This definition of a nearest integer differs from the Fortran definition used in pre-April 2009 versions of GIPSY. -----------------------------------------------------------------""" return floor(x+0.5) def prompthdu(key, mes): def askhdu(hdulist): """------------------------------------------------------------------------------- Helper function for Set class. The function prompts for a header data unit only if the file in the constructor was a FITS file. -------------------------------------------------------------------------------""" print hdulist.info() # Make sure the message is in the log file before the next keyword prompt stdout.flush() n = len(hdulist) if n > 1: while True: hnr = gipsy.userint(key, mes, default=1, defval=0) if hnr < n: break else: gipsy.reject(key, "Hdu number must be %d <= n < %d" %(0, n)) else: hnr = 0 return hnr return askhdu def promptalt(key, mes): def askalt(alts, header): """------------------------------------------------------------------------------- Helper function for Set class. The function prompts for a alternative header. A primary header is selected if a space is returned (' '). -------------------------------------------------------------------------------""" mes = "Select header ... " # Build a default message for c in alts: mes += "%s/"%c mes += "[Primary header]:" alth = '' while True: # Return only one uppercase character or ' ' alth = gipsy.usercharu(keyword=key, message=mes, default=1, defval=alth, width=1, nmax=1) # Is this selection allowed? if alth == '' or alth in alts: break else: gipsy.reject(key, "This header is not in the list!") return alth return askalt def MinimalHeader(header): result = {} naxis = header['NAXIS'] result['NAXIS'] = naxis for i in range(1,naxis+1): keyword = 'NAXIS%d' % i result[keyword] = header[keyword] keyword = 'CRPIX%d' % i try: result[keyword] = header[keyword] except: pass return result class WCSset(gipsy.Set): def __init__(self, setspec, create=False, write=False, gethdu=None, getalt=None, minimal=False): """--------------------------------------------------------------------------- Attributes: axperm_in: Axis permutation array for subset (FITS definition) axperm_out: Axis per. array for axes outside subset. wcsinfo: String with WCS information about subset axes transax_in: FITS axis number that corresponds to spectral axis in a subset or None transax_out: FITS axis number that corresponds to spectral axis outside a subset or None translations: List with allowed spectral translations for spectral axis in subset. Could be empty for incorrect FITS headers (e.g. missing rest frequency) projection: An object of class Projection from module wcs. Serves as a base object. subproj: A projection object that represents the entered subset. mixax: The axes number of the missing speatial axes or None. The application will be aborted if it cannot find a matching axis. projaxnum: The axes permutation array for the projection object associated with the subset. It can have more elements than the dimension of the subset when it contains the axis number of the missing spatial axis (e.g. in a XV map). Also, at least for two dimensional subsets the order need not to be the same as the order in GDS (i.e. fixed order of as found in header of set). spatialmap: If True, the subset is two dimensional and has two spatial axes slicepositions: Pixel positions on repeat axes to identify a slice. The positions are in order of the order of the repeat axes in the header. sliceaxnums: Axis numbers of the repeat axes. The numbers are 1 based and are ordered as the the order of the repeat axes in the header. Notes: The Set class defines an array with axis numbers in the given subset. Assume a (Ra, Dec, Freq) set and one defines a subset (Dec, Freq). The the (FITS) axis numbers for the entire set (i.e. 'setlevel') is (1,2,3). The axis permutation array for the subset is (2,3). The projection object associated with this Set defines the axis numbers (FITS numbering) of the longitude, latitude and spectral axis. So: Set.Projection.lonaxnum = 1 Set.Projection.lataxnum = 2 Set.Projection.specaxnum = 3 To test whether an axis is available in a subset one needs only to check if the required axis number is element of attribute axperm_in. On the other hand, this class creates a derived projection object that represents the subset as specified by a user. Then the axis numbers are redefined: Set.SubProjection.lonaxnum = None Set.SubProjection.lataxnum = 1 Set.SubProjection.specaxnum = 2 These numbers cannot be compared with the contents of axperm_in ------------------------------------------------------------------------------""" # Create a Set object using its init method gipsy.Set.__init__(self, setspec, create, write, gethdu, getalt) # We need one name for the entered data whether it is gds or a FITS file if self.fits: self.dataname = self.fitsname else: self.dataname = self.name # In the Kapteyn Package we have functionality which has the focus # on two dimensional images (maps). In GIPSY we allow for any dimension. # So we build our own projection object for the set input. try: if minimal: proj = wcs.Projection(MinimalHeader(self)) else: proj = wcs.Projection(self, alter=self.altwcs) except Exception, message: raise Exception, message # Next we do some analysis about the axes in this data set. # Do we have a missing spatial axis or are spectral translations possible? # First we convert the axis permutation array to the FITS standard (i.e. start with 1) self.axperm_in = [i + 1 for i in self.axperm('inside')] self.wcsinfo = [] transax_in = None # The axis number of the spectral axis in a subset for ax in self.axperm_in: i = ax - 1 # Axis numbers start with 1. if proj.types[i] == None: # Could be unknown type, or '-' at the wrong place # Then warn the user that he/she could get unexpected results info = "WCS Subset axis %s of unknown wcs type (units: %s). Axis is assumed linear!"%(proj.ctype[i], proj.units[i]) else: info = "WCS Subset axis %s of type '%s' (units: %s)"%(proj.ctype[i], proj.types[i], proj.units[i]) self.wcsinfo.append(info) if proj.types[i] == 'spectral': transax_in = i + 1 self.transax_in = transax_in self.axperm_out = [i + 1 for i in self.axperm('outside')] transax_out = None # The number of the spectral axis as repeat axis for ax in self.axperm_out: i = ax - 1 if proj.types[i] == 'spectral': transax_out = i + 1 self.transax_out = transax_out # Make Projection object an attribute self.projection = proj # It is possible that a user wants to inspect a data structure # with one spatial axis. For coordinate transformations we need # a matching spatial axis. Before we create a Projection object for # the subsets, we need to find a complete set of axis numbers. # Note that we don't change the dimensionality of the GDS structure. # The projection is only used for coordinate transformations. # If a missing spatial axis is found, the system can transform # between world coordinates and grids using the so called mixed # coordinate system from WCSLIB. newaxperm, mixax = self.findmatchingaxis() self.mixax = mixax # Store the missing axis number # In GIPSY's Python binding the axis numbers start at 0. # In WCSLIB the axis numbers start at 1. Then modify # 'newaxperm' for compatibility with WCSlib axnum = [i+1 for i in newaxperm] if len(axnum) == 0: # We have a problem here. The user probably specified a pixel (0 dimensional subset) self.subproj = None self.projaxnum = None else: self.subproj = self.projection.sub(axnum) # Modify the projection to support new axis number list self.subproj.allow_invalid = True # Do not raise exceptions after conversion problems self.projaxnum = axnum self.spatialmap = 'longitude' in self.subproj.types and 'latitude' in self.subproj.types # For a smooth integration with the maputils module from the Kapteyn # Package, we store the grids on the repeat axes for the slices as tuples. # Method grid() decodes a GDS coordinate word into a grid. In GIPSY # the grids are relative to the header value of CRPIX so in a # FITS environment one should take notice whether returned coordinates # are based on the grid- or the pixelsystem. slicepositions = [] sliceaxnums = [] setdim = self.naxis # self.ncoords includes hidden axes, which we don't want here subdim = self.ndims(subset=True) repeatdim = setdim - subdim if repeatdim > 0: for ax in self.axes(True, 'outside'): sliceaxnums.append(ax.axnum+1) # In gipsy module axnum is zero based for cword in self.subsets: spos = [] for ax in self.axes(True, 'outside'): grid = self.grid(ax.axnum, cword) pix = grid + int(nint(self.projection.crpix[ax.axnum])) spos.append(pix) # Here we have the following problem: The repeat axes are given by the grid() method # in order of the (user) input. For interaction with the maputils objects we # need slicepositions to be given in order of the axes in the header, so we need # to sort the slicepositions list. if repeatdim > 1: # For sorting tmp = zip(sliceaxnums, spos) tmp.sort() sortedaxnums, spos = zip(*tmp) slicepositions.append(spos) self.slicepositions = slicepositions self.sliceaxnums = sorted(sliceaxnums) def findmatchingaxis(self): """-------------------------------------------------------------------- Purpose: If the data, defined by the axes in 'axperm', has only one spatial axis, then coordinate transformations are impossible because a matching spatial axis is required. This routine checks whether a spatial axis is missing and returns the relevant information. Returns: ap- The original axis permutation array and if necessary extended with the axis number (0-based) of the missing spatial axis. mixax- The axis number of the missing spatial axis or None. -----------------------------------------------------------------------""" mixax = None setdim = self.naxis subdim = self.ndims(subset=True) axperm = self.axperm() settypes = self.projection.types # This is a list with axis types # ('longitude', 'latitude', 'spectral' etc.) ap = axperm[:subdim] # If a set has 2 or less axes, then there is no need to look for # a matching spatial axis because these are all the axes there are in a set. if setdim <= 2: return ap, mixax # Is a spatial axis missing in the sequence of subset axes subtypes = [] for i in axperm[:subdim]: # We need only the subset axis types subtypes.append(settypes[i]) lat = 'latitude' in subtypes lon = 'longitude' in subtypes if (lat and not lon) or (lon and not lat): # Then we need to find the missing axis if lat: missingaxis = 'longitude' else: missingaxis = 'latitude' try: # Have a look in the remaining types mixax = list(settypes).index(missingaxis) except ValueError: gipsy.error(4, "Cannot find a matching axis for the spatial axis!") else: # No need to find a missing spatial axis because this # data does not contain a spatial axis return ap, mixax # We are left with the last option i.e. we found a matching # spatial axis. We need to extend the axis permutation array ap.append(mixax) # Note that mixax numbering is 0 based. Not FITS! return ap, mixax def changeproj(self, skyout=None, spectra=None): #--------------------------------------------------------------------- """ Adjust sky system and/or spectral translation """ #--------------------------------------------------------------------- if spectra != None: self.subproj = self.projection.sub(self.projaxnum).spectra(spectra) if skyout != None: self.subproj.skyout = skyout def swapaxes(self): """-------------------------------------------------------------------- Purpose: Method to change the axes permutation array for two dimensional subsets only. Notes: Attribute projaxnum is a list with axis numbers in arbitrary order, which sets the projection object for the current subsets. Usually the order is the same as the order from the header of GDS set or FITS file, only for two dimensional subsets it makes sense to be able to swap axes. Note that this list can have more elements than the dimension of a subset because it can contain the number of the missing spatial axis. --------------------------------------------------------------------""" subdim = self.ndims(subset=True) if subdim != 2: raise Exception, "Cannot swap axes if subset dimension != 2" a = self.projaxnum[1] self.projaxnum[1] = self.projaxnum[0] self.projaxnum[0] = a # Recreate the projection object associated with the subset(s) self.subproj = self.projection.sub(self.projaxnum) # TODO add spectral translation and skyout. But, if this method is called # before any change of spectral translation or alternative sky definition # then there is no problem. def set_spectraltranslation(self, key=None, mes=None, default=1, inside=True): #--------------------------------------------------------------------- """ Set a spectral translation. Either for an position-velocity diagram or for a repeat axis with type 'spectral'. Method set_spectraltranslation Use: t,u = setobject.set_spectraltranslation( key, # Character string mes, # Character string default, # Integer inside, # Boolean) set_spectraltranslation Returns None, None if there is no possibility to return a spectral translation. Else, the spectral translation t as a string and also the unints are returned as a string. The output is minimal matched against the input (case insensitive). key Keyword. If nothing is entered then the default (SPECTRANS=) is used. mes Message for prompt default As in other input routines inside Boolean. If True then a translation is requested for one of the map axes and the projection object associated with the subset(s) is changed. Else a translation is requested for one of the repeat axes. Example: translation = gipsyset.set_spectraltranslation(key, mes, inside=True) """ #--------------------------------------------------------------------- # Conditions to ask a spectral translation: if inside: ok = not (self.transax_in is None) and len(self.subproj.altspec) > 0 if ok: altspecs = self.subproj.altspec transax = self.transax_in else: return None, None else: ok = not (self.transax_out is None) and len(self.projection.altspec) > 0 if ok: altspecs = self.projection.altspec transax = self.transax_out else: return None, None spectrans = None specunits = None translations, units = zip(*altspecs) # UNzip specaxname = self.projection.ctype[transax-1] ok = False if key == None: key = "SPECTRANS=" if mes == None: mes = "Spectral translation for axis %s or ? for help ... [native]:" % specaxname while not ok: st = gipsy.userchar(key, mes, defval="", default=default) if st == '?': if len(translations): gipsy.anyout("=========== Available spectral translations ===========", dev=1) n = 0 maxcol = 6 while n < len(translations): col = 0 trstr = "" while col < maxcol: if n < len(translations): if col > 0: trstr += ", " + translations[n] if units[n] != '': trstr += ' ('+units[n]+')' else: trstr += translations[n] if units[n] != '': trstr += ' ('+units[n]+')' n += 1 col += 1 gipsy.anyout(trstr, dev=1) gipsy.anyout("\n", dev=1) else: gipsy.anyout("No spectral translations found!", dev=1) ok = False gipsy.wkey(key) # First make empty gipsy.cancel(key) # Then cancel elif st != "": indx = positions.minmatch(st, translations) if indx == None: gipsy.reject(key, "Could not find a matching translation!" ) ok = False elif indx == -1: gipsy.reject(key, "Ambiguous input!") ok = False else: spectrans = translations[indx] specunits = units[indx] ok = True # Change projection object only if spectral axis is on of the subset axes. if inside: self.changeproj(spectra=spectrans) else: ok = True if ok: gipsy.wkey(key+st) # This prevents a semicolon in the keyword if first ? was used. else: default = 1 # You must be prompted return spectrans, specunits def set_skyout(self, key=None, mes=None, default=1): #--------------------------------------------------------------------- """ Set and return a sky definition for the axis labeling in spatial map. """ #--------------------------------------------------------------------- spatial = self.projection.lonaxnum in self.axperm_in or self.projection.lataxnum in self.axperm_in if not spatial: return None if key == None: key = "SKYOUT=" if mes == None: mes = "Enter sky definition (type ? for examples) ... [native]:" ok = False while not ok: skyout = gipsy.usertext(keyword=key, message=mes, defval=None, default=default) if skyout != None: if string_upper(skyout) == '?': gipsy.anyout("Sky systems: EQ=equatorial, EC=Ecliptic, GA=Galactic (II), SU=Supergalactic", dev=3) gipsy.anyout("Ref. systems: fk4 fk4_no_e,fk5, icrs", dev=3) gipsy.anyout(" FK4 mean place, old (Bessel-Newcomb) system", dev=3) gipsy.anyout(" FK4-NO-E mean place, old system but without e-terms", dev=3) gipsy.anyout(" FK5 mean place, new (IAU 1984) system", dev=3) gipsy.anyout(" ICRS International Celestial Reference System", dev=3) gipsy.anyout("Equinox, Epoch:", dev=3) s = """ B Besselian epoch. Example B 1950, b1950, B1983.5, -B1100 J Julian epoch. Example: j2000.7, J 2000, -j100.0 JD Julian date. This number of days (with decimals) that have elapsed since the initial epoch defined as noon Universal Time (UT) Monday, January 1, 4713 BC in the proleptic Julian calendar Example: JD2450123.7 MJD The Modified Julian Day (MJD) is the number of days that have elapsed since midnight at the beginning of Wednesday November 17, 1858. In terms of the Julian day: MJD = JD - 2400000.5 Example: mJD 24034, MJD50123.2 RJD The Reduced Julian Day (RJD): Julian date counted from nearly the same day as the MJD, but lacks the additional offset of 12 hours that MJD has. It therefore starts from the previous noon UT or TT, on Tuesday November 16, 1858. It is defined as: RJD = JD - 2400000 Example: rJD50123., Rjd 23433 F 1) DD/MM/YY Old FITS format Example: F29/11/57 2) YYYY-MM-DD FITS format Example: F2000-01-01 3) YYYY-MM-DDTHH:MM:SS FITS format with date and time. Example: F2002-04-04T09:42:42.1 """ gipsy.anyout(s, dev=3) gipsy.anyout("Example: SKYOUT=EQ FK4 J2000 B1950", dev=3) ok = False gipsy.wkey(key) gipsy.cancel(key) else: try: self.changeproj(skyout=skyout) ok = True except Exception, message: ok = False gipsy.reject(key, message) else: ok = True if ok: if skyout is None: gipsy.wkey(key) else: gipsy.wkey(key+skyout) else: default = 1 return skyout def printinfo(self): """-------------------------------------------------------------------- Purpose: Print a message in Hermes' log file The message contains information of the world coordinate system for all subset axes. For example: WCS Subset axis DEC--SIN of type 'latitude' (units: deg) WCS Subset axis FREQ of type 'spectral' (units: Hz) --------------------------------------------------------------------""" gipsy.anyout("\nWorld coordinate system for selected subset:") for axinfo in self.wcsinfo: # Print the WCS related axis information gipsy.anyout(axinfo) gipsy.anyout("\n") def wcspos(self, subnum, key="POS=", mes="", userpos=''): """-------------------------------------------------------------------- Purpose: Method for input of positions. Inputs: subnum- Index of current subset in subset arrray key- Keyword for prompt. If empty then the keyword becomes POS= mes- Message for prompt. If empty a default message is created. userpos- A string to replace the manual input. If not empty, then then no prompt will appear. Returns: grids- A NumPy array with grids world- A NumPy array with world coordinates in the system that is defined by the Projection object. This object can be modified by changing attribute *skyout* to set an alternative sky definition (sky system, reference system, equinox, epoch), or by setting a spectral translation with method *spectra()*, or by changing attribute *usedate* from False to True. Read more about *usedate* in the notes. pixels- A NumPy array with pixels. Pixels are FITS positions i.e. start with 1. They are needed if the result of the input needs to be converted to world coordinates with method toworld(). subsetunits- A list with units corresponding to the axes in the current subset. errmes- An error message if the input could not be processed into one or more positions. One should check the error message first before printing the output. Examples: Notes: 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 ('subnum'). If an error is made in a position then this routine will return with empty lists and a no empty error message. An object from class 'WCSset' inherits from class 'Set'. The most important difference is the addition of a Projection object from module 'wcs' of the Kapteyn Package. If your input has a sky reference system FK4 AND you transform to a non FK4 reference system (or to another sky system), then the date of observation could play a role in the transformation between grid/pixel- and world coordinates. The data of observation is taken into account if attribute 'usedate' of the Projection object is set to True (the default is False). If you transform from a non FK4 system to FK4, then the date of observation is included in the input of positions as in "{eq, B1950,fk4, J1983.5} 178.12830409 {} 53.93322241" independent of the value of attribute 'usedate'. -----------------------------------------------------------------------""" if self.subproj == None: # No subset available (i.e. a grid) return [], [], [], [],'', "No projection object available for this subset" # Create a default prompt if no string for 'mes' was given subdim = self.ndims(subset=True) # The subset dimension # The projection object associated with the subset contains the axes names in the right order. if mes == "": mes = "Enter position(s) in " for i in range(subdim): mes += self.subproj.ctype[i].split('-')[0] + ',' mes = mes.rstrip(' ,') mes += ' ... [Quit]:' subset = self.subsets[subnum] # Try to find the (FITS) pixel that corresponds # to the grid on the missing spatial axis, set by the CURRENT SUBSET. # Note that this pixel changes if the subset changes. # Example. If a data set has axes RA, DEC, FREQ and one entered # setname RA 0 1 2, then the subset is a DEC, FREQ plane. # The grids() method returns a list: grids = [0, None, None] # for the first subset, [1, None, None] for the second etc. # The number of elements is equal to the dimension of the data set. # The axes of the subset are set to None. The other axes are so # called 'repeat' axes the position of the subset corresponds to a grid # on these repeat axes. # The value of 'mixax' is always a number of an axis outside the subset. # and therefore one of the numbers from the list not equal to None. mixax = self.mixax crpix_mixax = None if mixax != None: # There is a missing spatial axis, number follows FITS definition ssgrids = self.grids(subset) # Find grids on axes outside subset mixgrid = ssgrids[mixax] # The grid on the missing spatial axis crpix_mixax = self.projection.crpix[self.mixax] mixpix = mixgrid mixpix += nint(crpix_mixax) # To get a pixel else: mixpix = None # Start to parse the input valid = False while not valid: errmess = '' if userpos != '': postxt = userpos else: postxt = gipsy.usertext(keyword=key, message=mes, defval="", default=1) if postxt == "": # Empty input --> return to calling environment return [], [], [], [], '' # Coordparser is a class that instantiates an object which either contains # one or more positions or an error message. It parses text representing a number # of coordinates. 'subdim' coordinates define 1 position. The text is parsed # according to the parser rules listed earlier. Each position is a list with # 'subdim' coordinates. A coordinate is a tuple consisting of a number, # a mode ('g' or 'w': grid or world coordinate), a sky system and a # spectral translation. The sky system and spectral translation modify the # projection object so that the transformations from pixels to world coordinates # v.v. take place in the modified system. The number could also be a # NumPy array with numbers read from a file. parsedpositions = positions.Coordparser(postxt, # Text containing positions as entered by user subdim, # The number of coordinates in 1 position self.subproj.units, # Units (for conversions) in order of subset axes self.subproj.types, # Axis types to distinguish spatial and spectral coords. self.subproj.crpix, # Crpix values for 'PC' (projection center) self.subproj.naxis, # Axis lengths for center 'AC' self.subproj.altspec,# List with allowed spectral translations self.subproj.source, # Header which was used to create object gipsygrids=True) # Used for internal conversion grids to pixels if len(parsedpositions.errmes) != 0: valid = False if userpos != '': return [], [], [], [], parsedpositions.errmes else: wor, pix, subsetunits, errmes = positions.dotrans(parsedpositions.positions, self.subproj, subdim, mixpix) if errmes != '': return [], [], [], [], errmes # Result is always pixels and we want an extra variable with grids grids = self.subproj.pixel2grid(pix) valid = True if not valid: gipsy.reject(key, parsedpositions.errmes) # Return the results as NumPy arrays because then in the calling environment # a user can easily extract one column. return grids, wor, pix, subsetunits, '' def setwcsbox(self, boxtxt, checklimits=True): """----------------------------------------------------------------------- Purpose: Set box in parameter boxtxt. If this fails, raise exception. Inputs: boxtxt- Box specification as a string checklimits- Boolean. If True, 'boxtxt' is not allowed to specify regions beyond the data regions. Returns: blo, bhi- Tuple with two lists. The first is the position of the lower boundary of the region. The second is the upper position of the region. For a two dimensional region one can interpret this as blo=(xlo,ylo) and bhi=(xhi,yhi) Examples: boxtxt = "-10 -10 20 30" boxtxt = "0 0 D 30 50" boxtxt = "2h10m30s 60d32m2.7s D 10 10" boxtxt = "2h10m30s 60d32m2.7s D 2 arcmin 3 arcmin" Notes: This routine is similar to 'wcsbox()', but then without the interaction with a user. So the parameter 'boxtxt' can be from an external source. If the parameter contains an invalid specification, there is no retry with a prompt for a new specification. Only an exception will be raised. -----------------------------------------------------------------------""" if self.subproj == None: # No subset available (i.e. a pixel) raise Exception, "No projection available" else: subproj = self.subproj subdim = self.ndims(subset=True) # Set default values for the limits defblo = [0.0]*subdim defbhi = [0.0]*subdim valid = False errmes = '?' # Get offsets to convert between (GIPSY) grids and (FITS) pixels offset = [0.0]*subdim for i in range(subdim): offset[i] = nint(self.subproj.crpix[i]) for i in range(subdim): d = nint(subproj.crpix[i]) defblo[i] = 1 - d defbhi[i] = subproj.naxis[i] - d if boxtxt == "": # Empty input --> return to calling environment # Parameter 'boxtxt' is empty. set the box to it limits. blo = defblo bhi = defbhi valid = True else: # Parameter 'boxtxt' is not empty? Then process it here. subset = self.subsets[0] blo = [0.0]*subdim bhi = [0.0]*subdim tokens = re_split('\s[dD]\s', boxtxt) # Split strings id 'D' is separator valid = True # At least the first one has to be a position # For a series of subsets, select the first subset for calculating positions. grids, world, pixels, units, errmes = self.wcspos(0, userpos=tokens[0]) valid = len(grids) == 2 and len(tokens) == 1 if valid: # Two positions. Sort and store the coordinates for i in range(subdim): blo[i] = min(grids[0][i], grids[1][i]) bhi[i] = max(grids[0][i], grids[1][i]) elif len(grids) == 0: valid = False errmes = "No positions found in: %s (%s)"%(tokens[0],errmes) elif len(grids) == 1 and len(tokens) == 1: valid = False errmes = "Only one position found" elif len(grids) > 2: # Always an lower and upper position of dimension 'subdim' valid = False errmes = "Too many positions for a box" elif len(grids) == 1 and len(tokens) == 2: # One position is found and the second token represents a size parsedpositions = positions.Coordparser(tokens[1], # Text containing positions as entered by user subdim, # The number of coordinates in 1 position self.subproj.units, # Units (for conversions) in order of subset axes self.subproj.types, # Axis types to distinguish spatial and spectral coords. self.subproj.crpix, # Crpix values for 'PC' (projection center) self.subproj.naxis, # Axis lengths for center 'AC' self.subproj.altspec,# List with allowed spectral translations self.subproj.source, # Header which was used to create object gipsygrids=False) # Sizes are always pixels if parsedpositions.errmes: valid = False errmes = parsedpositions.errmes else: valid = True p = parsedpositions.positions[0] for i in range(subdim): if p[i][1] == 'w': delta = subproj.cdelt[i] if delta != 0.0: size = abs(p[i][0][0] / delta) # Remember that p[i][0] is always a list else: valid = False errmes = "Cannot convert to grids because CDELT is 0" else: size = abs(p[i][0][0]) if valid: # Comparable to code in gdsbox() ci = nint(grids[0][i]) si = nint(size) halfsi = int(si)/2 hi = ci + halfsi lo = hi + 1 - si blo[i] = lo bhi[i] = hi if not valid: raise Exception, errmes if checklimits: # The caller required to check the limits. valid = True for i in range(subdim): valid = valid and blo[i] >= defblo[i] valid = valid and bhi[i] <= defbhi[i] if not valid: errmes = "Box outside limits of set" raise Exception, errmes # Update the box attributes for the Set object, then automatically the attributes # for the axes objects are correct because these are derived from these lists. self.axlimits = [[0.0, 0.0]]*subdim for i in range(subdim): self.axlimits[i] = [blo[i]+offset[i], bhi[i]+offset[i]] if subdim == 2 and self.projaxnum[0] > self.projaxnum[1]: # The axes are swapped, swap x and y of axes limits self.blo = [blo[1], blo[0]] self.bhi = [bhi[1], bhi[0]] else: self.blo = blo self.bhi = bhi return blo, bhi def wcsbox(self, key="BOX=", mes="", default=1, showdev=0, option=1, gui=False): """-------------------------------------------------------------------- Purpose: Function to set the grid limits of the axes in a subset. This function is similar to GIPSY's gdsbox() but it is based on WCSLIB and supports more projection- and sky systems. It needs the coordinate word of a subset to determine for which axes limits are needed. For these axes it creates a default prompt which can be overruled with keyword parameter 'mes'. A box is defined with two positions. The first position sets the lower limits of the axes and the second sets the upper limits. One position consists of n coordinates where n is equal to the dimensionality of the subset. Inputs: key- Keyword for prompt mes- Message for prompt. If empy then a default messafge is created. default- Set's the behaviour for defaults (as in the userxxxx routines). showdev- The destination of the informative messages. option- 0: No check on result 1: Check on upper limits for blo and bhi Returns: blo- bhi- Two lists with coordinates which set the axes limits of the input subset. The coordinates are sorted so that 'blo' represents the lower limits and 'bhi' the upper limits. Note that also the Set object (gset) is updated. Therefore you can retrieve the limits also by accessing attributes gset.blo, gset.bhi Attributes: blo, bhi- See above axlimits- Convenience attribute. It is a FITS representation of the box (i.e. not grids but pixels and first pixel is 1). Note that for 2-dim structures, the axes order could have been changed. This attributes follow the order as entered in the input, while attributes blo and bhi follow the GIPSY standard, i.e. the order is the same as the axes in the original data. Notes: 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). Here we follow the implementation in gdsbox. That is, that only the first subset is used for coordinate transformations. Examples: BOX=-10 -10 20 30 BOX= 0 0 D 30 50 BOX= 2h10m30s 60d32m2.7s D 10 10 BOX= 2h10m30s 60d32m2.7s D 2 arcmin 3 arcmin -----------------------------------------------------------------------""" if self.subproj == None: # No subset available (i.e. a pixel) return [], [], "No projection object available for this subset" else: subproj = self.subproj subdim = self.ndims(subset=True) # Set default values for the limits defblo = [0.0]*subdim defbhi = [0.0]*subdim # Create a default prompt if no string for 'mes' was given if mes == "": mes = "Enter BOX in " for i in range(subdim): mes += subproj.ctype[i].split('-')[0] + ',' mes = mes.rstrip(' ,') # The last comma is one to many blomes = '' bhimes = '' for i in range(subdim): d = nint(subproj.crpix[i]) defblo[i] = 1 - d blomes += "%d "%defblo[i] defbhi[i] = subproj.naxis[i] - d bhimes += "%d "%defbhi[i] blomes = blomes.rstrip() bhimes = bhimes.rstrip() mes += ' ... [%s, %s]:'% (blomes, bhimes) offset = [0.0]*subdim for i in range(subdim): offset[i] = nint(self.subproj.crpix[i]) # Prepare axis limits (box) subset = self.subsets[0] valid = False while not valid: errmess = '' boxtxt = gipsy.usertext(keyword=key, message=mes, defval="", default=default) if boxtxt == "": # Empty input --> return to calling environment blo = defblo bhi = defbhi valid = True else: blo = [0.0]*subdim bhi = [0.0]*subdim tokens = re_split('\s[dD]\s', boxtxt) valid = True # At least the first one has to be a position # For a series of subsets, select the first subset for calculating positions. grids, world, pixels, units, errmes = self.wcspos(0, userpos=tokens[0]) valid = len(grids) == 2 and len(tokens) == 1 if valid: # Two positions. Sort and store the coordinates for i in range(subdim): blo[i] = min(grids[0][i], grids[1][i]) bhi[i] = max(grids[0][i], grids[1][i]) elif len(grids) == 0: valid = False errmes = "No positions found" elif len(grids) == 1 and len(tokens) == 1: valid = False errmes = "Only one position found" elif len(grids) > 2: valid = False errmes = "Too many positions for a box" elif len(grids) == 1 and len(tokens) == 2: # One position is found and the second token represents a size # The trick is to calculate the second token as a position # wrt. (grid 0,0). Then, if necessary, divide by the # pixel size in CDELT to get the lengths in pixels wrt. (0,0) parsedpositions = positions.Coordparser(tokens[1], # Text containing positions as entered by user subdim, # The number of coordinates in 1 position self.subproj.units, # Units (for conversions) in order of subset axes self.subproj.types, # Axis types to distinguish spatial and spectral coords. self.subproj.crpix, # Crpix values for 'PC' (projection center) self.subproj.naxis, # Axis lengths for center 'AC' self.subproj.altspec,# List with allowed spectral translations self.subproj.source, # Header which was used to create object gipsygrids=False) # Size is same for GIPSY and FITS if parsedpositions.errmes: valid = False errmes = parsedpositions.errmes else: valid = True p = parsedpositions.positions[0] for i in range(subdim): if p[i][1] == 'w': delta = subproj.cdelt[i] if delta != 0.0: size = abs(p[i][0][0] / delta) # Remember that p[i][0] is always a list else: valid = False errmes = "Cannot convert to grids because cdelt is 0" else: size = abs(p[i][0][0]) if valid: # Comparable to code in gdsbox() ci = nint(grids[0][i]) si = nint(size) halfsi = int(si)/2 hi = ci + halfsi lo = hi + 1 - si blo[i] = lo bhi[i] = hi if option: # The caller required to check the limits. valid = True for i in range(subdim): valid = valid and blo[i] >= defblo[i] valid = valid and bhi[i] <= defbhi[i] if not valid: errmes = "Box outside limits of set" if not valid: gipsy.reject(key, errmes) if gui: return None, None # Write an informative message in the log/screen file (depends on 'showdev') header = "\nGrid range(s) for set %s :" % self.dataname gipsy.anyout(header, dev=showdev) for i in range(subdim): axname = subproj.ctype[i] # Long name from (alternative) header boxinfo = "%-18s from %6d to %6d" % (axname, blo[i], bhi[i]) gipsy.anyout(boxinfo, dev=showdev) gipsy.anyout(" ", dev=showdev) # Update the box attributes for the Set object, then automatically the attributes # for the axes objects are correct because these are derived from these lists. self.axlimits = [[0.0, 0.0]]*subdim #self.wcsblo = blo[:] #self.wcsbhi = bhi[:] for i in range(subdim): self.axlimits[i] = [blo[i]+offset[i], bhi[i]+offset[i]] if subdim == 2 and self.projaxnum[0] > self.projaxnum[1]: # The axes are swapped self.blo = [blo[1], blo[0]] self.bhi = [bhi[1], bhi[0]] else: self.blo = blo self.bhi = bhi return blo, bhi def listfiles(): #----------------------------------------------------------------------- """ List files related to GDS, FITS or zipped FITS in current directory """ #----------------------------------------------------------------------- li = [] filenames = listdir(".") for fname in filenames: if fnmatch(fname, '*.descr') or\ fnmatch(fname, '*.FITS') or\ fnmatch(fname, '*.fits') or\ fnmatch(fname, '*.fits.gz'): li.append(fname) li.sort() for fname in li: if fnmatch(fname, '*.descr'): try: basename, extension = path.splitext(fname) fname = basename + '.image' size = stat(fname)[6] / 1024.0 / 1024.0 sname = basename except: sname = '' else: sname = fname if sname: size = stat(fname)[6] / 1024.0 / 1024.0 gipsy.anyout("%-30s %6.2f Mb" % (sname, size), dev=3) def wcsinp(key="INSET=", mes="Enter data file:", create=False, write=False, hdukey="HDU=", hdumes="Enter number of Header Data Unit ... [0]:", altkey="ALTHEAD=", altmes=None, axkey="AXNAMES=", axmes=None, axsuffix=None, swapkey="AXSWAP=", swapmes=None, requireddim=0, listopt=True): """-------------------------------------------------------------------- Purpose: Routine for input of set/subsets, either from a GIPSY set or a FITS file. The parsing of a FITS file involves module PyFITS. Therefore we can also enter FITS data using a location on the web (url) or read data from a compressed file (.gz). If a FITS file is entered and the FITS file has multiple header Data Units (hdu's) then the user is prompted to give the number of the hdu that contains the wanted data. Inputs: key- Keyword for prompt mes- Message for prompt. This routine appends the string '[list files]' as the default action. When a user presses 'Enter' then a list with filenames is displayed. Keyword key appears again but now the default has been changed to abort the loop of asking a valid file to open and the string '[abort]' is appended. This last prompt appears until a valid name is entered or if you use the default (which aborts the loop). hdukey- Keyword for prompt of Header data unit hdumes- Message for prompt of Header data unit altkey- Keyword for prompt for a selection between alternate headers altmes- Message for 'altkey' prompt. axsuffix- A string to append to the axis names if a user needs to specify slices. For example if a 4-dim set is entered while a 2-dim subset is required, then you will be prompted to enter repeat axes and grids on the repeat axes (these are axes that do not belong to the subset). requireddim- Required (minimum) dimension of the specified subsets. If the dimension of the structure entered at the prompt is greater than the required dimension, then one should specify grids on the repeat axes. This function creates a prompt for each repeat axis. The default keyword can be appended by a user given suffix. Returns: An object of class WCSSet which inherits from Set. So it has all attributes from a object from class Set and some additional attributes related to WCS. Notes: If the constructor of the Set objects encounters multiple header data units, then a callback is triggered. The callback function deals with the user interaction. It returns a integer number. The same construction applies if an alternative header is encountered, Then a callback function returns a character for the alternative header or ' ' if the primary header must be selected. Examples: INSET=aurora freq 1:10 INSET=ngc1234 stokes 1 freq INSET=http://www.atnf.csiro.au/people/mcalabre/data/WCS/1904-66_ZPN.fits.gz INSET=ngcxxx.fits.gz velo 50 -----------------------------------------------------------------------""" validset = False swapaxes = False listed = False if listopt: mes2 = mes + " [list files]" else: mes2 = mes while not validset: datafile = gipsy.usertext(keyword=key, message=mes2, default=1, defval=None) if datafile == None: # Two options: Either the user wants a list (listopt = True) # or user wants nothing returned if listopt and not listed: listfiles() # User wants list with files listed = True gipsy.cancel(key) mes2 = mes + " [abort]" # Change default to abort input loop else: return None else: # If a user enters more than 1 string than this extra string # is part of a subset specification. We will use this property later. prespecified = len(datafile.split()) > 1 # Both a GDS set or a FITS file can be opened. A FITS file is # converted to a GDS set. try: # The prompthdu function is only active when a FITS file is detected # and the file has more than 1 header data unit (hdu) # The hdu number is stored in attribute 'hdunum' # The alternative header is stored in attribute 'altwcs' gipsyset = gipsy.Set(datafile, create=create, write=write, gethdu=prompthdu(hdukey, hdumes), getalt=promptalt(altkey, altmes)) validset = True except gipsy.GipsyError, message: validset = False errmes = str(message) #"Cannot open or process this set!" if validset: setdim = gipsyset.ndims(0) subdim = gipsyset.ndims(subset=True) if requireddim != 0: # Only if 0 the user is free to select validset = True if subdim == requireddim or subdim == setdim: validset = True else: validset = False errmes = "Wrong dimensionality" if requireddim > setdim: validset = False errmes = "Required dim. > dim. of set" if not validset: gipsy.reject(key, errmes) else: subsetspec = datafile # We always end up with a primary header because a selected alternative header # is migrated to a primary header if gipsyset.fits: name = gipsyset.fitsname else: name = gipsyset.name mes = "Set %s has %d axes:" % (name, setdim) gipsy.anyout(mes) for i in range(setdim): # Print some information about the axes crpix = gipsyset['CRPIX%d'%(i+1)] axlen = gipsyset['NAXIS%d'%(i+1)] axname = gipsyset['CTYPE%d'%(i+1)] # Follow the documented definition of the lower and upper axis limit: gipsy.anyout("%-20s from %6d to %6d"%(axname, 1-nint(crpix), axlen-nint(crpix))) # For users that work with a set with unknown axis names it is possible to # enter a set name only. Then prompt user to enter the axes of # the wanted data structure (subset). The keywords are the names of the # repeat axes (with a suffix, if given). if not prespecified: # User entered a name only if subdim > requireddim: # But we need to specify the subset axes # Prompt user to enter 'requireddim' names from 'setdim' axes if requireddim > 0: mes = "Enter %d names of " % requireddim else: mes = "Enter axes from " ax = gipsyset.axes(False) # Set method to get information about all axes for i in range(setdim): mes += ax[i].name if i < setdim - 1: mes += ',' # Build a default list with axis names and a message for the prompt deflist = " [" defval = [] if requireddim == 0: maxaxes = setdim else: maxaxes = requireddim for i in range(maxaxes): deflist += ax[i].name defval.append(ax[i].name) if i < maxaxes - 1: deflist += ',' deflist += ']:' mes += deflist # Ask user the name of the subset axes until a correct answer is entered valid = False nvalids = 0 validaxnums = [] while not valid: if requireddim == 0: default = 1 else: default = 4+1 axes = gipsy.userchar(axkey, mes, default=default, defval=defval, nmax=maxaxes) # 1=default, 4=exact for axselection in axes: for i in range(setdim): axnr = i + 1 axname = ax[i].name if axname.find(string_upper(axselection), 0, len(axselection)) > -1: # A valid axis name. Store the corresponding number nvalids += 1 validaxnums.append(axnr) # Make a sorted list with unique elements only. Then check if there are enough axes. sortedaxnums = validaxnums[:] uset = {} map(uset.__setitem__, validaxnums, []) newlist = uset.keys() # Is a copy, so validaxnums is unchanged! if requireddim in [0,2] and len(newlist) == 2: # In this situation (user specifies the axes of the subset), # one could have swapped the axes (Which the Set class does not # support). Then swap axes later. if validaxnums[0] > validaxnums[1]: swapaxes = True if requireddim > 0: valid = len(newlist) == requireddim else: valid = True if not valid: # No exceptions because we are in a loop gipsy.reject(axkey, "Incorrect input of image axes!") # At this point we have a set name, a list with axes, a hdu and an alt. header # so that we can make a reconstruction of the input in terms of the Set object # First retrieve grids on axes outside the subset if there are more # axes than the 'requireddim' axes. # For the keywords we use the axis names. Note that also for alternative headers # the axis names in the axes() method are correct because they were ectracted # from the selected header. if len(newlist) < setdim: grids = [] if gipsyset.fits: # One cannot use attribute gipsyset.name because for FITS files the url is stripped subsetspec = gipsyset.fitsname + " " else: subsetspec = gipsyset.name + " " axo = gipsyset.axes(False) for ax in axo: if ax.axnum+1 not in validaxnums: key = "%s" % ax.name if axsuffix != None: key += axsuffix # Add a user supplied suffix to distinguish from existing keywords # gipsy.anyout("name lo, hi, crpix %s %d %d %d"%(ax.name,ax.slo, ax.shi, ax.crpix)) unvalid = True while unvalid: defval = 0 if not (ax.slo <= defval <= ax.shi): defval = int((ax.shi+ax.slo)/2) prompt = "Enter grid position(s) on %s ... [%d:%d]:" % (ax.name, ax.slo, ax.shi) nmax = ax.shi - ax.slo + 1 # 'nmax' could be 1. Then a scalar is returned not a list. # With the next construction we force returning a list in # all situations. Also it allows a user to read the same subset # more than once. gr = gipsy.userint(key, prompt, defval=None, default=1, nmax=max(128,nmax)) outside = False # If numbers are entered then check them all if (gr != None): for g in gr: outside = not (g >= ax.slo and g <= ax.shi) if outside: break unvalid = outside if unvalid: gipsy.reject(key, "Grid not in range %d to %d! Try again." % (ax.slo, ax.shi)) else: # Copy the keyword contents to a string using the previous input gr = gipsy.usertext(key, "", defval="", default=1) subsetspec += ax.name + " " + gr + " " # Build the subset specification string. grids.append(gr) altwcs = gipsyset.altwcs # Store attributes before deletion hdunum = gipsyset.hdunum gipsyset.close() # Close set before deletion del gipsyset # Remove old one from memory first # Recreate the Set object but now as an object from class WCSset which inherits # from the Set class. Then it get specific attributes related to the world # coordinate system. stdout.flush() gipsyset = WCSset(subsetspec, create=create, write=write, gethdu=lambda x: hdunum, getalt=lambda x,y: altwcs) # There are more input options for which one possibly wants to be able # to swap axes. # 1) There was a pre-specification of the data set and the resulting # subset dimension is two. # 2) There was no pre specification and the set dimension of the # data set was 2 while the required dimension was also 2 if (prespecified and subdim == 2) or (not prespecified and requireddim == 2 and setdim == 2): if swapmes is None: ax = gipsyset.axes(level=gipsyset.subsets[0]) #, mode='inside') swapmes = "Do you want to swap %s and %s? ... Y/[N]"%(ax[0].name, ax[1].name) swapaxes = gipsy.userlog(swapkey, swapmes, defval=False, default=1) if swapaxes: # Could also have been set before (look for swapaxes) gipsyset.swapaxes() return gipsyset def dotest(): """-------------------------------------------------------------------- Purpose: Test routine for the parser in this module. It works only for 2 dim subsets. Input: Object from class WCSset -----------------------------------------------------------------------""" # First test WCSset. Create a minimal FITS file and test what we can do with it. 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= test123wcsRADECFREQ" gipsy.xeq(command) gipsyset = WCSset("test123wcsRADECFREQ freq 0") gipsyset.printinfo() # Create an Ascii file with dummy data for testing FILE command # The data in the Ascii file is composed of a fixed sequence of grids # that are transformed to their corresponding galactic coordinates. asciifile = "test123wcsRADECFREQ.txt" f = open(asciifile, 'w') gipsy.anyout("Content test file on disk: %s" % asciifile) s = "! Test file for Ascii data and the FILE command\n" f.write(s) gipsy.anyout(s) upos = "[0:9] [0:27:3]" oldsky = gipsyset.subproj.skyout gipsyset.subproj.skyout = "gal" gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) for i in range(10): s = "%f %.12f %f %.12f\n" % (gr[i][0], wor[i][0], gr[i][1], wor[i][1] ) f.write(s) gipsy.anyout(s) f.close() gipsyset.subproj.skyout = oldsky gipsy.anyout("\n\n=========== Testing numbers ==========") gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="3/4 3**2") gipsy.anyout("3/4=%f 3**2=%f" % (gr[0][0], gr[0][1])) gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="2.1234e+2 2.1234E-2") gipsy.anyout("2.1234e+2=%f 2.1234E-2=%f" % (gr[0][0], gr[0][1])) gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="(3*4)-5 1/5*(7-2)") gipsy.anyout("(3*4)-5=%f 1/5*(7-2)=%f" % (gr[0][0], gr[0][1])) # Numbers as described in http://www.astro.rug.nl/~gipsy/hermes/numbers.html#functions gipsy.anyout("\n") gipsy.anyout("=========== Testing constants ==========") gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="pi c_") gipsy.anyout("pi=%f c_=%f" % (gr[0][0], gr[0][1])) gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="h_ k_") gipsy.anyout("Too small for coordinate parsing: Planck h_=%.2g Boltzmann k_=%.2g" % (gr[0][0], gr[0][1])) gipsy.anyout("because crpix is added and crpix >> h_ or k_") gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="G_ s_") gipsy.anyout("Gravitation G_=%g Stefan-Boltzman s=%g" % (gr[0][0], gr[0][1])) gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos="M_ P_") gipsy.anyout("Solar mass M_=%f 1 parsec, P_=%f" % (gr[0][0], gr[0][1])) gipsy.anyout("\n") gipsy.anyout("=========== Testing special symbols and lists ==========") userpos = ["[log(10)]*4, [exp(2)]*4", "[log(10):log(100):2/4], 10**[0, 1, 2]", # yields 1.0 1.5 2.0 and 1 10 100 "[1:3]+100, range(3)", # yields 101 102 103 and 0 1 2 "sin([0, pi, pi/2]) [0, 1, 2]" ] for upos in userpos: gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos=upos) if errm: gipsy.anyout(errmes) else: gipsy.anyout("%s =\n %s" % (upos, array2string(gr, precision=1))) gipsy.anyout("\n") gipsy.anyout("=========== Testing grouping with quotes ==========") userpos = ["0 1 2 3 10 11 12 13", "'0, 1, 2, 3' '10, 11, 12, 13'", "'G_, 1, 2, 3' '10, 11, 12, 13'", "[G_*300000000, 1, 2, 3] [10, 11, 12, 13]", "{eq} '178.779200, 178.778200, 178.777200' {} '53.655000, 53.656000, 53.657000'", "{eq} [178.779200, 178.778200, 178.777200] {} [53.655000, 53.656000, 53.657000]", "'178.779200, 178.778200, 178.777200' deg '53.655000, 53.656000, 53.657000' deg", #'"178.779200 178.778200 178.777200" deg "53.655000 53.656000 53.657000" deg', ] for upos in userpos: gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos=upos) gipsy.anyout("%s =\n %s" % (upos, array2string(gr, precision=5))) gipsy.anyout("\n") gipsy.anyout("=========== Testing expressions ==========") userpos = ["abs(-10), sqrt(3)", "sin(radians(30)), degrees(asin(0.5))", "cos(radians(60)), degrees(acos(0.5))", "tan(radians(45)), degrees(atan(1.0))", "tanh(1.0), atan2(3, 4)", # Extra space in atan2, two arguments for function #"erf(0.1), erfc(0.1)", # Is a scipy.special function. Not available here "max([2,3]) min([4,10.3])", # no comma between coordinates. Note [] is required "sinc(0.0) sign(-4.4)", "randn(2), ranf(2)" ] for upos in userpos: gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos=upos) gipsy.anyout("%s = %f %f" % (upos, gr[0][0], gr[0][1])) gipsy.anyout("\n") gipsy.anyout("\n=========== Testing the native function readcol() for reading files ==========") gipsy.anyout("Read columns from an Ascii file on disk with the native readcol()") gipsy.anyout("function. Note the difference with GIPSY's file() function") userpos = ['readcol("test123wcsRADECFREQ.txt", 1) readcol("test123wcsRADECFREQ.txt", 3)', 'readcol("test123wcsRADECFREQ.txt", 1) readcol("test123wcsRADECFREQ.txt", 3)', 'readcol("test123wcsRADECFREQ.txt", col=1, rows=range(1,4)) [3:5]', 'readcol("test123wcsRADECFREQ.txt", 1, rows=range(4)) readcol("test123wcsRADECFREQ.txt", 3, rows=[0,1,2,3])', 'readcol("test123wcsRADECFREQ.txt", col=1, fromline=3, toline=5) readcol("test123wcsRADECFREQ.txt", col=3, fromline=3, toline=5)', 'readcol("test123wcsRADECFREQ.txt", 1, fromrow=3) readcol("test123wcsRADECFREQ.txt", 3, fromrow=3)', '{ga} readcol("test123wcsRADECFREQ.txt", 1) {} readcol("test123wcsRADECFREQ.txt", 3)' ] for upos in userpos: gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos=upos) if errm: gipsy.anyout("Error: %s"%errm) else: gipsy.anyout("%s = %s" % (upos, array2string(gr))) gipsy.anyout("\n") gipsy.anyout("=========== Testing sky systems for set: %s Output in native world coordinates ==========" % gipsyset.name) userpos = ["0 0", "eq 178.7792 eq 53.655", # e 10 will not work because e is not a symbol and an ambiguous sky system` "{eq} 178.7792 {} 53.655", "178.7792 deg 53.655 deg", "11h55m07.008s 53d39m18.0s", "11h55m07.008 53d39m18.0s", "178.779200d0m 53d39m18.0s", "{eq, B1950,fk4} 178.7792 {} 53.655", "{eq, B1950,fk4} 178.12830409 {} 53.93322241", "{eq, B1950,fk4, J1983.5} 178.12830409 {} 53.93322241", "ga 140.52382927 ga 61.50745891", "su 61.4767412, su 4.0520188", "su 61.4767412, {} 4.0520188", "ec 150.73844942 ec 47.22071243", "{} 178.7792 0.0", "0.0, {} 53.655", # "0.0 ga 61.50745891", are mixed and not native -> invalid input 'header("crpix1"), header("crpix2")', "ac", "pc" ] for upos in userpos: gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) #gipsy.anyout("errmes=%s"%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])) gipsy.anyout("\n") gipsy.anyout("=========== Testing sky systems for set: %s Output in {eq, B1950,fk4, J1983.5} coords.==========" % gipsyset.name) userpos = ["0 0", "eq 178.7792 eq 53.655", # e 10 will not work because e is not a symbol and an ambiguous sky system` "{eq} 178.7792 {} 53.655", "178.7792 deg 53.655 deg", "11h55m07.008s 53d39m18.0s", "{eq, B1950,fk4} 178.7792 {} 53.655", "{eq, B1950,fk4} 178.12830409 {} 53.93322241", "{eq, B1950,fk4, J1983.5} 178.12830409 {} 53.93322241", "ga 140.52382927 ga 61.50745891", "su 61.4767412, su 4.0520188", "ec 150.73844942 ec 47.22071243", "{} 178.7792 0.0", "0.0, {} 53.655", 'header("crpix1"), header("crpix2")' # "0.0 ga 61.50745891", are mixed and not native -> invalid input ] oldsky = gipsyset.subproj.skyout gipsyset.subproj.skyout = "{eq, B1950,fk4, J1983.5}" for upos in userpos: gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) 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])) gipsyset.subproj.skyout = oldsky # Restore # We re-open a set but now it wil be a XV map (DEC,FREQ) del gipsyset gipsyset = WCSset("test123wcsRADECFREQ RA 0") gipsyset.printinfo() gipsy.anyout("\n") gipsy.anyout("=========== Testing mixed systems for set: %s ==========" % gipsyset.name) gipsy.anyout("Notes: The header is a legacy header spectral type FREQ-OHEL") gipsy.anyout("This implies that the supplied velocity in the header (DRVAL3 or VELR) is an") gipsy.anyout("optical velocity (1.05000000e+06 m/s)") userpos = ["0 0", "{} 53.655 1.415418199417E+09 hz", "{} 53.655 1.415418199417E+03 Mhz", "{} 53.655 1.415418199417 Ghz", "{} 53.655 vopt 1.05000000e+06", "{} 53.655 , vopt 1050000 m/s", "0.0 , vopt 1050000 m/s", "10.0 , vopt 1050000 m/s", "{} 53.655 vrad 1.05000000e+06", "{} 53.655 FREQ 1.41541820e+09", ] for upos in userpos: gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) 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])) gipsy.anyout("\n") gipsy.anyout("=========== Testing mixed systems for set: %s [Output in VOPT-F2W] ==========" % gipsyset.name) gipsy.anyout("Notes: The header is a legacy header spectral type FREQ-OHEL") gipsy.anyout("This implies that the supplied velocity in the header (DRVAL3 or VELR) is an") gipsy.anyout("optical velocity (1050 km/s). Note that the spectral translation 'FREQ' is not") gipsy.anyout("equal to the native system (FREQ-OHEL) because FREQ means the barycentric frequency") gipsy.anyout("converted from the topocentric frequency.") userpos = ["0 0", "{} 53.655 freq 1.415418199417E+09 hz", "{} 53.655 fr 1.415418199417E+03 Mhz", "{} 53.655 fr 1.415418199417 Ghz", "{} 53.655 vopt 1.05000000e+06", "{} 53.655 , vopt 1050000 m/s", "0.0 , vopt 1050000 m/s", "10.0 , vopt 1050000 m/s", "{} 53.655 vrad 1.05000000e+06", "{} 53.655 FREQ 1.415418199417e+09", ] oldsubproj = gipsyset.subproj gipsyset.subproj = oldsubproj.spectra('VOPT-???') for upos in userpos: gr, wor, pix, units, errm = gipsyset.wcspos(0, userpos=upos) if errm != '': gipsy.anyout(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])) gipsyset.subproj = oldsubproj # Restore gipsy.anyout("\n=========== Testing errors and error messages ==========") gipsy.anyout("These are expressions with various errors") userpos = ['readcol("test123wcsRADECFREQ.txt, 0) readcol("test123wcsRADECFREQ.txt", 2)', 'readcol("test123wcsRADECFREQ.txt", 0, range(1:4)) 3:5', 'readcol("test123wcsRADECFREQ.txt", col=0, rows=range(3)) 3::5', 'readcol("test123wcsRADECFREQ.txt", 2, rows=[0,1,2,3])', 'readcol("test123wcsRADECFREQ.txt", 0, rows=range(1,5)) 3:5', 'readcol("testwcsRADECFREQ.txt", 0, rows=range(6)) readcol("test123wcsRADECFREQ.txt", 2, rows=[0,1,2,3,4,5])', 'readcol("test123wcsRADECFREQ.txt", 0, rowsslice(5,None)) readcol("test123wcsRADECFREQ.txt", 2, rowslice=(5,None))', 'readcol("test123wcsRADECFREQ.txt", 0, row=2) readcol("test123wcsRADECFREQ.txt", 2, row=2)', '{ga} readcol("test123wcsRADECFREQ.txt", 1) {} readcol("test123wcsRADECFREQ".txt, 3)' ] for upos in userpos: gr, wor, pix, units, errm= gipsyset.wcspos(0, userpos=upos) if errm: gipsy.anyout("%s"%errm) else: gipsy.anyout("%s = %s" % (upos, array2string(gr))) # Test the three main routines (start as application with ioutils.py) # These functions are: wcsinp, wcsbox and wcspos if __name__ == "__main__": gipsy.init() testing = gipsy.userlog("RUNTEST=", "Run test ... Y/[N]:", defval=False, default=1) if testing: dotest() gipsy.finis() dim = 2 dim = gipsy.userint("DIM=", "Required subset dimension ... [%d]:"%dim, defval=dim, default=1) gipsyset = wcsinp(requireddim=dim) if gipsyset == None: gipsy.finis() gipsyset.printinfo() key = "BOX=" blo, bhi = gipsyset.wcsbox(key, "", option=1) # Default output world coordinates in wcspos are in the native system # unless user specifies otherwise. key = "SKYOUT=" mes = None gipsyset.set_skyout(key, mes) # Same for the spectral coordinate key = "SPECTRANS=" mes = None gipsyset.set_spectraltranslation(key, mes, inside=True) # Loop for positions. Press enter to abort numsubsets = len(gipsyset.subsets) if numsubsets > 1: s = "I have %d subsets, but will display the results only for subset number 0!"%numsubsets gipsy.anyout(s) key = "POS=" while 1: grids, world, pix, units, errmes = gipsyset.wcspos(0, key, "") gipsy.cancel(key) if not len(grids) or errmes: if errmes: gipsy.reject(key, errmes) else: break else: gipsy.anyout("Grid coordinates: ") for g in grids: gipsy.anyout(str(g)) gipsy.anyout("World coordinates:") for i, w in enumerate(world): s = str(w) + ' ' + str(units) gipsy.anyout(s) gipsy.finis()