FITS file (Flexible Image Transport System) - Standard Astro Data Format

FITS (Flexible Image Transport System) is the data format most widely used within astronomy for transporting, analyzing, and archiving scientific data files. FITS is much more than just another image format (such as JPG or GIF) and is primarily designed to store scientific data sets consisting of multidimensional arrays (images) and 2-dimensional tables organized into rows and columns of information.



In [ ]:

    
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

import astropy.io.fits as fits



In [ ]:

    
plt.style.use('ggplot')
plt.rc('axes', grid=False)   # turn off the background grid for images

FITS files consist of at least two parts - `Header` and `Data`

A FITS file is comprised of segments called Header/Data Units (HDUs), where the first HDU is called the Primary HDU.
The primary data array can contain a 1-D spectrum, a 2-D image, or a 3-D data cube.
Every HDU consists of an ASCII formatted Header Unit.
Each header unit contains a sequence of fixed-length 80-character keyword (Cards)



In [ ]:

    
data_file = "./MyData/bsg01.fits"

my_fits_file = fits.open(data_file)

my_fits_file.info()



In [ ]:

    
image_data = my_fits_file[0].data
image_header = my_fits_file[0].header



In [ ]:

    
image_header



In [ ]:

    
print("The image has a shape [height,width] of {0}".format(image_data.shape))
print("The image is made up of data of type {0}".format(image_data.dtype))

FITS format preserves the full dynamic range of data



In [ ]:

    
print("The image has a maximum value of {0}".format(image_data.max()))
print("The image has a minimum value of {0}".format(image_data.min()))



In [ ]:

    
fig, ax = plt.subplots(1,2)

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(image_data,cmap=plt.cm.gray)

ax[1].hist(image_data.flatten(),bins=20);

You can use masks on images



In [ ]:

    
copy_data = np.copy(image_data)        # make a copy of the data to work with

cut_off = 40

mask = np.where(copy_data > cut_off)
copy_data[mask] = 60                   # You can not just throw data away, you have to set it to something.



In [ ]:

    
print("The cropped image has a maximum value of {0}".format(copy_data.max()))
print("The cropped image has a minimum value of {0}".format(copy_data.min()))

You can use specific bins for histograms



In [ ]:

    
my_bins = np.arange(0,110,5)

my_bins



In [ ]:

    
fig, ax = plt.subplots(1,2)

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(copy_data,cmap=plt.cm.gray)

ax[1].hist(image_data.flatten(),bins=my_bins)
ax[1].hist(copy_data.flatten(),bins=my_bins, alpha = 0.40);

You can add and subtract images



In [ ]:

    
another_image_file = "./MyData/bsg02.fits"

another_image_data = fits.getdata(another_image_file)     # a quick way to just get the data



In [ ]:

    
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].set_title("Original Image")
ax[1].set_title("Another Image")

ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(another_image_data, cmap=plt.cm.gray);

The two images above may look the same but they are not!

Subtracting the two images reveals the truth.



In [ ]:

    
real_image = image_data - another_image_data                 # Subtract the images pixel by pixel



In [ ]:

    
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].set_title("Original Image")
ax[1].set_title("Another Image")
ax[2].set_title("Real Image!")

ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(another_image_data, cmap=plt.cm.gray);
ax[2].imshow(real_image, cmap=plt.cm.gray);



In [ ]:

    
print("The real image has a maximum value of {0}".format(real_image.max()))
print("The real image has a minimum value of {0}".format(real_image.min()))

FITS Tables - An astronomical example

Stellar spectra data from the ESO Library of Stellar Spectra



In [ ]:

    
my_spectra_file = './MyData/Star_G5.fits'

my_spectra_fits = fits.open(my_spectra_file)

my_spectra_fits.info()



In [ ]:

    
spectra_data = my_spectra_fits[0].data
spectra_header = my_spectra_fits[0].header



In [ ]:

    
spectra_header[0:15]



In [ ]:

    
# The FITS header has the information to make an array of wavelengths

start_wavelength = spectra_header['CRVAL1']
num_of_points = spectra_header['NAXIS1']
width_of_points  = spectra_header['CDELT1']

end_wavelength = start_wavelength + (num_of_points * width_of_points)

wavelength_range = np.arange(start_wavelength,end_wavelength,width_of_points)



In [ ]:

    
fig, ax = plt.subplots(2,1)
fig.set_size_inches(11,8.5)

fig.tight_layout()

# Full spectra

ax[0].plot(wavelength_range, spectra_data, color='b')
ax[0].set_ylabel("Flux")
ax[0].set_xlabel("Wavelength [angstroms]")

