import matplotlib.pyplot as plt
from astropy.utils.data
import get_pkg_data_filename
from astropy.io
import fits
from astropy.wcs
import WCS
from astropy.visualization.wcsaxes.frame
import EllipticalFrame
image_file = get_pkg_data_filename('tutorials/FITS-images/HorseHead.fits')
image_header = fits.open(image_file)[0].header # extract header info
image_data = fits.getdata(image_file, ext = 0)
# use the WCS class to get coordinate info and projection axes to use
wcs = WCS(image_header)
# make the axes using the EllipticalFrame class
ax = plt.subplot(projection = wcs, frame_class = EllipticalFrame)
# add the image
im = ax.imshow(image_data, origin = 'lower')
# add a grid
overlay = ax.get_coords_overlay('fk5')
overlay.grid(color = 'white', ls = 'dotted')This example uses astropy.utils.data to download the file, astropy.io.fits to open the file, and matplotlib.pyplot to display the image.,This example opens an image stored in a FITS file and displays it to the screen.,Download the example FITS files used by this example:,Use astropy.io.fits.info() to display the structure of the file:
import matplotlib.pyplot as plt
from astropy.visualization
import astropy_mpl_style
plt.style.use(astropy_mpl_style)from astropy.utils.data
import get_pkg_data_filename
from astropy.io
import fits
image_file = get_pkg_data_filename('tutorials/FITS-images/HorseHead.fits')fits.info(image_file)
Filename: /home/docs / .astropy / cache / download / url / ff6e0b93871033c68022ca026a956d87 / contents
No.Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 161(891, 893) int16
1 er.mask 1 TableHDU 25 1600 R x 4 C[F6 .2, F6 .2, F6 .2, F6 .2]image_data = fits.getdata(image_file, ext = 0)
print(image_data.shape)
The image is plotted in terms of the pixel position on the \(x\) and \(y\) axes. Astropy contains a range of functions for plotting images in actual sky coordinates, overlaying coordinate grids, contours etc. (e.g. see the documentation for the astropy.wcs and astropy.visualization sub-packages).,For this particular file, the primary header is just a standard placeholder, which needs to be present but doesn’t convey any useful information. Other types of FITS data file may contain more extensive primary headers, e.g. containing important information about the observation (telescope, date, sky location) used to take an image or spectrum.,Returning to our list of HDUs, we see that HDU[1] has a more extensive header (with 67 cards) and a set of dimensions (927552 rows and 25 columns) and data formats corresponding to the table dimensions and the formats of the data in the columns therein.,Image data in FITS files takes the form of a 2-dimensional array where each item corresponds to a pixel value. For example, let’s look at a FITS image of the famous Horsehead nebula. You can find it in the Lesson data directory here). We’ll first open the file and look at its structure:
from astropy.io
import fits
gals = fits.open('gal_info_dr7_v5_2.fit')
gals.info()Filename: gal_info_dr7_v5_2.fit No.Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 4() 1 1 BinTableHDU 67 927552 R x 25 C[I, J, I, 5 I, E, E, 5 E, I, I, 19 A, 6 A, 21 A, E, E, I, E, E, E, E, E, E, 3 E, 3 E, 5 E, 12 A]
gals[0].headerSIMPLE = T / Dummy Created by MWRFITS v1 .6 a
BITPIX = 8 / Dummy primary header created by MWRFITS
NAXIS = 0 / No data is associated with this header
EXTEND = T / Extensions may(will!) be presentprint(gals[0].header['BITPIX'])
print(gals[0].header.comments['BITPIX'])8
Dummy primary header created by MWRFITSContextual clues are important when dealing with FITS files. Files that are explicitly single images almost always have the image stored in the 0th extension, while files that are explicitly table data tend to have the table stored in the 1st extension (with the 0th empty).,Next, we need to plot the image data. There are several operations that we almost always perform when plotting astronomical data, as well as several user-preferences for how we “by default” plot images before we begin tweaking things. If you spend any time as an astronomer, you will plot quite literally thousands of images — why set all these settings every time, when we can write a handy function to do it for us?,Often, for calibration purposes, we’d need to create apertures around all those sources in the image. We definitely don’t want to do that by hand! Instead, we’re going to use the sep package.,The FITS standard is pretty old, and may be retired soon, but almost all ground based telescopes still use it for their native image output, so it behooves us to know how to get image data out of FITS files.
import numpy as np def random_return(a, b, c, * args): '' ' A function which requires three inputs which are floats, accepts any number of additional inputs(of any kind), and returns one randomly chosen input. Parameters -- -- -- -- -- a: int description of this integer input b: int description of this integer input c: int description of this integer input * args: tuple any additional arguments get stored here Returns -- -- -- - choice The randomly selected input(type not specified) '' ' full_input_list = [a, b, c] + list(args) choice = np.random.choice(full_input_list) return choice
random_return(1, 5, 4, 6, 4, 21, 6)
6
help(random_return)
Help on
function random_return in module __main__:
random_return(a, b, c, * args)
A
function which requires three inputs which are floats, accepts any number of additional inputs(of any kind), and returns one randomly chosen input.
Parameters
-- -- -- -- --
a: int
description of this integer input
b: int
description of this integer input
c: int
description of this integer input *
args: tuple
any additional arguments get stored here
Returns
-- -- -- -
choice
The randomly selected input(type not specified)a: int, optional
Description of the thing.(
default 5)mosaics of multiple images (e.g. HI channel maps),plot that covers the entire sky (allsky plot),We have two examples of a mosaic of plots. First a mosaic is presented with an image and two position-velocity diagrams. The second is a classic examples which shows channel maps from an HI data cube at different velocities.,The base image has scheduled a function to interactively change its colors while the second image remains fixed. This enables you do compare the two images.
from kapteyn
import maputils
from matplotlib
import pyplot as plt
f = maputils.FITSimage("m101.fits")
fig = plt.figure()
frame = fig.add_subplot(1, 1, 1)
annim = f.Annotatedimage(frame)
annim.Image()
annim.Graticule()
annim.plot()
annim.interact_imagecolors()
plt.show()from kapteyn
import maputils
from matplotlib
import pyplot as plt
fitsobject = maputils.FITSimage('ngc6946.fits')
print("HEADER:\n")
print(fitsobject.str_header())
print("\nAXES INFO:\n")
print(fitsobject.str_axisinfo())
print("\nEXTENDED AXES INFO:\n")
print(fitsobject.str_axisinfo(long = True))
print("\nAXES INFO for image axes only:\n")
print(fitsobject.str_axisinfo(axnum = fitsobject.axperm))
print("\nAXES INFO for non existing axis:\n")
print(fitsobject.str_axisinfo(axnum = 4))
print("SPECTRAL INFO:\n")
fitsobject.set_imageaxes(axnr1 = 1, axnr2 = 3)
print(fitsobject.str_spectrans())HEADER:
SIMPLE = T / SIMPLE FITS FORMAT
BITPIX = -32 / NUMBER OF BITS PER PIXEL
NAXIS = 3 / NUMBER OF AXES
NAXIS1 = 100 / LENGTH OF AXIS
etc.
AXES INFO:
Axis 1: RA-- - NCP from pixel 1 to 100 {
crpix = 51 crval = -51.2821 cdelt = -0.007166(DEGREE)
} {
wcs type = longitude, wcs unit = deg
}
etc.
EXTENDED AXES INFO:
axisnr - Axis number: 1
axlen - Length of axis in pixels(NAXIS): 100
ctype - Type of axis(CTYPE): RA-- - NCP
axnamelong - Long axis name: RA-- - NCP
axname - Short axis name: RA
etc.
WCS INFO:
Current sky system: Equatorial
reference system: ICRS
Output sky system: Equatorial
Output reference system: ICRS
etc.
SPECTRAL INFO:
0 FREQ - V2F(Hz)
1 ENER - V2F(J)
2 WAVN - V2F(1 / m)
3 VOPT - V2W(m / s)
etc.from kapteyn
import maputils
print("Projection object from FITSimage:")
fitsobj = maputils.FITSimage("mclean.fits")
print("crvals:", fitsobj.convproj.crval)
fitsobj.set_imageaxes(1, 3)
print("crvals after axes specification:", fitsobj.convproj.crval)
fitsobj.set_spectrans("VOPT-???")
print("crvals after setting spectral translation:", fitsobj.convproj.crval)
print("Projection object from Annotatedimage:")
annim = fitsobj.Annotatedimage()
print("crvals:", annim.projection.crval)Projection object from FITSimage:
crvals: (178.7792, 53.655000000000001)
crvals after axes specification: (178.7792, 1415418199.4170001, 53.655000000000001)
crvals after setting spectral translation: (178.7792, 1050000.0000000042, 53.655000000000001)
Projection object from Annotatedimage:
crvals: (178.7792, 1050000.0000000042, 53.655000000000001)>>> fitsobject = maputils.FITSimage('alt2.fits', hdunr = 0, alter = 'A', memmap = 1) >>>
fitsobject = maputils.FITSimage(promptfie = maputils.prompt_fitsfile) >>>
fitsobject = maputils.FITSimage(promptfie = maputils.prompt_fitsfile, alter = 'A') >>>
fitsobject = maputils.FITSimage('NICMOSn4hk12010_mos.fits', hdunr = 'err')Sky Map Plotting (ligo.skymap.plot.allsky),Sky Map Plotting (ligo.skymap.plot.allsky) Examples Insets Complete Example , ligo.skymap v1.0.3.dev3+gba63039 »
import ligo.skymap.plot
from matplotlib
import pyplot as plt
ax = plt.axes(projection = 'astro hours mollweide')
ax.grid()import ligo.skymap.plot
from matplotlib
import pyplot as plt
ax = plt.axes(projection = 'geo aitoff')
ax.grid()import ligo.skymap.plot
from matplotlib
import pyplot as plt
ax = plt.axes(projection = 'astro zoom',
center = '5h -32d', radius = '5 deg', rotate = '20 deg')
ax.grid()import ligo.skymap.plot
from matplotlib
import pyplot as plt
ax = plt.axes(projection = 'geo globe', center = '-50d +23d')
ax.grid()import ligo.skymap.plot
from matplotlib
import pyplot as plt
fig = plt.figure(figsize = (9, 4), dpi = 100)
ax_globe = plt.axes(
[0.1, 0.1, 0.8, 0.8],
projection = 'astro degrees globe',
center = '120d +23d')
ax_zoom_rect = plt.axes(
[0.0, 0.2, 0.4, 0.4],
projection = 'astro degrees zoom',
center = '150d +30d',
radius = '9 deg')
ax_zoom_circle = plt.axes(
[0.55, 0.1, 0.6, 0.6],
projection = 'astro degrees zoom',
center = '120d +10d',
radius = '5 deg')
ax_globe.mark_inset_axes(ax_zoom_rect)
ax_globe.connect_inset_axes(ax_zoom_rect, 'upper left')
ax_globe.connect_inset_axes(ax_zoom_rect, 'lower right')
ax_globe.mark_inset_circle(ax_zoom_circle, '120d +10d', '4 deg')
ax_globe.connect_inset_circle(ax_zoom_circle, '120d +10d', '4 deg')
ax_globe.grid()
ax_zoom_rect.grid()
ax_zoom_circle.grid()
for ax in [ax_globe, ax_zoom_rect, ax_zoom_circle]:
ax.set_facecolor('none')
for key in ['ra', 'dec']:
ax.coords[key].set_auto_axislabel(False)from astropy.coordinates
import SkyCoord
from astropy.io
import fits
from astropy
import units as u
import ligo.skymap.plot
from matplotlib
import pyplot as plt
url = 'https://dcc.ligo.org/public/0146/G1701985/001/bayestar_no_virgo.fits.gz'
center = SkyCoord.from_name('NGC 4993')
fig = plt.figure(figsize = (4, 4), dpi = 100)
ax = plt.axes(
[0.05, 0.05, 0.9, 0.9],
projection = 'astro globe',
center = center)
ax_inset = plt.axes(
[0.59, 0.3, 0.4, 0.4],
projection = 'astro zoom',
center = center,
radius = 10 * u.deg)
for key in ['ra', 'dec']:
ax_inset.coords[key].set_ticklabel_visible(False)
ax_inset.coords[key].set_ticks_visible(False)
ax.grid()
ax.mark_inset_axes(ax_inset)
ax.connect_inset_axes(ax_inset, 'upper left')
ax.connect_inset_axes(ax_inset, 'lower left')
ax_inset.scalebar((0.1, 0.1), 5 * u.deg).label()
ax_inset.compass(0.9, 0.1, 0.2)
ax.imshow_hpx(url, cmap = 'cylon')
ax_inset.imshow_hpx(url, cmap = 'cylon')
ax_inset.plot(
center.ra.deg, center.dec.deg,
transform = ax_inset.get_transform('world'),
marker = ligo.skymap.plot.reticle(),
markersize = 30,
markeredgewidth = 3)The next display or contour will plot the image in square #2, leaving the image already plotted in #1 on the device. , Images can be placed side by side for comparison purposes using disp/left and disp/right/overlay. The currently active window is the last one displayed. , When the last viewport in the the configuration is reached, a new plotting device page will be opened and the numbering will begin again at 1. , The number following the dash (-) represents the viewport number in the file. If you display or contour at this point the image will be displayed in square #1 represented in the ASCII grid above. Using show again:
read / fits im1.fits!read the 1 st image
display / left!display the 1 st image
grid!overlay the sky grid
read / fits im2.fits!read the 2n d image
disp / right / overlay!display the two side by side
grid!overlay the sky grid