Python 2D array -- How to plug in x and retrieve y value?

Question

I have been looking for an answer since yesterday but no luck. So I have a 1D spectrum (.fits) file with flux value at each wavelength. I have converted them into a 2D array (x,y)=(wavelength, flux) and want to write a program which will return flux(y) at some assigned wavelengths(x). I have tried this:

#modules
import scipy
import numpy as np
import pyfits as pf

#Target Global Vaiables
hdulist_tg = pf.open('cutmask1-2.0001.fits')    
hdr_tg = hdulist_tg[0].header
flux_tg = hdulist_tg[0].data
crval_tg = hdr_tg['CRVAL1']             #Starting wavelength 
cdel_tg = hdr_tg['CDELT1']              #Wavelength axis width 
wave_tg = crval_tg + np.arange(3183)*cdel_tg        #Create an x-axis
wavelist = [6207,6315,6369,6438,6490,6565,6588]

wave_flux=[]
diff = 10
for wave in wave_tg:
    for flux in flux_tg:
        wave_flux.append((wave,flux))           

for item in wave_flux:
    wave = item[0]
    flux = item[1]

#Where I got my actual wavelength that exists in wave_tg
    diffmatch = np.abs(wave - wavelist[0])
    if diffmatch < diff:
        flux_wave = flux  
        diff = diffmatch 
        wavematch = wave

print wavelist[0],flux_wave,wavematch

but the program always return the same flux value even though the wavelength is different. Please help...

Answer 1

I would skip the creation of the two dimensional table altogether and just use interp:

fluxvalues = np.interp(wavelist, wave_tg, flux_tg)

For the file you posted, the code you posted doesn't work due to the hard-coded length of the wave_tg array. I would therefore recommend you rather use

wave_tg = crval_tg + np.arange(len(flux_tg))*cdel_tg

Also, for some reason it seems that the file you posted doesn't actually go up to the wavelengths you are looking up. You might need to check that you are calculating the corresponding wavelengths correctly or check that you are looking up the right wavelengths.

Answer 2

I've made some changes in your code:

• using numpy ot create wave_flux as a ndarray using np.hstack(), np.repeat() and np.tile()

• using fancy indexing to get the values matching your search

The resulting code is:

#modules
import scipy
import numpy as np
import pyfits as pf

#Target Global Vaiables
hdulist_tg = pf.open('cutmask1-2.0001.fits')
hdr_tg = hdulist_tg[0].header
flux_tg = hdulist_tg[0].data
crval_tg = hdr_tg['CRVAL1']             #Starting wavelength
cdel_tg = hdr_tg['CDELT1']              #Wavelength axis width
wave_tg = crval_tg + np.arange(3183)*cdel_tg        #Create an x-axis
wavelist = [6207,6315,6369,6438,6490,6565,6588]

wave_flux = np.vstack(( np.repeat(wave_tg, len(flux_tg)),
                        np.tile(flux_tg, len(wave_tg)) )).transpose()

wave_ref = wavelist[0]
diff = 10

print wave_flux[ np.abs(wave_flux[:,0]-wave_ref) < diff ]

Which will return a sub-group of wave_flux with the wave values in column 0 and flux values in column 1:

[[ 6197.10300138   500.21020508]
 [ 6197.10300138   523.24102783]
 [ 6197.10300138   510.6390686 ]
 ...,
 [ 6216.68436446   674.94732666]
 [ 6216.68436446   684.74255371]
 [ 6216.68436446   712.20098877]]