Sam's Astrophysics Coding Tutorial
This Notebook will cover the basics of using Python to analyze astrophysics data, particularly images.
#Our bookkeeping cell
import os #This package allows you to interact with your operating system
import numpy as np #standard math library: https://numpy.org/doc/stable/
import matplotlib.pyplot as plt #standard plotting library: https://matplotlib.org/stable/index.html
from astropy.io import fits #Astropy is a multi-purpose python package made by astronomers: https://docs.astropy.org/en/stable/index.html
from astropy.stats import * #Astropy is massive and we only want specific tools from it for now, so it's best to import only what we need.
from photutils.background import Background2D, MedianBackground #Photutils is another package used to manipulate images and find sources in our images.
from photutils.detection import DAOStarFinder #Photutils documentation: https://photutils.readthedocs.io/en/stable/
from photutils.aperture import CircularAperture
from photutils.aperture import aperture_photometry
import warnings
warnings.filterwarnings('ignore') #I turn off warnings because python complains about log values, but doesn't actually break. You may want to remove this
Astrophysics coding follows the same general steps most of the time
- Import your data and extract relevant information
- Clean your data
- Analyze your data
Lets make some code to handle a basic example, starting by importing and opening our data
Astrophysics data almost always comes in the form of fits files. Fits (flexible image transport system) files come with a data extension and a header extension. The data of a fits file is the image itself, sometimes along with other things like bad pixel masks. The header is a table of information about the image that is relevant to using data from it.
Included with this tutorial was an LCO image. We're going to open it.
#First define the image in python as a string.
imagefile = "Test_Image.fits" #Here I am assuming the image is in the same location on my computer as my python notebook. Use the below code for more flexibility:
# imagefile = r"C:\Users\Sam Whitebook\Test_Image.fits" #Replacing the path with wherever your image actually is. The r is necessary to define a "string literal" so python knows the \'s are for a path.
#Now we will use astropy to open the file. Note that sometimes LCO files come with multiple extensions for data and headers. You may need to index through your data to find your image.
image = fits.open(imagefile)[0].data
header = fits.open(imagefile)[0].header
The data of our image is a matrix of pixel values representing a brightness (photon count) at a pixel. We can show it in python with matplotlib imshow.
plt.figure(figsize=(10,10))
plt.imshow(np.log(image), #We usually take the log of an image to see the variance of the pixels better.
origin='lower', #This tells matplotlib where the [0,0] index is. For our images it is in the bottom.
cmap = "gray", #Use whatever colormap you like, but I find gray the easiest to visualize.
vmin = 6, vmax = 11) #This defines the data range we want to see. There are ways to automate this but I always do it manually.
plt.axis('off') #This removes the axes, which aren't really relevant for just viewing an image.
plt.show()
Our header is imported as an "astropy header" object, which is really just a fancy dictionary.
#The header has a bunch of 'cards' which are like dictionary keys. Calling header.cards lists all the possible cards. You may uncomment and call the function below, but it is a large output.
# print(header.cards)
#Lets call the Modified Julian Date for when this image was taken.
mjd = header['MJD-OBS']
print("The Modified Julian Date is", mjd)
The Modified Julian Date is 56370.08241256944
Ok, so now we know how to import our data, but as of now that's just a pretty picture. It still has background noise, and we don't know how many stars are in it, or how bright they are. Next we will do some basic photometry using the rest of the tools we have imported. Unfortunately many of these tools are what we call "black boxes" in that they take an input and give you an output, but attempting to understand how they do so is a difficult exercise that we often do not concern ourselves with.
#I will define a function to do our photometry for us. It appears long, but much of this code can be copy pasted for your own work.
#Most of these black box functions come from astropy stats and photutils. I will briefly explain what they do and what their inputs are, but not how they work.
def photometry(data, header):
sigma_clip = SigmaClip(sigma=3.0) #A sigma clip excludes values that are a certain number of standatd deviations, sigma, above the mean. We don't include these in our background because they're stars.
bkg_estimator = MedianBackground() #This specifies which method of estimating the background photutils uses. We almost always want to stick with median.
time = header["MJD-OBS"] #I am simply returning the modified julian date off the header for convenience.
mean, median, std = sigma_clipped_stats(data, sigma=3.0) #Here astropy just does some statistics for us. Use the same sigma as you choose for the clip. 3 is usually good.
bkg = Background2D(data, #Which image to make a background for
box_size=(50, 50), #Box size tells the function how much of the image to parse at once. Leave it at 50x50 usually.
filter_size=(3, 3), #3x3 is default for the filter size. I honestly don't know what a 'median filter' is, but this has always worked.
sigma_clip=sigma_clip, #Just use the sigma clip we defined earlier.
bkg_estimator=bkg_estimator) #Same thing for the background estimator method.
bkg_median = bkg.background_median #This pulls the median of the background values.
daofind = DAOStarFinder(fwhm=10, #This is the full width half maximum of the 'gaussian kernel' of the starfinder in units of pixels. Smaller values will find more sources, large values will find less.
threshold=5*std) #The threshold determines at what pixel value the starfinder will call something a source. We usually use a multiple of the standard deviation we found earlier. Again smaller for more.
sources = daofind(data - bkg_median) #We subtract the median we found earlier from the whole image.
positions = np.transpose((sources['xcentroid'], sources['ycentroid'])) #Positions are transposed by default.
apertures = CircularAperture(positions, r=15) #This defines an 'aperture' by radius around each source. The pixel values inside the aperture are summed to yield the photon count of the star.
phot_table = aperture_photometry(data - bkg_median, apertures) #We use the aperture photometry function to list the magnitudes and positions of the sources we found.
magnitudes = phot_table['aperture_sum']
positions = np.transpose([phot_table['xcenter'], phot_table['ycenter']])
return time, positions, magnitudes
That was a lot. Let
time, positions, magnitudes = photometry(image, header)
print("We found", len(magnitudes), "stars")
time, positions, magnitudes = photometry(image, header)
print("We found", len(magnitudes), "stars")
#Lets call a random star and see what information we found on it.
star_index = np.random.randint(0,283)
star_position = positions[star_index]
star_magnitude = magnitudes[star_index]
print("The star's pixel position is", star_position)
print("The star's photon count is", star_magnitude)
We found 283 stars The star's pixel position is [1462.00177082 104.93995699] The star's photon count is 13571.090080109187
What you would do with this information is dependent on your project, but for an example of a time series for exoplanet transits, you may wish to refer to a piece of code I have on my github for treating a transiting exoplanet (Actually the source of this image if you were curious.)
Link:https://github.com/Sewhitebook/Sewhitebook/blob/main/Transiting_Exoplanets.ipynb