From 3da3cb66bee0fa29b74c4fd1eade745b8e96ed82 Mon Sep 17 00:00:00 2001 From: zhangxin Date: Wed, 10 Apr 2024 16:50:31 +0800 Subject: [PATCH] add header info: id, obsflag, filename, shutter info... --- ObservationSim/Config/ChipOutput.py | 6 ++-- ObservationSim/Config/Header/ImageHeader.py | 31 +++++++++++------ ObservationSim/Config/Pointing.py | 34 +++++++++++++++---- ObservationSim/Config/_util.py | 14 +++++--- ObservationSim/ObservationSim.py | 4 ++- ObservationSim/sim_steps/__init__.py | 2 +- ObservationSim/sim_steps/add_objects.py | 3 ++ .../sim_steps/add_sky_background.py | 6 ++++ ObservationSim/sim_steps/prepare_headers.py | 18 ++++++++-- ObservationSim/sim_steps/readout_output.py | 12 ++++++- config/obs_config_SCI_WIDE_phot.yaml | 4 ++- 11 files changed, 102 insertions(+), 32 deletions(-) diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index 38862e2..4ad0762 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -3,7 +3,7 @@ import logging import ObservationSim.Config._util as _util class ChipOutput(object): - def __init__(self, config, focal_plane, chip, filt, imgKey0="", imgKey1="", imgKey2="", exptime=150., mjdTime="", ra_cen=None, dec_cen=None, pointing_type='SCIE', pointing_ID='0', subdir="./", prefix=""): + def __init__(self, config, focal_plane, chip, filt, imgKey0="", imgKey1="", imgKey2="", exptime=150., mjdTime="", ra_cen=None, dec_cen=None, pointing_type='SCI', pointing_ID='0',pointing_type_code='101', subdir="./", prefix=""): self.focal_plane = focal_plane self.chip = chip self.filt = filt @@ -21,7 +21,7 @@ class ChipOutput(object): self.dec_cen = config["obs_setting"]["dec_center"] self.chipLabel = focal_plane.getChipLabel(chip.chipID) - obs_id = _util.get_obs_id(img_type=pointing_type, project_cycle=config["project_cycle"], run_counter=config["run_counter"], pointing_num=pointing_ID) + obs_id = _util.get_obs_id(img_type=pointing_type, project_cycle=config["project_cycle"], run_counter=config["run_counter"], pointing_id=pointing_ID, pointing_type_code = pointing_type_code) # self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat" self.cat_name = "MSC_%s_chip_%s_filt_%s"%(obs_id, self.chipLabel, filt.filter_type) + ".cat" @@ -60,7 +60,7 @@ class ChipOutput(object): self.hdr += additional_column_names def create_output_file(self): - if self.pointing_type == 'SCIE': + if self.pointing_type == 'SCI': self.cat = open(os.path.join(self.subdir, self.cat_name), "w") self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name))) if not self.hdr.endswith("\n"): diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index cd022b5..3566a8d 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -343,7 +343,7 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r #TODO project_cycle is temporary， is not in header defined, delete in future -def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, time_pt = None, im_type = 'SCI', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, run_counter=0, chip_name="01"): +def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointing_id = '00000001', pointing_type_code='101', ra = 60, dec = -40, pixel_scale = 0.074, time_pt = None, im_type = 'SCI', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, run_counter=0, chip_name="01"): # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0 @@ -401,9 +401,10 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec # # OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0') # OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + str(int(run_counter)).rjust(2, '0') + pointNum.rjust(5,'0') - OBS_id = get_obs_id(img_type=im_type, project_cycle=project_cycle, run_counter=run_counter, pointing_num=pointNum) + OBS_id = get_obs_id(img_type=im_type, project_cycle=project_cycle, run_counter=run_counter, pointing_id=pointing_id, pointing_type_code = pointing_type_code) - h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7, '0') + # h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7, '0') + h_prim['OBJECT'] = pointing_id h_prim['OBSID'] = OBS_id # h_prim['TELFOCUS'] = 'f/14' h_prim['EXPTIME'] = exptime @@ -472,7 +473,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec return h_prim def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, pixel_scale = 0.074, pixel_size=1e-2, - extName='SCIE', row_num = None, col_num = None, xcen=None, ycen=None, timestamp = 1621915200,exptime = 150., readoutTime = 40.): + extName='SCIE', row_num = None, col_num = None, xcen=None, ycen=None, timestamp = 1621915200,exptime = 150., readoutTime = 40., t_shutter_open = 1.3, t_shutter_close = 1.3): e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/csst_msc_l0_ms.fits' f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst') @@ -547,19 +548,27 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p h_ext['PIXSCAL1'] = pixel_scale h_ext['PIXSCAL2'] = pixel_scale h_ext['EXPTIME'] = exptime - h_ext['DARKTIME'] = exptime + 2 + h_ext['DARKTIME'] = exptime + readoutTime datetime_obs = datetime.utcfromtimestamp(timestamp) tstart = Time(datetime_obs) + t_shutter_os = tstart + t_shutter_oe = Time(tstart.mjd + t_shutter_open / 86400., format="mjd") + t_shutter_co = Time(tstart.mjd + exptime / 86400., format="mjd") + t_shutter_ce = Time(tstart.mjd + (exptime + t_shutter_close) / 86400., format="mjd") + t_shutter_os1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_os.unix).timestamp(), 1)) + h_ext['SHTOPEN0'] = t_shutter_os1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] + t_shutter_oe1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_oe.unix).timestamp(), 1)) + h_ext['SHTOPEN1'] = t_shutter_oe1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] + t_shutter_co1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_co.unix).timestamp(), 1)) + h_ext['SHTCLOS0'] = t_shutter_co1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] + t_shutter_ce1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_shutter_ce.unix).timestamp(), 1)) + h_ext['SHTCLOS1'] = t_shutter_ce1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] + tstart_read = Time(tstart.mjd + exptime / 86400., format="mjd") tend_read = Time(tstart.mjd + (exptime + readoutTime) / 86400., format="mjd") - # t_s1 = str(tstart_read.datetime).split() - # h_ext['ROTIME0'] = t_s1[0]+'T'+t_s1[1] - # t_s2 = str(tend_read.datetime).split() - # h_ext['ROTIME1'] = t_s2[0] + 'T' + t_s2[1] - # tstart1=tstart.datetime.replace(microsecond=round(tstart.datetime.microsecond, -5)) - tstart1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tstart.unix).timestamp(), 1)) + tstart1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tstart_read.unix).timestamp(), 1)) h_ext['ROTIME0'] = tstart1.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5] # tend_read1 = tend_read.datetime.replace(microsecond=round(tend_read.datetime.microsecond, -5)) tend_read1 = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(tend_read.unix).timestamp(), 1)) diff --git a/ObservationSim/Config/Pointing.py b/ObservationSim/Config/Pointing.py index 78b15b9..a982ed3 100644 --- a/ObservationSim/Config/Pointing.py +++ b/ObservationSim/Config/Pointing.py @@ -6,7 +6,7 @@ from astropy.time import Time import ObservationSim.Instrument._util as _util class Pointing(object): - def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='SCIE', obs_config_file=None): + def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='SCI', pointing_type_code='101', pointing_id = '00000001', obs_config_file=None, t_shutter_open = 1.3, t_shutter_close = 1.3): self.id = id self.ra = ra self.dec = dec @@ -17,9 +17,26 @@ class Pointing(object): self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz self.exp_time = exp_time self.pointing_type = pointing_type + self.pointing_type_code = pointing_type_code + self.obs_id = pointing_id self.survey_field_type = 'WIDE' self.jdt = 0. self.obs_config_file = obs_config_file + self.t_shutter_open = t_shutter_open + self.t_shutter_close = t_shutter_close + if self.obs_config_file is not None: + with open(self.obs_config_file, "r") as stream: + try: + self.obs_param = yaml.safe_load(stream) + except yaml.YAMLError as exc: + print(exc) + if self.obs_param["obs_type"]: + self.pointing_type = self.obs_param["obs_type"] + if self.obs_param["obs_type_code"]: + self.pointing_type_code = self.obs_param["obs_type_code"] + if self.obs_param["obs_id"]: + self.obs_id = str(self.obs_param["obs_id"]) + def get_full_depth_exptime(self, filter_type): if self.survey_field_type == 'WIDE': @@ -44,7 +61,7 @@ class Pointing(object): return max(150., self.exp_time) # [TODO] for FGS - def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCIE'): + def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCI'): self.id = id col_len = len(columns) self.ra = float(columns[0]) @@ -74,12 +91,15 @@ class Pointing(object): if not self.obs_config_file: self.obs_config_file = str(columns[20]) + self.pointing_type_code = columns[21][0:3] + self.obs_id = columns[21][3:] + # Load the configuration file for this particular pointing - with open(self.obs_config_file, "r") as stream: - try: - self.obs_param = yaml.safe_load(stream) - except yaml.YAMLError as exc: - print(exc) + # with open(self.obs_config_file, "r") as stream: + # try: + # self.obs_param = yaml.safe_load(stream) + # except yaml.YAMLError as exc: + # print(exc) self.pointing_type = self.obs_param["obs_type"] else: self.timestamp = t diff --git a/ObservationSim/Config/_util.py b/ObservationSim/Config/_util.py index 2625ab7..7e49dee 100644 --- a/ObservationSim/Config/_util.py +++ b/ObservationSim/Config/_util.py @@ -1,10 +1,16 @@ -def get_obs_id(img_type='SCI', project_cycle=6, run_counter=0, pointing_num=0): +def get_obs_id(img_type='SCI', project_cycle=6, run_counter=0, pointing_id='00000001',pointing_type_code='101'): # obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'} - obs_type = {'SCIE': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'CAL': '01'} - obs_id = '1'+ obs_type[img_type] + str(int(project_cycle)).rjust(2, '0') + str(int(run_counter)) + str(pointing_num).rjust(8,'0') + # obs_type = {'SCIE': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'CAL': '01'} + obs_id = pointing_type_code + str(int(project_cycle)).rjust(2, '0') + str(int(run_counter)) + pointing_id return obs_id +# def get_obs_id(img_type='SCI', project_cycle=6, run_counter=0, pointing_num=0): +# # obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'} +# obs_type = {'SCIE': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'CAL': '01'} +# obs_id = '1'+ obs_type[img_type] + str(int(project_cycle)).rjust(2, '0') + str(int(run_counter)) + str(pointing_num).rjust(8,'0') +# return obs_id + def get_file_type(img_type='SCI'): - file_type = {'SCIE':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'} + file_type = {'SCI':'SCI', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'} return file_type[img_type] \ No newline at end of file diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 3960406..5960c3a 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -99,6 +99,7 @@ class Observation(object): sim_steps = SimSteps(overall_config=self.config, chip_output=chip_output, all_filters=self.all_filters) for step in pointing.obs_param["call_sequence"]: + print(step) if self.config["run_option"]["out_cat_only"]: if step != "scie_obs": continue @@ -321,7 +322,8 @@ class Observation(object): filt=filt, exptime=pointing.exp_time, pointing_type=pointing.pointing_type, - pointing_ID=pointing_ID, + pointing_ID=pointing.obs_id, + pointing_type_code = pointing.pointing_type_code, subdir=sub_img_dir, prefix=prefix, ra_cen=pointing.ra, diff --git a/ObservationSim/sim_steps/__init__.py b/ObservationSim/sim_steps/__init__.py index 0c928f6..65a8e21 100644 --- a/ObservationSim/sim_steps/__init__.py +++ b/ObservationSim/sim_steps/__init__.py @@ -6,7 +6,7 @@ class SimSteps: self.chip_output = chip_output self.all_filters = all_filters - from .prepare_headers import prepare_headers + from .prepare_headers import prepare_headers, updateHeaderInfo from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background from .add_objects import add_objects from .add_cosmic_rays import add_cosmic_rays diff --git a/ObservationSim/sim_steps/add_objects.py b/ObservationSim/sim_steps/add_objects.py index 94ef1c9..8c80e36 100644 --- a/ObservationSim/sim_steps/add_objects.py +++ b/ObservationSim/sim_steps/add_objects.py @@ -197,6 +197,9 @@ def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): if obs_param["shutter_effect"] == True: flat_normal = flat_normal * chip.shutter_img flat_normal = np.array(flat_normal, dtype='float32') + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True]) + else: + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','','']) chip.img *= flat_normal del flat_normal diff --git a/ObservationSim/sim_steps/add_sky_background.py b/ObservationSim/sim_steps/add_sky_background.py index 555aef7..52d13bd 100644 --- a/ObservationSim/sim_steps/add_sky_background.py +++ b/ObservationSim/sim_steps/add_sky_background.py @@ -10,6 +10,9 @@ def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param): if obs_param["shutter_effect"] == True: flat_normal = flat_normal * chip.shutter_img flat_normal = np.array(flat_normal, dtype='float32') + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True]) + else: + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','','']) if obs_param["enable_straylight_model"]: # Filter.sky_background, Filter.zodical_spec will be updated @@ -64,6 +67,9 @@ def add_sky_flat_calibration(self, chip, filt, tel, pointing, catalog, obs_param if self.overall_config["output_setting"]["shutter_output"] == True: # output 16-bit shutter effect image with pixel value <=65535 shutt_gsimg = galsim.ImageUS(chip.shutter_img*6E4) shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (self.chip_output.subdir, chip.chipID)) + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True]) + else: + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','','']) if chip.survey_type == "photometric": diff --git a/ObservationSim/sim_steps/prepare_headers.py b/ObservationSim/sim_steps/prepare_headers.py index 935add5..bfd88e9 100644 --- a/ObservationSim/sim_steps/prepare_headers.py +++ b/ObservationSim/sim_steps/prepare_headers.py @@ -4,7 +4,8 @@ def prepare_headers(self, chip, pointing): self.h_prim = generatePrimaryHeader( xlen=chip.npix_x, ylen=chip.npix_y, - pointNum = str(pointing.id), + pointing_id = pointing.obs_id, + pointing_type_code = pointing.pointing_type_code, ra=pointing.ra, dec=pointing.dec, pixel_scale=chip.pix_scale, @@ -34,5 +35,16 @@ def prepare_headers(self, chip, pointing): extName=pointing.pointing_type, timestamp = pointing.timestamp, exptime = pointing.exp_time, - readoutTime = chip.readout_time) - return self.h_prim, self.h_ext \ No newline at end of file + readoutTime = chip.readout_time, + t_shutter_open = pointing.t_shutter_open, + t_shutter_close = pointing.t_shutter_close) + + return self.h_prim, self.h_ext + +def updateHeaderInfo(self,header_flag='prim', keys = ['key'], values = [0]): + if header_flag == 'prim': + for key,value in zip(keys, values): + self.h_prim[key] = value + if header_flag == 'ext': + for key,value in zip(keys, values): + self.h_ext[key] = value diff --git a/ObservationSim/sim_steps/readout_output.py b/ObservationSim/sim_steps/readout_output.py index bd5a519..f893f2f 100644 --- a/ObservationSim/sim_steps/readout_output.py +++ b/ObservationSim/sim_steps/readout_output.py @@ -37,6 +37,7 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param) if not hasattr(self, 'h_ext'): _, _ = self.prepare_headers(chip=chip, pointing=pointing) + self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','','']) if obs_param["format_output"] == True: self.chip_output.Log_info(" Apply 1*16 format") @@ -48,10 +49,19 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param) chip.img.replaceNegative(replace_value=0) chip.img.quantize() chip.img = galsim.Image(chip.img.array, dtype=np.uint16) + fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits') + + f_name_size = 68 + if(len(self.h_prim['FILENAME'])>f_name_size): + self.updateHeaderInfo(header_flag='prim', keys = ['FILENAME'], values = [self.h_prim['FILENAME'][0:f_name_size]]) + + + hdu1 = fits.PrimaryHDU(header=self.h_prim) hdu1.add_checksum() hdu1.header.comments['CHECKSUM'] = 'HDU checksum' hdu1.header.comments['DATASUM'] = 'data unit checksum' + self.updateHeaderInfo(header_flag='ext', keys = ['DATASECT'], values = [str(chip.img.array.shape[1])+'x'+str(chip.img.array.shape[0])]) hdu2 = fits.ImageHDU(chip.img.array, header=self.h_ext) hdu2.add_checksum() hdu2.header.comments['XTENSION'] = 'extension type' @@ -59,6 +69,6 @@ def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param) hdu2.header.comments['DATASUM'] = 'data unit checksum' hdu2.header.comments["XTENSION"] = "image extension" hdu1 = fits.HDUList([hdu1, hdu2]) - fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits') hdu1.writeto(fname, output_verify='ignore', overwrite=True) + return chip, filt, tel, pointing diff --git a/config/obs_config_SCI_WIDE_phot.yaml b/config/obs_config_SCI_WIDE_phot.yaml index c6e81f5..45e1a7f 100644 --- a/config/obs_config_SCI_WIDE_phot.yaml +++ b/config/obs_config_SCI_WIDE_phot.yaml @@ -9,7 +9,9 @@ ############################################### # Observation type -obs_type: "SCIE" +obs_type: "SCI" +obs_type_code: "101" +obs_id: "00000001" # Define list of chips run_chips: [8] -- GitLab