Greg Ziemann sent me a short code that processes the spectral cube of IZw136 with a numpy-driven python code. Here it is:
This code: /home/sco/jimiL/red_Aug17_2020/gregz_aug17_2020_IZw136/gregz1.py
from astropy.io import fits
from astropy.table import Table
import numpy as np
# A simple python code would be (writing the spectrum to "spectrum.dat"
# as (wavelength and flux density in ergs/s/cm^2/s):
extract_radius = 2.25 #(in arcseconds)
sky_radius = 4.5 #(in arcseconds)
k = fits.open('eng_galaxy_1_LRS2B_cube.fits')
x = np.arange(k[0].header['NAXIS1'])*k[0].header['CDELT1'] + k[0].header['CRVAL1']
y = np.arange(k[0].header['NAXIS2'])*k[0].header['CDELT2'] + k[0].header['CRVAL2']
w = np.arange(k[0].header['NAXIS3'])*k[0].header['CDELT3'] + k[0].header['CRVAL3']
xgrid, ygrid = np.meshgrid(x, y)
data = k[0].data
d = np.sqrt(xgrid**2 + ygrid**2)
sel = d <= extract_radius
backsel = d>= sky_radius
background = np.nanmean(data[:, backsel], axis=1)
spectrum = np.nansum(data[:, sel] - background[:, np.newaxis], axis=1)
wave = w
Table([wave, spectrum], names=['wavelength', 'flux_density']).write('spectrum.dat', format='ascii.fixed_width_two_line')
There are a few things that I really do not understand, and GZ has suggested a
typical image processing problem, and he will explain how he would do it.
Hi Greg:
Thanks for the very kind offer of help. I have developed some test
codes where I read in test images as straight text or as FITS images
(using astropy's fitsio). I used this to play with some basic slicing and
indexing. I suppose a typical image processing problem I would want to
do in python (very similar to your example) is the following:
1) Read a 2D image into a numpy array (for example a 100x100 test
image with a "star" at x,y=66,70).
2) I test pixels around x,y=66,70 to find all those with a radius
of less that RAD = 5 pixels. I add up the flux contained in those
pixels. I.E. I intergate signal in fixed ciruclar apetrue
*** QUETSION: maybe in your "sel = d <= extract_radius" line
of you code, you are creating a pixel mask to flag
the aperture pixels for doing thsi?
3) I perform the same operation, but for a annulus around the circular
aperture. The annulus has an inner radius of 7 pixels and an outer
radius of 10 pixels. I integrate all of the flux in the annulus and
divide by the number of pixels to get a mean background flux per pixel.
4) I correct the circle aperture flux to subtract the background and
derive a total signal (in adu) for the circular aperture.
I think this is fairly clear. I would use an old guy's language to loop
in X,Y and gather arrays of background and aperture pixels. Actually, I'd
make a mask, like you, so that I could repeatedly survey an elliptical
annulus or aperture, yadda, yadda. Instead of a simple mask, I might make
masks with radius and angle values, so that I could compute azimutal
profiles in elliptical rings, etc.... But, I need to be able to address
pixels. What I am not understanding clearly is how you use one line to
mask pixels and integrate signals at the same time:
spectrum = np.nansum(data[:, sel] - background[:, np.newaxis], axis=1)
I know this is a very powerful way of shortening the command stream, but I
do not yet understand clearly the syntax.
Thanks for the offer of help. Sorry you have to hold my hand through
understanding what is probably a very basic thing.
Cheers,
Steve