Skip to content
GitLab
Commit 633fbbab authored by Fang Yuedong's avatar Fang Yuedong
Browse files

first test

parent 8d77fcd3
Hide whitespace changes
Inline Side-by-side
ObservationSim/Config/Header/ImageHeader.py
View file @ 633fbbab
......@@ -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')
......
ObservationSim/Config/Pointing.py
View file @ 633fbbab
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
ObservationSim/Config/_util.py
View file @ 633fbbab
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
......
ObservationSim/Instrument/Chip/ChipUtils.py
View file @ 633fbbab
......@@ -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)
......
ObservationSim/Instrument/FocalPlane.py
View file @ 633fbbab
......@@ -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)
......
ObservationSim/ObservationSim.py
View file @ 633fbbab
......@@ -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
ObservationSim/SimSteps.py 0 → 100644
View file @ 633fbbab
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
ObservationSim/_util.py
View file @ 633fbbab
......@@ -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"])
......
config/config_overall.yaml
View file @ 633fbbab
......@@ -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"
......
config/obs_config_SCI_WIDE.yaml → config/obs_config_SCI_WIDE_phot.yaml
View file @ 633fbbab
......@@ -11,15 +11,51 @@
# Observation type
obs_type: "SCIE"
# Define simulation sequence
# 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:
# Sky background simulation
sky_background:
# Accumulate fluxes from objects
scie_obs:
exptime: 150. # [s]
shutter_effect: YES
flat_fielding: YES
scie_obs:
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
run_sim.py
View file @ 633fbbab
......@@ -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"]:
......
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