diff --git a/README.md b/README.md new file mode 100644 index 0000000000000000000000000000000000000000..2113a2936c8986d85f5a96947a5ea54b0b343c6a --- /dev/null +++ b/README.md @@ -0,0 +1,45 @@ +# CCD 数据处理桌面应用 + +本项目是一个基于 PySide6 和 Qt6 的 CCD 数据处理桌面应用，集成了常用的天文 CCD 数据分析功能，支持图形化操作和日志输出。 + +## 主要功能 + +- **PTC/PRC**：进行 PTC/PRC 数据处理和自动绘图。 +- **Bias/Noise**：对单个 FITS 文件进行噪声分析，生成结果表格。 +- **Linearity**：线性度分析（功能待实现）。 +- **平场修正**：对单个 FITS 文件进行过扫区修正，生成修正后FITS文件。 +- **日志输出**：所有处理过程和结果均实时显示在界面右侧日志框。 + +## 操作方式 + +1. **启动程序** + ```bash + python main.py + ``` + +2. **PTC/PRC** + - 点击“PTC/PRC”按钮，弹出参数设置对话框。 + - 选择输入文件目录和输出目录，点击“开始处理”。 + - 程序会在输出目录下自动新建 `PTC_PRC_result_时间戳` 文件夹，所有结果和绘图保存在该文件夹中。 + +3. **Bias/Noise** + - 点击“Bias/Noise”按钮，弹出参数设置对话框。 + - 选择输入的 FITS 文件和输出目录（可选，默认为输入FITS文件所在目录）。 + - 生成噪声分析结果，名称为`原文件名_时间戳.tab`。 + +4. **Linearity** + - 目前为占位按钮，点击会弹窗提示“功能待实现”。 + +5. **Linearity** + - 点击“Bias/Noise”按钮，弹出参数设置对话框。 + - 选择输入的 FITS 文件和输出目录（可选，默认为输入FITS文件所在目录）。 + - 生成过扫区修正后的文件，名称为`原文件名_overscan.fits`。 + +6. **日志查看** + - 所有处理进度、结果和错误信息会实时显示在界面右侧的日志框中。 + +## 依赖环境 + +- Python 3.8+ +- PySide6 +- 其他依赖请参考 requirements.txt 或根据实际脚本需求安装（如 astropy、numpy、scipy、matplotlib 等） diff --git a/getNoise.py b/getNoise.py new file mode 100644 index 0000000000000000000000000000000000000000..b5c896d78d18f7f109bf71928ec4daeb3e6593b6 --- /dev/null +++ b/getNoise.py @@ -0,0 +1,84 @@ +#!/usr/bin/env python + +import os, sys +import time +from astropy.io import fits +import numpy as np +from scipy.stats import sigmaclip +from astropy.table import Table + +def write_log(log_msg): + + curTime = time.strftime("%Y-%m-%d %H:%M:%S", time.localtime()) + msg = f'[{curTime}] -- run getNoise.py and saved {log_msg}\n' + if os.path.exists('noise_log.txt'): + fp = open('noise_log.txt', 'a') + else: + fp = open('noise_log.txt', 'w') + + fp.write(msg) + fp.close() + +def remove_overscan(imgdata): + print('> doing overscan correction ...') + img = np.zeros((9232,9216), dtype=float) + + for i in range(4): + img[0:4616, i*1152:(i+1)*1152] = imgdata[0:4616,(i*1250+27):((i+1)*1250-71)] + img[4616:, i*1152:(i+1)*1152] = imgdata[4784:,(i*1250+27):((i+1)*1250-71)] + oc1 = np.mean(imgdata[0:4616,((i+1)*1250-71):((i+1)*1250)], axis=1).astype(float) + oc2 = np.mean(imgdata[4784:,((i+1)*1250-71):((i+1)*1250)], axis=1).astype(float) +# print('oc1.shape: {}'.format(oc1.shape)) + for j in range(4616): + img[j,i*1152:(i+1)*1152] = img[j,i*1152:(i+1)*1152] - oc1[j] + img[4616+j,i*1152:(i+1)*1152] = img[4616+j,i*1152:(i+1)*1152]- oc2[j] + for i in range(4,8): + img[0:4616, i*1152:(i+1)*1152] = imgdata[0:4616,(i*1250+71):((i+1)*1250-27)] + img[4616:, i*1152:(i+1)*1152] = imgdata[4784:,(i*1250+71):((i+1)*1250-27)] + oc1 = np.mean(imgdata[0:4616,(i*1250):(i*1250+71)], axis=1).astype(float) + oc2 = np.mean(imgdata[4784:,(i*1250):(i*1250+71)], axis=1).astype(float) + for j in range(4616): + img[j, i*1152:(i+1)*1152] = img[j, i*1152:(i+1)*1152] - oc1[j] + img[4616+j,i*1152:(i+1)*1152] = img[4616+j,i*1152:(i+1)*1152] - oc2[j] + + print('> overscan correction done!') + return img + +def main(): + fname = sys.argv[1] + imgdata = fits.getdata(fname).astype(float) + img = remove_overscan(imgdata) + + os = [] + noise_clip_4_4 = [] + noise_clip_3_3 = [] + for ch in range(16): +# print('> estimating rms for ch-{:02d}'.format(ch+1)) + os.append(ch+1) + if ch < 8: + ch_data = img[0:4616,ch*1152:(ch+1)*1152] + else: + ch_data = img[4616:, (15-ch)*1152:(16-ch)*1152] + + rms44 = np.std( sigmaclip(ch_data, 4, 4)[0] ) + rms33 = np.std( sigmaclip(ch_data, 3, 3)[0] ) + noise_clip_4_4.append( np.round(rms44,4) ) + noise_clip_3_3.append( np.round(rms33,4) ) + + tab = Table() + tab['OS'] = os + tab['Noise_Clip_4_4'] = noise_clip_4_4 + tab['Noise_Clip_3_3'] = noise_clip_3_3 + print(tab) + curTime = time.strftime("%Y-%m-%dT%H:%M:%S", time.localtime()) + + tab_name = fname.replace('.fits','') + '_{}.tab'.format(curTime) + tab.write(tab_name, format='ipac', overwrite=True) + write_log(tab_name) + +# print('\n# Channel noise(clip,4,4) noise(clip,3,3)\n') +# for i in range(16): +# print('> OS-{:02d} {:6.3f} {:6.3f}'.format(i, noise_clip_4_4[i], noise_clip_3_3[i])) + +if __name__ == '__main__': + main() diff --git a/main.py b/main.py new file mode 100644 index 0000000000000000000000000000000000000000..b0048e1f625cb2702d140d041979d58f7f7ccf6f --- /dev/null +++ b/main.py @@ -0,0 +1,438 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +数据处理GUI应用程序 +基于PySide6和Qt6的桌面应用 +""" + +import sys +from PySide6.QtWidgets import (QApplication, QMainWindow, QWidget, QVBoxLayout, + QHBoxLayout, QLabel, QLineEdit, QPushButton, + QFileDialog, QMessageBox, QRadioButton, QButtonGroup, QTextEdit, QDialog) +from PySide6.QtCore import Qt, QThread, Signal +from PySide6.QtGui import QFont, QIcon +import datetime +import threading +import os +from processPTC import processPTC +import importlib.util +import subprocess +from astropy.io import fits +import numpy as np +from getNoise import remove_overscan + + +class DataProcessingGUI(QMainWindow): + """数据处理GUI主窗口""" + + def __init__(self): + super().__init__() + self.init_ui() + + def init_ui(self): + """初始化用户界面""" + # 设置窗口标题和大小 + self.setWindowTitle("CCD数据处理工具") + self.setGeometry(100, 100, 900, 500) + + # 创建中央部件 + central_widget = QWidget() + self.setCentralWidget(central_widget) + + # ====== 新增：主布局为左右结构 ====== + main_h_layout = QHBoxLayout(central_widget) + # 左侧为原有内容 + left_widget = QWidget() + main_layout = QVBoxLayout(left_widget) + main_layout.setSpacing(20) + main_layout.setContentsMargins(30, 30, 30, 30) + + # 添加标题 + title_label = QLabel("CCD数据处理工具") + title_label.setAlignment(Qt.AlignCenter) + title_font = QFont() + title_font.setPointSize(18) + title_font.setBold(True) + title_label.setFont(title_font) + main_layout.addWidget(title_label) + + # 删除文件路径和MONITOR路径输入区域 + + # 添加一些间距 + main_layout.addSpacing(20) + + # 单选按钮区域 + radio_layout = QHBoxLayout() + radio_label = QLabel("选择CCD类型:") + radio_label.setMinimumWidth(100) + + # 创建单选按钮组 + self.radio_group = QButtonGroup() + self.radio_290_00 = QRadioButton("290-99") + self.radio_47_20 = QRadioButton("47-20") + + # 将单选按钮添加到按钮组 + self.radio_group.addButton(self.radio_290_00) + self.radio_group.addButton(self.radio_47_20) + + # 默认选择第一个选项 + self.radio_290_00.setChecked(True) + + radio_layout.addWidget(radio_label) + radio_layout.addWidget(self.radio_290_00) + radio_layout.addWidget(self.radio_47_20) + radio_layout.addStretch() # 添加弹性空间 + + main_layout.addLayout(radio_layout) + + # 添加一些间距 + main_layout.addSpacing(10) + + # 处理按钮区域，分两行布局 + buttons_layout1 = QHBoxLayout() + buttons_layout1.setSpacing(15) + buttons_layout2 = QHBoxLayout() + buttons_layout2.setSpacing(15) + # 第一行按钮 + self.process_button1 = QPushButton("PTC/PRC") + self.process_button1.setMinimumHeight(40) + self.process_button1.clicked.connect(self.process_action1) + buttons_layout1.addWidget(self.process_button1) + self.process_button2 = QPushButton("Bias/Noise") + self.process_button2.setMinimumHeight(40) + self.process_button2.clicked.connect(self.process_action2) + buttons_layout1.addWidget(self.process_button2) + # 第二行按钮 + self.process_button3 = QPushButton("Linearity") + self.process_button3.setMinimumHeight(40) + self.process_button3.clicked.connect(self.process_action3) + buttons_layout2.addWidget(self.process_button3) + self.process_button4 = QPushButton("过扫区修正") + self.process_button4.setMinimumHeight(40) + self.process_button4.clicked.connect(self.process_action4) + buttons_layout2.addWidget(self.process_button4) + main_layout.addLayout(buttons_layout1) + main_layout.addLayout(buttons_layout2) + + # 添加弹性空间 + main_layout.addStretch() + # 右侧为日志框 + self.log_text_edit = QTextEdit() + self.log_text_edit.setReadOnly(True) + self.log_text_edit.setMinimumWidth(250) + self.log_text_edit.setPlaceholderText("日志输出...") + # 添加到主布局 + main_h_layout.addWidget(left_widget, stretch=3) + main_h_layout.addWidget(self.log_text_edit, stretch=2) + + def browse_path(self, line_edit): + """浏览文件夹对话框""" + folder_path = QFileDialog.getExistingDirectory( + self, + "选择文件夹", + "", + QFileDialog.ShowDirsOnly + ) + if folder_path: + line_edit.setText(folder_path) + + def get_selected_type(self): + """获取当前选中的类型""" + if self.radio_290_00.isChecked(): + return "290-99" + elif self.radio_47_20.isChecked(): + return "47-20" + else: + return "未选择" + + def process_action1(self): + """PTC/PRC处理按钮动作：弹出参数设置对话框""" + dlg = PTCDialog(self) + if dlg.exec() == QDialog.Accepted: + fits_dir, out_base_dir = dlg.get_dirs() + if not fits_dir or not out_base_dir: + QMessageBox.warning(self, "警告", "请输入输入和输出目录！") + return + now_str = datetime.datetime.now().strftime("%Y%m%d_%H%M%S") + out_dir = os.path.join(out_base_dir, f"PTC_PRC_result_{now_str}") + self.log_text_edit.append(f"开始PTC/PRC处理：\n输入目录: {fits_dir}\n输出目录: {out_dir}") + self.ptc_worker = PTCWorker(fits_dir, out_dir) + self.ptc_worker.finished_signal.connect(self.log_text_edit.append) + self.ptc_worker.result_signal.connect(self.run_plot_prc) + self.ptc_worker.start() + + def process_action2(self): + """Bias/Noise按钮动作：弹出参数设置对话框""" + dlg = NoiseDialog(self) + if dlg.exec() == QDialog.Accepted: + fits_path, out_dir = dlg.get_paths() + if not fits_path: + QMessageBox.warning(self, "警告", "请输入FITS文件路径！") + return + # 如果输出目录为空，默认在输入文件同目录下 + if not out_dir: + out_dir = os.path.dirname(fits_path) + self.log_text_edit.append(f"开始噪声分析：\n输入文件: {fits_path}\n输出目录: {out_dir}") + self.noise_worker = NoiseWorker(fits_path, out_dir) + self.noise_worker.log_signal.connect(self.log_text_edit.append) + self.noise_worker.finished_signal.connect(self.log_text_edit.append) + self.noise_worker.start() + + def process_action3(self): + """线性度按钮动作（待实现）""" + QMessageBox.information(self, "信息", "线性度按钮被点击，功能待实现。") + + def process_action4(self): + """过扫区修正按钮动作""" + dlg = OverscanDialog(self) + if dlg.exec() == QDialog.Accepted: + file_path, out_dir = dlg.get_paths() + if not file_path: + QMessageBox.warning(self, "警告", "请输入FITS文件路径！") + return + import os + if not out_dir: + out_dir = os.path.dirname(file_path) + base_name = os.path.basename(file_path) + name, ext = os.path.splitext(base_name) + out_file = os.path.join(out_dir, f"{name}_overscan{ext}") + # 读取FITS数据 + try: + with fits.open(file_path) as hdul: + imgdata = hdul[0].data.astype(float) + header = hdul[0].header + except Exception as e: + QMessageBox.critical(self, "错误", f"读取FITS文件失败：{e}") + return + # 过扫区修正 + try: + self.log_text_edit.append(f"开始修正过扫区...") + img_corr = remove_overscan(imgdata) + except Exception as e: + QMessageBox.critical(self, "错误", f"过扫区修正失败：{e}") + return + # 保存修正后的数据 + try: + fits.writeto(out_file, img_corr, header, overwrite=True) + self.log_text_edit.append(f"修正后的FITS文件已保存：\n{out_file}") + except Exception as e: + QMessageBox.critical(self, "错误", f"保存FITS文件失败：{e}") + + def run_plot_prc(self, out_dir): + try: + import importlib.util + plot_prc_path = os.path.join(os.path.dirname(__file__), 'plot_prc.py') + spec = importlib.util.spec_from_file_location("plot_prc", plot_prc_path) + plot_prc = importlib.util.module_from_spec(spec) + spec.loader.exec_module(plot_prc) + plot_prc.plotPRC(out_dir) + self.log_text_edit.append("绘图完成，所有结果已保存到输出文件夹！") + except Exception as e: + self.log_text_edit.append(f"绘图出错: {e}") + + +class PTCWorker(QThread): + finished_signal = Signal(str) + result_signal = Signal(str) + + def __init__(self, fits_dir, out_dir): + super().__init__() + self.fits_dir = fits_dir + self.out_dir = out_dir + + def run(self): + try: + script_path = os.path.join(os.path.dirname(__file__), 'processPTC.py') + cmd = [sys.executable, script_path, self.fits_dir, self.out_dir] + process = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True) + for line in process.stdout: + self.finished_signal.emit(line.rstrip()) + process.wait() + if process.returncode == 0: + self.result_signal.emit(self.out_dir) + self.finished_signal.emit("PTC/PRC数据处理完成！") + else: + self.finished_signal.emit(f"PTC/PRC处理出错，返回码: {process.returncode}") + except Exception as e: + self.finished_signal.emit(f"处理出错: {e}") + + +class NoiseWorker(QThread): + log_signal = Signal(str) + finished_signal = Signal(str) + def __init__(self, fits_path, out_dir): + super().__init__() + self.fits_path = fits_path + self.out_dir = out_dir + def run(self): + try: + script_path = os.path.join(os.path.dirname(__file__), 'getNoise.py') + # 传递输出目录参数给getNoise.py（如需支持） + cmd = [sys.executable, script_path, self.fits_path, self.out_dir] + process = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True) + for line in process.stdout: + self.log_signal.emit(line.rstrip()) + process.wait() + if process.returncode == 0: + self.finished_signal.emit("噪声分析完成！") + else: + self.finished_signal.emit(f"噪声分析出错，返回码: {process.returncode}") + except Exception as e: + self.finished_signal.emit(f"噪声分析出错: {e}") + + +class PTCDialog(QDialog): + def __init__(self, parent=None): + super().__init__(parent) + self.setWindowTitle("PTC/PRC参数设置") + self.setMinimumWidth(400) + layout = QVBoxLayout() + # 输入目录 + in_layout = QHBoxLayout() + in_label = QLabel("输入目录:") + self.in_edit = QLineEdit() + in_browse = QPushButton("浏览") + in_browse.clicked.connect(self.browse_in) + in_layout.addWidget(in_label) + in_layout.addWidget(self.in_edit) + in_layout.addWidget(in_browse) + layout.addLayout(in_layout) + # 输出目录 + out_layout = QHBoxLayout() + out_label = QLabel("输出目录:") + self.out_edit = QLineEdit() + out_browse = QPushButton("浏览") + out_browse.clicked.connect(self.browse_out) + out_layout.addWidget(out_label) + out_layout.addWidget(self.out_edit) + out_layout.addWidget(out_browse) + layout.addLayout(out_layout) + # 开始处理按钮 + btn_layout = QHBoxLayout() + btn_layout.addStretch() + self.start_btn = QPushButton("开始处理") + self.start_btn.clicked.connect(self.accept) + btn_layout.addWidget(self.start_btn) + layout.addLayout(btn_layout) + self.setLayout(layout) + def browse_in(self): + dir_path = QFileDialog.getExistingDirectory(self, "选择输入文件目录") + if dir_path: + self.in_edit.setText(dir_path) + def browse_out(self): + dir_path = QFileDialog.getExistingDirectory(self, "选择输出目录") + if dir_path: + self.out_edit.setText(dir_path) + def get_dirs(self): + return self.in_edit.text().strip(), self.out_edit.text().strip() + + +class NoiseDialog(QDialog): + def __init__(self, parent=None): + super().__init__(parent) + self.setWindowTitle("Bias/Noise参数设置") + self.setMinimumWidth(400) + layout = QVBoxLayout() + # 输入文件 + in_layout = QHBoxLayout() + in_label = QLabel("输入FITS文件:") + self.in_edit = QLineEdit() + in_browse = QPushButton("浏览") + in_browse.clicked.connect(self.browse_in) + in_layout.addWidget(in_label) + in_layout.addWidget(self.in_edit) + in_layout.addWidget(in_browse) + layout.addLayout(in_layout) + # 输出目录 + out_layout = QHBoxLayout() + out_label = QLabel("输出目录:") + self.out_edit = QLineEdit() + out_browse = QPushButton("浏览") + out_browse.clicked.connect(self.browse_out) + out_layout.addWidget(out_label) + out_layout.addWidget(self.out_edit) + out_layout.addWidget(out_browse) + layout.addLayout(out_layout) + # 开始处理按钮 + btn_layout = QHBoxLayout() + btn_layout.addStretch() + self.start_btn = QPushButton("开始处理") + self.start_btn.clicked.connect(self.accept) + btn_layout.addWidget(self.start_btn) + layout.addLayout(btn_layout) + self.setLayout(layout) + def browse_in(self): + file_path, _ = QFileDialog.getOpenFileName(self, "选择FITS文件", "", "FITS Files (*.fits);;All Files (*)") + if file_path: + self.in_edit.setText(file_path) + def browse_out(self): + dir_path = QFileDialog.getExistingDirectory(self, "选择输出目录") + if dir_path: + self.out_edit.setText(dir_path) + def get_paths(self): + return self.in_edit.text().strip(), self.out_edit.text().strip() + + +class OverscanDialog(QDialog): + def __init__(self, parent=None): + super().__init__(parent) + self.setWindowTitle("过扫区修正参数设置") + self.setMinimumWidth(400) + layout = QVBoxLayout() + # 输入文件 + in_layout = QHBoxLayout() + in_label = QLabel("输入FITS文件:") + self.in_edit = QLineEdit() + in_browse = QPushButton("浏览") + in_browse.clicked.connect(self.browse_in) + in_layout.addWidget(in_label) + in_layout.addWidget(self.in_edit) + in_layout.addWidget(in_browse) + layout.addLayout(in_layout) + # 输出目录 + out_layout = QHBoxLayout() + out_label = QLabel("输出目录:") + self.out_edit = QLineEdit() + out_browse = QPushButton("浏览") + out_browse.clicked.connect(self.browse_out) + out_layout.addWidget(out_label) + out_layout.addWidget(self.out_edit) + out_layout.addWidget(out_browse) + layout.addLayout(out_layout) + # 开始处理按钮 + btn_layout = QHBoxLayout() + btn_layout.addStretch() + self.start_btn = QPushButton("开始处理") + self.start_btn.clicked.connect(self.accept) + btn_layout.addWidget(self.start_btn) + layout.addLayout(btn_layout) + self.setLayout(layout) + def browse_in(self): + file_path, _ = QFileDialog.getOpenFileName(self, "选择FITS文件", "", "FITS Files (*.fits *.fit)") + if file_path: + self.in_edit.setText(file_path) + def browse_out(self): + dir_path = QFileDialog.getExistingDirectory(self, "选择输出目录") + if dir_path: + self.out_edit.setText(dir_path) + def get_paths(self): + return self.in_edit.text().strip(), self.out_edit.text().strip() + + +def main(): + """主函数""" + app = QApplication(sys.argv) + + # 设置应用程序样式 + app.setStyle('Fusion') + + # 创建并显示主窗口 + window = DataProcessingGUI() + window.show() + + # 运行应用程序 + sys.exit(app.exec()) + + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/plot_prc.py b/plot_prc.py new file mode 100755 index 0000000000000000000000000000000000000000..9a2bb499ad7db615cb9070fe854ac6506123aa7b --- /dev/null +++ b/plot_prc.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python + +import os, sys +import numpy as np +from astropy.table import Table +import matplotlib.pylab as plt +from scipy import odr + +def plotPRC(out_dir): + tab = Table.read(f'{out_dir}/ptcdata.tab', format='ipac') + + model = odr.polynomial(1) + + fig = plt.figure(figsize=(20,40)) + + for ch in range(16): + row = ch // 4 + pos1 = row*8 + ch - row*4 + 1 + pos2 = pos1 + 4 + + # print(f'pos1: {pos1}, pos2: {pos2}') + + plt.subplot(8,4,pos1) + + index = tab['channel'] == ch+1 + + t_exp = tab[index]['t_exp'] + flux = tab[index]['flux_4'] + + # fit prc and calculate residuals + data = odr.Data(t_exp, flux, we=1/(t_exp**2), wd=1/(flux**2) ) + # data = odr.Data(t_exp, flux, we=1/(t_exp**2) ) + # data = odr.Data(t_exp, flux) + fitter = odr.ODR(data, model, ifixb=[1,1], beta0=[2500,800]) + res = fitter.run() + poly = np.poly1d(res.beta[::-1]) + nonlin = (flux-poly(t_exp)) / poly(t_exp) * 100 + + plt.scatter(t_exp, flux/1e3, zorder=1, alpha=0.85, color='m', edgecolors='none', s=30) + # plt.plot(t_exp, flux/1e3, 'g--o', alpha=0.75) + plt.plot(t_exp, poly(t_exp)/1e3, 'b--', alpha=0.75) + plt.xlabel('Expourse Time (s)') + plt.ylabel('Flux (kDN)') + plt.ylim([-1,6.1e1]) + # plt.title('channel {}'.format(ch+1)) + plt.text(5,55,'channel {}'.format(ch+1), color='limegreen', fontsize=14) + plt.grid() + + + plt.subplot(8,4,pos2) + print(res.beta) + plt.scatter(t_exp, nonlin, alpha=0.75, color='darkorange', edgecolors='none', s=30) + plt.hlines(y=1.8, xmin=t_exp.min(), xmax=t_exp.max(), linestyle=':', color='dimgrey') + plt.hlines(y=-1.8, xmin=t_exp.min(), xmax=t_exp.max(), linestyle=':', color='dimgrey') + plt.plot(t_exp, np.zeros(len(t_exp)), '--', color='grey') + plt.ylim([-3.25,3.25]) + plt.xlabel('Expourse Time (s)') + plt.ylabel('PRC Nonlinearity (%)') + # plt.title('channel {}'.format(ch+1)) + plt.text(5,2.5,'channel {}'.format(ch+1), color='brown', fontsize=14) + plt.grid() + + plt.tight_layout() + plt.savefig(f'{out_dir}/prc.pdf') + plt.close() + +if __name__ == '__main__': + plotPRC(sys.argv[1]) diff --git a/processPTC.py b/processPTC.py new file mode 100755 index 0000000000000000000000000000000000000000..983b25b216a167c82da98c7114268a507af3b7f2 --- /dev/null +++ b/processPTC.py @@ -0,0 +1,327 @@ +#!/usr/bin/env python + +import glob, os, sys +from tqdm import trange +from astropy.table import Table +from astropy.io import fits +from scipy.stats import sigmaclip +from scipy import odr +import numpy as np +from matplotlib import pyplot as plt + +def remove_overscan(imgdata, ch): + # print('> removing 71 cols overscan ...') + img = np.zeros((4616,1152), dtype=float) + +# for i in range(4): +# img[0:4616, i*1152:(i+1)*1152] = imgdata[0:4616,(i*1250+27):((i+1)*1250-71)] +# img[4616:, i*1152:(i+1)*1152] = imgdata[4784:,(i*1250+27):((i+1)*1250-71)] +# oc1 = np.median(imgdata[0:4616,((i+1)*1250-71):((i+1)*1250)], axis=1).astype(float) +# oc2 = np.median(imgdata[4784:,((i+1)*1250-71):((i+1)*1250)], axis=1).astype(float) +# # print('oc1.shape: {}'.format(oc1.shape)) +# for j in range(4616): +# img[j,i*1152:(i+1)*1152] = img[j,i*1152:(i+1)*1152] - oc1[j] +# img[4616+j,i*1152:(i+1)*1152] = img[4616+j,i*1152:(i+1)*1152]- oc2[j] +# for i in range(4,8): +# img[0:4616, i*1152:(i+1)*1152] = imgdata[0:4616,(i*1250+71):((i+1)*1250-27)] +# img[4616:, i*1152:(i+1)*1152] = imgdata[4784:,(i*1250+71):((i+1)*1250-27)] +# oc1 = np.median(imgdata[0:4616,(i*1250):(i*1250+71)], axis=1).astype(float) +# oc2 = np.median(imgdata[4784:,(i*1250):(i*1250+71)], axis=1).astype(float) +# for j in range(4616): +# img[j, i*1152:(i+1)*1152] = img[j, i*1152:(i+1)*1152] - oc1[j] +# img[4616+j,i*1152:(i+1)*1152] = img[4616+j,i*1152:(i+1)*1152] - oc2[j] + + if ch < 4: + img = imgdata[0:4616,(ch*1250+27):((ch+1)*1250-71)] + oc1 = np.median(imgdata[0:4616,((ch+1)*1250-71):((ch+1)*1250)], axis=1).astype(float) + for j in range(4616): + img[j,:] = img[j,:] - oc1[j] + elif ch < 8: + img = imgdata[0:4616,(ch*1250+71):((ch+1)*1250-27)] + oc1 = np.median(imgdata[0:4616,(ch*1250):(ch*1250+71)], axis=1).astype(float) + for j in range(4616): + img[j, :] = img[j,:] - oc1[j] + elif ch < 12: + img = imgdata[4784:,((15-ch)*1250+71):((16-ch)*1250-27)] + oc2 = np.median(imgdata[4784:,((15-ch)*1250):((15-ch)*1250+71)], axis=1).astype(float) + for j in range(4616): + img[j,:] = img[j,:] - oc2[j] + else: + img = imgdata[4784:,((15-ch)*1250+27):((16-ch)*1250-71)] + oc2 = np.median(imgdata[4784:,((16-ch)*1250-71):((16-ch)*1250)], axis=1).astype(float) + for j in range(4616): + img[j,:] = img[j,:]- oc2[j] + + return img + +def bindata(image, n): + ny, nx = image.shape + nx2 = nx // n + x_res = nx - nx2 * n + x1 = x_res // 2 + x2 = x_res - x1 + ny2 = ny // n + y_res = ny - ny2 * n + y1 = y_res // 2 + y2 = y_res - y1 + data = image[y1:ny-y2, x1:nx-x2] + # data = data.reshape((ny2, n, nx2, n)).sum(axis=1).sum(axis=-1) + data = data.reshape((ny2, n, nx2, n)).mean(axis=1).mean(axis=-1) + return data + +#blist = glob.glob("fits/bias/*.fits") + +#bias = None +#bcnt = 0 +#for i in range(len(blist)): +# if i==0: +# bias = fits.getdata(blist[i]).astype(float) +# else: +# bias += fits.getdata(blist[i]).astype(float) + +# bcnt += 1 + +#bias = bias / bcnt + + +def processPTC(fits_dir, out_dir): + + if os.path.exists(out_dir) is False: + os.system(f'mkdir {out_dir}') + + ftable = f'{out_dir}/flat_filelist.tab' + fdata = f'{out_dir}/ptcdata.tab' + + if os.path.exists(ftable) is False: + flist = glob.glob(f'{fits_dir}/*ptc*.fits') + t_exp = [] + seq = [] + for f in flist: + name = os.path.basename(f) + tmp = name.split('_') +# print(tmp) + t_exp.append(float(tmp[-4].replace('s', ''))) + seq.append(int(tmp[-2].replace('.fits', ''))) + tab = Table() + tab['name'] = flist + tab['t_exp'] = t_exp + tab['seq'] = seq + tab.sort(['t_exp', 'seq']) + tab.write(ftable, format='ipac', overwrite=True) + +# print('#########') +# sys.exit(0) + + clip_n = 4 + tab = Table.read(ftable, format='ipac') + if os.path.exists(fdata) is False: + t_exp = np.unique(tab['t_exp']) + t_exp = t_exp[t_exp < 35] + t = [] + chans = [] + flux_1 = [] + var_1 = [] + flux_2 = [] + var_2 = [] + flux_4 = [] + var_4 = [] + flux_8 = [] + var_8 = [] + for chan in range(16): + for i in trange(len(t_exp), desc='channel {}'.format(chan+1)): + index = tab['t_exp'] == t_exp[i] + f1, f2 = tab['name'][index][0], tab['name'][index][1] + with fits.open(f1) as img1, fits.open(f2) as img2: + t.append(t_exp[i]) + chans.append(chan+1) + # m1 = img1[0].data.astype(float)[:4617, chan*1250+72:chan*1250+1224] + # m2 = img2[0].data.astype(float)[:4617, chan*1250+72:chan*1250+1224] + + # if chan < 4: + # m1 = img1[0].data.astype(float)[2:4617, chan*1250+30:chan*1250+30+1170] - bias[2:4617, chan*1250+30:chan*1250+30+1170]*0 + # m2 = img2[0].data.astype(float)[2:4617, chan*1250+30:chan*1250+30+1170] - bias[2:4617, chan*1250+30:chan*1250+30+1170]*0 + # elif chan >= 4 and chan < 8: + # m1 = img1[0].data.astype(float)[2:4617, chan*1250+72:chan*1250+72+1170] - bias[2:4617, chan*1250+72:chan*1250+72+1170]*0 + # m2 = img2[0].data.astype(float)[2:4617, chan*1250+72:chan*1250+72+1170] - bias[2:4617, chan*1250+72:chan*1250+72+1170]*0 + # elif chan >= 8 and chan < 12: + # m1 = img1[0].data.astype(float)[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170] - bias[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170]*0 + # m2 = img2[0].data.astype(float)[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170] - bias[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170]*0 + # else: + # m1 = img1[0].data.astype(float)[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170] - bias[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170]*0 + # m2 = img2[0].data.astype(float)[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170] - bias[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170]*0 + +# if chan < 4: +# m1 = img1[0].data.astype(float)[2:4617, chan*1250+30:chan*1250+30+1170] #- bias[2:4617, chan*1250+30:chan*1250+30+1170]*0 +# m2 = img2[0].data.astype(float)[2:4617, chan*1250+30:chan*1250+30+1170] #- bias[2:4617, chan*1250+30:chan*1250+30+1170]*0 +# elif chan >= 4 and chan < 8: +# m1 = img1[0].data.astype(float)[2:4617, chan*1250+72:chan*1250+72+1170] #- bias[2:4617, chan*1250+72:chan*1250+72+1170]*0 +# m2 = img2[0].data.astype(float)[2:4617, chan*1250+72:chan*1250+72+1170] #- bias[2:4617, chan*1250+72:chan*1250+72+1170]*0 +# elif chan >= 8 and chan < 12: +# m1 = img1[0].data.astype(float)[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170] #- bias[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170]*0 +# m2 = img2[0].data.astype(float)[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170] #- bias[4785:, (15-chan)*1250+72:(15-chan)*1250+72+1170]*0 +# else: +# m1 = img1[0].data.astype(float)[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170] #- bias[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170]*0 +# m2 = img2[0].data.astype(float)[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170] #- bias[4785:, (15-chan)*1250+30:(15-chan)*1250+30+1170]*0 + + m1 = remove_overscan(img1[0].data, chan) + m2 = remove_overscan(img2[0].data, chan) + + flux_1.append(np.median(m1 + m2) / 2) + var_1.append(sigmaclip(m1 - m2, clip_n, clip_n)[0].var() / 2) + + mm1 = bindata(m1, 2) + mm2 = bindata(m2, 2) + flux_2.append(np.median(mm1 + mm2) / 2) + var_2.append(sigmaclip(mm1 - mm2, clip_n, clip_n)[0].var() / 2 * 4) + + mm1 = bindata(m1, 4) + mm2 = bindata(m2, 4) + flux_4.append(np.median(mm1 + mm2) / 2) + var_4.append(sigmaclip(mm1 - mm2, clip_n, clip_n)[0].var() / 2 * 16) + + mm1 = bindata(m1, 8) + mm2 = bindata(m2, 8) + flux_8.append(np.median(mm1 + mm2) / 2) + var_8.append(sigmaclip(mm1 - mm2, clip_n, clip_n)[0].var() / 2 * 64) + + cat = Table() + cat['t_exp'] = t + cat['channel'] = chans + cat['flux_1'] = flux_1 + cat['var_1'] = var_1 + cat['flux_2'] = flux_2 + cat['var_2'] = var_2 + cat['flux_4'] = flux_4 + cat['var_4'] = var_4 + cat['flux_8'] = flux_8 + cat['var_8'] = var_8 + cat.sort(['channel', 't_exp']) + cat.write(fdata, format='ipac', overwrite=True) + + + cat = Table.read(fdata, format='ipac') + for chan in np.unique(cat['channel']): + index = cat['channel'] == chan + ind = np.argmin(cat['t_exp'][index]) + cat['flux_1'][index] = cat['flux_1'][index] - cat['flux_1'][index][ind] + cat['var_1'][index] = cat['var_1'][index] - cat['var_1'][index][ind] + cat['flux_2'][index] = cat['flux_2'][index] - cat['flux_2'][index][ind] + cat['var_2'][index] = cat['var_2'][index] - cat['var_2'][index][ind] + cat['flux_4'][index] = cat['flux_4'][index] - cat['flux_4'][index][ind] + cat['var_4'][index] = cat['var_4'][index] - cat['var_4'][index][ind] + cat['flux_8'][index] = cat['flux_8'][index] - cat['flux_8'][index][ind] + cat['var_8'][index] = cat['var_8'][index] - cat['var_8'][index][ind] + index = cat['t_exp'] > 0 + cat = cat[index] + + model = odr.polynomial(1) + plt.figure(figsize=(20, 40)) + output_gains = [] + output_chans = [] + output_bins = [] + for chan in range(16): + index = cat['channel'] == chan+1 + bins = ['1', '2', '4', '8'] + colors = ['red', 'orange', 'blue', 'green'] + handlers = [] + txts = [] + row = chan // 4 + pos1 = row * 8 + chan - row * 4 + 1 + pos2 = pos1 + 4 + for i in range(4): + if i != 2: + pass + else: + x = cat['flux_'+bins[i]][index] + y = cat['var_'+bins[i]][index] + print('x = ', x) + print('y = ', y) + idx = x > 0 + x = x[idx] + y = y[idx] + fw = x[np.argmax(y)] + ind = (x > fw * 0.1) * (x < fw * 0.8) + data = odr.Data(x[ind], y[ind], we=1/(x[ind]**2), wd=1/(y[ind]**2)) + fitter = odr.ODR(data, model, ifixb=[0, 1], beta0=[0, 1]) + res = fitter.run() + poly = np.poly1d(res.beta[::-1]) + nonlin = (y - poly(x)) / poly(x) * 100 + + plt.subplot(8, 4, pos1) + plt.scatter(x[~ind], y[~ind], color=colors[i], zorder=1, alpha=0.25, edgecolors='none', s=20) + s = plt.scatter(x[ind], y[ind], color=colors[i], zorder=2, edgecolors='none', s=20) + l, = plt.plot(x, poly(x), color=colors[i], zorder=3, ls='--') + handlers.append(s) + txts.append('bin {}: gain = {:.3f} e-/DN'.format(bins[i], 1/res.beta[1])) + + output_gains.append(1/res.beta[1]) + output_chans.append(chan+1) + output_bins.append(bins[i]) + + plt.subplot(8, 4, pos2) + plt.scatter(x[~ind], nonlin[~ind], color=colors[i], zorder=2, alpha=0.25, edgecolors='none', s=20) + plt.scatter(x[ind], nonlin[ind], color=colors[i], zorder=3, edgecolors='none', s=20) + plt.plot(x, np.zeros(len(x)), color='grey', zorder=1) + + +# x = cat['flux_'+bins[i]][index] +# y = cat['var_'+bins[i]][index] +# print('x = ', x) +# print('y = ', y) +# idx = x > 0 +# x = x[idx] +# y = y[idx] +# fw = x[np.argmax(y)] +# ind = (x > fw * 0.1) * (x < fw * 0.8) +# data = odr.Data(x[ind], y[ind], we=1/(x[ind]**2), wd=1/(y[ind]**2)) +# fitter = odr.ODR(data, model, ifixb=[0, 1], beta0=[0, 1]) +# res = fitter.run() +# poly = np.poly1d(res.beta[::-1]) +# nonlin = (y - poly(x)) / poly(x) * 100 + +# plt.subplot(8, 4, pos1) +# plt.scatter(x[~ind], y[~ind], color=colors[i], zorder=1, alpha=0.25, edgecolors='none', s=20) +# s = plt.scatter(x[ind], y[ind], color=colors[i], zorder=2, edgecolors='none', s=20) +# l, = plt.plot(x, poly(x), color=colors[i], zorder=3, ls='--') +# handlers.append(s) +# txts.append('bin {}: gain = {:.3f} e-/DN'.format(bins[i], 1/res.beta[1])) + +# output_gains.append(1/res.beta[1]) +# output_chans.append(chan+1) +# output_bins.append(bins[i]) + +# plt.subplot(8, 4, pos2) +# plt.scatter(x[~ind], nonlin[~ind], color=colors[i], zorder=2, alpha=0.25, edgecolors='none', s=20) +# plt.scatter(x[ind], nonlin[ind], color=colors[i], zorder=3, edgecolors='none', s=20) +# plt.plot(x, np.zeros(len(x)), color='grey', zorder=1) + + plt.subplot(8, 4, pos1) + plt.xlabel('flux (DN)') + plt.ylabel('variance (DN^2)') + plt.title('channel {}'.format(chan+1)) + # plt.text(4e4,6e3,'channel {}'.format(chan+1), color='limegreen', fontsize=14) + plt.legend(handlers, txts) + plt.grid(True) + + plt.subplot(8, 4, pos2) + # plt.xlabel('flux (ADU)') + plt.xlabel('flux (DN)') + plt.ylabel('non-linearity (%)') + plt.title('channel {}'.format(chan+1)) + # plt.text(4e4,8,'channel {}'.format(chan+1), color='brown', fontsize=14) + plt.ylim([-10, 10]) + plt.grid(True) + + plt.tight_layout() + plt.savefig(f'{out_dir}/ptc.pdf') + plt.close() + tab = Table() + tab['channel'] = output_chans + tab['bin'] = output_bins + tab['gain'] = output_gains + tab.write(f'{out_dir}/gain.tab', format='ipac', overwrite=True) + + +if __name__ == '__main__': + processPTC(sys.argv[1], sys.argv[2]) +