# Just the visible range with the hydrogen Balmer lines

ax[1].set_xlim(4000,7000)
ax[1].set_ylim(0.6,1.2)
ax[1].plot(wavelength_range, spectra_data, color='b')
ax[1].set_ylabel("Flux")
ax[1].set_xlabel("Wavelength [angstroms]")

H_Balmer = [6563,4861,4341,4102,3970,3889,3835,3646]

ax[1].vlines(H_Balmer,0,2, color='r', linewidth=3, alpha = 0.25);

Stellar spectral classes



In [ ]:

    
import glob

star_list = glob.glob('./MyData/Star_*5.fits')

star_list



In [ ]:

    
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,5)

fig.tight_layout()

# Full spectra

ax.set_ylabel("Flux")
ax.set_xlabel("Wavelength [angstroms]")
ax.set_ylim(0.0, 1.05)

for file in star_list:
        
    spectra = fits.getdata(file)
    
    spectra_normalized = spectra / spectra.max()
    
    ax.plot(wavelength_range, spectra_normalized, label=file)

ax.legend(loc=0,shadow=True);

FITS Images - An astronomical example



In [ ]:

    
star_file = "./MyData/star_field.fits"

star_fits_file = fits.open(star_file)

star_fits_file.info()



In [ ]:

    
star_data = star_fits_file[0].data
star_header = star_fits_file[0].header



In [ ]:

    
plt.imshow(star_data, cmap=plt.cm.gray)

Notice the origin is in the upper left corner (the image is upside down)



In [ ]:

    
plt.imshow(star_data, origin='lower', cmap=plt.cm.gray)

Better, the origin is in the lower left corner

World Coordinate System `wcs`



In [ ]:

    
from astropy.wcs import WCS



In [ ]:

    
my_wcs = WCS(star_header)



In [ ]:

    
my_wcs



In [ ]:

    
fig = plt.figure()
ax = fig.add_subplot(111, projection=my_wcs)

fig.set_size_inches(6,6)

fig.tight_layout()

ax.grid(color='y', ls='-')
ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)')

ax.imshow(star_data, origin='lower', cmap=plt.cm.gray);

RGB - Pseudocolor Images



In [ ]:

    
from astropy.visualization import make_lupton_rgb



In [ ]:

    
red_img = fits.getdata("./MyData/m51_IR.fits").astype(float)
green_img = fits.getdata("./MyData/m51_red.fits").astype(float)
blue_img= fits.getdata("./MyData/m51_blue.fits").astype(float)



In [ ]:

    
clean_r =  red_img - np.median(red_img)
clean_g =  green_img - np.median(green_img)
clean_b =  blue_img - np.median(blue_img)



In [ ]:

    
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].set_title("R")
ax[1].set_title("G")
ax[2].set_title("B")

ax[0].imshow(clean_r, origin='lower', cmap=plt.cm.gray)
ax[1].imshow(clean_g, origin='lower', cmap=plt.cm.gray);
ax[2].imshow(clean_b, origin='lower', cmap=plt.cm.gray);



In [ ]:

    
image = make_lupton_rgb(clean_r, clean_g, clean_b, stretch = 5000)



In [ ]:

    
fig, ax = plt.subplots(1,1)
fig.set_size_inches(6,6)

fig.tight_layout()

ax.imshow(image, origin='lower', cmap=plt.cm.gray);

Sigma Clipping



In [ ]:

    
from astropy import stats



In [ ]:

    
my_data = np.array([1, 5, 6, 8, 100, 5, 3, 2, 4, 5])

my_data.mean(), my_data.std()



In [ ]:

    
filtered_data = stats.sigma_clip(my_data, sigma=2, iters=5) 

filtered_data

Rejected Data



In [ ]:

    
my_data[filtered_data.mask]

Accepted Data



In [ ]:

    
my_data[~filtered_data.mask]

Sigma clipping an image



In [ ]:

    
star_data.mean(), star_data.std(), star_data.max(), star_data.min()



In [ ]:

    
clip_star = stats.sigma_clip(star_data, sigma=8, iters=5)



In [ ]:

    
copy_data = np.copy(star_data)        # make a copy of the data to work with

copy_data[~clip_star.mask] = star_data.min()



In [ ]:

    
fig, ax = plt.subplots(1,2)

fig.set_size_inches(12,6)

fig.tight_layout()

ax[0].imshow(star_data, origin='lower', cmap=plt.cm.gray)
ax[1].imshow(copy_data, origin='lower', cmap=plt.cm.gray);