当前位置：首页 > article >正文

获得两类相关点之间的线性关系

article 2024/9/24 15:17:07

例如我是有两幅图，但是他们的单位不一样（一个温度，一个央斯基）原则上他们只是相差一个因子。可以利用一些线性拟合的关系来获得那个因子。其实代码和那个计算谱指数差不多，只是做了一些容易犯错误的修改（斜率的计算上）

from astropy.io import fits as pf
from astropy.wcs import WCS
import sys
from copy import deepcopy
from astropy.coordinates import SkyCoord
from reproject import reproject_interp
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u
import math
from statistics import mode
from astropy.convolution import convolve, Gaussian2DKernel
from astropy.modeling import Fittable1DModel, Parameter, models, fitting

#例如data和data_cite 记得先平滑到一个分辨率上
cal_spectral(data[mark1==1].flatten()[::1], data_cite[mark1==1].flatten()[::1], 1.248, 1.248)

def cal_alpha(data1,data2):
    n = data1.size
    x = np.zeros(n) 
    y = np.zeros(n)
    kk = 0
    for i in range(n):
        if np.isnan(data1[i]) or np.isnan(data2[i]): continue
        x[kk] = data1[i]
        y[kk] = data2[i]
        kk += 1

    if kk<10: return np.nan,np.nan,np.nan,np.nan

    p_init = models.Polynomial1D(1)
    fit_p = fitting.LinearLSQFitter()
    p = fit_p(p_init, x[:kk], y[:kk])
    a, b = p.c0.value, p.c1.value
    return a, b

def cal_spectral(data1, data2, freq1, freq2):

    a1, b1 = cal_alpha(data1, data2)
    a2, b2 = cal_alpha(data2, data1)
    print(b1, 1/b2)
    alpha_1 = (b1 + 1/b2) / 2.
    d_alpha_1 = abs(1/b2 - b1)
    a1_p = a1 * 1.
    b1_p = b1 * 1.
    a2_p = -a2/b2
    b2_p = 1./b2
    #画图
    plt.rcParams.update({'font.size': 25})
    fig = plt.figure(figsize=(15,11))
    ax = fig.add_subplot(1, 1, 1)
    ax.plot(data1, data2, 'o')
    ax.plot(data1, a1_p + b1_p * data1)
    ax.plot(data1, a2_p + b2_p * data1)
    ss = r'$-$' + '%4.2f' % abs(alpha_1) + r'$\pm$' + '%4.2f' % d_alpha_1
    ax.text(0.1, 0.9, ss, horizontalalignment='left', verticalalignment='center', transform=ax.transAxes)
    ax.set_xlabel('T (K, at %5.3f GHz)' % freq1)
    ax.set_ylabel('T(K, Extrapolating E. M. Berkhuijsen observ-\n-ational data to a frequency of 1.248 GHz)')
    ax.xaxis.labelpad = 15
    ax.yaxis.labelpad = 8
    plt.savefig('Test_Temperature calibration.png', bbox_inches = 'tight')
    plt.show()