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