Skip to content
GitLab
Commit 8df06b27 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

Merge branch 'sim_scheduler' into develop

parents 81e2570f 93270bbf
Hide whitespace changes
Inline Side-by-side
ObservationSim/ObservationSim.py
View file @ 8df06b27
......@@ -2,51 +2,50 @@ import os
import numpy as np
import mpi4py.MPI as MPI
import galsim
import logging
import psutil
import gc
from astropy.io import fits
from datetime import datetime
import traceback
from ObservationSim.Config import config_dir, ChipOutput
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Config import ChipOutput
from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
from ObservationSim.Instrument.Chip import Effects
from ObservationSim.Straylight import calculateSkyMap_split_g
from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp
from ObservationSim._util import get_shear_field, makeSubDir_PointingList
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
from ObservationSim.sim_steps 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, wcs_fp=None):
# 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_img, _ = chip_utils.get_flat(img=chip.img, seed=int(self.config["random_seeds"]["seed_flat"]))
if chip.chipID > 30:
chip.shutter_img = np.ones_like(chip.img.array)
else:
chip.shutter_img = 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):
......@@ -60,16 +59,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":
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)
......@@ -98,368 +87,112 @@ 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)
# 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
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 == 'MS':
# 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.exp_time, 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
# Prepare necessary chip properties for simulation
chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing)
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 = 40.)
chip_wcs = galsim.FitsWCS(header=h_ext)
for j in range(self.nobj):
# (DEBUG)
# if j >= 10:
# break
obj = self.cat.objs[j]
# 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()]))
# Initialize SimSteps
sim_steps = SimSteps(overall_config=self.config, chip_output=chip_output, all_filters=self.all_filters)
# 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()
for step in pointing.obs_param["call_sequence"]:
if self.config["run_option"]["out_cat_only"]:
if step != "scie_obs":
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.real_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:
# 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()
del psf_model
del self.cat
gc.collect()
chip_output.Log_info("Starting simulation step: %s, calling function: %s"%(step, SIM_STEP_TYPES[step]))
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=self.tel,
pointing=pointing,
catalog=self.Catalog,
obs_param=obs_param)
chip_output.Log_info("Finished simulation step: %s"%(step))
except Exception as e:
traceback.print_exc()
chip_output.Log_error(e)
chip_output.Log_error("Failed simulation on step: %s"%(step))
break
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) ))
# Detector Effects
# ===========================================================
# 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 pointing.pointing_type == 'MS':
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],
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=40.)
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) ))
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
pointing_ID = pointing.obs_id
pointing.make_output_pointing_dir(overall_config=self.config, copy_obs_config=True)
self.focal_plane = FocalPlane(chip_list=pointing.obs_param["run_chips"])
# Make Chip & Filter lists
self.chip_list = []
self.filter_list = []
self.all_filters = []
for i in range(self.focal_plane.nchips):
chipID = i + 1
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)
chip = run_chips[ichip]
filt = run_filts[ichip]
# chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid))
chip_output = ChipOutput(
config=self.config,
focal_plane=self.focal_plane,
chip=chip,
filt=filt,
exptime=pointing.exp_time,
pointing_type=pointing.pointing_type,
pointing_ID=pointing_ID,
subdir=sub_img_dir,
prefix=prefix)
chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid))
config = self.config,
chip = chip,
filt = filt,
pointing = pointing
)
chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(int(pointing_ID), chip.chipID, pid))
self.run_one_chip(
chip=chip,
filt=filt,
chip_output=chip_output,
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/PSF/PSFInterp.py
View file @ 8df06b27
......@@ -15,7 +15,6 @@ import h5py
from ObservationSim.PSF.PSFModel import PSFModel
LOG_DEBUG = False #***#
NPSF = 900 #***# 30*30
PixSizeInMicrons = 5. #***# in microns
......@@ -205,8 +204,9 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=Tru
###define PSFInterp###
class PSFInterp(PSFModel):
def __init__(self, chip, npsf=NPSF, PSF_data=None, PSF_data_file=None, PSF_data_prefix="", sigSpin=0, psfRa=0.15, HocBuild=False):
if LOG_DEBUG:
def __init__(self, chip, npsf=NPSF, PSF_data=None, PSF_data_file=None, PSF_data_prefix="", sigSpin=0, psfRa=0.15, HocBuild=False, LOG_DEBUG=False):
self.LOG_DEBUG = LOG_DEBUG
if self.LOG_DEBUG:
print('===================================================')
print('DEBUG: psf module for csstSim ' \
+time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
......@@ -225,7 +225,7 @@ class PSFInterp(PSFModel):
self.npsf = npsf
self.PSF_data = self._loadPSF(self.iccd, PSF_data_file, PSF_data_prefix)
if LOG_DEBUG:
if self.LOG_DEBUG:
print('nwave-{:} on ccd-{:}::'.format(self.nwave, self.iccd), flush=True)
print('self.PSF_data ... ok', flush=True)
print('Preparing self.[psfMat,cen_col,cen_row] for psfMaker ... ', end='', flush=True)
......@@ -252,7 +252,7 @@ class PSFInterp(PSFModel):
self.hoc.append(hoc)
self.hoclist.append(hoclist)
if LOG_DEBUG:
if self.LOG_DEBUG:
print('ok', flush=True)
......@@ -292,7 +292,7 @@ class PSFInterp(PSFModel):
psfSet.append(psfWave)
fq.close()
if LOG_DEBUG:
if self.LOG_DEBUG:
print('psfSet has been loaded:', flush=True)
print('psfSet[iwave][ipsf][keys]:', psfSet[0][0].keys(), flush=True)
return psfSet
......@@ -342,6 +342,7 @@ class PSFInterp(PSFModel):
assert(self.hoc != 0), 'hoclist should be built correctly!'
imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True)
'''
############TEST: START
TestGaussian = False
if TestGaussian:
......@@ -349,6 +350,7 @@ class PSFInterp(PSFModel):
#pointing_pa = -23.433333
imPSF= gsx.shear(g1=0.8, g2=0.).rotate(0.*galsim.degrees).drawImage(nx = 256, ny=256, scale=pixSize).array
############TEST: END
'''
if galsimGSObject:
imPSFt = np.zeros([257,257])
......@@ -370,6 +372,7 @@ class PSFInterp(PSFModel):
return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
return imPSF
'''
def PSFspin(self, x, y):
"""
The PSF profile at a given image position relative to the axis center
......@@ -392,7 +395,4 @@ class PSFInterp(PSFModel):
qr = np.sqrt((1.0+ell)/(1.0-ell))
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return self.psf.shear(PSFshear), PSFshear
if __name__ == '__main__':
pass
'''
ObservationSim/PSF/PSFInterpSLS.py 0 → 100644
View file @ 8df06b27
'''
PSF interpolation for CSST-Sim
NOTE: [iccd, iwave, ipsf] are counted from 1 to n, but [tccd, twave, tpsf] are counted from 0 to n-1
'''
import sys
import time
import copy
import numpy as np
import scipy.spatial as spatial
import galsim
import h5py
from ObservationSim.PSF.PSFModel import PSFModel
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
import os
from astropy.io import fits
from astropy.modeling.models import Gaussian2D
from scipy import signal
LOG_DEBUG = False #***#
NPSF = 900 #***# 30*30
PixSizeInMicrons = 5. #***# in microns
###find neighbors-KDtree###
# def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
# """
# find nearest neighbors by 2D-KDTree
#
# Parameters:
# tx, ty (float, float): a given position
# px, py (numpy.array, numpy.array): position data for tree
# dr (float-optional): distance
# dn (int-optional): nearest-N
# OnlyDistance (bool-optional): only use distance to find neighbors. Default: True
#
# Returns:
# dataq (numpy.array): index
# """
# datax = px
# datay = py
# tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
#
# dataq=[]
# rr = dr
# if OnlyDistance == True:
# dataq = tree.query_ball_point([tx, ty], rr)
# if OnlyDistance == False:
# while len(dataq) < dn:
# dataq = tree.query_ball_point([tx, ty], rr)
# rr += dr
# dd = np.hypot(datax[dataq]-tx, datay[dataq]-ty)
# ddSortindx = np.argsort(dd)
# dataq = np.array(dataq)[ddSortindx[0:dn]]
# return dataq
#
# ###find neighbors-hoclist###
# def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy):
# if np.max(partx) > nhocx*dhocx:
# print('ERROR')
# sys.exit()
# if np.max(party) > nhocy*dhocy:
# print('ERROR')
# sys.exit()
#
# npart = partx.size
# hoclist= np.zeros(npart, dtype=np.int32)-1
# hoc = np.zeros([nhocy, nhocx], dtype=np.int32)-1
# for ipart in range(npart):
# ix = int(partx[ipart]/dhocx)
# iy = int(party[ipart]/dhocy)
# hoclist[ipart] = hoc[iy, ix]
# hoc[iy, ix] = ipart
# return hoc, hoclist
#
# def hocFind(px, py, dhocx, dhocy, hoc, hoclist):
# ix = int(px/dhocx)
# iy = int(py/dhocy)
#
# neigh=[]
# it = hoc[iy, ix]
# while it != -1:
# neigh.append(it)
# it = hoclist[it]
# return neigh
#
# def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None):
# nhocy = nhocx = 20
#
# pxMin = np.min(px)
# pxMax = np.max(px)
# pyMin = np.min(py)
# pyMax = np.max(py)
#
# dhocx = (pxMax - pxMin)/(nhocx-1)
# dhocy = (pyMax - pyMin)/(nhocy-1)
# partx = px - pxMin +dhocx/2
# party = py - pyMin +dhocy/2
#
# if hoc is None:
# hoc, hoclist = hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy)
# return hoc, hoclist
#
# if hoc is not None:
# tx = tx - pxMin +dhocx/2
# ty = ty - pyMin +dhocy/2
# itx = int(tx/dhocx)
# ity = int(ty/dhocy)
#
# ps = [-1, 0, 1]
# neigh=[]
# for ii in range(3):
# for jj in range(3):
# ix = itx + ps[ii]
# iy = ity + ps[jj]
# if ix < 0:
# continue
# if iy < 0:
# continue
# if ix > nhocx-1:
# continue
# if iy > nhocy-1:
# continue
#
# #neightt = myUtil.hocFind(ppx, ppy, dhocx, dhocy, hoc, hoclist)
# it = hoc[iy, ix]
# while it != -1:
# neigh.append(it)
# it = hoclist[it]
# #neigh.append(neightt)
# #ll = [i for k in neigh for i in k]
# if dn != -1:
# ptx = np.array(partx[neigh])
# pty = np.array(party[neigh])
# dd = np.hypot(ptx-tx, pty-ty)
# idx = np.argsort(dd)
# neigh= np.array(neigh)[idx[0:dn]]
# return neigh
#
#
# ###PSF-IDW###
# def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False):
# """
# psf interpolation by IDW
#
# Parameters:
# px, py (float, float): position of the target
# PSFMat (numpy.array): image
# cen_col, cen_row (numpy.array, numpy.array): potions of the psf centers
# IDWindex (int-optional): the power index of IDW
# OnlyNeighbors (bool-optional): only neighbors are used for psf interpolation
#
# Returns:
# psfMaker (numpy.array)
# """
#
# minimum_psf_weight = 1e-8
# ref_col = px
# ref_row = py
#
# ngy, ngx = PSFMat[0, :, :].shape
# npsf = PSFMat[:, :, :].shape[0]
# psfWeight = np.zeros([npsf])
#
# if OnlyNeighbors == True:
# if hoc is None:
# neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=4, OnlyDistance=False)
# if hoc is not None:
# neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist)
#
# neighFlag = np.zeros(npsf)
# neighFlag[neigh] = 1
#
# for ipsf in range(npsf):
# if OnlyNeighbors == True:
# if neighFlag[ipsf] != 1:
# continue
#
# dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
# if IDWindex == 1:
# psfWeight[ipsf] = dist
# if IDWindex == 2:
# psfWeight[ipsf] = dist**2
# if IDWindex == 3:
# psfWeight[ipsf] = dist**3
# if IDWindex == 4:
# psfWeight[ipsf] = dist**4
# psfWeight[ipsf] = max(psfWeight[ipsf], minimum_psf_weight)
# psfWeight[ipsf] = 1./psfWeight[ipsf]
# psfWeight /= np.sum(psfWeight)
#
# psfMaker = np.zeros([ngy, ngx], dtype=np.float32)
# for ipsf in range(npsf):
# if OnlyNeighbors == True:
# if neighFlag[ipsf] != 1:
# continue
#
# iPSFMat = PSFMat[ipsf, :, :].copy()
# ipsfWeight = psfWeight[ipsf]
#
# psfMaker += iPSFMat * ipsfWeight
# psfMaker /= np.nansum(psfMaker)
#
# return psfMaker
###define PSFInterp###
class PSFInterpSLS(PSFModel):
def __init__(self, chip, filt,PSF_data_prefix="", sigSpin=0, psfRa=0.15, pix_size = 0.005):
if LOG_DEBUG:
print('===================================================')
print('DEBUG: psf module for csstSim ' \
+time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
print('===================================================')
self.sigSpin = sigSpin
self.sigGauss = psfRa
self.grating_ids = chip_utils.getChipSLSGratingID(chip.chipID)
_,self.grating_type = chip.getChipFilter(chipID=chip.chipID)
self.data_folder = PSF_data_prefix
self.getPSFDataFromFile(filt)
self.pixsize = pix_size # um
def getPSFDataFromFile(self, filt):
gratingInwavelist = {'GU':0,'GV':1,'GI':2}
grating_orders = ['0','1']
waveListFn = self.data_folder + '/wavelist.dat'
wavelists = np.loadtxt(waveListFn)
self.waveList = wavelists[:,gratingInwavelist[self.grating_type]]
bandranges = np.zeros([4,2])
midBand = (self.waveList[0:3] + self.waveList[1:4])/2.*10000.
bandranges[0,0] = filt.blue_limit
bandranges[1:4,0] = midBand
bandranges[0:3, 1] = midBand
bandranges[3,1] = filt.red_limit
self.bandranges = bandranges
self.grating1_data = {}
g_folder = self.data_folder + '/' + self.grating_ids[0] + '/'
for g_order in grating_orders:
g_folder_order = g_folder + 'PSF_Order_' + g_order + '/'
grating_order_data = {}
for bandi in [1,2,3,4]:
subBand_data = {}
subBand_data['bandrange'] = bandranges[bandi-1]
final_folder = g_folder_order + str(bandi) + '/'
print(final_folder)
pca_fs = os.listdir(final_folder)
for fname in pca_fs:
if ('_PCs.fits' in fname) and (fname[0] != '.'):
fname_ = final_folder + fname
hdu = fits.open(fname_)
subBand_data['band_data'] = hdu
grating_order_data['band'+str(bandi)] = subBand_data
self.grating1_data['order'+g_order] = grating_order_data
self.grating2_data = {}
g_folder = self.data_folder + '/' + self.grating_ids[1] + '/'
for g_order in grating_orders:
g_folder_order = g_folder + 'PSF_Order_' + g_order + '/'
grating_order_data = {}
for bandi in [1, 2, 3, 4]:
subBand_data = {}
subBand_data['bandrange'] = bandranges[bandi - 1]
final_folder = g_folder_order + str(bandi) + '/'
print(final_folder)
pca_fs = os.listdir(final_folder)
for fname in pca_fs:
if ('_PCs.fits' in fname) and (fname[0] != '.'):
fname_ = final_folder + fname
hdu = fits.open(fname_)
subBand_data['band_data'] = hdu
grating_order_data['band' + str(bandi)] = subBand_data
self.grating2_data['order' + g_order] = grating_order_data
#
#
#
# def _getPSFwave(self, iccd, PSF_data_file, PSF_data_prefix):
# # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r')
# fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r')
# nwave = len(fq.keys())
# fq.close()
# return nwave
#
#
# def _loadPSF(self, iccd, PSF_data_file, PSF_data_prefix):
# psfSet = []
# # fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r')
# fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r')
# for ii in range(self.nwave):
# iwave = ii+1
# psfWave = []
#
# fq_iwave = fq['w_{:}'.format(iwave)]
# for jj in range(self.npsf):
# ipsf = jj+1
# psfInfo = {}
# psfInfo['wavelength']= fq_iwave['wavelength'][()]
#
# fq_iwave_ipsf = fq_iwave['psf_{:}'.format(ipsf)]
# psfInfo['pixsize'] = PixSizeInMicrons
# psfInfo['field_x'] = fq_iwave_ipsf['field_x'][()]
# psfInfo['field_y'] = fq_iwave_ipsf['field_y'][()]
# psfInfo['image_x'] = fq_iwave_ipsf['image_x'][()]
# psfInfo['image_y'] = fq_iwave_ipsf['image_y'][()]
# psfInfo['centroid_x']= fq_iwave_ipsf['cx'][()]
# psfInfo['centroid_y']= fq_iwave_ipsf['cy'][()]
# psfInfo['psfMat'] = fq_iwave_ipsf['psfMat'][()]
#
# psfWave.append(psfInfo)
# psfSet.append(psfWave)
# fq.close()
#
# if LOG_DEBUG:
# print('psfSet has been loaded:', flush=True)
# print('psfSet[iwave][ipsf][keys]:', psfSet[0][0].keys(), flush=True)
# return psfSet
#
#
# def _findWave(self, bandpass):
# if isinstance(bandpass,int):
# twave = bandpass
# return twave
#
# for twave in range(self.nwave):
# bandwave = self.PSF_data[twave][0]['wavelength']
# if bandpass.blue_limit < bandwave and bandwave < bandpass.red_limit:
# return twave
# return -1
#
#
def convolveWithGauss(self, img=None, sigma=1):
offset = int(np.ceil(sigma * 3))
g_size = 2 * offset + 1
m_cen = int(g_size / 2)
print('-----',g_size)
g_PSF_ = Gaussian2D(1, m_cen, m_cen, sigma, sigma)
yp, xp = np.mgrid[0:g_size, 0:g_size]
g_PSF = g_PSF_(xp, yp)
psf = g_PSF / g_PSF.sum()
convImg = signal.fftconvolve(img, psf, mode='full', axes=None)
convImg = convImg/np.sum(convImg)
return convImg
def get_PSF(self, chip, pos_img_local = [1000,1000], bandNo = 1, galsimGSObject=True, folding_threshold=5.e-3, g_order = 'A', grating_split_pos=3685):
"""
Get the PSF at a given image position
Parameters:
chip: A 'Chip' object representing the chip we want to extract PSF from.
pos_img: A 'galsim.Position' object representing the image position.
bandpass: A 'galsim.Bandpass' object representing the wavelength range.
pixSize: The pixels size of psf matrix
findNeighMode: 'treeFind' or 'hoclistFind'
Returns:
PSF: A 'galsim.GSObject'.
"""
order_IDs = {'A': '1', 'B': '0' ,'C': '0', 'D': '0', 'E': '0'}
contam_order_sigma = {'C':0.28032344707964174,'D':0.39900182912061344,'E':1.1988309797685412} #arcsec
x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
# print(pos_img.x - x_start)
pos_img_x = pos_img_local[0] + x_start
pos_img_y = pos_img_local[1] + y_start
pos_img = galsim.PositionD(pos_img_x, pos_img_y)
if pos_img_local[0] < grating_split_pos:
psf_data = self.grating1_data
else:
psf_data = self.grating2_data
grating_order = order_IDs[g_order]
# if grating_order in ['-2','-1','2']:
# grating_order = '1'
# if grating_order in ['0', '1']:
psf_order = psf_data['order'+grating_order]
psf_order_b = psf_order['band'+str(bandNo)]
psf_b_dat = psf_order_b['band_data']
pos_p = psf_b_dat[1].data
pc_coeff = psf_b_dat[2].data
pcs = psf_b_dat[0].data
# print(max(pos_p[:,0]), min(pos_p[:,0]),max(pos_p[:,1]), min(pos_p[:,1]))
# print(chip.x_cen, chip.y_cen)
# print(pos_p)
px = pos_img.x*chip.pix_size
py = pos_img.y*chip.pix_size
dist2=(pos_p[:,1] - px)*(pos_p[:,1] - px) + (pos_p[:,0] - py)*(pos_p[:,0] - py)
temp_sort_dist = np.zeros([dist2.shape[0],2])
temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0],1)
temp_sort_dist[:, 1] = dist2
# print(temp_sort_dist)
dits2_sortlist = sorted(temp_sort_dist, key=lambda x:x[1])
# print(dits2_sortlist)
nearest4p = np.zeros([4,2])
pc_coeff_4p = np.zeros([pc_coeff.data.shape[0],4])
for i in np.arange(4):
smaller_ids = int(dits2_sortlist[i][0])
nearest4p[i, 0] = pos_p[smaller_ids, 1]
nearest4p[i, 1] = pos_p[smaller_ids, 0]
pc_coeff_4p[:,i] = pc_coeff[:,smaller_ids]
idw_dist = 1/(np.sqrt((px-nearest4p[:,0]) * (px-nearest4p[:,0]) + (py-nearest4p[:,1]) * (py-nearest4p[:,1])))
coeff_int = np.zeros(pc_coeff.data.shape[0])
for i in np.arange(4):
coeff_int = coeff_int + pc_coeff_4p[:,i]*idw_dist[i]
coeff_int = coeff_int / np.sum(coeff_int)
npc = 10
m_size = int(pcs.shape[0]**0.5)
PSF_int = np.dot(pcs[:,0:npc],coeff_int[0:npc]).reshape(m_size,m_size)
# PSF_int = PSF_int/np.sum(PSF_int)
PSF_int_trans = np.flipud(np.fliplr(PSF_int))
PSF_int_trans = np.fliplr(PSF_int_trans.T)
# PSF_int_trans = np.abs(PSF_int_trans)
# ids_szero = PSF_int_trans<0
# PSF_int_trans[ids_szero] = 0
# print(PSF_int_trans[ids_szero].shape[0],PSF_int_trans.shape)
PSF_int_trans = PSF_int_trans/np.sum(PSF_int_trans)
# from astropy.io import fits
# fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)
# if g_order in ['C','D','E']:
# g_simgma = contam_order_sigma[g_order]/pixel_size_arc
# PSF_int_trans = self.convolveWithGauss(PSF_int_trans,g_simgma)
# n_m_size = int(m_size/2)
#
# n_PSF_int = np.zeros([n_m_size, n_m_size])
#
# for i in np.arange(n_m_size):
# for j in np.arange(n_m_size):
# n_PSF_int[i,j] = np.sum(PSF_int[2*i:2*i+2, 2*j:2*j+2])
#
# n_PSF_int = n_PSF_int/np.sum(n_PSF_int)
# chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
# chip.img.wcs = galsim.wcs.AffineTransform
if galsimGSObject:
# imPSFt = np.zeros([257,257])
# imPSFt[0:256, 0:256] = imPSF
# # imPSFt[120:130, 0:256] = 1.
pixel_size_arc = np.rad2deg(self.pixsize * 1e-3 / 28) * 3600
img = galsim.ImageF(PSF_int_trans, scale=pixel_size_arc)
gsp = galsim.GSParams(folding_threshold=folding_threshold)
############TEST: START
# Use sheared PSF to test the PSF orientation
# self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.)
############TEST: END
self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
# if g_order in ['C','D','E']:
# add_psf = galsim.Gaussian(sigma=contam_order_sigma[g_order], flux=1.0)
# self.psf = galsim.Convolve(self.psf, add_psf)
wcs = chip.img.wcs.local(pos_img)
scale = galsim.PixelScale(0.074)
self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img))
# return self.PSFspin(x=px/0.01, y=py/0.01)
return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
return PSF_int_trans, PSF_int
# pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize
#
# # assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID'
# twave = self._findWave(bandpass)
# if twave == -1:
# print("!!!PSF bandpass does not match.")
# exit()
# PSFMat = self.psfMat[twave]
# cen_col= self.cen_col[twave]
# cen_row= self.cen_row[twave]
#
# px = (pos_img.x - chip.cen_pix_x)*0.01
# py = (pos_img.y - chip.cen_pix_y)*0.01
# if findNeighMode == 'treeFind':
# imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, PSFCentroidWgt=True)
# if findNeighMode == 'hoclistFind':
# assert(self.hoc != 0), 'hoclist should be built correctly!'
# imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True)
#
# ############TEST: START
# TestGaussian = False
# if TestGaussian:
# gsx = galsim.Gaussian(sigma=0.04)
# #pointing_pa = -23.433333
# imPSF= gsx.shear(g1=0.8, g2=0.).rotate(0.*galsim.degrees).drawImage(nx = 256, ny=256, scale=pixSize).array
# ############TEST: END
#
# if galsimGSObject:
# imPSFt = np.zeros([257,257])
# imPSFt[0:256, 0:256] = imPSF
# # imPSFt[120:130, 0:256] = 1.
#
# img = galsim.ImageF(imPSFt, scale=pixSize)
# gsp = galsim.GSParams(folding_threshold=folding_threshold)
# ############TEST: START
# # Use sheared PSF to test the PSF orientation
# # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.)
# ############TEST: END
# self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
# wcs = chip.img.wcs.local(pos_img)
# scale = galsim.PixelScale(0.074)
# self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img))
#
# # return self.PSFspin(x=px/0.01, y=py/0.01)
# return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
# return imPSF
#
# def PSFspin(self, x, y):
# """
# The PSF profile at a given image position relative to the axis center
#
# Parameters:
# theta : spin angles in a given exposure in unit of [arcsecond]
# dx, dy: relative position to the axis center in unit of [pixels]
#
# Return:
# Spinned PSF: g1, g2 and axis ratio 'a/b'
# """
# a2Rad = np.pi/(60.0*60.0*180.0)
#
# ff = self.sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
# rc = np.sqrt(x*x + y*y)
# cpix = rc*(self.sigSpin*a2Rad)
#
# beta = (np.arctan2(y,x) + np.pi/2)
# ell = cpix**2/(2.0*ff**2+cpix**2)
# qr = np.sqrt((1.0+ell)/(1.0-ell))
# PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
# return self.psf.shear(PSFshear), PSFshear
from ObservationSim.Instrument import Filter, FilterParam, Chip
import yaml
if __name__ == '__main__':
configfn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/config/config_C6_dev.yaml'
with open(configfn, "r") as stream:
try:
config = yaml.safe_load(stream)
for key, value in config.items():
print (key + " : " + str(value))
except yaml.YAMLError as exc:
print(exc)
chip = Chip(chipID=1,config=config)
filter_id, filter_type = chip.getChipFilter()
filt = Filter(filter_id=filter_id,
filter_type=filter_type,
filter_param=FilterParam())
psf_i = PSFInterpSLS(chip, filt,PSF_data_prefix="/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/")
pos_img = galsim.PositionD(x=25155, y=-22060)
psf_im = psf_i.get_PSF(chip, pos_img = pos_img, g_order = '1')
ObservationSim/PSF/__init__.py
View file @ 8df06b27
......@@ -2,4 +2,5 @@ from .PSFModel import PSFModel
from .PSFGauss import PSFGauss
# from .PSFInterp.PSFInterp import PSFInterp
from .PSFInterp import PSFInterp
from .PSFInterpSLS import PSFInterpSLS
from .FieldDistortion import FieldDistortion
\ No newline at end of file
ObservationSim/_util.py
View file @ 8df06b27
......@@ -13,7 +13,7 @@ def parse_args():
'''
parser = argparse.ArgumentParser()
parser.add_argument('--config_file', type=str, required=True, help='.yaml config file for simulation settings.')
parser.add_argument('--catalog', type=str, required=True, help='name of the catalog interface class to be loaded.')
parser.add_argument('--catalog', type=str, help='name of the catalog interface class to be loaded.')
parser.add_argument('-c', '--config_dir', type=str, help='Directory that houses the .yaml config file.')
parser.add_argument('-d', '--data_dir', type=str, help='Directory that houses the input data.')
parser.add_argument('-w', '--work_dir', type=str, help='The path for output.')
......@@ -29,47 +29,45 @@ def generate_pointing_list(config, pointing_filename=None, data_dir=None):
# Calculate starting time(s) for CAL exposures
# NOTE: temporary implementation
t = datetime.timestamp(t0)
ncal = config['obs_setting']['np_cal']
ipoint = 0
for i in range(ncal):
pointing = Pointing(
id = ipoint,
ra=config["obs_setting"]["ra_center"],
dec=config["obs_setting"]["dec_center"],
img_pa=config["obs_setting"]["image_rot"],
timestamp=t,
pointing_type='CAL')
t += 3 * delta_t * 60. # 3 calibration exposures for each pointing
pointing_list.append(pointing)
ipoint += 1
run_pointings = config['obs_setting']['run_pointings']
if config["obs_setting"]["exp_time"]:
exp_time = config["obs_setting"]["exp_time"]
if "obs_config_file" in config['obs_setting']:
obs_config_file = config['obs_setting']["obs_config_file"]
else:
exp_time = 150.
obs_config_file = None
if pointing_filename and data_dir:
pointing_file = os.path.join(data_dir, pointing_filename)
# if pointing_filename and data_dir:
if pointing_filename:
if data_dir:
pointing_file = os.path.join(data_dir, pointing_filename)
else:
pointing_file = pointing_filename
f = open(pointing_file, 'r')
for _ in range(1):
header = f.readline()
# for _ in range(1):
# header = f.readline()
iline = 0
for line in f:
if len(line.strip()) == 0 or line[0]=='#':
continue
if run_pointings and iline not in run_pointings:
iline += 1
ipoint += 1
continue
line = line.strip()
columns = line.split()
pointing = Pointing(exp_time=exp_time)
pointing.read_pointing_columns(columns=columns, id=ipoint, t=t)
pointing = Pointing(obs_config_file=obs_config_file)
pointing.read_pointing_columns(columns=columns, id=ipoint)
t += delta_t * 60.
pointing_list.append(pointing)
iline += 1
ipoint += 1
f.close()
else:
if config["obs_setting"]["exp_time"]:
exp_time = config["obs_setting"]["exp_time"]
else:
exp_time = 150.
pointing = Pointing(
id=ipoint,
ra=config["obs_setting"]["ra_center"],
......@@ -77,7 +75,8 @@ def generate_pointing_list(config, pointing_filename=None, data_dir=None):
img_pa=config["obs_setting"]["image_rot"],
timestamp=t,
exp_time=exp_time,
pointing_type='MS'
pointing_type='SCI',
obs_config_file=obs_config_file
)
t += delta_t * 60.
pointing_list.append(pointing)
......@@ -96,46 +95,16 @@ def make_run_dirs(work_dir, run_name, pointing_list):
os.makedirs(imgDir, exist_ok=True)
except OSError:
pass
prefix = "MSC_"
for pointing in pointing_list:
fname=prefix + str(pointing.id).rjust(7, '0')
subImgDir = os.path.join(imgDir, fname)
if not os.path.exists(subImgDir):
try:
os.makedirs(subImgDir, exist_ok=True)
except OSError:
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
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
def makeSubDir_PointingList(path_dict, config, pointing_ID=0):
def make_output_pointing_dir(path_dict, config, pointing_ID=0):
imgDir = os.path.join(path_dict["work_dir"], config["run_name"])
if not os.path.exists(imgDir):
try:
os.makedirs(imgDir, exist_ok=True)
except OSError:
pass
prefix = "MSC_" + str(pointing_ID).rjust(7, '0')
prefix = "MSC_" + str(pointing_ID).rjust(8, '0')
subImgdir = os.path.join(imgDir, prefix)
if not os.path.exists(subImgdir):
try:
......
ObservationSim/sim_steps/__init__.py 0 → 100644
View file @ 8df06b27
import os
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
from .prepare_headers import prepare_headers, updateHeaderInfo
from .add_sky_background import add_sky_background_sci, add_sky_flat_calibration, add_sky_background
from .add_objects import add_objects
from .add_cosmic_rays import add_cosmic_rays
from .add_pattern_noise import apply_PRNU, add_poisson_and_dark, add_detector_defects, add_nonlinearity, add_blooming, add_bias
from .add_brighter_fatter_CTE import add_brighter_fatter, apply_CTE
from .readout_output import add_prescan_overscan, add_readout_noise, apply_gain, quantization_and_output
from .add_LED_flat import add_LED_Flat
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",
"led_calib_model":"add_LED_Flat",
"sky_flatField":"add_sky_flat_calibration",
}
\ No newline at end of file
ObservationSim/sim_steps/add_LED_flat.py 0 → 100644
View file @ 8df06b27
import numpy as np
from ObservationSim.MockObject import FlatLED
import galsim
from astropy.time import Time
from datetime import datetime, timezone
import gc
def add_LED_Flat(self, chip, filt, tel, pointing, catalog, obs_param):
if not hasattr(self, 'h_ext'):
_, _ = self.prepare_headers(chip=chip, pointing=pointing)
chip_wcs = galsim.FitsWCS(header = self.h_ext)
pf_map = np.zeros_like(chip.img.array)
if obs_param["LED_TYPE"] is not None:
if len(obs_param["LED_TYPE"]) != 0:
print("LED OPEN--------")
led_obj = FlatLED(chip, filt)
led_flat, ledstat, letts = led_obj.drawObj_LEDFlat(led_type_list=obs_param["LED_TYPE"], exp_t_list=obs_param["LED_TIME"])
pf_map = led_flat
self.updateHeaderInfo(header_flag='ext', keys = ['LEDSTAT'], values = [ledstat])
self.updateHeaderInfo(header_flag='ext', keys = ['LEDT01','LEDT02','LEDT03','LEDT04','LEDT05','LEDT06','LEDT07','LEDT08','LEDT09','LEDT10','LEDT11','LEDT12','LEDT13','LEDT14'], values = letts)
if obs_param["shutter_effect"] == True:
pf_map = pf_map * chip.shutter_img
pf_map = np.array(pf_map, dtype='float32')
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
else:
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [True,'','','',''])
chip.img = chip.img + pf_map
# renew header info
datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
t_obs = Time(datetime_obs)
##ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光
t_obs_renew = Time(t_obs.mjd - (2.) / 86400., format="mjd")
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time :
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [pointing.exp_time])
gc.collect()
return chip, filt, tel, pointing
\ No newline at end of file
ObservationSim/sim_steps/add_brighter_fatter_CTE.py 0 → 100644
View file @ 8df06b27
import numpy as np
import galsim
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from ObservationSim.Instrument.Chip.libCTI.CTI_modeling import CTI_sim
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 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 = chip.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
\ No newline at end of file
ObservationSim/sim_steps/add_cosmic_rays.py 0 → 100644
View file @ 8df06b27
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
def add_cosmic_rays(self, chip, filt, tel, pointing, catalog, obs_param):
self.chip_output.Log_info(" Adding Cosmic-Ray")
# Get exposure time
if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
exptime = obs_param["exptime"]
else:
exptime = pointing.exp_time
chip.img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(
img=chip.img,
chip=chip,
exptime=exptime,
seed=self.overall_config["random_seeds"]["seed_CR"]+pointing.id*30+chip.chipID)
# Save cosmic ray image
if (obs_param) and ("save_cosmic_img" in obs_param) and (obs_param["save_cosmic_img"] is not None):
if obs_param["save_cosmic_img"]:
chip_utils.output_fits_image(
chip=chip,
pointing=pointing,
img=crmap_gsimg,
output_dir=self.chip_output.subdir,
img_type='CRS',
img_type_code=pointing.pointing_type_code,
project_cycle=self.overall_config["project_cycle"],
run_counter=self.overall_config["run_counter"]
)
return chip, filt, tel, pointing
\ No newline at end of file
ObservationSim/sim_steps/add_objects.py 0 → 100644
View file @ 8df06b27
import os
import gc
import psutil
import traceback
import numpy as np
import galsim
from ObservationSim._util import get_shear_field
from ObservationSim.PSF import PSFGauss, FieldDistortion, PSFInterp, PSFInterpSLS
from astropy.time import Time
from datetime import datetime, timezone
def add_objects(self, chip, filt, tel, pointing, catalog, obs_param):
# Get exposure time
if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
exptime = obs_param["exptime"]
else:
exptime = pointing.exp_time
# Load catalogues
if catalog is None:
self.chip_output.Log_error("Catalog interface class must be specified for SCIE-OBS")
raise ValueError("Catalog interface class must be specified for SCIE-OBS")
cat = catalog(config=self.overall_config, chip=chip, pointing=pointing, chip_output=self.chip_output, filt=filt)
# Prepare output file(s) for this chip
# [NOTE] Headers of output .cat file may be updated by Catalog instance
# this should be called after the creation of Catalog instance
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(cat.objs)
missed_obj = 0
bright_obj = 0
dim_obj = 0
for j in range(nobj):
# # [DEBUG] [TODO]
# if j >= 10:
# break
obj = cat.objs[j]
# load and convert SED; also caculate object's magnitude in all CSST bands
try:
sed_data = cat.load_sed(obj)
norm_filt = cat.load_norm_filt(obj)
obj.sed, obj.param["mag_%s"%filt.filter_type.lower()], obj.param["flux_%s"%filt.filter_type.lower()] = 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()] = 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()
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=exptime,
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=exptime,
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')
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
else:
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [True,'','','',''])
chip.img *= flat_normal
del flat_normal
# renew header info
datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
t_obs = Time(datetime_obs)
##ccd刷新2s,等待0.s，开始曝光
t_obs_renew = Time(t_obs.mjd - (2.+0.) / 86400., format="mjd")
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time : 曝光时间+刷新后等带时间0.s+关快门后读出前等待0.s
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.+0.+pointing.exp_time])
return chip, filt, tel, pointing
\ No newline at end of file
ObservationSim/sim_steps/add_pattern_noise.py 0 → 100644
View file @ 8df06b27
from numpy.random import Generator, PCG64
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from ObservationSim.Instrument.Chip import Effects
def apply_PRNU(self, chip, filt, tel, pointing, catalog, obs_param):
chip.img *= chip.prnu_img
if self.overall_config["output_setting"]["prnu_output"] == True:
chip.prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (self.chip_output.subdir, str(chip.chipID).rjust(2, '0')))
return chip, filt, tel, pointing
def add_poisson_and_dark(self, chip, filt, tel, pointing, catalog, obs_param):
# Add dark current & Poisson noise
# Get exposure time
if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
exptime = obs_param["exptime"]
else:
exptime = pointing.exp_time
if obs_param["add_dark"] == True:
chip.img, _ = chip_utils.add_poisson(img=chip.img,
chip=chip,
exptime=pointing.exp_time,
poisson_noise=chip.poisson_noise,
InputDark=None)
else:
chip.img, _ = chip_utils.add_poisson(img=chip.img,
chip=self,
exptime=exptime,
poisson_noise=chip.poisson_noise,
dark_noise=0.)
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 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
ObservationSim/sim_steps/add_sky_background.py 0 → 100644
View file @ 8df06b27
import numpy as np
import galsim
from ObservationSim.Straylight import calculateSkyMap_split_g
from ObservationSim.Instrument import FilterParam
from astropy.time import Time
from datetime import datetime, timezone
def add_sky_background_sci(self, chip, filt, tel, pointing, catalog, obs_param):
# Get exposure time
if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
exptime = obs_param["exptime"]
else:
exptime = pointing.exp_time
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')
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
else:
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [False,'','','',''])
if obs_param["enable_straylight_model"]:
# Filter.sky_background, Filter.zodical_spec will be updated
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 = 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)*exptime
# sky_map = sky_map * tel.pupil_area * obs_param["exptime"]
chip.img += sky_map
return chip, filt, tel, pointing
def add_sky_flat_calibration(self, chip, filt, tel, pointing, catalog, obs_param):
if not hasattr(self, 'h_ext'):
_, _ = self.prepare_headers(chip=chip, pointing=pointing)
chip_wcs = galsim.FitsWCS(header = self.h_ext)
# Get exposure time
if (obs_param) and ("exptime" in obs_param) and (obs_param["exptime"] is not None):
exptime = obs_param["exptime"]
else:
exptime = pointing.exp_time
skyback_level = obs_param["flat_level"]
filter_param = FilterParam()
sky_level_filt = obs_param["flat_level_filt"]
norm_scaler = skyback_level/exptime / filter_param.param[sky_level_filt][5]
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 self.overall_config["output_setting"]["shutter_output"] == True: # output 16-bit shutter effect image with pixel value <=65535
shutt_gsimg = galsim.ImageUS(chip.shutter_img*6E4)
shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (self.chip_output.subdir, str(chip.chipID).rjust(2, '0')))
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT'], values = [True])
else:
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1'], values = [True,'','','',''])
if chip.survey_type == "photometric":
sky_map = flat_normal * np.ones_like(chip.img.array) * norm_scaler * filter_param.param[chip.filter_type][5] / tel.pupil_area * exptime
elif chip.survey_type == "spectroscopic":
# 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')
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)
sky_map = sky_map * norm_scaler * exptime
chip.img += sky_map
return chip, filt, tel, pointing
def add_sky_background(self, chip, filt, tel, pointing, catalog, obs_param):
if not hasattr(self, 'h_ext'):
_, _ = self.prepare_headers(chip=chip, pointing=pointing)
chip_wcs = galsim.FitsWCS(header = self.h_ext)
if "flat_level" not in obs_param or "flat_level_filt" not in obs_param:
chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
else:
if obs_param.get('flat_level') is None or obs_param.get('flat_level_filt')is None:
chip, filt, tel, pointing = self.add_sky_background_sci(chip, filt, tel, pointing, catalog, obs_param)
else:
chip, filt, tel, pointing = self.add_sky_flat_calibration(chip, filt, tel, pointing, catalog, obs_param)
# renew header info
datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
t_obs = Time(datetime_obs)
##ccd刷新2s,等待0.5s，开始曝光
t_obs_renew = Time(t_obs.mjd - (2.+0.5) / 86400., format="mjd")
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
#dark time : 曝光时间+刷新后等带时间0.5s+关闭快门时间1.5s+管快门后读出前等待0.5s
self.updateHeaderInfo(header_flag='ext', keys = ['DARKTIME'], values = [0.5+1.5+0.5+pointing.exp_time])
return chip, filt, tel, pointing
\ No newline at end of file
ObservationSim/sim_steps/prepare_headers.py 0 → 100644
View file @ 8df06b27
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
def prepare_headers(self, chip, pointing):
self.h_prim = generatePrimaryHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
pointing_id = pointing.obs_id,
pointing_type_code = pointing.pointing_type_code,
ra=pointing.ra,
dec=pointing.dec,
pixel_scale=chip.pix_scale,
time_pt = pointing.timestamp,
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,
t_shutter_open = pointing.t_shutter_open,
t_shutter_close = pointing.t_shutter_close)
return self.h_prim, self.h_ext
def updateHeaderInfo(self,header_flag='prim', keys = ['key'], values = [0]):
if header_flag == 'prim':
for key,value in zip(keys, values):
self.h_prim[key] = value
if header_flag == 'ext':
for key,value in zip(keys, values):
self.h_ext[key] = value
ObservationSim/sim_steps/readout_output.py 0 → 100644
View file @ 8df06b27
import os
import galsim
import numpy as np
from astropy.io import fits
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from ObservationSim.Instrument.Chip import Effects
from astropy.time import Time
from datetime import datetime, timezone
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)
if obs_param["add_dark"] == True:
ny = int(chip.npix_y/2)
base_dark = (ny-1)*(chip.readout_time/ny)*chip.dark_noise
chip.img.array[(chip.prescan_y+ny):-(chip.prescan_y+ny),:] = base_dark
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 = chip.nsecy,
nsecx=chip.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):
if not hasattr(self, 'h_ext'):
_, _ = self.prepare_headers(chip=chip, pointing=pointing)
self.updateHeaderInfo(header_flag='ext', keys = ['SHTSTAT','SHTOPEN0','SHTOPEN1','SHTCLOS0','SHTCLOS1','EXPTIME'], values = [False,'','','','',0.0])
# renew header info
datetime_obs = datetime.utcfromtimestamp(pointing.timestamp)
datetime_obs = datetime_obs.replace(tzinfo=timezone.utc)
t_obs = Time(datetime_obs)
##ccd刷新2s,等待0.5s，开灯后等待0.5s,开始曝光
t_obs_renew = Time(t_obs.mjd - 2. / 86400., format="mjd")
t_obs_utc = datetime.utcfromtimestamp(np.round(datetime.utcfromtimestamp(t_obs_renew.unix).replace(tzinfo=timezone.utc).timestamp(), 1))
self.updateHeaderInfo(header_flag='prim', keys = ['DATE-OBS'], values = [t_obs_utc.strftime("%Y-%m-%dT%H:%M:%S.%f")[:-5]])
gains1 = list(chip.gain_channel[0:8])
gains2 = list(chip.gain_channel[8:])
gains2.reverse()
gains = np.append(gains1,gains2)
self.updateHeaderInfo(header_flag='ext', keys = ['GAIN01','GAIN02','GAIN03','GAIN04','GAIN05','GAIN06','GAIN07','GAIN08','GAIN09','GAIN10','GAIN11','GAIN12','GAIN13','GAIN14','GAIN15','GAIN16'], values = gains)
if obs_param["format_output"] == True:
self.chip_output.Log_info(" Apply 1*16 format")
chip.img = chip_utils.formatOutput(GSImage=chip.img)
chip.nsecy = 1
chip.nsecx = 16
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)
fname = os.path.join(self.chip_output.subdir, self.h_prim['FILENAME'] + '.fits')
f_name_size = 68
if(len(self.h_prim['FILENAME'])>f_name_size):
self.updateHeaderInfo(header_flag='prim', keys = ['FILENAME'], values = [self.h_prim['FILENAME'][0:f_name_size]])
hdu1 = fits.PrimaryHDU(header=self.h_prim)
hdu1.add_checksum()
hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
hdu1.header.comments['DATASUM'] = 'data unit checksum'
self.updateHeaderInfo(header_flag='ext', keys = ['DATASECT'], values = [str(chip.img.array.shape[1])+'x'+str(chip.img.array.shape[0])])
hdu2 = fits.ImageHDU(chip.img.array, header=self.h_ext)
hdu2.add_checksum()
hdu2.header.comments['XTENSION'] = 'extension type'
hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
hdu2.header.comments['DATASUM'] = 'data unit checksum'
hdu2.header.comments["XTENSION"] = "image extension"
hdu1 = fits.HDUList([hdu1, hdu2])
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
return chip, filt, tel, pointing
README.md
View file @ 8df06b27
# CSST主巡天仿真软件
## 重要更新或问题修复：
* v3.0版本相关内容：（待补充）
* 2023.07.29: 更新至v2.1版本，内容包括：
* 加入杂散光模块，config文件中添加开关：enable_straylight_model: True/False
* 调整WCS定义：x：E, y：-N
......@@ -22,12 +23,11 @@
## 使用方法及相关说明
* 软件安装和使用方法可参考：[软件说明文档](https://kdocs.cn/l/ckrtHibdV3OS)。目前还较为粗糙，持续更新中
* 软件安装和使用方法可参考：[软件说明文档](https://kdocs.cn/l/cjyiU0SXGyn2)。目前还较为粗糙，持续更新中
<!---
* 仿真数据说明文档：[C6数据说明文档](https://www.kdocs.cn/l/cuRmPf71Ly7v)
-->
* 在线版使用说明：[CSST网页接口使用说明](https://kdocs.cn/l/cl81flAXe1YQ)
-->
* 软件安装及问题反馈表：[问题反馈表](https://f.kdocs.cn/w/XYcBkE63/)
## 相关数据
......
config/config_50sqdeg.yaml
View file @ 8df06b27
......@@ -9,13 +9,13 @@
# Base diretories and naming setup
# Can add some of the command-line arguments here as well;
# OK to pass either way or both, as long as they are consistent
work_dir: "/share/home/fangyuedong/csst-simulation/workplace/"
work_dir: "/share/home/weichengliang/CSST_git/test_new_sim/outputs/"
data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
run_name: "rotate_0"
run_name: "testRun0"
# Whether to use MPI
run_option:
use_mpi: NO
use_mpi: YES
# NOTE: "n_threads" paramters is currently not used in the backend
# simulation codes. It should be implemented later in the web frontend
# in order to config the number of threads to request from NAOC cluster
......@@ -45,10 +45,10 @@ catalog_options:
# AGN_SED_WAVE: "wave_ross13.npy"
# Only simulate stars?
star_only: NO
star_only: YES
# Only simulate galaxies?
galaxy_only: YES
galaxy_only: NO
# rotate galaxy ellipticity
rotateEll: 0. # [degree]
......@@ -85,7 +85,7 @@ obs_setting:
# if you just want to run default pointing:
# - pointing_dir: null
# - pointing_file: null
pointing_dir: "/share/home/fangyuedong/50sqdeg_pointings/"
pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/"
pointing_file: "pointing_50_combined.dat"
# Number of calibration pointings
......@@ -163,6 +163,7 @@ ins_effects:
add_dark: ON # Whether to add dark noise
add_readout: ON # Whether to add read-out (Gaussian) noise
add_bias: ON # Whether to add bias-level to images
add_prescan: OFF
bias_16channel: ON # Whether to add different biases for 16 channels
gain_16channel: ON # Whether to make different gains for 16 channels
shutter_effect: ON # Whether to add shutter effect
......@@ -171,12 +172,13 @@ ins_effects:
non_linear: ON # Whether to add non-linearity
cosmic_ray: ON # Whether to add cosmic-ray
cray_differ: ON # Whether to generate different cosmic ray maps CAL and MS output
cte_trail: ON # Whether to simulate CTE trails
cte_trail: OFF # Whether to simulate CTE trails
saturbloom: ON # Whether to simulate Saturation & Blooming
add_badcolumns: ON # Whether to add bad columns
add_hotpixels: ON # Whether to add hot pixels
add_deadpixels: ON # Whether to add dead(dark) pixels
bright_fatter: ON # Whether to simulate Brighter-Fatter (also diffusion) effect
format_output: OFF
# Values:
# default values have been defined individually for each chip in:
......@@ -218,4 +220,4 @@ random_seeds:
seed_badcolumns: 20240309 # Seed for bad columns
seed_defective: 20210304 # Seed for defective (bad) pixels
seed_readout: 20210601 # Seed for read-out gaussian noise
...
\ No newline at end of file
...
config/config_C6.yaml
View file @ 8df06b27
......@@ -9,9 +9,12 @@
# Base diretories and naming setup
# Can add some of the command-line arguments here as well;
# OK to pass either way or both, as long as they are consistent
work_dir: "/share/home/fangyuedong/csst-simulation/workplace/"
work_dir: "/share/home/fangyuedong/new_sim/workplace/"
# work_dir: "/share/C6_new_sim_2sq"
data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
run_name: "profile_C6"
run_name: "C6_new_sim_2sq_run1"
project_cycle: 6
run_counter: 1
# Whether to use MPI
run_option:
......@@ -44,7 +47,7 @@ catalog_options:
AGN_SED_WAVE: "wave_ross13.npy"
# Only simulate stars?
star_only: NO
star_only: YES
# Only simulate galaxies?
galaxy_only: NO
......@@ -112,7 +115,8 @@ obs_setting:
cut_in_band: "z"
# saturation magnitude margin
mag_sat_margin: -2.5
# mag_sat_margin: -2.5
mag_sat_margin: -15.
# limiting magnitude margin
mag_lim_margin: +1.0
......@@ -135,7 +139,7 @@ psf_setting:
# path to PSF data
# NOTE: only valid for "Interp" PSF
psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"
psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/"
###############################################
# Shear setting
###############################################
......@@ -217,4 +221,4 @@ random_seeds:
seed_badcolumns: 20240309 # Seed for bad columns
seed_defective: 20210304 # Seed for defective (bad) pixels
seed_readout: 20210601 # Seed for read-out gaussian noise
...
\ No newline at end of file
...
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