Bo Zhang mailto:bozhang@nao.cas.cn
In this notebook, we describe how to use hrs package to reduce HRS data.
To install hrs package, type pip install hrs in your terminal.
In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
import ccdproc
import hrs
# from astropy.table import Table, Column, Row
os.chdir('/town/HRS/20161110test')
# %matplotlib
In [ ]:
ls
In [ ]:
logt = hrs.LogTable.read_log('./20161110.txt')
In [ ]:
logt.pprint(1000, 1000)
In [ ]:
ind_bias = logt['obj']=='bias'
bias = ccdproc.combine(','.join(logt['filename'][ind_bias]), unit='adu', method='median')
In [ ]:
ind_flat1_1 = (logt['obj']=='flat')*(logt['exp_time']==8)
ind_flat1_2 = (logt['obj']=='flat')*(logt['exp_time']==9)
ind_flat2_1 = (logt['obj']=='flat')*(logt['exp_time']==32)
ind_flat2_2 = (logt['obj']=='flat')*(logt['exp_time']==36)
ind_flat3 = (logt['obj']=='flat')*(logt['exp_time']==144)
flat1_1 = ccdproc.combine(','.join(logt['filename'][ind_flat1_1]), unit='adu', method='median')
flat1_2 = ccdproc.combine(','.join(logt['filename'][ind_flat1_2]), unit='adu', method='median')
flat2_1 = ccdproc.combine(','.join(logt['filename'][ind_flat2_1]), unit='adu', method='median')
flat2_2 = ccdproc.combine(','.join(logt['filename'][ind_flat2_2]), unit='adu', method='median')
flat3 = ccdproc.combine(','.join(logt['filename'][ind_flat3]), unit='adu', method='median')
flat1 = ccdproc.combine([flat1_1.multiply(9.).divide(8.), flat1_2], method='average')
flat2 = ccdproc.combine([flat2_1.multiply(36.).divide(32.), flat2_2], method='average')
flat1_bias = ccdproc.subtract_bias(flat1, bias)
flat2_bias = ccdproc.subtract_bias(flat2, bias)
flat3_bias = ccdproc.subtract_bias(flat3, bias)
flat1_bias_trim = ccdproc.trim_image(flat1_bias[:, :4096])
flat2_bias_trim = ccdproc.trim_image(flat2_bias[:, :4096])
flat3_bias_trim = ccdproc.trim_image(flat3_bias[:, :4096])
In [ ]:
flat_list = [flat1_bias_trim, flat2_bias_trim, flat3_bias_trim]
# find & combine & group apertures
ap_comb = hrs.combine_apertures(flat_list, n_jobs=10)
cheb_coefs, ap_uorder_interp = hrs.group_apertures(ap_comb, start_col=2100, order_dist=10)
# combine flat
flat_comb, flat_origin = hrs.combine_flat(flat_list, ap_uorder_interp, sat_count=45000, p=95)
flat_comb = ccdproc.CCDData(flat_comb, unit='adu')
# scattered light substraction
flat_comb_sl = hrs.substract_scattered_light(flat_comb, ap_uorder_interp, ap_width=10, shrink=.85)
flat1d = hrs.extract_1dspec(flat_comb_sl, ap_uorder_interp, ap_width=7)[0]
In [ ]:
ind_thar1 = (logt['obj']=='thar')*(logt['exp_time']==30)
ind_thar2 = (logt['obj']=='thar')*(logt['exp_time']==60)
ind_thar3 = (logt['obj']=='thar')*(logt['exp_time']==120)
thar1 = ccdproc.combine(','.join(logt['filename'][ind_thar1]), unit='adu', method='average')
thar2 = ccdproc.combine(','.join(logt['filename'][ind_thar2]), unit='adu', method='average')
thar3 = ccdproc.CCDData.read(','.join(logt['filename'][ind_thar3]), unit='adu', method='average')
thar1_bias = ccdproc.subtract_bias(thar1, bias)
thar2_bias = ccdproc.subtract_bias(thar2, bias)
thar3_bias = ccdproc.subtract_bias(thar3, bias)
thar1_bias_trim = ccdproc.trim_image(thar1_bias[:, :4096])
thar2_bias_trim = ccdproc.trim_image(thar2_bias[:, :4096])
thar3_bias_trim = ccdproc.trim_image(thar3_bias[:, :4096])
In [ ]:
""" combine thar """
thar_list = [thar1_bias_trim, thar2_bias_trim, thar3_bias_trim]
thar_comb, thar_origin = hrs.combine_flat(thar_list, ap_uorder_interp, sat_count=45000, p=100)
thar_comb = ccdproc.CCDData(thar_comb, unit='adu')
""" extract 1d thar """
#thar_comb_sl = substract_scattered_light(thar_comb, ap_uorder_interp, ap_width=10, shrink=.5)
thar1d = hrs.extract_1dspec(thar_comb, ap_uorder_interp, ap_width=7)[0]
thar1d_scaled = thar1d/flat1d
In [ ]:
# laod template thar
thar_temp_path = '/home/cham/PycharmProjects/hrs/hrs/calibration/thar_template/thar_temp_w20160120022t.fits'
wave_temp, thar_temp, order_temp = hrs.load_thar_temp(thar_temp_path)
# fix thar
thar1d_fixed = hrs.fix_thar_sat_neg(thar1d, arm=30, sat_count=500000)
thar_temp_fixed = hrs.fix_thar_sat_neg(thar_temp, arm=30, sat_count=500000)
""" initial estimation of wavelength """
shift, corr2d = hrs.thar_corr2d(thar1d_fixed, thar_temp, ytrim=(30, 70), y_shiftmax = 3, x_shiftmax=20)
print ("@Cham: the shift is ", shift)
wave_init = hrs.interpolate_wavelength_shift(wave_temp, shift, thar_temp, thar1d_fixed)
order_init = hrs.interpolate_order_shift(order_temp, shift, thar1d_fixed)
In [ ]:
thar_list = np.loadtxt('/home/cham/PycharmProjects/hrs/hrs/calibration/iraf/thar.dat')
poly_order = (3, 5)
# refine thar positions
lc_coord, lc_order, lc_thar, popt, pcov = hrs.calibration.refine_thar_positions(
wave_init, order_init, thar1d_fixed, thar_list, n_jobs=20, verbose=10)
# fit grating function
x_mini_lsq, ind_good_thar, scaler_coord, scaler_order, scaler_ml = hrs.calibration.fit_grating_equation(
lc_coord, lc_order, lc_thar, popt, pcov, poly_order=(3, 5), max_dev_threshold=1., iter=False)
# recover the wavelength (grid)
grid_coord = np.repeat(np.arange(wave_init.shape[1]).reshape(1,-1), 104, axis=0)
grid_order = order_init
wave_fitted = hrs.calibration.grating_equation_predict(
grid_coord, grid_order, x_mini_lsq, poly_order, scaler_coord, scaler_order, scaler_ml)
# recover the wavelength (thar)
lc_fitted = hrs.grating_equation_predict(
lc_coord, lc_order, x_mini_lsq, poly_order, scaler_coord, scaler_order, scaler_ml)
lc_fitted_diff = lc_fitted - lc_thar
In [25]:
""" residual check zoom in """
%matplotlib inline
fig = plt.figure(figsize=(16,8));
ax = fig.add_subplot(221)
ax.plot(lc_order,lc_fitted-lc_thar, 'b+', alpha=.8)
ax.plot(lc_order[ind_good_thar],lc_fitted[ind_good_thar]-lc_thar[ind_good_thar], 'r+', alpha=.8)
plt.xlabel("Order")
plt.axis([61, 161, -0.1, 0.1])
ax = fig.add_subplot(222)
ax.plot(lc_coord,lc_fitted-lc_thar, 'b+', alpha=.8)
ax.plot(lc_coord[ind_good_thar],lc_fitted[ind_good_thar]-lc_thar[ind_good_thar], 'r+', alpha=.8)
plt.xlabel("CCD X Coordinate")
plt.xticks(np.arange(5)*1024)
plt.axis([0, 4096, -0.1, 0.1])
ax = fig.add_subplot(212)
ax.plot(lc_thar,lc_fitted-lc_thar, 'b+', alpha=.8)
ax.plot(lc_thar[ind_good_thar],lc_fitted[ind_good_thar]-lc_thar[ind_good_thar], 'r+', alpha=.8)
plt.xlabel("True ThAr wavelength")
plt.axis([3700, 10000, -0.1, 0.1])
fig.tight_layout()
fig.show()
In [ ]:
"""
Master variables:
bias
ap_uorder_interp
flat1d
wave_fitted
"""
# 1. read
sci = ccdproc.CCDData.read(logt['filename'][41], unit='adu')
# 2. BIAS correction
sci_bias = ccdproc.subtract_bias(sci, bias)
# 3. trim
sci_bias_trim = ccdproc.trim_image(sci_bias[:, :4096])
# 4. scattered light correction
sci_bias_trim_sl = hrs.substract_scattered_light(sci_bias_trim, ap_uorder_interp, ap_width=10, shrink=.85)
# 5. extract 1D spectra
spec1d = hrs.extract_1dspec(sci_bias_trim_sl, ap_uorder_interp, ap_width=7)[0]
# 6. de-blaze
spec1d_db = spec1d/flat1d
# 7. roughly scale to 1.
spec1d_nm = spec1d_db/np.percentile(spec1d_db, 95, axis=1)[:, None]
In [32]:
plt.figure(figsize=(16, 5))
plt.plot(wave_fitted.T, spec1d_nm.T);
plt.xlim(5900, 6000)
plt.ylim(-0.1, 1.2)
Out[32]:
In [ ]: