Examples

In general, these examples show the Python equivalent of the corresponding example in the "How To..." section of the documentation for the native C version of AST, SUN/211, which should be consulted for further information on each example.

...Read a WCS Calibration from a Dataset

import pyfits
import starlink.Grf as Grf
import starlink.Ast as Ast
import starlink.Atl as Atl
import matplotlib.pyplot

#  Open the FITS file using (for instance) pyfits. A list of the HDUs
#  in the FITS file is returned.
hdu_list = pyfits.open( 'test.fit' )

#  Create a FitsChan and fill it with the FITS header cards in the
#  primary HDU (element zero of the list returned by pyfits.open).
fitschan = Ast.FitsChan( Atl.PyFITSAdapter(hdu_list[ 0 ]) )

#  Note which encoding has been used for the WCS information. This is only
#  necessary if you later wish to write modified headers out using the
#  same encoding. See section "...Write a Modified WCS Calibration to a
#  Dataset" below.
encoding = fitschan.Encoding

#  Read WCS information from the FitsChan. Additionally, this removes all
#  WCS information from the FitsChan.
wcsinfo = fitschan.read()

...Validate WCS Information

#  Check that the FITS header contained WCS in a form that can be
#  understood by AST.
if wcsinfo == None:
   print( "Failed to read WCS information from test.fit" )

#  Check that the object read from the FitsChan is of the expected class
#  (Ast.FrameSet).
elif not isinstance( wcsinfo, Ast.FrameSet ):
   print( "A "+wcsinfo.__class__.__name__+" was read from test.fit - "
          "was expecting an Ast.FrameSet" )

#  This particular code example can only handle WCS with 2 pixel axes
#  (given by Nin) and 2 WCS axes (given by Nout).
elif wcsinfo.Nin != 2 or wcsinfo.Nout != 2:
   print( "The world coordinate system is not 2-dimensional" )

#  Proceed if the WCS information was read OK.
else:
   ...

...Display AST Data

#  Display the internal representation of the AST FrameSet
   print( wcsinfo )

#  or alternatively...
   wcsinfo.show()

...Convert Between Pixel and World Coordinates

#  Store the x and y pixel coordinates at 3 points in the FITS image
   xpixel = [ 120., 150., 180. ]
   ypixel = [ 175., 150., 125. ]

#  Convert to world coordinates.
   ( xworld, yworld ) = wcsinfo.tran( [ xpixel, ypixel ] )

#  Convert back to pixel coordinates.
   ( xpix, ypix ) = wcsinfo.tran( [ xworld, yworld ], False )

...Test if a WCS is a Celestial Coordinate System

#  Obtain a pointer to the current Frame and determine if it is a SkyFrame.
#  The frame index argument could have been omitted in this case since it
#  defaults to Ast.CURRENT.
   frame = wcsinfo.getframe( Ast.CURRENT )
   issky = isinstance( frame, Ast.SkyFrame )

#  Or equivalently.
   issky = frame.isaskyframe()

...Format Coordinates for Display

# Loop round all the positions to be formatted
   for (x,y) in zip( xworld, yworld):

#  Get the formatted value for the first WCS axis, and append the units
#  string. Then do the same for the second WCS axis.
      xtext = wcsinfo.format( 1, x ) + " " + wcsinfo.Unit_1
      ytext = wcsinfo.format( 2, y ) + " " + wcsinfo.Unit_2

#  Print both values.
      print( "Position = "+xtext+", "+ytext )

...Display Coordinates as they are Transformed

#  Tell the Mapping to report each position as it is transformed.
   wcsinfo.Report = True

#  Convert to world coordinates.
   ( xworld, yworld ) = wcsinfo.tran( [ xpixel, ypixel ] )

#  Convert back to pixel coordinates.
   ( xpix, ypix ) = wcsinfo.tran( [ xworld, yworld ], False )

...Read Coordinates Entered by a User

#  Obtain the number of coordinate axes (if not already known).
#  The values of Naxes and Nout are equal for a FrameSet.
   naxes = wcsinfo.Naxes

#  Loop to read each line of input text, in this case from the
#  standard input stream.
   pos = []
   while True:
      text = input( "Coords: " )
      if text == "":
         break

#  Attempt to read a coordinate for each axis.
      start = 0
      end = len( text )
      for iaxis in range( 1, naxes + 1 ):
         ( n, axis_value ) = wcsinfo.unformat( iaxis, text[ start: ] )

#  If nothing was read and this is not the first axis or the end-of-string,
#  try stepping over a separator and reading again.
         if n == 0 and ( iaxis > 1 ) and ( start < end ):
            start += 1
            ( n, axis_value ) = wcsinfo.unformat( iaxis, text[ start: ] )

#  Quit if nothing was read, otherwise move on to the next coordinate. */
         if n == 0 :
            break
         start += n

#  Store axis value
         pos.append( axis_value )

#  Error in input data at character text[n]. */
      if start < end or n == 0:
         break

#  Coordinates were read OK.
      else:
         print( pos )
         pos = []

