How Do You Extract A Point Spread Function From A Fits Image?

I'm new to python coding and I have a .fits image of a star and I want to extract the PSF from it. How would i go about doing this? After that how would I fit a Gaussian curve to i

Solution 1:

What exactly do you mean with "extracting the PSF from it"? Do you want find its position? Do you want to cut out a box around it? What is it that you in the end want to do with it?

Assuming that you just have an image image.fits with a PSF in it, you can use astropy.modeling to fit a 2D Gaussian to your PSF, after figuring out where its center in the provided image is.

import numpy as np
import matplotlib.pyplot as plt
from astropy.modeling import models, fitting
from astropy.io import fits

# Load the data and find center of PSF
image = fits.getdata('path/image.fits')
cents = np.where(image == np.max(image))
xc = int(cents[1])
yc = int(cents[0])

# Cut out smaller box around PSF
bb = 30
box = image[yc-bb:yc+bb,xc-bb:xc+bb]
yp, xp = box.shape

# Generate grid of same size like box to put the fit on
y, x, = np.mgrid[:yp, :xp]
# Declare what function you want to fit to your data
f_init = models.Gaussian2D()
# Declare what fitting function you want to use
fit_f = fitting.LevMarLSQFitter()

# Fit the model to your data (box)
f = fit_f(f_init, x, y, box)

# Plot the data with the best-fit model
plt.figure(figsize=(8, 2.5))
plt.subplot(1, 3, 1)
plt.imshow(box)
plt.title("Data")
plt.subplot(1, 3, 2)
plt.imshow(f(x, y))
plt.title("Model")
plt.subplot(1, 3, 3)
plt.imshow(box - f(x, y))
plt.title("Residual")
plt.show()

Image: 2D Gaussian fit to PSF