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 as plt # 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()
# The axes within a WCS FrameSet obtained from an external source may be
# in any arbitrary order. This order may be determined by examining the
# attributes (Domain, Label, IsLonAxis, IsLatAxis, etc) of the current
# Frame in the FrameSet. If your application code for some reason requires
# the axes to be in a specific order, this may be arranged using the
# findframe() method. You need to provide it with a Frame to act as a
# template that defines the nature and order of the axes required by
# your application code. For instance, if your code is restricted to
# 2-dimensional WCS FrameSets that represent (RA,Dec) axes, in that
# order, then you create a suitable SkyFrame with these properties to
# act as a template, and call findframe on the original WCS FrameSet to
# create a modified WCS FrameSet with the same properties as the
# template (if possible). For the above example, a default SkyFrame is
# an appropriate template since the default SkyFrame represents ICRS
# (RA,Dec) axes. So if "wcsinfo" is the original WCS FrameSet we can
# do...
newwcs = wcsinfo.findframe( Ast.SkyFrame() )
# It may be that the "wcsinfo" FrameSet did not contain any celestial
# axes, or was not 2-dimensional. in which case the above invocation of
# findframe will have failed. Check for this.
if newwcs:
# Demonstrate that the above has succeeded by printing the labels of the
# two axes in the current Frame.
print( "Axis 1: {0}".format( newwcs.Label_1 ) )
print( "Axis 2: {0}".format( newwcs.Label_2 ) )
# 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 figure, 12x9 inches in size.
dx=12.0
dy=9.0
fig = plt.figure( figsize=(dx,dy) )
fig_aspect_ratio = dy/dx
# 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
# Set up the bounding box of the image in pixel coordinates, and get
# the aspect ratio of th eimage.
bbox = (0.5, 0.5, naxis1 + 0.5, naxis2 + 0.5)
fits_aspect_ratio = ( bbox[3] - bbox[1] )/( bbox[2] - bbox[0] )
# Set up the bounding box of the image as fractional offsets within the
# figure. The hx and hy variables hold the horizontal and vertical half
# widths of the image, as fractions of the width and height of the figure.
# Shrink the image area by a factor of 0.7 to leave room for annotated axes.
if fig_aspect_ratio > fits_aspect_ratio :
hx = 0.5
hy = 0.5*fits_aspect_ratio/fig_aspect_ratio
else:
hx = 0.5*fig_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 )
# Add an Axes structure to the figure and display the image within it,
# scaled between data values zero and 100. Suppress the axes as we will
# be using AST to create axes.
ax_image = fig.add_axes( [ gbox[0], gbox[1], gbox[2] - gbox[0],
gbox[3] - gbox[1] ], zorder=1 )
ax_image.xaxis.set_visible( False )
ax_image.yaxis.set_visible( False )
ax_image.imshow( hdu_list[0].data, vmin=0, vmax=200, cmap=plt.cm.gist_heat,
origin='lower', aspect='auto')
# Add another Axes structure to the figure to hold the annotated axes
# produced by AST. It is displayed on top of the previous Axes
# structure. Make it transparent so that the image will show through.
ax_plot = fig.add_axes( [ 0, 0, 1, 1 ], zorder=2 )
ax_plot.xaxis.set_visible(False)
ax_plot.yaxis.set_visible(False)
ax_plot.patch.set_alpha(0.0)
# Create a drawing object that knows how to draw primitives (lines,
# marks and strings) into this second Axes structure.
grf = Grf.grf_matplotlib( ax_plot )
# 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 )
# Set some nice attributes for the plot.
plot.Colour_Border = "lime"
plot.Colour_Ticks = "lime"
plot.MajTickLen = 0.02
plot.Width_TextLab = 2
plot.Width_NumLab = 2
plot.Width_Title = 2
# And finally, draw the annotated WCS axes.
plot.grid()
# Make the matplotlib plotting area visible
plt.show()
from astropy.io import fits
import starlink.Ast as Ast
import numpy as np
# This example creates a FITS binary table containing a description of
# a polygonal area of sky using the IVOA Multi-Order Coverage (MOC)
# recommendation. AST provides several other ways to define the area
# of sky to be described by the MOC. For instance, other forms of
# 2-dimensional shape may be used (circles, ellipses, boxes, etc),
# or a contour within a pixel array may be used (see astAddPixelMask),
# or individual HEALPix cells may be specified (see astAddCell).
# Define the polygon vertices. This example uses galactic coordinates,
# although they may be defined in any celestial coordinate system
# supported by AST. They will be converted to ICRS automatically before
# creating the MOC. The values in this list are in degrees. Note,
# reversing the order of the vertices will cause the Region to describe
# the area *outside* the polygon.
glon_list = [260.8126, 252.9279, 266.8518, 288.8583, 299.2654]
glat_list = [-58.5782, -46.3082, -40.9669, -46.8499, -60.9320]
vertices = [glon_list, glat_list]
# Create an AST polygon describing this region in galactic coordinates.
# Convert the latitude and longitude values to radians, as required by
# AST.
polygon = Ast.Polygon( Ast.SkyFrame("System=galactic"),
np.radians(vertices) )
# Create an empty Moc, indicating that it should be accurate to 60
# arc-seconds (this corresponds to a MOC order of 12).
moc = Ast.Moc( "MaxRes=60")
# Add the polygon into the Moc.
moc.addregion( polygon )
# Get a list of integers that describe the MOC in a form suitable for
# storing in a FITS binary table.
coldata = moc.getmocdata()
# Get a FitsChan holding the FITS header values required for the binary
# table.
fc = moc.getmocheader()
# Use astropy.io.fits to create the binary table in file "astmoc.fits",
# copying in all the headers from "fc". (Note, the docs say that the
# "name" value is optional when creating a Column, but omitting it causes
# an exception to be raised in BinTableHDU.from_columns).
col = fits.Column( name=fc['TTYPE1'], format=fc['TFORM1'], array=coldata )
header = fits.Header()
for cardimage in fc:
header.append( fits.Card.fromstring( cardimage ) )
table = fits.BinTableHDU.from_columns( [col], header )
table.writeto( 'astmoc.fits', overwrite=True )
from astropy.io import fits import starlink.Ast as Ast import starlink.Atl as Atl import numpy as np # This example reads an IVOA Multi-Order Coverage (MOC) map from a FITS # binary table held in file "moc.fits" and tests if it overlaps the # region of sky covered by an image held in FITS file "image.fits". # Use astropy.io.fits to read the integer values defining a MOC from # a fits file called "moc.fits". For simplicity, it is assumed in this # example that the MOC is in the first extension. A real application # would use a more sophisticated method to find the MOC. hdul = fits.open( "moc.fits" ) data = hdul[1].data.field(0) # Create an empty AST Moc object, leaving the resolution unspecified so # that it is set to the resolution of the MOC importred from the FITS # file. moc = Ast.Moc() # Add the data array read from "moc.fits" into the Moc. moc.addmocdata( data ) # Now read the WCS from a FITS image held in file "image.fits". This is # described in more detail in an earlier pyast example. hdu_list = fits.open( 'image.fits' ) fc = Ast.FitsChan( Atl.PyFITSAdapter(hdu_list[ 0 ]) ) wcsinfo = fc.read() # Create an AST Region describing the sky area covered by the pixel # array. First create a Box in pixel coords (the base Frame in the # "wcsinfo" FrameSet ). Then map it into the sky coordinate system # defined by the current Frame of the FrameSet. A FrameSet can act as # both a Mapping and a Frame - the first argument to "pixbox.mapregion" # uses "wcsinfo" as a Mapping and the second one uses it as a Frame. lbnd = [ 0.5, 0.5 ] ubnd = [ fc['NAXIS1'], fc['NAXIS2'] ] pixbox = Ast.Box( wcsinfo.getframe( Ast.BASE ), 1, lbnd, ubnd ) wcsbox = pixbox.mapregion( wcsinfo, wcsinfo ) # See if the wcsbox overlaps the moc, and report. overlap = moc.overlap( wcsbox ) if overlap == 2: print( "The MOC is completely inside the image") elif overlap == 3: print( "The image is completely inside the MOC") elif overlap == 4: print( "The MOC and the image overlap partially") elif overlap == 5: print( "The MOC and image cover nearly identical regions") else: print( "The MOC and image do not overlap")
import starlink.Ast as Ast
# Read the contents of the text file as a single string.
with open('moc.txt', 'r') as file:
data = file.read()
# Create an empty Moc.
moc = Ast.Moc()
# Add the data read from the text file into the empty Moc.
is_json = moc.addmocstring( data )
if is_json:
print("Read a JSON-encoded MOC from file 'moc.txt'")
else:
print("Read a string-encoded MOC from file 'moc.txt'")
# After establishing a Plot as shown in the earlier example, "Display a
# Graphical Coordinate Grid Over a Greyscale Image", the outline of the
# MOC can be drawn as follows.
plot.regionoutline( moc )
import starlink.Ast as Ast # Create an empty MOC, specifying its maximum order as 6 (about 1 degree) moc = Ast.Moc() moc.MaxOrder = 6 # Add in the four order 5 cells that touch the ICRS origin (RA=0,Dec=0) moc.addcell( 5, 4522 ) moc.addcell( 5, 4864 ) moc.addcell( 5, 4693 ) moc.addcell( 5, 4351 ) # Remove the four order 6 cells that touch the ICRS origin, thus leaving # a hole in the coverage at (RA=0,Dec=0). moc.addcell( 6, 17407, Ast.XOR ) moc.addcell( 6, 18090, Ast.XOR ) moc.addcell( 6, 18773, Ast.XOR ) moc.addcell( 6, 19456, Ast.XOR ) # Convert the Moc to a string using the "string-encoding" defined in the # IVOA's MOC standard. print( moc.getmocstring() )