...Modify a WCS Calibration

#  The MatrixMap may be constructed directly from the matrix "m", which
#  in this particular case describes a 30 degree rotation.
   m = [ [ 0.866, 0.5 ], [-0.5, 0.866] ]
   matrixmap = Ast.MatrixMap( m )

#  Here we take the simpler option of using a ShiftMap rather than a
#  WinMap. The "z" list specifies the shifts.
   z = [ -1.0, 1.0 ]
   shiftmap = Ast.ShiftMap( z )

#  Join the two Mappings together, so that they are applied one after
#  the other.
   newmap = Ast.CmpMap( matrixmap, shiftmap )

#  If necessary, make a copy of the WCS calibration, since we are
#  about to alter it.
   wcsinfo2 = wcsinfo.copy()

#  Re-map the base Frame so that it refers to the new data grid
#  instead of the old one.
   wcsinfo2.remapframe( Ast.BASE, newmap )

...Write a Modified WCS Calibration to a Dataset

#  Create an empty FitsChan to contain the output WCS information. Use a
#  PyFITSAdapter (defined in the Atl module) to provide the sink function
#  that FitsChan uses to write FITS header cards to a PyFITS HDU. Note,
#  doing it this way will mean the FITS file does not inherit the non-WCS
#  keywords from the input FITS file. If this is a problem, then the
#  FitsChan that was used for reading the WCS should also be used for
#  writing the WCS. This will retain non-WCS keywords in the output FITS
#  file.
   fitschan2 = Ast.FitsChan( None, Atl.PyFITSAdapter(hdu_list[ 0 ]) )

#  Ensure the FitsChan will create FITS keywords that use the same WCS
#  encoding that was used by the original FITS file.
   fitschan2.Encoding = encoding

#  Attempt to write the output WCS information to the FitsChan (i.e.
#  convert the FrameSet to a set of FITS header cards stored in the FitsChan).
   if fitschan2.write( wcsinfo2 ) == 0 :

#  If this didn't work (the WCS FrameSet has become too complex to be
#  represented using the original WCS encoding), then use the native AST
#  encoding instead.
      fitschan2.Encoding = "NATIVE"
      fitschan2.write( wcsinfo2 )

#  The cards are currently stored in the FitsChan. They will be written
#  to the PyFITS HDU when the FitsChan is deleted. But in Python, object
#  deletion can happen at unexpected times, so to be in complete control,
#  we ensure the cards are written out now by invoking the writefits method
#  now.
   fitschan2.writefits()

...Display a Graphical Coordinate Grid

#  Create a matplotlib plotting region. Ensure that the matplotlib axis
#  annotations are not drawn.
   dx=12.0
   dy=9.0
   fig = matplotlib.pyplot.figure( figsize=(dx,dy) )
   ax = fig.add_axes( ( 0, 0, 1, 1 ) )
   ax.xaxis.set_visible( False )
   ax.yaxis.set_visible( False )

#  We map the entire FITS pixel grid onto this box. So get the bounds
#  of the pixel grid from the FitsChan. If the NAXIS1/2 keywords are
#  not in the FitsChan, use defaults of 500.
   if "NAXIS1" in fitschan:
      naxis1 = fitschan["NAXIS1"]
   else:
      naxis1 = 500

   if "NAXIS2" in fitschan:
      naxis2 = fitschan["NAXIS2"]
   else:
      naxis2 = 500
   bbox = ( 0.5, 0.5, naxis1 + 0.5, naxis2 + 0.5 )

#  Set the bounds (in matplotlib data coordinates) of the largest rectangle
#  that can be drawn on the matplotlib plotting area that has the same
#  aspect ratio as the FITS array. Shrink the box to leave room for axis
#  annotation.
   fits_aspect_ratio = ( bbox[3] - bbox[1] )/( bbox[2] - bbox[0] )
   grf_aspect_ratio = dy/dx
   if grf_aspect_ratio > fits_aspect_ratio :
      hx = 0.5
      hy = 0.5*fits_aspect_ratio/grf_aspect_ratio
   else:
      hx = 0.5*grf_aspect_ratio/fits_aspect_ratio
      hy = 0.5

   hx *= 0.7
   hy *= 0.7

   gbox = ( 0.5 - hx, 0.5 - hy, 0.5 + hx, 0.5 + hy )

#  Create a drawing object that knows how to draw primitives (lines,
#  marks and strings) into the matplotlib plotting region.
   grf = Grf.grf_matplotlib( ax )

#  Create the AST Plot, using the above object to draw the primitives. The
#  Plot is based on the FrameSet that describes the WCS read from the FITS
#  headers, so the plot knows how to convert from WCS coords to pixel
#  coords, and then to matplotlib data coords.
   plot = Ast.Plot( wcsinfo, gbox, bbox, grf )

#  And finally, draw the annotated WCS axes and any coordinate grid
#  requested in the plotting attributes.
   plot.grid()

#  Make the matplotlib plotting area visible
   matplotlib.pyplot.show()