1

从昨天开始，我一直在寻找答案，但没有运气。所以我有一个一维光谱（.fits）文件，每个波长都有通量值。我已将它们转换为二维数组（x，y）=（波长，通量）并想编写一个程序，该程序将在某些指定的波长（x）处返回通量（y）。我试过这个：

#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

但即使波长不同，程序也总是返回相同的通量值。请帮忙...

4

2 回答 2

1

我对您的代码进行了一些更改：

  • 使用 numpy ot createwave_flux作为ndarrayusingnp.hstack()和np.repeat()np.tile()

  • 使用精美的索引来获取与您的搜索匹配的值

结果代码是：

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

它将返回一个子组，wave_flux其中列中的波值和列中的0通量值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]]
于 2013-07-02T07:30:58.380 回答
1

我会完全跳过二维表的创建，只使用interp：

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

对于您发布的文件，由于 wave_tg 数组的硬编码长度，您发布的代码不起作用。因此，我建议您宁愿使用

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

此外，由于某种原因，您发布的文件似乎实际上并没有达到您正在查找的波长。您可能需要检查您是否正确计算了相应的波长或检查您是否正在查找正确的波长。

于 2013-07-02T10:19:45.037 回答