From 633fbbabf587c15a6ef61efdfca9fdf6b6cedb46 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Thu, 18 Jan 2024 02:40:29 +0800 Subject: [PATCH] first test --- ObservationSim/Config/Header/ImageHeader.py | 7 +- ObservationSim/Config/Pointing.py | 14 +- ObservationSim/Config/_util.py | 5 +- ObservationSim/Instrument/Chip/ChipUtils.py | 2 +- ObservationSim/Instrument/FocalPlane.py | 13 +- ObservationSim/ObservationSim.py | 782 ++++++++++---------- ObservationSim/SimSteps.py | 456 ++++++++++++ ObservationSim/_util.py | 38 +- config/config_overall.yaml | 3 - config/obs_config_SCI_WIDE.yaml | 25 - config/obs_config_SCI_WIDE_phot.yaml | 61 ++ run_sim.py | 8 +- 12 files changed, 979 insertions(+), 435 deletions(-) create mode 100644 ObservationSim/SimSteps.py delete mode 100644 config/obs_config_SCI_WIDE.yaml create mode 100644 config/obs_config_SCI_WIDE_phot.yaml diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 171e7ca..764f241 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -342,7 +342,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, date='200930', time_obs='120000', im_type = 'MS', 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, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'SCIE', 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 @@ -391,7 +391,8 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec # # Define file types # file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'} # h_prim['FILETYPE'] = file_type[im_type] - h_prim['FILETYPE'] = get_file_type(img_type=im_type) + # h_prim['FILETYPE'] = get_file_type(img_type=im_type) + h_prim['FILETYPE'] = im_type co = coord.SkyCoord(ra, dec, unit='deg') @@ -450,7 +451,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='SCI', 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.): e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header' f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst') diff --git a/ObservationSim/Config/Pointing.py b/ObservationSim/Config/Pointing.py index dccf918..3db8a98 100644 --- a/ObservationSim/Config/Pointing.py +++ b/ObservationSim/Config/Pointing.py @@ -1,5 +1,6 @@ import numpy as np import galsim +import yaml from astropy.time import Time import ObservationSim.Instrument._util as _util @@ -42,13 +43,13 @@ class Pointing(object): return max(150., self.exp_time) # [TODO] for FGS - def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCI'): + def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCIE'): self.id = id col_len = len(columns) self.ra = float(columns[0]) self.dec = float(columns[1]) self.img_pa = float(columns[4]) * galsim.degrees - self.pointing_type = pointing_type + # self.pointing_type = pointing_type if col_len > 5: jdt = np.double(columns[5]) t_temp = Time(jdt, format='jd') @@ -68,5 +69,14 @@ class Pointing(object): # [TODO] Can also define other survey types if is_deep != -1.0: self.survey_field_type = "DEEP" + + # Load the configuration file for this particular pointing + self.obs_config_file = "/share/home/fangyuedong/20231211/csst-simulation/config/obs_config_SCI_WIDE_phot.yaml" + 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_config_file["obs_type"] else: self.timestamp = t diff --git a/ObservationSim/Config/_util.py b/ObservationSim/Config/_util.py index 964d5a2..b7beed9 100644 --- a/ObservationSim/Config/_util.py +++ b/ObservationSim/Config/_util.py @@ -1,8 +1,7 @@ -def get_obs_id(img_type='SCI', project_cycle=6, run_counter=0, pointing_num=0): +def get_obs_id(img_type='SCIE', 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 = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99', 'CAL': '01'} - # obs_id = '1'+ obs_type[img_type] + str(int(project_cycle)) + str(int(run_counter)).rjust(2, '0') + str(pointing_num).rjust(5,'0') + 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 diff --git a/ObservationSim/Instrument/Chip/ChipUtils.py b/ObservationSim/Instrument/Chip/ChipUtils.py index 8eb95f3..4b7aa65 100644 --- a/ObservationSim/Instrument/Chip/ChipUtils.py +++ b/ObservationSim/Instrument/Chip/ChipUtils.py @@ -87,7 +87,7 @@ def generateHeader(chip, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime pixel_size=chip.pix_size, xcen=chip.x_cen, ycen=chip.y_cen, - extName='SCI', + extName='SCIE', timestamp = timestamp, exptime = exptime, readoutTime = chip.readout_time) diff --git a/ObservationSim/Instrument/FocalPlane.py b/ObservationSim/Instrument/FocalPlane.py index ddca78f..70ccbef 100755 --- a/ObservationSim/Instrument/FocalPlane.py +++ b/ObservationSim/Instrument/FocalPlane.py @@ -2,17 +2,24 @@ import galsim import numpy as np class FocalPlane(object): - def __init__(self, config=None, survey_type='Photometric', bad_chips=None): + def __init__(self, config=None, chip_list=None, survey_type='Photometric', bad_chips=None): """Get the focal plane layout """ self.nchips = 42 + self.ignore_chips = [] + if bad_chips == None: self.bad_chips = [] else: self.bad_chips = bad_chips + for chip_id in bad_chips: + self.ignore_chips.append(chip_id) - self.ignore_chips = [] - if survey_type == 'Photometric': + if chip_list is not None: + for i in range(42): + if not (i+1 in chip_list): + self.ignore_chips.append(i+1) + elif survey_type == 'Photometric': for i in range(5): self.ignore_chips.append(i+1) self.ignore_chips.append(i+26) diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 4ca354d..9151075 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -2,7 +2,6 @@ import os import numpy as np import mpi4py.MPI as MPI import galsim -import logging import psutil import gc from astropy.io import fits @@ -14,40 +13,44 @@ from ObservationSim.Config import config_dir, ChipOutput from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip from ObservationSim.Instrument.Chip import Effects +from ObservationSim.Instrument.Chip import ChipUtils as chip_utils from ObservationSim.Straylight import calculateSkyMap_split_g from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS from ObservationSim._util import get_shear_field, makeSubDir_PointingList from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position from ObservationSim.MockObject import FlatLED +from ObservationSim.SimSteps import SimSteps, SIM_STEP_TYPES + class Observation(object): def __init__(self, config, Catalog, work_dir=None, data_dir=None): self.path_dict = config_dir(config=config, work_dir=work_dir, data_dir=data_dir) self.config = config self.tel = Telescope() - self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) self.filter_param = FilterParam() - self.chip_list = [] - self.filter_list = [] - self.all_filter = [] self.Catalog = Catalog - # Construct chips & filters: - for i in range(self.focal_plane.nchips): - chipID = i + 1 + def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing): + # Get WCS for the focal plane + if wcs_fp == None: + wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale) - # Make Chip & Filter lists - chip = Chip( - chipID=chipID, - config=self.config) - filter_id, filter_type = chip.getChipFilter() - filt = Filter(filter_id=filter_id, - filter_type=filter_type, - filter_param=self.filter_param) - if not self.focal_plane.isIgnored(chipID=chipID): - self.chip_list.append(chip) - self.filter_list.append(filt) - self.all_filter.append(filt) + # Create chip Image + chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + chip.img.wcs = wcs_fp + + # Get random generators for this chip + chip.rng_poisson, chip.poisson_noise = chip_utils.get_poisson( + seed=int(self.config["random_seeds"]["seed_poisson"]) + pointing.id*30 + chip.chipID, sky_level=0.) + + # Get flat, shutter, and PRNU images + _, chip.flat_normal = chip_utils.get_flat(img=chip.img, seed=int(self.config["random_seeds"]["seed_flat"])) + chip.shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735, dt=1E-3) + chip.prnu_img = Effects.PRNU_Img(xsize=chip.npix_x, ysize=chip.npix_y, sigma=0.01, + seed=int(self.config["random_seeds"]["seed_prnu"]+chip.chipID)) + + return chip def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, cat_dir=None, sed_dir=None): @@ -61,19 +64,6 @@ class Observation(object): chip_output.Log_info('Chip : %d' % chip.chipID) chip_output.Log_info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::') - if self.config["psf_setting"]["psf_model"] == "Gauss": - psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"]) - elif self.config["psf_setting"]["psf_model"] == "Interp": - if chip.survey_type == "spectroscopic": - psf_model = PSFInterpSLS(chip, filt,PSF_data_prefix=self.path_dict["psf_sls_dir"]) - else: - psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.path_dict["psf_dir"]) - else: - chip_output.Log_error("unrecognized PSF model type!!", flush=True) - - # Figure out shear fields - self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config) - # Apply astrometric simulation for pointing if self.config["obs_setting"]["enable_astrometric_model"]: dt = datetime.utcfromtimestamp(pointing.timestamp) @@ -102,249 +92,272 @@ class Observation(object): ra_cen = pointing.ra dec_cen = pointing.dec - # Get WCS for the focal plane - if wcs_fp == None: - wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale) + # # Get WCS for the focal plane + # if wcs_fp == None: + # wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale) - # Create chip Image - chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) - chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) - chip.img.wcs = wcs_fp + # # Create chip Image + # chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) + # chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + # chip.img.wcs = wcs_fp + chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing) - if self.config["obs_setting"]["enable_straylight_model"]: - filt.setFilterStrayLightPixel(jtime = pointing.jdt, sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]), pointing_radec = np.array([pointing.ra,pointing.dec]), sun_pos = np.array([pointing.sun_x,pointing.sun_y,pointing.sun_z])) + # Load catalogues + self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output, filt=filt) - chip_output.Log_info("========================sky pix========================") - chip_output.Log_info(filt.sky_background) + # Initialize SimSteps + sim_steps = SimSteps(overall_config=self.config, chip_output=chip_output, all_filters=self.all_filters) - if chip.survey_type == "photometric": - sky_map = None - elif chip.survey_type == "spectroscopic": - # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits') - flat_normal = np.ones_like(chip.img.array) - if self.config["ins_effects"]["flat_fielding"] == True: - chip_output.Log_info("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID) - msg = str(chip.img.bounds) - chip_output.Log_info(msg) - flat_img = Effects.MakeFlatSmooth( - chip.img.bounds, - int(self.config["random_seeds"]["seed_flat"])) - flat_normal = flat_normal * flat_img.array / np.mean(flat_img.array) - if self.config["ins_effects"]["shutter_effect"] == True: - chip_output.Log_info("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID) - shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735, - dt=1E-3) # shutter effect normalized image for this chip - flat_normal = flat_normal*shuttimg - flat_normal = np.array(flat_normal,dtype='float32') - sky_map = calculateSkyMap_split_g( - skyMap=flat_normal, - blueLimit=filt.blue_limit, - redLimit=filt.red_limit, - conf=chip.sls_conf, - pixelSize=chip.pix_scale, - isAlongY=0, - flat_cube=chip.flat_cube, zoldial_spec = filt.zodical_spec) - sky_map = sky_map+filt.sky_background - del flat_normal - - if pointing.pointing_type == 'SCI': - # Load catalogues and templates - self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output, filt=filt) - chip_output.create_output_file() - self.nobj = len(self.cat.objs) - - for ifilt in range(len(self.all_filter)): - temp_filter = self.all_filter[ifilt] - # Update the limiting magnitude using exposure time in pointing - temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip) - - # Select cutting band filter for saturation/limiting magnitude - if temp_filter.filter_type.lower() == self.config["obs_setting"]["cut_in_band"].lower(): - cut_filter = temp_filter - - if self.config["ins_effects"]["field_dist"] == True: - self.fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) - else: - self.fd_model = None - - # Loop over objects - missed_obj = 0 - bright_obj = 0 - dim_obj = 0 + for step in pointing.obs_param["call_sequence"]: + obs_param = pointing.obs_param["call_sequence"][step] + step_name = SIM_STEP_TYPES[step] + try: + step_func = getattr(sim_steps, step_name) + chip, filt, tel, pointing = step_func( + chip=chip, + filt=filt, + tel=tel, + pointing=pointing, + catalog=self.cat, + obs_param=obs_param) + except Exception as e: + traceback.print_exc() + chip_output.Log_error(e) + continue + + # if self.config["obs_setting"]["enable_straylight_model"]: + # filt.setFilterStrayLightPixel(jtime = pointing.jdt, sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]), pointing_radec = np.array([pointing.ra,pointing.dec]), sun_pos = np.array([pointing.sun_x,pointing.sun_y,pointing.sun_z])) + # chip_output.Log_info("========================sky pix========================") + # chip_output.Log_info(filt.sky_background) + + # if chip.survey_type == "photometric": + # sky_map = None + # elif chip.survey_type == "spectroscopic": + # # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits') + # flat_normal = np.ones_like(chip.img.array) + # if self.config["ins_effects"]["flat_fielding"] == True: + # chip_output.Log_info("SLS flat preprocess,CHIP %d : Creating and applying Flat-Fielding"%chip.chipID) + # msg = str(chip.img.bounds) + # chip_output.Log_info(msg) + # flat_img = Effects.MakeFlatSmooth( + # chip.img.bounds, + # int(self.config["random_seeds"]["seed_flat"])) + # flat_normal = flat_normal * flat_img.array / np.mean(flat_img.array) + # if self.config["ins_effects"]["shutter_effect"] == True: + # chip_output.Log_info("SLS flat preprocess,CHIP %d : Apply shutter effect"%chip.chipID) + # shuttimg = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735, + # dt=1E-3) # shutter effect normalized image for this chip + # flat_normal = flat_normal*shuttimg + # flat_normal = np.array(flat_normal,dtype='float32') + # sky_map = calculateSkyMap_split_g( + # skyMap=flat_normal, + # blueLimit=filt.blue_limit, + # redLimit=filt.red_limit, + # conf=chip.sls_conf, + # pixelSize=chip.pix_scale, + # isAlongY=0, + # flat_cube=chip.flat_cube, zoldial_spec = filt.zodical_spec) + # sky_map = sky_map+filt.sky_background + # del flat_normal + + # if pointing.pointing_type == 'SCI': + # # Load catalogues and templates + # self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output, filt=filt) + # chip_output.create_output_file() + # self.nobj = len(self.cat.objs) + + # for ifilt in range(len(self.all_filter)): + # temp_filter = self.all_filter[ifilt] + # # Update the limiting magnitude using exposure time in pointing + # temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip) + + # # Select cutting band filter for saturation/limiting magnitude + # if temp_filter.filter_type.lower() == self.config["obs_setting"]["cut_in_band"].lower(): + # cut_filter = temp_filter + + # if self.config["ins_effects"]["field_dist"] == True: + # self.fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) + # else: + # self.fd_model = None + + # # Loop over objects + # missed_obj = 0 + # bright_obj = 0 + # dim_obj = 0 - h_ext = generateExtensionHeader( - chip=chip, - xlen=chip.npix_x, - ylen=chip.npix_y, - ra=pointing.ra, - dec=pointing.dec, - pa=pointing.img_pa.deg, - gain=chip.gain, - readout=chip.read_noise, - dark=chip.dark_noise, - saturation=90000, - pixel_scale=chip.pix_scale, - pixel_size=chip.pix_size, - xcen=chip.x_cen, - ycen=chip.y_cen, - extName='SCI', - timestamp = pointing.timestamp, - exptime = pointing.exp_time, - readoutTime = chip.readout_time) + # h_ext = generateExtensionHeader( + # chip=chip, + # xlen=chip.npix_x, + # ylen=chip.npix_y, + # ra=pointing.ra, + # dec=pointing.dec, + # pa=pointing.img_pa.deg, + # gain=chip.gain, + # readout=chip.read_noise, + # dark=chip.dark_noise, + # saturation=90000, + # pixel_scale=chip.pix_scale, + # pixel_size=chip.pix_size, + # xcen=chip.x_cen, + # ycen=chip.y_cen, + # extName='SCI', + # timestamp = pointing.timestamp, + # exptime = pointing.exp_time, + # readoutTime = chip.readout_time) - chip_wcs = galsim.FitsWCS(header=h_ext) + # chip_wcs = galsim.FitsWCS(header=h_ext) - for j in range(self.nobj): + # for j in range(self.nobj): - # # (DEBUG) - # if j >= 10: - # break - - obj = self.cat.objs[j] - - # (DEBUG) - # if obj.getMagFilter(filt)>20: - # continue - - # load and convert SED; also caculate object's magnitude in all CSST bands - try: - sed_data = self.cat.load_sed(obj) - norm_filt = self.cat.load_norm_filt(obj) - obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = self.cat.convert_sed( - mag=obj.param["mag_use_normal"], - sed=sed_data, - target_filt=filt, - norm_filt=norm_filt, - ) - _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = self.cat.convert_sed( - mag=obj.param["mag_use_normal"], - sed=sed_data, - target_filt=cut_filter, - norm_filt=norm_filt, - ) - except Exception as e: - traceback.print_exc() - chip_output.Log_error(e) - continue + # # # (DEBUG) + # # if j >= 10: + # # break + + # obj = self.cat.objs[j] + + # # (DEBUG) + # # if obj.getMagFilter(filt)>20: + # # continue + + # # load and convert SED; also caculate object's magnitude in all CSST bands + # try: + # sed_data = self.cat.load_sed(obj) + # norm_filt = self.cat.load_norm_filt(obj) + # obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = self.cat.convert_sed( + # mag=obj.param["mag_use_normal"], + # sed=sed_data, + # target_filt=filt, + # norm_filt=norm_filt, + # ) + # _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = self.cat.convert_sed( + # mag=obj.param["mag_use_normal"], + # sed=sed_data, + # target_filt=cut_filter, + # norm_filt=norm_filt, + # ) + # except Exception as e: + # traceback.print_exc() + # chip_output.Log_error(e) + # continue - # [TODO] Testing - # chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()])) - - # Exclude very bright/dim objects (for now) - if cut_filter.is_too_bright( - mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()], - margin=self.config["obs_setting"]["mag_sat_margin"]): - chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()])) - bright_obj += 1 - obj.unload_SED() - continue - if filt.is_too_dim( - mag=obj.getMagFilter(filt), - margin=self.config["obs_setting"]["mag_lim_margin"]): - chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt))) - dim_obj += 1 - obj.unload_SED() - continue - - # Get corresponding shear values - if self.config["shear_setting"]["shear_type"] == "constant": - if obj.type == 'star': - obj.g1, obj.g2 = 0., 0. - else: - obj.g1, obj.g2 = self.g1_field, self.g2_field - elif self.config["shear_setting"]["shear_type"] == "catalog": - pass - else: - chip_output.Log_error("Unknown shear input") - raise ValueError("Unknown shear input") - - # Get position of object on the focal plane - pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=h_ext) - - # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane - # Otherwise they will be considered missed objects - # if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)): - if pos_img.x == -1 or pos_img.y == -1: - chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f'%(obj.ra, obj.dec, obj.ra_orig, obj.dec_orig)) - chip_output.Log_error("Objected missed: %s"%(obj.id)) - missed_obj += 1 - obj.unload_SED() - continue - - # Draw object & update output catalog - try: - if self.config["run_option"]["out_cat_only"]: - isUpdated = True - obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs) - pos_shear = 0. - elif chip.survey_type == "photometric" and not self.config["run_option"]["out_cat_only"]: - isUpdated, pos_shear = obj.drawObj_multiband( - tel=self.tel, - pos_img=pos_img, - psf_model=psf_model, - bandpass_list=filt.bandpass_sub_list, - filt=filt, - chip=chip, - g1=obj.g1, - g2=obj.g2, - exptime=pointing.exp_time, - fd_shear=fd_shear) - - elif chip.survey_type == "spectroscopic" and not self.config["run_option"]["out_cat_only"]: - isUpdated, pos_shear = obj.drawObj_slitless( - tel=self.tel, - pos_img=pos_img, - psf_model=psf_model, - bandpass_list=filt.bandpass_sub_list, - filt=filt, - chip=chip, - g1=obj.g1, - g2=obj.g2, - exptime=pointing.exp_time, - normFilter=norm_filt, - fd_shear=fd_shear) - - if isUpdated == 1 and self.config["run_option"]["out_psf"]: - obj.drawObj_PSF( - tel=self.tel, - pos_img=pos_img, - psf_model=psf_model, - bandpass_list=filt.bandpass_sub_list, - filt=filt, - chip=chip, - g1=obj.g1, - g2=obj.g2, - exptime=pointing.exp_time, - fd_shear=fd_shear, - chip_output=chip_output) - - if isUpdated == 1: - # TODO: add up stats - chip_output.cat_add_obj(obj, pos_img, pos_shear) - pass - elif isUpdated == 0: - missed_obj += 1 - chip_output.Log_error("Objected missed: %s"%(obj.id)) - else: - chip_output.Log_error("Draw error, object omitted: %s"%(obj.id)) - continue - except Exception as e: - traceback.print_exc() - chip_output.Log_error(e) + # # [TODO] Testing + # # chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()])) + + # # Exclude very bright/dim objects (for now) + # if cut_filter.is_too_bright( + # mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()], + # margin=self.config["obs_setting"]["mag_sat_margin"]): + # chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()])) + # bright_obj += 1 + # obj.unload_SED() + # continue + # if filt.is_too_dim( + # mag=obj.getMagFilter(filt), + # margin=self.config["obs_setting"]["mag_lim_margin"]): + # chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt))) + # dim_obj += 1 + # obj.unload_SED() + # continue + + # # Get corresponding shear values + # if self.config["shear_setting"]["shear_type"] == "constant": + # if obj.type == 'star': + # obj.g1, obj.g2 = 0., 0. + # else: + # obj.g1, obj.g2 = self.g1_field, self.g2_field + # elif self.config["shear_setting"]["shear_type"] == "catalog": + # pass + # else: + # chip_output.Log_error("Unknown shear input") + # raise ValueError("Unknown shear input") + + # # Get position of object on the focal plane + # pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=h_ext) + + # # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane + # # Otherwise they will be considered missed objects + # # if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)): + # if pos_img.x == -1 or pos_img.y == -1: + # chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f'%(obj.ra, obj.dec, obj.ra_orig, obj.dec_orig)) + # chip_output.Log_error("Objected missed: %s"%(obj.id)) + # missed_obj += 1 + # obj.unload_SED() + # continue + + # # Draw object & update output catalog + # try: + # if self.config["run_option"]["out_cat_only"]: + # isUpdated = True + # obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs) + # pos_shear = 0. + # elif chip.survey_type == "photometric" and not self.config["run_option"]["out_cat_only"]: + # isUpdated, pos_shear = obj.drawObj_multiband( + # tel=self.tel, + # pos_img=pos_img, + # psf_model=psf_model, + # bandpass_list=filt.bandpass_sub_list, + # filt=filt, + # chip=chip, + # g1=obj.g1, + # g2=obj.g2, + # exptime=pointing.exp_time, + # fd_shear=fd_shear) + + # elif chip.survey_type == "spectroscopic" and not self.config["run_option"]["out_cat_only"]: + # isUpdated, pos_shear = obj.drawObj_slitless( + # tel=self.tel, + # pos_img=pos_img, + # psf_model=psf_model, + # bandpass_list=filt.bandpass_sub_list, + # filt=filt, + # chip=chip, + # g1=obj.g1, + # g2=obj.g2, + # exptime=pointing.exp_time, + # normFilter=norm_filt, + # fd_shear=fd_shear) + + # if isUpdated == 1 and self.config["run_option"]["out_psf"]: + # obj.drawObj_PSF( + # tel=self.tel, + # pos_img=pos_img, + # psf_model=psf_model, + # bandpass_list=filt.bandpass_sub_list, + # filt=filt, + # chip=chip, + # g1=obj.g1, + # g2=obj.g2, + # exptime=pointing.exp_time, + # fd_shear=fd_shear, + # chip_output=chip_output) + + # if isUpdated == 1: + # # TODO: add up stats + # chip_output.cat_add_obj(obj, pos_img, pos_shear) + # pass + # elif isUpdated == 0: + # missed_obj += 1 + # chip_output.Log_error("Objected missed: %s"%(obj.id)) + # else: + # chip_output.Log_error("Draw error, object omitted: %s"%(obj.id)) + # continue + # except Exception as e: + # traceback.print_exc() + # chip_output.Log_error(e) - # # [C6 TEST] - # chip_output.Log_info("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) - # chip_output.Log_info('draw object %s'%obj.id) - # chip_output.Log_info('mag = %.3f'%obj.param['mag_use_normal']) - - # Unload SED: - obj.unload_SED() - del obj - gc.collect() + # # # [C6 TEST] + # # chip_output.Log_info("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) + # # chip_output.Log_info('draw object %s'%obj.id) + # # chip_output.Log_info('mag = %.3f'%obj.param['mag_use_normal']) + + # # Unload SED: + # obj.unload_SED() + # del obj + # gc.collect() - del psf_model - del self.cat - gc.collect() + # del psf_model + # del self.cat + # gc.collect() chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) @@ -352,79 +365,79 @@ class Observation(object): # =========================================================== # whether to output zero, dark, flat calibration images. - if not self.config["run_option"]["out_cat_only"]: - chip.img = chip.addEffects( - config=self.config, - img=chip.img, - chip_output=chip_output, - filt=filt, - ra_cen=pointing.ra, - dec_cen=pointing.dec, - img_rot=pointing.img_pa, - exptime=pointing.exp_time, - pointing_ID=pointing.id, - timestamp_obs=pointing.timestamp, - pointing_type=pointing.pointing_type, - sky_map=sky_map, tel = self.tel, - logger=chip_output.logger) + # if not self.config["run_option"]["out_cat_only"]: + # chip.img = chip.addEffects( + # config=self.config, + # img=chip.img, + # chip_output=chip_output, + # filt=filt, + # ra_cen=pointing.ra, + # dec_cen=pointing.dec, + # img_rot=pointing.img_pa, + # exptime=pointing.exp_time, + # pointing_ID=pointing.id, + # timestamp_obs=pointing.timestamp, + # pointing_type=pointing.pointing_type, + # sky_map=sky_map, tel = self.tel, + # logger=chip_output.logger) - if pointing.pointing_type == 'SCI': - datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) - date_obs = datetime_obs.strftime("%y%m%d") - time_obs = datetime_obs.strftime("%H%M%S") - h_prim = generatePrimaryHeader( - xlen=chip.npix_x, - ylen=chip.npix_y, - pointNum = str(pointing.id), - ra=pointing.ra, - dec=pointing.dec, - pixel_scale=chip.pix_scale, - date=date_obs, - time_obs=time_obs, - exptime=pointing.exp_time, - im_type='SCI', - sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], - sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz], - project_cycle=self.config["project_cycle"], - run_counter=self.config["run_counter"], - chip_name=str(chip.chipID).rjust(2, '0')) - h_ext = generateExtensionHeader( - chip=chip, - xlen=chip.npix_x, - ylen=chip.npix_y, - ra=pointing.ra, - dec=pointing.dec, - pa=pointing.img_pa.deg, - gain=chip.gain, - readout=chip.read_noise, - dark=chip.dark_noise, - saturation=90000, - pixel_scale=chip.pix_scale, - pixel_size=chip.pix_size, - xcen=chip.x_cen, - ycen=chip.y_cen, - extName='SCI', - timestamp=pointing.timestamp, - exptime=pointing.exp_time, - readoutTime=chip.readout_time) + # if pointing.pointing_type == 'SCIE': + # datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) + # date_obs = datetime_obs.strftime("%y%m%d") + # time_obs = datetime_obs.strftime("%H%M%S") + # h_prim = generatePrimaryHeader( + # xlen=chip.npix_x, + # ylen=chip.npix_y, + # pointNum = str(pointing.id), + # ra=pointing.ra, + # dec=pointing.dec, + # pixel_scale=chip.pix_scale, + # date=date_obs, + # time_obs=time_obs, + # exptime=pointing.exp_time, + # im_type='SCI', + # sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], + # sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz], + # project_cycle=self.config["project_cycle"], + # run_counter=self.config["run_counter"], + # chip_name=str(chip.chipID).rjust(2, '0')) + # h_ext = generateExtensionHeader( + # chip=chip, + # xlen=chip.npix_x, + # ylen=chip.npix_y, + # ra=pointing.ra, + # dec=pointing.dec, + # pa=pointing.img_pa.deg, + # gain=chip.gain, + # readout=chip.read_noise, + # dark=chip.dark_noise, + # saturation=90000, + # pixel_scale=chip.pix_scale, + # pixel_size=chip.pix_size, + # xcen=chip.x_cen, + # ycen=chip.y_cen, + # extName='SCI', + # timestamp=pointing.timestamp, + # exptime=pointing.exp_time, + # readoutTime=chip.readout_time) - chip.img = galsim.Image(chip.img.array, dtype=np.uint16) - hdu1 = fits.PrimaryHDU(header=h_prim) - hdu1.add_checksum() - hdu1.header.comments['CHECKSUM'] = 'HDU checksum' - hdu1.header.comments['DATASUM'] = 'data unit checksum' - hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) - hdu2.add_checksum() - hdu2.header.comments['XTENSION'] = 'extension type' - hdu2.header.comments['CHECKSUM'] = 'HDU checksum' - hdu2.header.comments['DATASUM'] = 'data unit checksum' - hdu1 = fits.HDUList([hdu1, hdu2]) - fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits') - hdu1.writeto(fname, output_verify='ignore', overwrite=True) - - chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, self.nobj)) - chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, self.nobj)) - chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, self.nobj)) + # chip.img = galsim.Image(chip.img.array, dtype=np.uint16) + # hdu1 = fits.PrimaryHDU(header=h_prim) + # hdu1.add_checksum() + # hdu1.header.comments['CHECKSUM'] = 'HDU checksum' + # hdu1.header.comments['DATASUM'] = 'data unit checksum' + # hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) + # hdu2.add_checksum() + # hdu2.header.comments['XTENSION'] = 'extension type' + # hdu2.header.comments['CHECKSUM'] = 'HDU checksum' + # hdu2.header.comments['DATASUM'] = 'data unit checksum' + # hdu1 = fits.HDUList([hdu1, hdu2]) + # fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits') + # hdu1.writeto(fname, output_verify='ignore', overwrite=True) + + # chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, self.nobj)) + # chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, self.nobj)) + # chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, self.nobj)) del chip.img chip_output.Log_info("check running:2: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) @@ -561,38 +574,57 @@ class Observation(object): chip_output.Log_info("check running:2: pointing-%d chip-%d pid-%d memory-%6.2fGB" % ( pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024))) - def runExposure_MPI_PointingList(self, pointing_list,chips=None, use_mpi=False): + def runExposure_MPI_PointingList(self, pointing_list, chips=None, use_mpi=False): if use_mpi: comm = MPI.COMM_WORLD ind_thread = comm.Get_rank() num_thread = comm.Get_size() - if chips is None: - nchips_per_fp = len(self.chip_list) - run_chips = self.chip_list - run_filts = self.filter_list - else: - # Only run a particular set of chips - run_chips = [] - run_filts = [] - nchips_per_fp = len(chips) - for ichip in range(len(self.chip_list)): - chip = self.chip_list[ichip] - filt = self.filter_list[ichip] - if chip.chipID in chips: - run_chips.append(chip) - run_filts.append(filt) - - + process_counter = 0 for ipoint in range(len(pointing_list)): + # Construct chips & filters: + pointing = pointing_list[ipoint] + pointing_ID = pointing.id + self.focal_plane = FocalPlane(chip_list=pointing.obs_param["run_chips"]) + # Make Chip & Filter lists + for i in range(self.focal_plane.nchips): + chipID = i + 1 + self.chip_list = [] + self.filter_list = [] + self.all_filters = [] + chip = Chip(chipID=chipID, config=self.config) + filter_id, filter_type = chip.getChipFilter() + filt = Filter( + filter_id=filter_id, + filter_type=filter_type, + filter_param=self.filter_param) + if not self.focal_plane.isIgnored(chipID=chipID): + self.chip_list.append(chip) + self.filter_list.append(filt) + self.all_filters.append(filt) + + if chips is None: + # Run all chips defined in configuration of this pointing + run_chips = self.chip_list + run_filts = self.filter_list + nchips_per_fp = len(self.chip_list) + else: + # Only run a particular set of chips (defined in the overall config file) + run_chips = [] + run_filts = [] + for ichip in range(len(self.chip_list)): + chip = self.chip_list[ichip] + filt = self.filter_list[ichip] + if chip.chipID in chips: + run_chips.append(chip) + run_filts.append(filt) + nchips_per_fp = len(chips) + for ichip in range(nchips_per_fp): - i = ipoint*nchips_per_fp + ichip - pointing = pointing_list[ipoint] - pointing_ID = pointing.id + i_process = process_counter + ichip if use_mpi: - if i % num_thread != ind_thread: + if i_process % num_thread != ind_thread: continue - pid = os.getpid() sub_img_dir, prefix = makeSubDir_PointingList(path_dict=self.path_dict, config=self.config, pointing_ID=pointing_ID) @@ -611,20 +643,26 @@ class Observation(object): subdir=sub_img_dir, prefix=prefix) chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid)) - if self.config["obs_setting"]["survey_type"] == "CALIBRATION": - self.run_one_chip_calibration(chip=chip, - filt=filt, - chip_output=chip_output, - pointing=pointing, - skyback_level = self.config["obs_setting"]["FLAT_LEVEL"], - sky_level_filt = self.config["obs_setting"]["FLAT_LEVEL_FIL"]) - else: - self.run_one_chip( - chip=chip, - filt=filt, - chip_output=chip_output, - pointing=pointing) + self.run_one_chip( + chip=chip, + filt=filt, + chip_output=chip_output, + pointing=pointing) + # if self.config["obs_setting"]["survey_type"] == "CALIBRATION": + # self.run_one_chip_calibration(chip=chip, + # filt=filt, + # chip_output=chip_output, + # pointing=pointing, + # skyback_level = self.config["obs_setting"]["FLAT_LEVEL"], + # sky_level_filt = self.config["obs_setting"]["FLAT_LEVEL_FIL"]) + # else: + # self.run_one_chip( + # chip=chip, + # filt=filt, + # chip_output=chip_output, + # pointing=pointing) chip_output.Log_info("finished running chip#%d..."%(chip.chipID)) for handler in chip_output.logger.handlers[:]: chip_output.logger.removeHandler(handler) gc.collect() + process_counter += nchips_per_fp diff --git a/ObservationSim/SimSteps.py b/ObservationSim/SimSteps.py new file mode 100644 index 0000000..b111993 --- /dev/null +++ b/ObservationSim/SimSteps.py @@ -0,0 +1,456 @@ +import os +import galsim +import traceback +import gc +import psutil +import numpy as np +from astropy.io import fits +from datetime import datetime +from numpy.random import Generator, PCG64 + +from ObservationSim._util import get_shear_field +from ObservationSim.Straylight import calculateSkyMap_split_g +from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader +from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS +from ObservationSim.Instrument.Chip import ChipUtils as chip_utils +from ObservationSim.Instrument.Chip import Effects +from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim + +class SimSteps: + def __init__(self, overall_config, chip_output, all_filters): + self.overall_config = overall_config + self.chip_output = chip_output + self.all_filters = all_filters + + def prepare_headers(self, chip, pointing): + datetime_obs = datetime.utcfromtimestamp(pointing.timestamp) + date_obs = datetime_obs.strftime("%y%m%d") + time_obs = datetime_obs.strftime("%H%M%S") + self.h_prim = generatePrimaryHeader( + xlen=chip.npix_x, + ylen=chip.npix_y, + pointNum = str(pointing.id), + ra=pointing.ra, + dec=pointing.dec, + pixel_scale=chip.pix_scale, + date=date_obs, + time_obs=time_obs, + exptime=pointing.exp_time, + im_type=pointing.pointing_type, + sat_pos=[pointing.sat_x, pointing.sat_y, pointing.sat_z], + sat_vel=[pointing.sat_vx, pointing.sat_vy, pointing.sat_vz], + project_cycle=self.overall_config["project_cycle"], + run_counter=self.overall_config["run_counter"], + chip_name=str(chip.chipID).rjust(2, '0')) + self.h_ext = generateExtensionHeader( + chip=chip, + xlen=chip.npix_x, + ylen=chip.npix_y, + ra=pointing.ra, + dec=pointing.dec, + pa=pointing.img_pa.deg, + gain=chip.gain, + readout=chip.read_noise, + dark=chip.dark_noise, + saturation=90000, + pixel_scale=chip.pix_scale, + pixel_size=chip.pix_size, + xcen=chip.x_cen, + ycen=chip.y_cen, + extName=pointing.pointing_type, + timestamp = pointing.timestamp, + exptime = pointing.exp_time, + readoutTime = chip.readout_time) + return self.h_prim, self.h_ext + + def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param): + flat_normal = np.ones_like(chip.img.array) + if obs_param["flat_fielding"] == True: + flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array) + if obs_param["shutter_effect"] == True: + flat_normal = flat_normal * chip.shutter_img + flat_normal = np.array(flat_normal, dtype='float32') + + if obs_param["enable_straylight_model"]: + # Filter.sky_background, Filter.zodical_spec will be updated + filt.setFilterStrayLightPixel( + jtime = pointing.jdt, + sat_pos = np.array([pointing.sat_x, pointing.sat_y, pointing.sat_z]), + pointing_radec = np.array([pointing.ra,pointing.dec]), + sun_pos = np.array([pointing.sun_x, pointing.sun_y, pointing.sun_z])) + self.chip_output.Log_info("================================================") + self.chip_output.Log_info("sky background + stray light pixel flux value: %.5f"%(filt.sky_background)) + + if chip.survey_type == "photometric": + sky_map = filt.getSkyNoise(exptime = obs_param["exptime"]) + sky_map = sky_map * np.ones_like(chip.img.array) * flat_normal + sky_map = galsim.Image(array=sky_map) + else: + # chip.loadSLSFLATCUBE(flat_fn='flat_cube.fits') + sky_map = calculateSkyMap_split_g( + skyMap=flat_normal, + blueLimit=filt.blue_limit, + redLimit=filt.red_limit, + conf=chip.sls_conf, + pixelSize=chip.pix_scale, + isAlongY=0, + flat_cube=chip.flat_cube, + zoldial_spec = filt.zodical_spec) + sky_map = sky_map + filt.sky_background + + sky_map = sky_map * tel.pupil_area * obs_param["exptime"] + chip.img += sky_map + return chip, filt, tel, pointing + + def add_objects(self, chip, filt, tel, pointing, catalog, obs_param): + + # Prepare output file(s) for this chip + self.chip_output.create_output_file() + + # Prepare the PSF model + if self.overall_config["psf_setting"]["psf_model"] == "Gauss": + psf_model = PSFGauss(chip=chip, psfRa=self.overall_config["psf_setting"]["psf_rcont"]) + elif self.overall_config["psf_setting"]["psf_model"] == "Interp": + if chip.survey_type == "spectroscopic": + psf_model = PSFInterpSLS(chip, filt, PSF_data_prefix=self.overall_config["psf_setting"]["psf_sls_dir"]) + else: + psf_model = PSFInterp(chip=chip, npsf=chip.n_psf_samples, PSF_data_file=self.overall_config["psf_setting"]["psf_pho_dir"]) + else: + self.chip_output.Log_error("unrecognized PSF model type!!", flush=True) + + # Apply field distortion model + if obs_param["field_dist"] == True: + fd_model = FieldDistortion(chip=chip, img_rot=pointing.img_pa.deg) + else: + fd_model = None + + # Update limiting magnitudes for all filters based on the exposure time + # Get the filter which will be used for magnitude cut + for ifilt in range(len(self.all_filters)): + temp_filter = self.all_filters[ifilt] + temp_filter.update_limit_saturation_mags(exptime=pointing.get_full_depth_exptime(temp_filter.filter_type), chip=chip) + if temp_filter.filter_type.lower() == self.overall_config["obs_setting"]["cut_in_band"].lower(): + cut_filter = temp_filter + + # Read in shear values from configuration file if the constant shear type is used + if self.overall_config["shear_setting"]["shear_type"] == "constant": + g1_field, g2_field, _ = get_shear_field(config=self.overall_config) + self.chip_output.Log_info("Use constant shear: g1=%.5f, g2=%.5f"%(g1_field, g2_field)) + + # Get chip WCS + if not hasattr(self, 'h_ext'): + _, _ = self.prepare_headers(chip=chip, pointing=pointing) + chip_wcs = galsim.FitsWCS(header = self.h_ext) + + # Loop over objects + nobj = len(catalog.objs) + missed_obj = 0 + bright_obj = 0 + dim_obj = 0 + for j in range(nobj): + # # [DEBUG] [TODO] + # if j >= 10: + # break + obj = catalog.objs[j] + + # load and convert SED; also caculate object's magnitude in all CSST bands + try: + sed_data = catalog.load_sed(obj) + norm_filt = catalog.load_norm_filt(obj) + obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = catalog.convert_sed( + mag=obj.param["mag_use_normal"], + sed=sed_data, + target_filt=filt, + norm_filt=norm_filt, + ) + _, obj.param["mag_%s"%cut_filter.filter_type.lower()], obj.param["flux_%s"%cut_filter.filter_type.lower()] = catalog.convert_sed( + mag=obj.param["mag_use_normal"], + sed=sed_data, + target_filt=cut_filter, + norm_filt=norm_filt, + ) + except Exception as e: + traceback.print_exc() + self.chip_output.Log_error(e) + continue + + # [TODO] Testing + # self.chip_output.Log_info("mag_%s = %.3f"%(filt.filter_type.lower(), obj.param["mag_%s"%filt.filter_type.lower()])) + + # Exclude very bright/dim objects (for now) + if cut_filter.is_too_bright( + mag=obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()], + margin=self.overall_config["obs_setting"]["mag_sat_margin"]): + self.chip_output.Log_info("obj %s too birght!! mag_%s = %.3f"%(obj.id, cut_filter.filter_type, obj.param["mag_%s"%self.overall_config["obs_setting"]["cut_in_band"].lower()])) + bright_obj += 1 + obj.unload_SED() + continue + if filt.is_too_dim( + mag=obj.getMagFilter(filt), + margin=self.overall_config["obs_setting"]["mag_lim_margin"]): + self.chip_output.Log_info("obj %s too dim!! mag_%s = %.3f"%(obj.id, filt.filter_type, obj.getMagFilter(filt))) + dim_obj += 1 + obj.unload_SED() + continue + + # Get corresponding shear values + if self.overall_config["shear_setting"]["shear_type"] == "constant": + if obj.type == 'star': + obj.g1, obj.g2 = 0., 0. + else: + # Figure out shear fields from overall configuration shear setting + obj.g1, obj.g2 = g1_field, g2_field + elif self.overall_config["shear_setting"]["shear_type"] == "catalog": + pass + else: + self.chip_output.Log_error("Unknown shear input") + raise ValueError("Unknown shear input") + + # Get position of object on the focal plane + pos_img, _, _, _, fd_shear = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=fd_model, chip=chip, verbose=False, chip_wcs=chip_wcs, img_header=self.h_ext) + + # [TODO] For now, only consider objects which their centers (after field distortion) are projected within the focal plane + # Otherwise they will be considered missed objects + # if pos_img.x == -1 or pos_img.y == -1 or (not chip.isContainObj(x_image=pos_img.x, y_image=pos_img.y, margin=0.)): + if pos_img.x == -1 or pos_img.y == -1: + self.chip_output.Log_info('obj_ra = %.6f, obj_dec = %.6f, obj_ra_orig = %.6f, obj_dec_orig = %.6f'%(obj.ra, obj.dec, obj.ra_orig, obj.dec_orig)) + self.chip_output.Log_error("Objected missed: %s"%(obj.id)) + missed_obj += 1 + obj.unload_SED() + continue + + # Draw object & update output catalog + try: + if self.overall_config["run_option"]["out_cat_only"]: + isUpdated = True + obj.real_pos = obj.getRealPos(chip.img, global_x=obj.posImg.x, global_y=obj.posImg.y, img_real_wcs=obj.chip_wcs) + pos_shear = 0. + elif chip.survey_type == "photometric" and not self.overall_config["run_option"]["out_cat_only"]: + isUpdated, pos_shear = obj.drawObj_multiband( + tel=tel, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=obj.g1, + g2=obj.g2, + exptime=pointing.exp_time, + fd_shear=fd_shear) + + elif chip.survey_type == "spectroscopic" and not self.overall_config["run_option"]["out_cat_only"]: + isUpdated, pos_shear = obj.drawObj_slitless( + tel=tel, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=obj.g1, + g2=obj.g2, + exptime=pointing.exp_time, + normFilter=norm_filt, + fd_shear=fd_shear) + + if isUpdated == 1: + # TODO: add up stats + self.chip_output.cat_add_obj(obj, pos_img, pos_shear) + pass + elif isUpdated == 0: + missed_obj += 1 + self.chip_output.Log_error("Objected missed: %s"%(obj.id)) + else: + self.chip_output.Log_error("Draw error, object omitted: %s"%(obj.id)) + continue + except Exception as e: + traceback.print_exc() + self.chip_output.Log_error(e) + + # Unload SED: + obj.unload_SED() + del obj + gc.collect() + del psf_model + gc.collect() + + self.chip_output.Log_info("Running checkpoint #1 (Object rendering finished): pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) )) + + self.chip_output.Log_info("# objects that are too bright %d out of %d"%(bright_obj, nobj)) + self.chip_output.Log_info("# objects that are too dim %d out of %d"%(dim_obj, nobj)) + self.chip_output.Log_info("# objects that are missed %d out of %d"%(missed_obj, nobj)) + + # Apply flat fielding (with shutter effects) + flat_normal = np.ones_like(chip.img.array) + if obs_param["flat_fielding"] == True: + flat_normal = flat_normal * chip.flat_img.array / np.mean(chip.flat_img.array) + if obs_param["shutter_effect"] == True: + flat_normal = flat_normal * chip.shutter_img + flat_normal = np.array(flat_normal, dtype='float32') + chip.img *= flat_normal + del flat_normal + + return chip, filt, tel, pointing + + def add_cosmic_rays(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info(msg=" Adding Cosmic-Ray", logger=self.logger) + chip.img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays( + img=chip.img, + chip=chip, + exptime=pointing.exptime, + seed=self.overall_config["random_seeds"]["seed_CR"]+pointing.id*30+chip.chipID) + # [TODO] output cosmic ray image + return chip, filt, tel, pointing + + def apply_PRNU(self, chip, filt, tel, pointing, catalog, obs_param): + chip.img *= chip.prnu_img + return chip, filt, tel, pointing + + def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param): + # Add dark current & Poisson noise + InputDark = False + if obs_param["add_dark"] == True: + if InputDark: + chip.img = chip_utils.add_inputdark(img=chip.img, chip=chip, exptime=pointing.exptime) + else: + chip.img, _ = chip_utils.add_poisson(img=chip.img, chip=chip, exptime=pointing.exptime, poisson_noise=chip.poisson_noise) + else: + chip.img, _ = chip_utils.add_poisson(img=chip.img, chip=self, exptime=pointing.exptime, poisson_noise=chip.poisson_noise, dark_noise=0.) + return chip, filt, tel, pointing + + def add_brighter_fatter(self, chip, filt, tel, pointing, catalog, obs_param): + chip.img = chip_utils.add_brighter_fatter(img=chip.img) + return chip, filt, tel, pointing + + def add_detector_defects(self, chip, filt, tel, pointing, catalog, obs_param): + # Add Hot Pixels or/and Dead Pixels + rgbadpix = Generator(PCG64(int(self.overall_config["random_seeds"]["seed_defective"]+chip.chipID))) + badfraction = 5E-5*(rgbadpix.random()*0.5+0.7) + chip.img = Effects.DefectivePixels( + chip.img, + IfHotPix=obs_param["hot_pixels"], + IfDeadPix=obs_param["dead_pixels"], + fraction=badfraction, + seed=self.overall_config["random_seeds"]["seed_defective"]+chip.chipID, biaslevel=0) + # Apply Bad columns + if obs_param["bad_columns"] == True: + chip.img = Effects.BadColumns(chip.img, seed=self.overall_config["random_seeds"]["seed_badcolumns"], chipid=chip.chipID) + return chip, filt, tel, pointing + + def add_nonlinearity(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info(" Applying Non-Linearity on the chip image") + chip.img = Effects.NonLinearity(GSImage=chip.img, beta1=5.e-7, beta2=0) + return chip, filt, tel, pointing + + def add_blooming(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info(" Applying CCD Saturation & Blooming") + chip.img = Effects.SaturBloom(GSImage=chip.img, nsect_x=1, nsect_y=1, fullwell=int(chip.full_well)) + return chip, filt, tel, pointing + + def apply_CTE(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info(" Apply CTE Effect") + ### 2*8 -> 1*16 img-layout + img = chip_utils.formatOutput(GSImage=chip.img) + chip.nsecy = 1 + chip.nsecx = 16 + + img_arr = img.array + ny, nx = img_arr.shape + dx = int(nx/chip.nsecx) + dy = int(ny/chip.nsecy) + newimg = galsim.Image(nx, int(ny+chip.overscan_y), init_value=0) + for ichannel in range(16): + print('\n***add CTI effects: pointing-{:} chip-{:} channel-{:}***'.format(pointing.id, chip.chipID, ichannel+1)) + noverscan, nsp, nmax = self.overscan_y, 3, 10 + beta, w, c = 0.478, 84700, 0 + t = np.array([0.74, 7.7, 37],dtype=np.float32) + rho_trap = np.array([0.6, 1.6, 1.4],dtype=np.float32) + trap_seeds = np.array([0, 1000, 10000],dtype=np.int32) + ichannel + chip.chipID*16 + release_seed = 50 + ichannel + pointing.id*30 + chip.chipID*16 + newimg.array[:, 0+ichannel*dx:dx+ichannel*dx] = CTI_sim(img_arr[:, 0+ichannel*dx:dx+ichannel*dx],dx,dy,noverscan,nsp,nmax,beta,w,c,t,rho_trap,trap_seeds,release_seed) + newimg.wcs = img.wcs + del img + img = newimg + + ### 1*16 -> 2*8 img-layout + chip.img = chip_utils.formatRevert(GSImage=img) + chip.nsecy = 2 + chip.nsecx = 8 + + # [TODO] make overscan_y == 0 + chip.overscan_y = 0 + return chip, filt, tel, pointing + + def add_prescan_overscan(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info("Apply pre/over-scan") + chip.img = chip_utils.AddPreScan(GSImage=chip.img, + pre1=chip.prescan_x, + pre2=chip.prescan_y, + over1=chip.overscan_x, + over2=chip.overscan_y) + return chip, filt, tel, pointing + + def add_bias(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info(" Adding Bias level and 16-channel non-uniformity") + if obs_param["bias_16channel"] == True: + chip.img = Effects.AddBiasNonUniform16(chip.img, + bias_level=float(chip.bias_level), + nsecy = chip.nsecy, nsecx=chip.nsecx, + seed=self.overall_config["random_seeds"]["seed_biasNonUniform"]+chip.chipID) + elif obs_param["bias_16channel"] == False: + chip.img += self.bias_level + return chip, filt, tel, pointing + + def add_readout_noise(self, chip, filt, tel, pointing, catalog, obs_param): + seed = int(self.overall_config["random_seeds"]["seed_readout"]) + pointing.id*30 + chip.chipID + rng_readout = galsim.BaseDeviate(seed) + readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=chip.read_noise) + chip.img.addNoise(readout_noise) + return chip, filt, tel, pointing + + def apply_gain(self, chip, filt, tel, pointing, catalog, obs_param): + self.chip_output.Log_info(" Applying Gain") + if obs_param["gain_16channel"] == True: + chip.img, chip.gain_channel = Effects.ApplyGainNonUniform16( + chip.img, gain=chip.gain, + nsecy = self.nsecy, nsecx=self.nsecx, + seed=self.overall_config["random_seeds"]["seed_gainNonUniform"]+chip.chipID) + elif obs_param["gain_16channel"] == False: + chip.img /= chip.gain + return chip, filt, tel, pointing + + def quantization_and_output(self, chip, filt, tel, pointing, catalog, obs_param): + chip.img.array[chip.img.array > 65535] = 65535 + chip.img.replaceNegative(replace_value=0) + chip.img.quantize() + chip.img = galsim.Image(chip.img.array, dtype=np.uint16) + hdu1 = fits.PrimaryHDU(header=self.h_prim) + hdu1.add_checksum() + hdu1.header.comments['CHECKSUM'] = 'HDU checksum' + hdu1.header.comments['DATASUM'] = 'data unit checksum' + hdu2 = fits.ImageHDU(chip.img.array, header=self.h_ext) + hdu2.add_checksum() + hdu2.header.comments['XTENSION'] = 'extension type' + hdu2.header.comments['CHECKSUM'] = 'HDU checksum' + hdu2.header.comments['DATASUM'] = 'data unit checksum' + 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 + +SIM_STEP_TYPES = { + "scie_obs": "add_objects", + "sky_background": "add_sky_background", + "cosmic_rays": "add_cosmic_rays", + "PRNU_effect": "apply_PRNU", + "poisson_and_dark": "add_poisson_and_dark", + "bright_fatter": "add_brighter_fatter", + "detector_defects": "add_detector_defects", + "nonlinearity": "add_nonlinearity", + "blooming": "add_blooming", + "CTE_effect": "apply_CTE", + "prescan_overscan": "add_prescan_overscan", + "bias": "add_bias", + "readout_noise": "add_readout_noise", + "gain": "apply_gain", + "quantization_and_output": "quantization_and_output" +} \ No newline at end of file diff --git a/ObservationSim/_util.py b/ObservationSim/_util.py index 67ac785..d7059fa 100755 --- a/ObservationSim/_util.py +++ b/ObservationSim/_util.py @@ -96,7 +96,7 @@ def make_run_dirs(work_dir, run_name, pointing_list): os.makedirs(imgDir, exist_ok=True) except OSError: pass - prefix = "MSC_" + # prefix = "MSC_" # for pointing in pointing_list: # fname=prefix + str(pointing.id).rjust(7, '0') # subImgDir = os.path.join(imgDir, fname) @@ -107,26 +107,26 @@ def make_run_dirs(work_dir, run_name, pointing_list): # pass return imgDir -def imgName(tt=0): - ut = datetime.utcnow() - eye, emo, eda, eho, emi, ese = str(ut.year), str(ut.month), str(ut.day), str(ut.hour), str(ut.minute), str(ut.second) - emse = str(ut.microsecond) - if int(emo)<10: emo = "0%s"%emo - if int(eda)<10: eda = "0%s"%eda - if int(eho)<10: eho = "0%s"%eho - if int(emi)<10: emi = "0%s"%emi - if int(ese)<10: ese = "0%s"%ese +# def imgName(tt=0): +# ut = datetime.utcnow() +# eye, emo, eda, eho, emi, ese = str(ut.year), str(ut.month), str(ut.day), str(ut.hour), str(ut.minute), str(ut.second) +# emse = str(ut.microsecond) +# if int(emo)<10: emo = "0%s"%emo +# if int(eda)<10: eda = "0%s"%eda +# if int(eho)<10: eho = "0%s"%eho +# if int(emi)<10: emi = "0%s"%emi +# if int(ese)<10: ese = "0%s"%ese - if tt==0: - namekey = "CSST%s%s%sT%s%s%s"%(eye,emo,eda,eho,emi,ese) - elif tt==1: - namekey = "%s-%s-%sT%s:%s:%s.%s"%(eye,emo,eda,eho,emi,ese,emse) - elif tt==2: - namekey = "%s%s%s%s%s%s"%(eye,emo,eda,eho,emi,ese) - else: - raise ValueError("!!! Give a right 'tt' value.") +# if tt==0: +# namekey = "CSST%s%s%sT%s%s%s"%(eye,emo,eda,eho,emi,ese) +# elif tt==1: +# namekey = "%s-%s-%sT%s:%s:%s.%s"%(eye,emo,eda,eho,emi,ese,emse) +# elif tt==2: +# namekey = "%s%s%s%s%s%s"%(eye,emo,eda,eho,emi,ese) +# else: +# raise ValueError("!!! Give a right 'tt' value.") - return namekey +# return namekey def makeSubDir_PointingList(path_dict, config, pointing_ID=0): imgDir = os.path.join(path_dict["work_dir"], config["run_name"]) diff --git a/config/config_overall.yaml b/config/config_overall.yaml index 5cf6fd8..b1d4d0d 100644 --- a/config/config_overall.yaml +++ b/config/config_overall.yaml @@ -73,9 +73,6 @@ obs_setting: # Whether to enable astrometric modeling enable_astrometric_model: True - # Whether to enable straylight model - enable_straylight_model: True - # Cut by saturation magnitude in which band? cut_in_band: "z" diff --git a/config/obs_config_SCI_WIDE.yaml b/config/obs_config_SCI_WIDE.yaml deleted file mode 100644 index 346f7f5..0000000 --- a/config/obs_config_SCI_WIDE.yaml +++ /dev/null @@ -1,25 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# For single exposure type: -# SCI-WIDE -# CSST-Sim Group, 2024/01/08 -# -############################################### - -# Observation type -obs_type: "SCIE" - -# Define simulation sequence -call_sequence: - # Sky background simulation - sky_background: - exptime: 150. # [s] - shutter_effect: YES - flat_fielding: YES - scie_obs: - exptime: 150. # [s] - shutter_effect: YES - flat_fielding: YES -... \ No newline at end of file diff --git a/config/obs_config_SCI_WIDE_phot.yaml b/config/obs_config_SCI_WIDE_phot.yaml new file mode 100644 index 0000000..2002fb4 --- /dev/null +++ b/config/obs_config_SCI_WIDE_phot.yaml @@ -0,0 +1,61 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# For single exposure type: +# SCI-WIDE +# CSST-Sim Group, 2024/01/08 +# +############################################### + +# Observation type +obs_type: "SCIE" + +# Define list of chips +run_chips: [6, 7, 8, 9, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 22, 23, 24, 25] + +# Define observation sequence +call_sequence: + # Accumulate fluxes from objects + scie_obs: + exptime: 150. # [s] + shutter_effect: YES + flat_fielding: YES + field_dist: YES + # Accumulate fluxes from sky background + sky_background: + exptime: 150. # [s] + shutter_effect: YES + flat_fielding: YES + enable_straylight_model: True + # Apply PRNU to accumulated photons + PRNU_effect: {} + # Accumulate photons caused by cosmic rays + cosmic_rays: {} + # Add Poission noise and dark current + poisson_and_dark: + add_dark: YES + # Simulate brighter fatter effects + bright_fatter: {} + # Add detector defects: hot/warm pixels, bad columns + detector_defects: + hot_pixels: YES + dead_pixels: YES + bad_columns: YES + # Apply response nonlinearity + nonlinearity: {} + # Apply CCD Saturation & Blooming + blooming: {} + # Run CTE simulation + CTE_effect: {} + # Add prescan and overscan + prescan_overscan: {} + # Add bias + bias: + bias_16channel: YES + # Add readout noise + readout_noise: {} + # Apply gain + gain: + gain_16channel: YES +... \ No newline at end of file diff --git a/run_sim.py b/run_sim.py index b0a3428..c9d1ade 100755 --- a/run_sim.py +++ b/run_sim.py @@ -50,10 +50,10 @@ def run_sim(): config['work_dir'] = args.work_dir # Some default values - if "bias_16channel" not in config["ins_effects"]: - config["ins_effects"]["bias_16channel"] = False - if "gain_16channel" not in config["ins_effects"]: - config["ins_effects"]["gain_16channel"] = False + # if "bias_16channel" not in config["ins_effects"]: + # config["ins_effects"]["bias_16channel"] = False + # if "gain_16channel" not in config["ins_effects"]: + # config["ins_effects"]["gain_16channel"] = False if "mag_sat_margin" not in config["obs_setting"]: config["obs_setting"]["mag_sat_margin"] = -2.5 if "mag_lim_margin" not in config["obs_setting"]: -- GitLab