starlink.Atl
module, which provide wrappers
for some of the code fragments shown below.
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.
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()
# 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 the internal representation of the AST FrameSet print( wcsinfo ) # or alternatively... wcsinfo.show()
# 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 )
# 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()
# 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 )
# 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 )
# 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 = []
# 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 )
# 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()
# 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()