Commit 3da3cb66 authored by Zhang Xin's avatar Zhang Xin
Browse files

add header info: id, obsflag, filename, shutter info...

parent 349c1960
......@@ -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"):
......
......@@ -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))
......
......@@ -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
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
......@@ -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,
......
......@@ -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
......
......@@ -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
......
......@@ -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":
......
......@@ -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
......@@ -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
......@@ -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]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment