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...

could you provide a link to the `cutmask1-2.0001.fits` file? — Saullo G. P. Castro, Jul 02 '13 at 06:58

score 1 · Answer 1 · answered Jul 02 '13 at 07:30

1

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]]

answered Jul 02 '13 at 07:30

Saullo G. P. Castro

56,802
26
179
234

Thank you so much. I have a question -- How could one wavelength value return so many fluxes? – N.K. Jul 02 '13 at 09:11
In the way you are building `wave_flux`, each value in `wave_tg` will have all fluxes in `flux_tg`. The sub-group returns all waves matching your criterion. I am not technically familiar with your topic, so I don't know if, physically, one wavelength can have more than one flux... – Saullo G. P. Castro Jul 02 '13 at 09:16
Or could it be that not all the decimal number is printed in column 0? By the way, your code is real neat. Nice and clean. Thank you again! – N.K. Jul 02 '13 at 09:17
all the decimal number in column 0 are waves and in column 1 are fluxes... the code also runs much faster because the Python `for` loops were removed and it is more memory efficient because you are using an np.ndarray instead of a list for `wave_flux` – Saullo G. P. Castro Jul 02 '13 at 09:19
Another question, what is the different between wave_flux = np.vstack(( np.repeat(wave_tg, len(flux_tg)), np.tile(flux_tg, len(wave_tg)) )).transpose() and wave_flux = np.vstack(( np.repeat(wave_tg, len(flux_tg)), np.repeat(flux_tg, len(wave_tg)) )).transpose() Thank you. – N.K. Jul 02 '13 at 10:18
the `np.repeat()` will repeat each element a given number of times while the `np.tile()` will repeat the whole array a given number of times... – Saullo G. P. Castro Jul 02 '13 at 10:21
Then shouldn't I be using np.repeat for both wave_tg and flux_tg instead? It does yield a different result and each wavelength does not return multiple result of fluxes anymore (Which is really how it's supposed to be ) – N.K. Jul 02 '13 at 10:32
The way posted in the answer to create `wave_flux` produces exactly the same result as your two-level for loop – Saullo G. P. Castro Jul 02 '13 at 10:38

score 1 · Accepted Answer · answered Jul 02 '13 at 10:19

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.

Now that you mentioned, you maybe right. Because the file I posted is the one I only cut out the wavelength portion I need so the axis was shorten. Thank you! — N.K., Jul 02 '13 at 10:42

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

2 Answers2