diff --git a/BKGscript.py b/BKGscript.py new file mode 100644 index 0000000000000000000000000000000000000000..baab3e5382630b7bfa95ec32f356eb223727eb19 --- /dev/null +++ b/BKGscript.py @@ -0,0 +1,120 @@ +from concurrent import futures +from glob import glob +from astropy.io import fits +import os +import numpy as np +from time import time +from concurrent.futures import ThreadPoolExecutor + +A, B = 2, 2 + + +def split(data: np.ndarray, i: int) -> np.ndarray: + h, w = data.shape + x, y = h//A, w//B + n, m = i//B, i % B + data = data[n*x:(n+1)*x, + m*y:(m+1)*y].copy() + return data + + +def concatenate(data_list): + data_u = np.concatenate(data_list[:B], axis=1) + data_d = np.concatenate(data_list[B:], axis=1) + data = np.concatenate([data_u, data_d], axis=0) + return data + + +def array_combine(ndarray, mode="mean") -> np.ndarray: + """ Function to combine 3-D data array + + Parameters + ---------- + ndarray: array, input data cube (3D) + model: mean, median, sum, mean_clip, median_clip, default is mean + """ + if mode == "median": + array = np.median(ndarray, axis=0) + elif mode == "median_clip": + ndarray = np.sort(ndarray, axis=0)[1:-1] + array = np.median(ndarray, axis=0) + elif mode == "sum": + array = np.sum(ndarray, axis=0) + elif mode == "mean": + array = np.mean(ndarray, axis=0) + elif mode == "mean_clip": + ndarray = np.sort(ndarray, axis=0)[1:-1] + array = np.mean(ndarray, axis=0) + return array + + +def load_bias(path: str, i: int) -> np.ndarray: + with fits.open(path) as hdul: + du = hdul[1].data + du = du.astype(int) + du = split(du, i) + return du + + +def load_dark(path: str, i: int, bias: np.ndarray) -> np.ndarray: + with fits.open(path) as hdul: + du = hdul[1].data + hu = hdul[0].header + du = du.astype(int) + du = du - bias + du = du / hu["EXPTIME"] + du = split(du, i) + return du + + +def load_flat(path: str, i: int, bias: np.ndarray, dark: np.ndarray) -> np.ndarray: + with fits.open(path) as hdul: + du = hdul[1].data + hu = hdul[0].header + du = du.astype(int) + du = du - bias - dark * hu["EXPTIME"] + du = du / hu["EXPTIME"] + du = du / np.median(du) + du = split(du, i) + return du + + +def save_fits(func, mode: str, path_list, save_path, *args) -> np.ndarray: + ch_list = [] + with ThreadPoolExecutor() as tpe: + for i in range(A*B): + futures = [tpe.submit(func, path, i, *args) for path in path_list] + du_list = [future.result() for future in futures] + du = array_combine(du_list, mode) + ch_list.append(du) + du = concatenate(ch_list) + du_output = os.path.basename(path_list[0]) + du_output = du_output.replace("raw", "combine") + du_output = os.path.join(save_path, du_output) + du = du.astype(np.float32) + du_fits = fits.HDUList([fits.PrimaryHDU(data=du)]) + du_fits.writeto(du_output, overwrite=True) + return du + + +def main(input_path : str, save_path : str, number_list : list): + for number in number_list: + print(number, end=' ') + bias_path = glob(os.path.join(input_path, "MSC*/*CLB*_" + number + '_*')) + bias = save_fits(load_bias, "median", bias_path, save_path) + print('bias finish', end=' ') + dark_path = glob(os.path.join(input_path, "MSC*/*CLD*_" + number + '_*')) + dark = save_fits(load_dark, "median", dark_path, save_path, bias) + print('dark finish', end=' ') + flat_path = glob(os.path.join(input_path, "MSC*/*CLF*_" + number + '_*')) + flat = save_fits(load_flat, "median", flat_path, save_path, bias, dark) + print('flat finish') + + +if __name__ == "__main__": + input_path = "/data/test20211012/ref/" + save_path = "/data/test20211012/ref/" + number_list = ['06', '07', '08', '09', '11', '12', '13', '14', + '15', '16', '17', '18', '19', '20', '22', '23', '24', '25'] + main(input_path, save_path, number_list) + diff --git a/CSST_2021-12-30_CCD23_epoch20.pth b/CSST_2021-12-30_CCD23_epoch20.pth new file mode 100644 index 0000000000000000000000000000000000000000..c1757cfe8977bea0a7f7ad0f9499c0031501e685 Binary files /dev/null and b/CSST_2021-12-30_CCD23_epoch20.pth differ diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000000000000000000000000000000000000..62b47a27e9d7661fafda92df9311d566dc88ca42 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,21 @@ +FROM hub.cstcloud.cn/csst/python:3.8 + +RUN pip install --no-deps ccdproc && pip install deepCR + +RUN apt-get update && apt-get install -y vim + +RUN mkdir -p /home/csstpipeline/code/csst-msc-instrument-effect/csst_msc_instrument_effect +# COPY run.sh /app/bin/ +COPY instrument_effect_wrapper.py setup.py requirements.txt \ + /home/csstpipeline/code/csst-msc-instrument-effect/ + +COPY csst_msc_instrument_effect/* \ + /home/csstpipeline/code/csst-msc-instrument-effect/csst_msc_instrument_effect/ + +COPY MSC_crmask.ini /home/csstpipeline/csst-msc-instrument-effect-master/ + +RUN cd /home/csstpipeline/code/csst-msc-instrument-effect/ \ + && python setup.py install + +COPY cleanup.sh /app/bin/ +COPY csst_msc_iec.sh /app/bin/run.sh diff --git a/MSC_crmask.ini b/MSC_crmask.ini new file mode 100644 index 0000000000000000000000000000000000000000..5cf881d7a36787bce2340cc5760c0d6b9ec789c1 --- /dev/null +++ b/MSC_crmask.ini @@ -0,0 +1,75 @@ +[global] +;which method used for cosmic ray detection +;model=deepCR + +update_flag=True +;whether save CR cleaned image +save_flag=True +;flag image suffix +flag_suffix=flg +;CR masked image suffix +data_suffix=crclean +;whether append CR masked image to the output FITS file +append_flag=False +;mask image suffix +mask_suffix=crmask +;whether use GPU +gpu_flag=False +;whether fill appropriate value for CR contaminated pixels +fill_flag=True +;filling method +fill_method=meanmask +;torch thread +torch_thread=1 + +[deepCR] +threshold=0.5 +inpaint=True +binary=True +patch=256 +segment=True +parallel=True +;Cosmic ray clean training model +clean_model=/home/csstpipeline/csst-msc-instrument-effect-main/CSST_2021-12-30_CCD23_epoch20.pth +;Inpaint model +inpaint_model=ACS-WFC-F606W-2-32 +n_jobs=4 +;train parameters +;ignore pixel, e.g., bad pixel, saturation, npy or fits +;ignore=ignore.fits +;sky=sky.fits +;aug_sky_min=0 +;aug_sky_max=0 +;name=xxxx +hidden=50 +gpu=False +epoch=50 +batch=16 +lr=0.005 +auto_lr_decay=True +lr_decay_patience=4 +lr_decay_factor=0.1 +save_after=1 +plot_every=10 +use_tqdm=False +use_tqdm_notebook=False +plot_roc=False +directory=./data/ + +[lacosmic] +sigclip=4.5 +sigfrac=0.3 +objlim=5.0 +gain=1.0 +readnoise=6.5 +satlevel=65535.0 +pssl=0.0 +niter=4 +sepmed=True +cleantype=meanmask +fsmode=median +psfmodel=gauss +psffwhm=2.5 +psfsize=7 +psfbeta=4.765 +gain_apply=True diff --git a/Makefile b/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..2df9082e8845a86799271df3915879ffbe37e6d1 --- /dev/null +++ b/Makefile @@ -0,0 +1,13 @@ +IMAGE=hub.cstcloud.cn/csst/detector-effect-correction + +build: + docker build --network=host -t $(IMAGE) . + +push: + docker push $(IMAGE) + +clean: + docker rmi $(IMAGE) + +rmi: + for H in h0 h1 h2 ; do echo $$H ; ssh $$H docker rmi $(IMAGE) ; done diff --git a/README.md b/README.md new file mode 100644 index 0000000000000000000000000000000000000000..43da1a42ea6f1c905279a71dec64a7806f53db3b --- /dev/null +++ b/README.md @@ -0,0 +1,34 @@ +# csst_msc_instrument_effect 仪器效应改正 + +用于改正CSST的0级数据中的仪器效应问题，如本底、暗流、平场，并标记出坏像元、热像元、暖像元、溢出像元、宇宙线。也会运算出权重weight。 + +## 在命令行运行 + +```bash +sh csst_msc_iec.sh [id] +``` + +## 使用python运行 +```python +from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection + + +iec = InstrumentEffectCorrection( + data_path=data_filename, # 科学图fits文件路径 + bias_path=bias_filename, # 本底fits文件路径(可以是路径列表) + dark_path=dark_filename, # 暗场fits文件路径(可以是路径列表) + flat_path=flat_filename, # 平场fits文件路径(可以是路径列表) + cray_path=cray_filename, # 模拟宇宙线文件路径(因为模拟数据的问题, 暂时添加, 之后会去除) + output_path=output_path, # 输出路径 + config_path=config_path, # crmask配置文件路径, 指的就是MSC_crmask.ini +) +iec.run() +``` +生成的文件会有img, flag, weight, header四个 + +## 使用BKGscript.py + +运行BKGscript.py 更改main函数内input_path与save_path与number_list三个变量 +input_path为参考源文件路径 +save_path为输出路径 +number_list为需要进行合并的参考文件编号 diff --git a/cleanup.sh b/cleanup.sh new file mode 100644 index 0000000000000000000000000000000000000000..251ed0f97e0685e0deff49c6832bbc42c11e7459 --- /dev/null +++ b/cleanup.sh @@ -0,0 +1,4 @@ +#!/bin/bash + +# 当前0级数据处理完成，告知下阶段继续处理 +touch /next-stage/$1 diff --git a/csst_msc_iec.sh b/csst_msc_iec.sh new file mode 100644 index 0000000000000000000000000000000000000000..20612ea0dfce265c80ac90173874769b51eb027a --- /dev/null +++ b/csst_msc_iec.sh @@ -0,0 +1,4 @@ +#!/bin/bash + +cd / +python /home/csstpipeline/csst-msc-instrument-effect-main/instrument_effect_wrapper.py $1 diff --git a/csst_msc_instrument_effect/__init__.py b/csst_msc_instrument_effect/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..5bd6efae077658e3d69005e31cbb3a23f0814b01 --- /dev/null +++ b/csst_msc_instrument_effect/__init__.py @@ -0,0 +1 @@ +# init file \ No newline at end of file diff --git a/csst_msc_instrument_effect/msc_crmask.py b/csst_msc_instrument_effect/msc_crmask.py new file mode 100644 index 0000000000000000000000000000000000000000..478076e3826298c4a263d389dbe13e64a4ee82f2 --- /dev/null +++ b/csst_msc_instrument_effect/msc_crmask.py @@ -0,0 +1,660 @@ +#!/usr/bin/env python3 + +############################################################## +#Author Hu Yi (NAOC, huyi.naoc@gmail.com), Zhang Yajie (TJU)# +############################################################## + +#Thanks to Zhang Keming who is the author of deepCR. +#Thanks to the authors of astropy and ccdproc. +#Software package dependencies, numpy, astropy, ccdproc and deepCR +#Installation dependencies on Ubuntu 20.04 +#apt install python3-numpy python3-astropy python3-ccdproc python3-pip +#python3 -m pip install deepCR +#Version 0.1 + +import configparser +import numpy as np +from astropy.io import fits as pyfits +from astropy import units as u +from astropy.nddata import CCDData +from pathlib import Path +import ccdproc + +from deepCR import deepCR +from deepCR import train +import torch + +import datetime +import sys + + +__all__ = ['CRMask'] + +class CRMask: + def __init__(self, obj, flag = None, mask = None, sky = None, save_flag = True, update_flag = True, save_name = None, flag_suffix = 'flg', clean_suffix = 'crclean', append_flag = False, mask_suffix = 'crmask', model = 'deepCR', fill_flag = True, fill_method = 'inpainting', gpu_flag = False, config_path = 'MSC_crmask.ini', **kwargs): + """ + Instantiation of CRMask with specified model configuration. + Parameters + ---------- + obj : string, Path, astropy.io.fits.HDUList, numpy.ndarray, astropy.nddata.CCDData or list of string + if model is ``deepCR``, ``lacosmic``, obj is input image to be cosmic-ray masked + if model is ``deepCR_train``, obj is input training and validating images + flag : (optional) string, Path, astropy.io.fits.HDUList, numpy.ndarray or astropy.nddata.CCDData + flag image(s), default is None + mask : (optional) string or list of string + mask image(s), default is None, necessary when model is ``deepCR_train``, otherwise, let it alone + sky : (optional) string or list of string + sky image(s), default is None, optional when model is ``deepCR_train``, otherwise, let it alone + save_flag : (optional) boolean + whether save CR mask (and cleaned) image, default is True + update_flag : (optional) boolean + whether update flag image, default is True + save_name : (optional) string + output mask, cleaned and flag filename. default is None. If save_name is None, use the filename as the input. And if save_flag is True, save_name is None, and obj is a numpy.ndarray or a astropy.nddata.CCDData, cr_mask will raise a ValueError exception + flag_suffix : (optional) string + suffix name of flag file, if flag is a numpy.ndarray or a astropy.nddata.CCDData, default is ``flg`` + mask_suffix : (optional) string + suffix name of mask file, default is ``crmask`` + clean_suffix : (optional) string + suffix name of cleaned file, default is ``crclean`` + model : (optional) string + model type, can be ``deepCR``, ``lacosmic``, ``deepCR_train``, default is ``deepCR`` + fill_flag : (optional) boolean + whehter generate cleaned image, default is True + fill_method : (optional) string + fill method for CR contaminated pixel, can be ``inpainting``, ``meanmask``, ``meanmed``, default is ``inpainting`` + gpu_flag : (optional) boolean + whether use GPU, default is False + config_path : (optional) string + configuration file path, default is ``../conf/MSC_crmask.ini`` + """ + if model == 'deepCR_train': + self.image_sets = obj + self.mask_sets = mask + self.ignore_sets = flag + self.sky_sets = sky + self.gpu_flag = gpu_flag + if config_path != None: + config = configparser.ConfigParser() + config.read(config_path) + if config.has_option('global', 'gpu_flag'): + self.append_flag = config.getboolean('global', 'gpu_flag') + self.config = config + return + + if isinstance(obj, str) or isinstance(obj, Path): + self.hdulist = pyfits.open(obj) + elif isinstance(obj, pyfits.HDUList): + self.hdulist = obj + elif isinstance(obj, np.ndarray) or isinstance(obj, CCDData): + hdulist = pyfits.HDUList() + primary = pyfits.PrimaryHDU() + header = primary.header + header['NEXTEND'] = 1 + hdulist.append(primary) + image = pyfits.ImageHDU() + header = image.header + header['BITPIX'] = -32 + header['NAXIS'] = 2 + header.insert('NAXIS', ('NAXIS1', obj.shape[0]), after = True) + header.insert('NAXIS1', ('NAXIS2', obj.shape[1]), after = True) + header['EXTNAME'] = 'SCI' + header['EXTVER'] = 1 + if isinstance(obj, CCDData): + image.data = np.asarray(CCDData) + else: + image.data = obj + hdulist.append(image) + self.hdulist = hdulist + else: + raise TypeError('For cosmic ray masking and cleaning, obj must be a string, or a pathlib.Path object, or a astropy.io.fits.HDUList object!') + if flag != None: + if isinstance(flag, str) or isinstance(flag, Path): + flag_hdulist = pyfits.open(flag, mode = 'update') + elif isinstance(flag, pyfits.HDUList): + flag_hdulist = flag + elif isinstance(flag, np.ndarray) or isinstance(flag, CCDData): + flag_hdulist = pyfits.HDUList() + flag_primary = pyfits.PrimaryHDU(header = self.hdulist[0].header) + flag_hdulist.append(flag_hdulist) + if isinstance(obj, CCDData): + flag_data = np.asarray(CCDData) + flag_image = pyfits.ImageHDU(data=flag_data, header=self.hdulist[1].header) + else: + flag_image = pyfits.ImageHDU(data=flag, header=self.hdulist[1].header) + flag_hdulist.append(flag_image) + else: + raise TypeError('For cosmic ray masking and cleaning, mask must be a string, or a pathlib.Path object, or a astropy.io.fits.HDUList object!') + self.flag = flag_hdulist + else: + self.flag = None + + self.update_flag = update_flag + self.save_flag = save_flag + self.save_name = save_name + self.flag_suffix = flag_suffix + self.clean_suffix = clean_suffix + self.append_flag = append_flag + self.mask_suffix = mask_suffix + self.gpu_flag = gpu_flag + self.fill_method = fill_method + self.model = model + + if config_path != None: + config = configparser.ConfigParser() + config.read(config_path) + + if config.has_option('global', 'save_flag'): + self.save_flag = config.getboolean('global', 'save_flag') + + if config.has_option('global', 'clean_suffix'): + self.clean_suffix = config.get('global', 'clean_suffix') + + if config.has_option('global', 'append_flag'): + self.append_flag = config.getboolean('global', 'append_flag') + + if config.has_option('global', 'flag_suffix'): + self.flag_suffix = config.get('global', 'flag_suffix') + + if config.has_option('global', 'mask_suffix'): + self.mask_suffix = config.get('global', 'mask_suffix') + + if config.has_option('global', 'gpu_flag'): + self.append_flag = config.getboolean('global', 'gpu_flag') + + if config.has_option('global', 'fill_flag'): + self.fill_flag = config.getboolean('global', 'fill_flag') + + if config.has_option('global', 'fill_method'): + self.fill_method = config.get('global', 'fill_method') + + if config.has_option('global', 'model'): + self.model = config.get('global', 'model') + + if config.has_option('global', 'torch_thread'): + self.torch_thread = config.get('global', 'torch_thread') + torch.set_num_threads(self.torch_thread) + self.config = config + + def cr_mask_lacosmic(self): + + config = self.config + + if config.has_option('lacosmic', 'sigclip'): + sigclip = config.getfloat('lacosmic', 'sigclip') + else: + sigclip = 4.0 + + if config.has_option('lacosmic', 'sigfrac'): + sigfrac = config.getfloat('lacosmic', 'sigfrac') + else: + sigfrac = 0.3 + + if config.has_option('lacosmic', 'objlim'): + objlim = config.getfloat('lacosmic', 'objlim') + else: + objlim = 5.0 + + if config.has_option('lacosmic', 'gain'): + gain = config.getfloat('lacosmic', 'gain') + else: + gain = 1.0 + + if config.has_option('lacosmic', 'readnoise'): + readnoise = config.getfloat('lacosmic', 'readnoise') + else: + readnoise = 6.5 + + if config.has_option('lacosmic', 'satlevel'): + satlevel = config.getfloat('lacosmic', 'satlevel') + else: + satlevel = 65535.0 + + if config.has_option('lacosmic', 'pssl'): + pssl = config.getfloat('lacosmic', 'pssl') + else: + pssl = 0.0 + + if config.has_option('lacosmic', 'niter'): + niter = config.getint('lacosmic', 'niter') + else: + niter = 4 + + if config.has_option('lacosmic', 'sepmed'): + sepmed = config.getboolean('lacosmic', 'sepmed') + else: + sepmed = True + + + if config.has_option('lacosmic', 'cleantype'): + cleantype = config.get('lacosmic', 'cleantype') + else: + cleantype = self.fill_method + + if config.has_option('lacosmic', 'fsmode'): + fsmode = config.get('lacosmic', 'fsmode') + else: + fsmode = 'median' + + if config.has_option('lacosmic', 'psfmodel'): + psfmodel = config.get('lacosmic', 'psfmodel') + else: + psfmodel = 'gauss' + + if config.has_option('lacosmic', 'psffwhm'): + psffwhm = config.getfloat('lacosmic', 'psffwhm') + else: + psffwhm = 2.5 + + if config.has_option('lacosmic', 'psfsize'): + psfsize = config.getint('lacosmic', 'psfsize') + else: + psfsize = 7 + + if config.has_option('lacosmic', 'psfk'): + psfk = config.get('lacosmic', 'psfk') + else: + psfk = None + + if config.has_option('lacosmic', 'psfbeta'): + psfbeta = config.getfloat('lacosmic', 'psfbeta') + else: + psfbeta = 4.765 + + if config.has_option('lacosmic', 'gain_apply'): + gain_apply = config.getboolean('lacosmic', 'gain_apply') + else: + gain_apply = True + + data = self.hdulist[1].data + + start_time = datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S.%f") + cleaned, masked = ccdproc.cosmicray_lacosmic(data, sigclip, sigfrac, objlim, gain, readnoise, satlevel, pssl, niter, sepmed, cleantype, fsmode, psfmodel, psffwhm, psfsize, psfk, psfbeta, False, gain_apply) + end_time = datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S.%f") + masked_hdulist = pyfits.HDUList() + masked_primary = pyfits.PrimaryHDU(header=self.hdulist[0].header) + masked_hdulist.append(masked_primary) + masked = masked.astype(np.uint8) + masked_image = pyfits.ImageHDU(data=masked, header=self.hdulist[1].header) + masked_hdulist.append(masked_image) + + #TODO + #add history keywords here + mask = None + if self.flag != None: + self.flag[1].data |= (masked<<4) + self.flag[1].header.add_history('Use ccdproc.cosimicray_lacosmicfor cosmic ray detecting') + value = 'CRMask start at {0}'.format(start_time) + self.flag[1].header.add_history(value) + value = 'CRMask end at {0}'.format(end_time) + self.flag[1].header.add_history(value) + + masked_hdulist[1].header.add_history('Use ccdproc.cosimicray_lacosmic for cosmic ray detecting') + value = 'CRMask start at {0}'.format(start_time) + masked_hdulist[1].header.add_history(value) + value = 'CRMask end at {0}'.format(end_time) + masked_hdulist[1].header.add_history(value) + + if self.fill_flag: + cleaned_hdulist = pyfits.HDUList() + cleaned_primary = pyfits.PrimaryHDU(header=self.hdulist[0].header) + cleaned_hdulist.append(cleaned_primary) + cleaned_image = pyfits.ImageHDU(data=cleaned, header=self.hdulist[1].header) + cleaned_hdulist.append(cleaned_image) + cleaned_hdulist[1].header.add_history('Use ccdproc.cosimicray_lacosmic for cosmic ray detecting') + value = 'CRMask start at {0}'.format(start_time) + cleaned_hdulist[1].header.add_history(value) + value = 'CRMask end at {0}'.format(end_time) + cleaned_hdulist[1].header.add_history(value) + + return masked_hdulist, cleaned_hdulist + else: + return masked_hdulist + + def cr_mask_deepCR(self): + + config = self.config + + if config.has_option('deepCR', 'threshold'): + threshold = config.getfloat('deepCR', 'threshold') + else: + threshold = 0.5 + + if config.has_option('deepCR', 'inpaint'): + inpaint = config.getboolean('deepCR', 'inpaint') + else: + inpaint = True + + if config.has_option('deepCR', 'binary'): + binary = config.getboolean('deepCR', 'binary') + else: + binary = True + + if config.has_option('deepCR', 'patch'): + patch = config.getint('deepCR', 'patch') + else: + patch = 256 + + if config.has_option('deepCR', 'segment'): + segment = config.getboolean('deepCR', 'segment') + else: + segment = True + + if config.has_option('deepCR', 'parallel'): + parallel = config.getboolean('deepCR', 'parallel') + else: + parallel = True + + if config.has_option('deepCR', 'clean_model'): + clean_model = config.get('deepCR', 'clean_model') + else: + clean_model = 'ACS-WFC-F606W-2-32' + + if config.has_option('deepCR', 'inpaint_model'): + inpaint_model = config.get('deepCR', 'inpaint_model') + else: + inpaint_model = 'ACS-WFC-F606W-2-32' + + if config.has_option('deepCR', 'n_jobs'): + n_jobs = config.getint('deepCR', 'n_jobs') + else: + n_jobs = -1 + + if self.gpu_flag: + model = deepCR(clean_model, inpaint_model, device = 'GPU', hidden=50) + else: + model = deepCR(clean_model, inpaint_model, device = 'CPU', hidden=50) + + data = self.hdulist[1].data + + masked_hdulist = pyfits.HDUList() + masked_primary = pyfits.PrimaryHDU(header=self.hdulist[0].header) + masked_hdulist.append(masked_primary) + if inpaint: + cleaned_hdulist = pyfits.HDUList() + cleaned_primary = pyfits.PrimaryHDU(header=self.hdulist[0].header) + cleaned_hdulist.append(cleaned_primary) + + start_time = datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S.%f") + if binary: + if inpaint: + masked, cleaned = model.clean(data, threshold = threshold, inpaint = inpaint, segment = segment, patch = patch, parallel = parallel, n_jobs = n_jobs) + else: + masked = model.clean(data, threshold = threshold, inpaint = inpaint, segment = segment, patch = patch, parallel = parallel, n_jobs = n_jobs) + else: + if inpaint: + masked, cleaned = model.clean(data, threshold = threshold, inpaint = inpaint, binary = False, segment = segment, patch = patch, parallel = parallel, n_jobs = n_jobs) + else: + masked = model.clean(data, threshold = threshold, inpaint = inpaint, binary = False, segment = segment, patch = patch, parallel = parallel, n_jobs = n_jobs) + end_time = datetime.datetime.now().strftime("%Y-%m-%dT%H:%M:%S.%f") + + masked = masked.astype(np.uint8) + masked_image = pyfits.ImageHDU(data=masked, header=self.hdulist[1].header) + masked_hdulist.append(masked_image) + + #TODO + #add history keywords here + if self.flag != None: + self.flag[1].data |= (masked<<4) + self.flag[1].header.add_history('Use deepCR for cosmic ray detecting') + value = 'CRMask start at {0}'.format(start_time) + self.flag[1].header.add_history(value) + value = 'CRMask end at {0}'.format(end_time) + self.flag[1].header.add_history(value) + masked_hdulist[1].header.add_history('Use deepCR for cosmic ray detecting') + value = 'CRMask start at {0}'.format(start_time) + masked_hdulist[1].header.add_history(value) + value = 'CRMask end at {0}'.format(end_time) + masked_hdulist[1].header.add_history(value) + if inpaint: + cleaned_image = pyfits.ImageHDU(data=cleaned, header=self.hdulist[1].header) + cleaned_hdulist.append(cleaned_image) + cleaned_hdulist[1].header.add_history('Use deepCR for cosmic ray detecting') + value = 'CRMask start at {0}'.format(start_time) + cleaned_hdulist[1].header.add_history(value) + value = 'CRMask end at {0}'.format(end_time) + cleaned_hdulist[1].header.add_history(value) + return masked_hdulist,cleaned_hdulist + else: + return masked_hdulist + + #return a hdulist of a masked image, and a hdulist of a cleaned image if fill_flag is set. + def cr_mask(self): + #here do cr mask task + if self.model == 'lacosmic': + if self.fill_flag: + masked, cleaned = CRMask.cr_mask_lacosmic(self) + else: + masked = CRMask.cr_mask_lacosmic(self) + elif self.model == 'deepCR': + if self.fill_flag: + masked, cleaned = CRMask.cr_mask_deepCR(self) + else: + masked = CRMask.cr_mask_deepCR(self) + + else: + raise ValueError('Cosmic ray model are not supported!') + + #save the result to a file. + output = None + if self.save_flag: + if self.save_name != None: + output = self.save_name + else: + output = self.hdulist.filename() + #append to the input + if self.append_flag: + self.hdulist.append(masked[1]) + if self.fill_flag: + self.hdulist.append(cleaned[1]) + if output != None: + self.hdulist.writeto(output) + else: + #TODO + #raise an exception. + pass + else: + if output != None: + if output.split('.')[-1] == 'fits': + mask_output = output[0:-4] + self.mask_suffix + '.fits' + clean_output = output[0:-4] + self.clean_suffix + '.fits' + else: + mask_output = output + '.' + self.mask_suffix + '.fits' + clean_output = output + '.' + self.clean_suffix + '.fits' + masked.writeto(mask_output) + if self.fill_flag: + cleaned.writeto(clean_output) + else: + #TODO + #raise an exception here + pass + + #update flag file. + if self.update_flag and self.flag != None: + flag_output = self.flag.filename() + if flag_output == None: + if output != None: + if output.split('.')[-1] == 'fits': + flag_output = output[0:-4] + self.flag_suffix + '.fits' + else: + flag_output = output + '.' + self.flag_suffix + '.fits' + else: + #raise an exception here. + pass + self.flag.writeto(flag_output) + + if self.fill_flag: + return (masked, cleaned) + else: + return masked + + def cr_train_deepCR_image_to_ndarray(self, image_sets, patch): + + if isinstance(image_sets, str): + if image_sets[:-4] == '.npy': + input_image = np.load(image_sets) + elif image_sets[:-4] == '.fits': + image = pyfits.open(image_sets) + data = image[1].data + x, y = data.shape + m = x // patch + n = y // patch + start_x = np.mod(x, patch) // 2 + start_y = np.mod(y, patch) // 2 + input_image = data[start_x:start_x + m * patch][ start_y:start_y + m * patch].reshape[m*n][patch][patch] + else: + #TODO + #raise an exception + pass + if isinstance(image_sets, list): + input_list = [] + for image_file in image_sets: + if isinstance(image_sets, str): + if image_file[:-4] == '.npy': + input_image = np.load(image_file) + elif image_file[:-4] == '.fits': + image = pyfits.open(image_file) + data = image[1].data + x, y = data.shape + m = x // patch + n = y // patch + start_x = np.mod(x, patch) // 2 + start_y = np.mod(y, patch) // 2 + input_image = data[start_x:start_x + m * patch][start_y:start_y + m * patch].reshape[m*n][patch][patch] + else: + #TODO + #raise an exception + pass + input_list.append(input_image) + else: + #TODO + #raise an exception + pass + input_image = np.concatenate(input_list) + return input_image + + def cr_train_deepCR_prepare_data(self, patch): + if self.image_sets != None: + self.training_image = CRMask.cr_train_deepCR_image_to_ndarray(self, self.image_sets, patch) + else: + raise ValueError('training image should not be None') + + if self.mask_sets != None: + self.training_mask = CRMask.cr_train_deepCR_image_to_ndarray(self, self.mask_sets, patch) + else: + raise ValueError('mask image should not be None') + + if self.flag != None: + self.training_ignore = CRMask.cr_train_deepCR_image_to_ndarray(self, self.flag, patch) + else: + self.training_ignore = None + + if self.sky != None: + self.training_sky = CRMask.cr_train_deepCR_image_to_ndarray(self, self.sky, patch) + else: + self.training_sky = None + + + def cr_train_deepCR(self): + config = self.config + + if config.has_option('deepCR', 'aug_sky'): + aug_sky = np.array(config.get('deepCR', 'aug_sky').split()).astype(np.float32) + aug_sky0 = aug_sky[0] + aug_sky1 = aug_sky[1] + else: + aug_sky0 = 0. + aug_sky1 = 0. + + if config.has_option('deepCR', 'model_name'): + model_name = config.get('deepCR', 'model_name') + else: + model_name = 'model' + + if config.has_option('deepCR', 'hidden'): + hidden = config.getint('deepCR', 'hidden') + else: + hidden = 50 + + if config.has_option('deepCR', 'epoch'): + epoch = config.getint('deepCR', 'epoch') + else: + epoch = 50 + + if config.has_option('deepCR', 'patch'): + patch = config.getint('deepCR', 'patch') + else: + patch = 256 + + if config.has_option('deepCR', 'batch'): + batch = config.getint('deepCR', 'batch') + else: + batch = 16 + + if config.has_option('deepCR', 'lr'): + lr = config.getfloat('deepCR', 'lr') + else: + lr = 0.005 + + if config.has_option('deepCR', 'auto_lr_decay'): + auto_lr_decay = config.getboolean('deepCR', 'auto_lr_decay') + else: + auto_lr_decay = True + + if config.has_option('deepCR', 'lr_decay_patience'): + lr_decay_patience = config.getint('deepCR', 'lr_decay_patience') + else: + lr_decay_patience = 4 + + if config.has_option('deepCR', 'lr_decay_factor'): + lr_decay_factor = config.getfloat('deepCR', 'lr_decay_factor') + else: + lr_decay_factor = 0.1 + + if config.has_option('deepCR', 'save_after'): + save_after = config.getfloat('deepCR', 'save_after') + else: + save_after = 100000.0 + + if config.has_option('deepCR', 'plot_every'): + plot_every = config.getint('deepCR', 'plot_every') + else: + plot_every = 4 + + if config.has_option('deepCR', 'verbose'): + verbose = config.getboolean('deepCR', 'verbose') + else: + verbose = True + + if config.has_option('deepCR', 'use_tqdm'): + use_tqdm = config.getboolean('deepCR', 'use_tqdm') + else: + use_tqdm = False + + if config.has_option('deepCR', 'use_tqdm_notebook'): + use_tqdm_notebook = config.getboolean('deepCR', 'use_tqdm_notebook') + else: + use_tqdm_notebook = False + + if config.has_option('deepCR', 'directory'): + directory = config.get('deepCR', 'directory') + else: + directory = './' + + CRMask.cr_train_deepCR_prepare_data(self, patch) + trainer = train(self.training_image, self.training_mask, self.training_ignore, self.training_sky, [aug_sky0, aug_sky1], model_name, hidden, self.gpu_flag, epoch, batch, lr, auto_lr_decay, lr_decay_patience, lr_decay_patience, save_after, plot_every, verbose, use_tqdm, use_tqdm_notebook, directory) + trainer.train() + + def cr_train(self): + if self.model == 'deepCR_train': + CRMask.cr_train_deepCR(self) + else: + raise ValueError('Unsupported training model') + +#to run the code +#./MSC_crmask xxxx.fits deepCR for deepCR +#./MSC_crmask xxxx.fits lacosmic for lacosmic + +if __name__ == '__main__': + crobj = CRMask(sys.argv[1], model = sys.argv[2]) + crobj.cr_mask() diff --git a/csst_msc_instrument_effect/msc_iec.py b/csst_msc_instrument_effect/msc_iec.py new file mode 100644 index 0000000000000000000000000000000000000000..624c819b62854f8e14e771f0d068f6223f49e61d --- /dev/null +++ b/csst_msc_instrument_effect/msc_iec.py @@ -0,0 +1,276 @@ +import os +from typing import List +from astropy.io.fits import header +from datetime import datetime + +import numpy as np +from astropy.io import fits +from deepCR import deepCR +import torch + +from csst_msc_instrument_effect.msc_crmask import CRMask + + +class InstrumentEffectCorrection: + def __init__( + self, + data_input, + bias_input, + dark_input, + flat_input, + output_path, + config_path, + cray_input=None, + ) -> None: + self.data_input = data_input + self.bias_input = bias_input + self.dark_input = dark_input + self.flat_input = flat_input + self.cray_input = cray_input + self.output = output_path + self.config_path = config_path + # RDNOISE + self.RDNOISE = "RDNOISE1" + # GAIN + self.GAIN = "GAIN1" + # EXPTIME + self.EXPTIME = "EXPTIME" + # output + if not os.path.exists(self.output): + os.mkdir(self.output) + + def set_data(self): + # data + with fits.open(self.data_input) as hdul: + self.primary_header = hdul[0].header + self.data_raw = hdul[1].data.astype(int) + self.image_header = hdul[1].header + + def set_bias(self, mode): + # bias + def func(path): + du = fits.getdata(path) + return du.astype(int) + + if isinstance(self.bias_input, list) and len(self.bias_input) != 1: + bias = np.array([*map(func, self.bias_input)]) + self.bias, var = array_combine(bias, mode) + elif isinstance(self.bias_input, list) and len(self.bias_input) == 1: + self.bias = func(self.bias_input[0]) + elif isinstance(self.bias_input, np.ndarray) and len(self.bias_input.shape) == 2: + self.bias = self.bias_input + else: + self.bias = func(self.bias_input) + + def set_dark(self, mode): + # dark + def func(path): + with fits.open(path) as hdul: + du = hdul[1].data + hu = hdul[0].header + du = du.astype(int) + du = du - self.bias + du = du - self.cray / self.image_header[self.GAIN] # tmp + return du / hu[self.EXPTIME] + + if isinstance(self.dark_input, list) and len(self.dark_input) != 1: + dark = np.array([*map(func, self.dark_input)]) + self.dark, var = array_combine(dark, mode) + elif isinstance(self.dark_input, list) and len(self.dark_input) == 1: + self.dark = func(self.dark_input[0]) + elif isinstance(self.dark_input, np.ndarray) and len(self.dark_input.shape) == 2: + self.dark = self.dark_input + else: + self.dark = func(self.dark_input) + + def set_flat(self, mode): + # flat + def func(path): + with fits.open(path) as hdul: + du = hdul[1].data + hu = hdul[0].header + du = du.astype(int) + du = du - self.bias + du = du / hu[self.EXPTIME] + du = du - self.dark + return du / np.median(du) + + if isinstance(self.flat_input, list) and len(self.flat_input) != 1: + flat = np.array([*map(func, self.flat_input)]) + self.flat, var = array_combine(flat, mode) + elif isinstance(self.flat_input, list) and len(self.flat_input) == 1: + self.flat = func(self.flat_input[0]) + elif isinstance(self.flat_input, np.ndarray) and len(self.flat_input.shape) == 2: + self.flat = self.flat_input + else: + self.flat = func(self.flat_input) + + def set_cray(self): + # cray + if not self.cray_input == None: + self.cray = fits.getdata(self.cray_input).astype(int) + else: + self.cray = 0 + + def fix_data(self,): + self.data_fix0 = np.divide( + self.data_raw - self.bias - self.dark * self.primary_header[self.EXPTIME], + self.flat, + out=np.zeros_like(self.data_raw, float), + where=(self.flat != 0), + ) + + def set_flag(self,): + flag = np.zeros_like(self.data_raw, dtype=np.uint16) + # 00000001: 坏像元 + # 因为探测器本身原因造成不可用于有效科学研究的像元, 像元响应小于0.5中值 大于1.5中值 + med = np.median(self.flat) + flg = (self.flat < 0.5 * med) | (1.5 * med < self.flat) + flag = flag | (flg * 1) + # 00000010: 热像元 + # 因为探测器本身原因造成的影响科学研究结果的像元. 像元在150秒积分时间内, 暗流计数大于探测器平均读出噪声的平方. + dark = self.dark.copy() * self.primary_header[self.EXPTIME] + dark[dark < 0] = 0 + flg = 1 * self.image_header[self.RDNOISE] ** 2 <= dark # 不确定是否包含 暂定包含 + flag = flag | (flg * 2) + # 00000100: 暖像元 + # 因为探测器本身原因造成的影响科学研究结果的像元. 像元的150秒积分时间内, 暗流计数大于0.5被读出噪声的平方, 但小于1被读出噪声的平方. + flg = (0.5 * self.image_header[self.RDNOISE] ** 2 < dark) & ( + dark < 1 * self.image_header[self.RDNOISE] ** 2 + ) + flag = flag | (flg * 4) + # 00001000: 饱和溢出像元 + # 饱和像元及流量溢出污染的像元. + flg = self.data_raw == 65535 + flag = flag | (flg * 8) + # 00010000: 宇宙线像元 + # 宇宙线污染的像元 + crobj = CRMask(self.data_fix0, config_path=self.config_path) + flag_fits, data_fits = crobj.cr_mask() + flag_fits = flag_fits[1].data.astype(np.uint16) + flag = flag | (flag_fits * 16) + # 00100000: 卫星或者人造移动天体轨迹污染的像元. + pass + # 01000000: 鬼像污染的像元. + pass + # 10000000: 散射光污染的像元(包括dragon breath) + pass + self.flag = flag + # self.data_fix1 = data_fits[1].data + + def set_weight(self): + data = self.data_fix0.copy() + data[self.data_fix0 < 0] = 0 + weight_raw = 1. / (self.image_header[self.GAIN] * data + self.image_header[self.RDNOISE] ** 2) + bias_weight = np.std(self.bias) + weight = 1. / (1. / weight_raw + 1. / bias_weight) * (self.primary_header[self.EXPTIME] ** 2) + weight[self.flag > 0] = 0 + self.weight = weight + + def save(self): + data_filename = self.data_input + data_basename0 = os.path.basename(data_filename).replace("L0", "L1") + data_basename = os.path.basename(data_basename0).replace(".fits", "_img.fits") + data_output = os.path.join(self.output, data_basename) + data_header = self.image_header.copy() + data_header["EXTNAME"] = "img" + data_header[self.GAIN] = ( + data_header[self.GAIN] * self.primary_header[self.EXPTIME] + ) + data_header["BUNIT"] = "e/s" + data_fits = fits.HDUList( + [ + fits.PrimaryHDU(header=self.primary_header), + fits.ImageHDU( + data=self.data_fix0.astype(np.float32) + / self.primary_header[self.EXPTIME], + header=data_header, + ), + ] + ) + data_fits.writeto(data_output, overwrite=True) + self.data_output = data_output + + flag_output = data_output.replace("img", "flg") + flag_header = self.image_header.copy() + flag_header["EXTNAME"] = "flg" + flag_fits = fits.HDUList( + [ + fits.PrimaryHDU(header=self.primary_header), + fits.ImageHDU(data=self.flag.astype(np.uint16), header=flag_header), + ] + ) + flag_fits.writeto(flag_output, overwrite=True) + self.flag_output = flag_output + + weight_output = data_output.replace("img", "wht") + weight_header = self.image_header.copy() + weight_header["EXTNAME"] = "wht" + weight_fits = fits.HDUList( + [ + fits.PrimaryHDU(header=self.primary_header), + fits.ImageHDU( + data=self.weight.astype(np.float32), header=weight_header + ), + ] + ) + weight_fits.writeto(weight_output, overwrite=True) + self.weight_output = weight_output + + header_name = data_basename.replace(".fits", ".head") + header_output = os.path.join(self.output, header_name) + hu = self.image_header + hu["INSTRU_S"] = (0, "instrument effect status") + hu["INST_TOL"] = (1234, "instrument effect operation time") + hu["INSTRU_V"] = ("0.1", "instrument effect version") + hu["INSTRU_P"] = ("test.conf", "instrument effect config file") + header_fits = fits.HDUList([fits.PrimaryHDU(header=hu)]) + header_fits.writeto(header_output, overwrite=True) + self.header_output = header_output + + # bias_output = data_output.replace("img", "bias") + # bias_fits = fits.HDUList([fits.PrimaryHDU(data=self.bias)]) + # bias_fits.writeto(bias_output, overwrite=True) + # dark_output = data_output.replace("img", "dark") + # dark_fits = fits.HDUList([fits.PrimaryHDU(data=self.dark)]) + # dark_fits.writeto(dark_output, overwrite=True) + # flat_output = data_output.replace("img", "flat") + # flat_fits = fits.HDUList([fits.PrimaryHDU(data=self.flat)]) + # flat_fits.writeto(flat_output, overwrite=True) + + def run(self): + self.set_data() + self.set_cray() + self.set_bias(mode="median_clip") + self.set_dark(mode="mean_clip") + self.set_flat(mode="mean_clip") + self.fix_data() + self.set_flag() + self.set_weight() + self.save() + + +def array_combine(ndarray, mode="mean"): + """ Function to combine 3-D data array + + Parameters + ---------- + ndarray: array, input data cube (3D) + model: mean, median, sum, mean_clip, median_clip, default is mean + """ + var = np.std(ndarray, axis=0) + ndarray_sort = np.sort(ndarray, axis=0) ### sort the flux for each pixel + ndarray_cut = ndarray_sort[1:-1] ### cut the min and max flux for each pixel + + if mode == "median": + array = np.median(ndarray_sort, axis=0) + elif mode == "median_clip": + array = np.median(ndarray_cut, axis=0) + elif mode == "sum": + array = np.sum(ndarray_sort, axis=0) + elif mode == "mean": + array = np.mean(ndarray_sort, axis=0) + elif mode == "mean_clip": + array = np.mean(ndarray_cut, axis=0) + return array, var + diff --git a/demo.py b/demo.py new file mode 100644 index 0000000000000000000000000000000000000000..e303bab85175fd9600a6f445f244babacaa4a169 --- /dev/null +++ b/demo.py @@ -0,0 +1,14 @@ +from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection + + + +iec = InstrumentEffectCorrection( + data_input="../data/test/MSC_MS_210525160000_100000022_24_raw.fits",#"../data/test/MSC_MS_210525123000_100000001_07_raw.fits",#"/data/L0/MSC_MS_210525121500_100000001_08_raw.fits", + bias_input="/data/BKG/MSC_CLB_210525120000_100000000_24_raw.fits", + dark_input="/data/BKG/MSC_CLD_210525121000_100000000_24_raw.fits", + flat_input="/data/BKG/MSC_CLF_210525120500_100000000_24_raw.fits", + cray_input="/data/BKG/MSC_CRD_210525121000_100000000_24_raw.fits", + output_path="/data/test/", + config_path="MSC_crmask.ini", +) +iec.run() diff --git a/instrument_effect_wrapper.py b/instrument_effect_wrapper.py new file mode 100644 index 0000000000000000000000000000000000000000..e54668d75e8d69cea1b85855514e3b08e570870d --- /dev/null +++ b/instrument_effect_wrapper.py @@ -0,0 +1,60 @@ +import sys +from datetime import datetime + +from csst_dfs_api.facility.level0 import Level0DataApi +from csst_dfs_api.facility.level0prc import Level0PrcApi +from csst_dfs_api.facility.calmerge import CalMergeApi +from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection + +l0Api = Level0DataApi() +prcApi = Level0PrcApi() +calApi = CalMergeApi() + +# ----------------------------------------------------- +# get data +l0_id = sys.argv[1] + +rec = l0Api.get(level0_id=l0_id) +print(rec) +print("Find Level0 Data :" + rec.data.file_path) + +bias_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="bias") +dark_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="dark") +flat_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="flat") +cray_rec = calApi.get_latest_by_l0(level0_id=l0_id, ref_type="cray") # temp + +print("Find Level0 Bias :" + bias_rec.data.file_path) +print("Find Level0 Dark :" + dark_rec.data.file_path) +print("Find Level0 Flat :" + flat_rec.data.file_path) +# ------------------------------------------------------ +# to-do data processing here +iec = InstrumentEffectCorrection( + data_input=rec.data.file_path, + bias_input=bias_rec.data.file_path, + dark_input=dark_rec.data.file_path, + flat_input=flat_rec.data.file_path, + output_path="/home/csstpipeline/data/L05_test/", + config_path="/home/csstpipeline/csst-msc-instrument-effect-master/MSC_crmask.ini", + cray_input=cray_rec.data.file_path, # temp +) +iec.run() + +params_file_path = iec.config_path +result_file_path = ', '.join([ + iec.data_output, + iec.flag_output, + iec.weight_output, + iec.header_output, +]) +dt = datetime.now() +# ------------------------------------------------------ +# write result +rec = prcApi.write( + level0_id=l0_id, + pipeline_id="P1", + prc_module="MSC-IE", # csst-msc-instrument-effect + params_file_path=params_file_path, + prc_status=0, + prc_time=dt.strftime('%Y-%m-%d %H:%M:%S'), + result_file_path=result_file_path) +print('Write Result:', rec) diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000000000000000000000000000000000000..d947126f613d6f72b8df03c9e9819cd465c161be --- /dev/null +++ b/requirements.txt @@ -0,0 +1,3 @@ +astropy +ccdproc +deepCR \ No newline at end of file diff --git a/run.py b/run.py new file mode 100644 index 0000000000000000000000000000000000000000..5403f207667fb0023a92c5f8743e762211641064 --- /dev/null +++ b/run.py @@ -0,0 +1,62 @@ +import os +from glob import glob +import argparse +from astropy.io.fits import getdata, getheader +from time import time + +from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection + + +def run(data_path, output, ref, conf): + header = getheader(data_path) + number = header['DETECTOR'][-2:] + # number = os.path.basename(data_path).split('_')[4] # 原来版本 + number_list = ['06', '07', '08', '09', '11', '12', '13', '14', '15', + '16', '17', '18', '19', '20', '22', '23', '24', '25'] + if number not in number_list: + print(number, 'not in') + return + bias = getdata(glob(os.path.join(ref, "*CLB*_" + number + '_*'))[0]) + dark = getdata(glob(os.path.join(ref, "*CLD*_" + number + '_*'))[0]) + flat = getdata(glob(os.path.join(ref, "*CLF*_" + number + '_*'))[0]) + iec = InstrumentEffectCorrection( + data_input=data_path, + bias_input=bias, + dark_input=dark, + flat_input=flat, + output_path=output, + config_path=conf + ) + iec.run() + print(number, 'done') + + +def main(): + parser = argparse.ArgumentParser() + parser.add_argument( + '-i', '--input', help='fits file path or folder path', required=True) + parser.add_argument( + '-o', '--output', help='output path', default='/data/L05/') + parser.add_argument( + '-r', '--ref', help='ref path', default='/data/ref/') + parser.add_argument( + '-c', '--conf', help='config path', default="/home/csstpipeline/csst-msc-instrument-effect-main/MSC_crmask.ini") + args = parser.parse_args() + if os.path.isdir(args.input): + for path in glob(os.path.join(args.input.strip(), '*_MS_SCI_*.fits')): + #for path in glob(os.path.join(args.input.strip(), '*_MS_*.fits')): + run(path, args.output, args.ref, args.conf) + else: + print('run test') + begin = time() + run(args.input.strip(), args.output.strip(), args.ref.strip(), args.conf.strip()) + print(time() - begin) + + +if __name__ == '__main__': + # python run.py -i "/data/test20211012/30s/MSC_MS_210525220000_100000020_13_raw.fits" -o "/data/test20211012/output/150s/" + # python run.py -i "/data/test20211012/150s/" -o '/data/test20211012/output/150s/' + # python run.py -i "/data/cali_20211012/L0/150s/" -o '/data/cali_20211012/output/' + # python run.py -i '/data/test20211012/30s/' -o '/data/test20211012/output/30s/' + # python run.py -i "/data/sim_data/20220413/MSC_0000100/CSST_MSC_MS_CRS_20270810081950_20270810082220_100000100_13_L0_1.fits" -o "/home/csstpipeline/test/" + main() diff --git a/run.sh b/run.sh new file mode 100644 index 0000000000000000000000000000000000000000..2a18d686aeed41c6829c97308b25ed1ab48a2e33 --- /dev/null +++ b/run.sh @@ -0,0 +1,13 @@ +# 运行: sh /home/csstpipeline/csst-msc-instrument-effect-main/run.sh [输入文件路径] [输出文件路径] [背景文件路径] [配置文件路径] +# 结果默认保存在/data/cali_20211012/output/ +# 背景文件路径默认在/data/test20211012/ref/ +# 配置文件路径默认在/home/csstpipeline/csst-msc-instrument-effect-main/MSC_crmask.ini +input=${1} +output=${2:-'/data/cali_20211012/output/'} +ref=${3:-'/data/test20211012/ref/'} +conf=${4:-"/home/csstpipeline/csst-msc-instrument-effect-main/MSC_crmask.ini"} +echo ${input} +echo ${output} +echo ${ref} +echo ${conf} +python /home/csstpipeline/csst-msc-instrument-effect-main/run.py -i ${input} -o ${output} -r ${ref} -c ${conf} diff --git a/run_18.py b/run_18.py new file mode 100644 index 0000000000000000000000000000000000000000..41fd4be99ab5389364e86b97f8889ac35ef503a2 --- /dev/null +++ b/run_18.py @@ -0,0 +1,41 @@ +import os +os.environ['KMP_DUPLICATE_LIB_OK'] = 'TRUE' +from glob import glob +from astropy.io.fits import getdata, getheader +from concurrent.futures import ThreadPoolExecutor +from time import time + +from csst_msc_instrument_effect.msc_iec import InstrumentEffectCorrection + + +def run(data_path, output, ref, conf): + header = getheader(data_path) + number = header['DETECTOR'][-2:] + number_list = ['06', '07', '08', '09', '11', '12', '13', '14', '15', + '16', '17', '18', '19', '20', '22', '23', '24', '25'] + if number not in number_list: + print(number, 'not in') + return + bias = getdata(glob(os.path.join(ref, "*CLB*_" + number + '_*'))[0]) + dark = getdata(glob(os.path.join(ref, "*CLD*_" + number + '_*'))[0]) + flat = getdata(glob(os.path.join(ref, "*CLF*_" + number + '_*'))[0]) + iec = InstrumentEffectCorrection( + data_input=data_path, + bias_input=bias, + dark_input=dark, + flat_input=flat, + output_path=output, + config_path=conf + ) + iec.run() + print(number, 'done') + + +if __name__ == '__main__': + input_path = 'D:/Desktop/data/' + output_path = 'D:/Desktop/data/output/' + ref_path = r'D:/Desktop/data/ref/' + conf_path = 'D:/Desktop/test/csst-msc-instrument-effect/CSST_2021-12-30_CCD23_epoch20.pth' + with ThreadPoolExecutor(18) as tpe: + for path in glob(os.path.join(input_path, '*_MS_SCI_*.fits')): + tpe.submit(run, path, output_path, ref_path, conf_path) diff --git a/setup.py b/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..042793904b64a1e7e547932cea11918c713517e0 --- /dev/null +++ b/setup.py @@ -0,0 +1,14 @@ +from setuptools import setup + + +setup( + name='csst_msc_instrument_effect', + author='xiezhou', + version='0.1', + packages=['csst_msc_instrument_effect'], + entry_points={ + 'console_scripts': [ + 'csst-msc-iec=csst_msc_instrument_effect.msc_iec:main' + ] + } +)