2017-10-25 6 views
2

2次元のサーフェスをデータに合わせようとしています。具体的には、IRAFのFITCOORDSと同じように、ピクセル座標を波長座標にマッピングする関数を探したいと思います。 Pythonの3.6にastropy 2.0.2numpy 1.13.3データポイントのいくつかがNaNである場合の2次元関数のフィット方法は?

import numpy as np 
from astropy.modeling.models import Chebyshev2D 
from astropy.modeling.fitting import LevMarLSQFitter 
#%% 
test = np.array([[7473, 7040, 6613, 6183, 5753, 5321, 4888], 
    [7474, 7042, 6616, 6186, np.nan, 5325, 4893], 
    [7476, 7044, 6619, 6189, 5759, 5328, 4897], 
    [7479, 7047, np.nan, 6192, 5762, 5331, 4900]]) 
grid_pix, grid_wave = np.mgrid[:4, :7] 
fitter = LevMarLSQFitter() 
c2_init = Chebyshev2D(x_degree=3, y_degree=3) 
c2_fit = fitter(c2_init, grid_wave, grid_pix, test) 
print(c2_fit) 

ResultI:

は一例として、私は次のスニペットでtest配列へのフィット感を見つけたい

Model: Chebyshev2D 
Inputs: ('x', 'y') 
Outputs: ('z',) 
Model set size: 1 
X-Degree: 3 
Y-Degree: 3 
Parameters: 
    c0_0 c1_0 c2_0 c3_0 c0_1 c1_1 c2_1 c3_1 c0_2 c1_2 c2_2 c3_2 c0_3 c1_3 c2_3 c3_3 
    ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- 
    0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 
WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting] 

それは明らかだフィッティングこと働いたことがない

np.nanを特定の値に変更すると、当てはまりがうまくいきます(たとえば、手動でnp.nanを0,1などに変更します)。

妥当な結果を得るにはどうすればよいですか?フィッターにnp.nanの値を無視させることはできますか?

+0

ダミー値で置き換えるのではなく、nanを削除するだけですか? – Julien

+0

@Julien配列から 'np.nan'を削除するにはどうしたらいいですか? –

+0

'np.where'たとえば – Julien

答えて

2

データ内のnanと、グリッドの対応する「インデックス」を単に削除することができます。 boolean indexingと例えば:

notnans = np.isfinite(test) # array containing True for finite values and False for nans/infs 
c2_fit = fitter(c2_init, grid_wave[notnans], grid_pix[notnans], test[notnans]) 
print(c2_fit) 

それはまだパラメータに線形に関する警告を出力しますが、値は間違いなく非ゼロである:

Model: Chebyshev2D 
Inputs: ('x', 'y') 
Outputs: ('z',) 
Model set size: 1 
X-Degree: 3 
Y-Degree: 3 
Parameters: 
     c0_0   c1_0   c2_0  ...  c2_3   c3_3  
    ------------- -------------- -------------- ... --------------- ---------------- 
    7473.01546325 -431.633443323 0.471726190475 ... 0.0229037267082 -0.0012077294686 

ここでのトリックはxyとあなたのdataは必要がないということです2D配列である場合は、2Dグリッドを「表す」限り、(ブールインデックスによって返されるような)1D配列にすることができます。

NaNを含む「大きな領域」がある場合には、このアプローチは、フィッターがそこに適合することができるため、十分ではない可能性があります。その場合は、あなたはastropy.convolution.convolveを使用して、これらの領域の上に補間することができ、その後、convolveの結果であなたのdataのNaNを交換する:いくつかのまばらに分布NaNのためしかし

from astropy.convolution import convolve 
# Just for illustration I used a 5x5 mean filter here, the size must be adjusted 
# depending on the size of all-nan-regions in your data. 
mean_convolved = convolve(test, np.ones((5, 5)), boundary='extend') 
# Replacing NaNs in test with the mean_convolved values: 
nan_mask = np.isnan(test) 
test[nan_mask] = mean_convolved[nan_mask] 

# Now pass test to the fitter: 
c2_fit = fitter(c2_init, grid_wave, grid_pix, test) 

コンボリューションは必要ありません。おそらく、2つのアプローチを比較し、どちらがより「信頼できる」結果をもたらすかを確認する必要があります。欠損値はフィッティングの実際の問題になります。

+0

ああ!これはとても単純で正確に私が探していたものです!ありがとうございます:D –

関連する問題