Loading Catalog/NGPCatalog.py +28 −12 Original line number Diff line number Diff line Loading @@ -23,15 +23,13 @@ except ImportError: NSIDE = 128 class NGPCatalog(CatalogBase): def __init__(self, config, chip, pointing, **kwargs): def __init__(self, config, chip, pointing, chip_output, **kwargs): super().__init__() self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"]) self.seed_Av = config["random_seeds"]["seed_Av"] if "logger" in kwargs: self.logger = kwargs["logger"] else: self.logger = None self.chip_output = chip_output self.logger = chip_output.logger with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: self.normF_star = Table.read(str(filter_path)) Loading @@ -58,9 +56,17 @@ class NGPCatalog(CatalogBase): else: self.rotation = 0. # Update output .cat header with catalog specific output columns self._add_output_columns_header() self._get_healpix_list() self._load() def _add_output_columns_header(self): self.add_hdr = " model_tag teff logg feh" self.add_fmt = " %10s %8.4f %8.4f %8.4f" self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr) def _get_healpix_list(self): self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax Loading Loading @@ -133,10 +139,8 @@ class NGPCatalog(CatalogBase): param['dec'] = dec_arr[igals] param['ra_orig'] = gals['ra_true'][igals] param['dec_orig'] = gals['dec_true'][igals] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] if param['mag_use_normal'] >= 24: if param['mag_use_normal'] >= 26.5: continue param['z'] = gals['redshift_true'][igals] param['model_tag'] = 'None' Loading Loading @@ -173,15 +177,22 @@ class NGPCatalog(CatalogBase): param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im param['star'] = 2 # Quasar # NOTE: this cut cannot be put before the SED type has been assigned if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue self.ids += 1 # param['id'] = self.ids param['id'] = gals['galaxyID'][igals] if param['star'] == 0: obj = Galaxy(param, self.rotation, logger=self.logger) self.objs.append(obj) if param['star'] == 2: obj = Quasar(param, logger=self.logger) # Need to deal with additional output columns obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.) self.objs.append(obj) def _load_stars(self, stars, pix_id=None): Loading Loading @@ -234,7 +245,7 @@ class NGPCatalog(CatalogBase): if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = stars['app_sdss_g'][istars] if param['mag_use_normal'] >= 26: if param['mag_use_normal'] >= 26.5: continue self.ids += 1 # param['id'] = self.ids Loading @@ -247,6 +258,10 @@ class NGPCatalog(CatalogBase): param['z'] = 0.0 param['star'] = 1 # Star obj = Star(param, logger=self.logger) # Append additional output columns to the .cat file obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh']) self.objs.append(obj) def _load(self, **kwargs): Loading Loading @@ -302,7 +317,8 @@ class NGPCatalog(CatalogBase): raise ValueError("Object type not known") speci = interpolate.interp1d(wave, flux) # lamb = np.arange(2500, 10001 + 0.5, 0.5) lamb = np.arange(2400, 11001 + 0.5, 0.5) # lamb = np.arange(2400, 11001 + 0.5, 0.5) lamb = np.arange(2000, 18001 + 0.5, 0.5) y = speci(lamb) # erg/s/cm2/A --> photo/s/m2/A all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 Loading ObservationSim/Config/ChipOutput.py +20 −55 Original line number Diff line number Diff line Loading @@ -11,6 +11,7 @@ class ChipOutput(object): self.imgKey2 = imgKey2 self.exptime = exptime self.mjdTime = mjdTime self.pointing_type = pointing_type if (ra_cen is not None) and (dec_cen is not None): self.ra_cen = ra_cen self.dec_cen = dec_cen Loading @@ -23,95 +24,59 @@ class ChipOutput(object): # self.cat_name = 'MSC_' + config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat" self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat" self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat" self.subdir = subdir # Setup logger for each chip logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log" logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log" self.logger = logging.getLogger() fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8') fh.setLevel(logging.DEBUG) self.logger.setLevel(logging.DEBUG) formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') logging.getLogger('numba').setLevel(logging.WARNING) # formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s') fh.setFormatter(formatter) self.logger.addHandler(fh) hdr1 = "# obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type " hdr2 = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge g1 g2 " hdr3 = "sed_type av redden " hdr4 = "pm_ra pm_dec RV parallax\n" hdr4 = "pm_ra pm_dec RV parallax" fmt1 = "%10d %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s " fmt2 = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f " fmt3 = "%2d %8.4f %8.4f " fmt4 = "%15.8f %15.8f %15.8f %15.8f\n" fmt4 = "%15.8f %15.8f %15.8f %15.8f" self.hdr = hdr1 + hdr2 + hdr3 + hdr4 self.fmt = fmt1 + fmt2 + fmt3 + fmt4 self.logger.info("pointing_type = %s\n"%(pointing_type)) if pointing_type == 'MS': def update_ouptut_header(self, additional_column_names=""): self.hdr += additional_column_names def create_output_file(self): if self.pointing_type == 'MS': self.cat = open(os.path.join(self.subdir, self.cat_name), "w") self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name))) if not self.hdr.endswith("\n"): self.hdr += "\n" self.cat.write(self.hdr) # def updateHDR(self, hdr): # hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise, "comment":"read noise in e-/pixel"}, # {"name":"DARK", "value":self.chip.dark_noise, "comment":"Dark noise (e-/pixel/s)"}, # {"name":"EXPTIME", "value":self.exptime, "comment":"exposure time in second"}, # {"name":"GAIN", "value":self.chip.gain, "comment":"CCD gain in e-/ADU"}, # {"name":"SATURATE","value":65535.0, "comment":"saturation level"}, # {"name":"CCDCHIP", "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"}, # {"name":"FILTER", "value":self.filt.filter_type, "comment":"filter name"}, # {"name":"MJD-OBS", "value":self.mjdTime, "comment":"Modified Julian Date (MJD) of observation"}, # {"name":"DATE-OBS","value":self.imgKey1, "comment":"Date of observation"}, # {"name":"EQUINOX", "value":2000.0}, # {"name":"RADECSYS","value":"ICRS"}, # {"name":"RA", "value":self.ra_cen, "comment":"telescope pointing center"}, # {"name":"DEC", "value":self.dec_cen, "comment":"telescope pointing center"}, # {"name":"OBJECT", "value":"CSS-OS"}, # {"name":"WCSDIM", "value":2.0, "comment":"WCS Dimensionality"}, # {"name":"EXTNAME", "value":"IM1", "comment":"Extension name"}, # {"name":"BSCALE", "value":1.0}, # {"name":"BZERO", "value":0.0}, # {"name":"OBSID", "value":self.imgKey0, "comment":"Observation ID"}, # {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"}, # {"name":"RSPEED", "value":10.0, "comment":"Read speed"}, # {"name":"CHIPTEMP","value":-100.0, "comment":"Chip temperature"}, # {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"}, # {"name":"CCDSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, # {"name":"NSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, # {"name":"AUTHOR", "value":"CSST-Sim Group"}, # {"name":"GROUP", "value":"Weak Lensing Working Group for CSST"}] # for item in hdrNew: # hdr.add_record(item) # return hdr # def cat_add_obj(self, obj, pos_img, snr, pos_shear, g1, g2): def cat_add_obj(self, obj, pos_img, pos_shear): ximg = pos_img.x - self.chip.bound.xmin + 1.0 yimg = pos_img.y - self.chip.bound.ymin + 1.0 # if obj.type == 'galaxy': # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, # obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0) # elif obj.type == "quasar": # line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, # obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0) # else: # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh']) # print( # obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, # obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2, # obj.sed_type, obj.av, obj.redden, # obj.pmra, obj.pmdec, obj.rv, obj.parallax) line = self.fmt%( obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2, obj.sed_type, obj.av, obj.redden, obj.pmra, obj.pmdec, obj.rv, obj.parallax) line += obj.additional_output_str if not line.endswith("\n"): line += "\n" self.cat.write(line) ObservationSim/Instrument/Chip/Chip.py +16 −14 Original line number Diff line number Diff line Loading @@ -33,6 +33,7 @@ class Chip(FocalPlane): self.dark_exptime = float(config["ins_effects"]['dark_exptime']) self.flat_exptime = float(config["ins_effects"]['flat_exptime']) self.readout_time = float(config["ins_effects"]['readout_time']) self.full_well = int(config["ins_effects"]["full_well"]) self.logger = logger self.direct_Img_sls = direct_Img_sls Loading Loading @@ -96,16 +97,17 @@ class Chip(FocalPlane): if filter_type in ['g', 'r', 'GV']: filename = 'Astro_MB.txt' if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt' # Mirror efficiency: if filter_type == 'nuv': mirror_eff = 0.54 if filter_type == 'u': mirror_eff = 0.68 if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8 if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right # if filter_type == 'nuv': mirror_eff = 0.54 # if filter_type == 'u': mirror_eff = 0.68 # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8 # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right # path = os.path.join(self.ccdEffCurve_dir, filename) # table = Table.read(path, format='ascii') with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path: table = Table.read(ccd_path, format='ascii') throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear') # throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear') throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear') bandpass = galsim.Bandpass(throughput, wave_type='nm') return bandpass Loading Loading @@ -303,7 +305,6 @@ class Chip(FocalPlane): hdu1.writeto(fname, output_verify='ignore', overwrite=True) def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None, logger=None): SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"]) SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"]) SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"]) Loading Loading @@ -340,7 +341,7 @@ class Chip(FocalPlane): img.array[img.array > fullwell] = fullwell img = img - sky_map - self.dark_noise * (self.exptime + 0.5 * self.readout_time) return img return img/self.exptime if config["ins_effects"]["add_hotpixels"] == True: Loading Loading @@ -533,7 +534,7 @@ class Chip(FocalPlane): # Output images for calibration pointing ###################################################################################### # Bias output if config["output_setting"]["bias_output"] == True and pointing_type=='CAL': if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type=='CAL': if self.logger is not None: self.logger.info(" Output N frame Bias files") else: Loading @@ -546,6 +547,7 @@ class Chip(FocalPlane): ncombine=1, read_noise=self.read_noise, gain=1, seed=SeedBiasNonuni+self.chipID, logger=self.logger) # Readout noise for Biases is not generated with random seeds. So readout noise for bias images can't be reproduced. if config["ins_effects"]["cosmic_ray"] == True: if config["ins_effects"]["cray_differ"] == True: cr_map, cr_event_num = effects.produceCR_Map( Loading Loading @@ -601,7 +603,7 @@ class Chip(FocalPlane): del BiasCombImg # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file if config["output_setting"]["flat_output"] == True and pointing_type=='CAL': if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type=='CAL': if self.logger is not None: self.logger.info(" Output N frame Flat-Field files") else: Loading Loading @@ -673,7 +675,7 @@ class Chip(FocalPlane): # Add Read-out Noise if config["ins_effects"]["add_readout"] == True: seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 3 rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) FlatCombImg.addNoise(readout_noise) Loading Loading @@ -710,7 +712,7 @@ class Chip(FocalPlane): del flat_img # Export Dark current images if config["output_setting"]["dark_output"] == True and pointing_type=='CAL': if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type=='CAL': if self.logger is not None: self.logger.info(" Output N frame Dark Current files") else: Loading Loading @@ -795,7 +797,7 @@ class Chip(FocalPlane): # Add Read-out Noise if config["ins_effects"]["add_readout"] == True: seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 2 rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) DarkCombImg.addNoise(readout_noise) Loading ObservationSim/Instrument/Chip/Effects.py +1 −1 Original line number Diff line number Diff line Loading @@ -197,7 +197,7 @@ def MakeFlatSmooth(GSBounds, seed): Flty, Fltx = np.mgrid[ymin:(ymax+1), xmin:(xmax+1)] rg = Generator(PCG64(int(seed))) p1,p2,bg=rg.poisson(1000, 3) Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20 Fltz = 0.6*1e-7*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20 FlatImg = galsim.ImageF(Fltz) return FlatImg Loading ObservationSim/Instrument/Filter.py +26 −6 Original line number Diff line number Diff line Loading @@ -3,7 +3,7 @@ import pylab as pl import os import numpy as np from ObservationSim.Instrument._util import photonEnergy from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag from ObservationSim.Instrument.FilterParam import FilterParam try: Loading @@ -13,7 +13,7 @@ except ImportError: import importlib_resources as pkg_resources class Filter(object): def __init__(self, filter_id, filter_type, filter_param, ccd_bandpass): def __init__(self, filter_id, filter_type, filter_param, ccd_bandpass=None): self.filter_id = filter_id self.filter_type = filter_type self.ccd_bandpass = ccd_bandpass Loading @@ -36,14 +36,15 @@ class Filter(object): self.efficiency = filter_param.param[filter_type][4] self.sky_background = filter_param.param[filter_type][5] self.mag_saturation = filter_param.param[filter_type][6] self.mag_dim = filter_param.param[filter_type][7] self.mag_limiting = filter_param.param[filter_type][7] # self.filter_dir = filter_param.filter_dir def is_too_bright(self, mag): return mag <= self.mag_saturation - 1.0 return mag <= self.mag_saturation - 2.5 # return mag <= 14.0 def is_too_dim(self, mag): return mag >= self.mag_dim + 1.0 return mag >= self.mag_limiting + 1.0 def _get_bandpasses(self, filter_dir=None, unit='A'): if self.filter_id < 7: # Photometric Loading @@ -51,8 +52,10 @@ class Filter(object): # filter_file = os.path.join(filter_dir, self.filter_type+".dat") # bandpass_full = galsim.Bandpass(filter_file, wave_type=unit) with pkg_resources.path('ObservationSim.Instrument.data.filters', self.filter_type.lower() + '.txt') as filter_file: self.filter_bandpass = galsim.Bandpass(str(filter_file), wave_type=unit) with pkg_resources.path('ObservationSim.Instrument.data.throughputs', self.filter_type.lower() + '_throughput.txt') as filter_file: bandpass_full = galsim.Bandpass(str(filter_file), wave_type=unit) bandpass_full = bandpass_full * self.ccd_bandpass # bandpass_full = bandpass_full * self.ccd_bandpass # Get sub-bandpasses bandpass_sub_list = [] Loading Loading @@ -95,3 +98,20 @@ class Filter(object): def getSkyNoise(self, exptime, gain=1.): return self.sky_background * exptime / gain def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None): if self.filter_type in ["GI", "GV", "GU"]: return if chip is not None: pix_scale = chip.pix_scale read_noise = chip.read_noise dark_noise = chip.dark_noise full_well = chip.full_well else: pix_scale = 0.074 read_noise = 5.0 dark_noise = 0.02 full_well = 90000 throughput_file = self.filter_type.lower() + '_throughput.txt' self.mag_limiting, self.mag_saturation = calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well) Loading
Catalog/NGPCatalog.py +28 −12 Original line number Diff line number Diff line Loading @@ -23,15 +23,13 @@ except ImportError: NSIDE = 128 class NGPCatalog(CatalogBase): def __init__(self, config, chip, pointing, **kwargs): def __init__(self, config, chip, pointing, chip_output, **kwargs): super().__init__() self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"]) self.seed_Av = config["random_seeds"]["seed_Av"] if "logger" in kwargs: self.logger = kwargs["logger"] else: self.logger = None self.chip_output = chip_output self.logger = chip_output.logger with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: self.normF_star = Table.read(str(filter_path)) Loading @@ -58,9 +56,17 @@ class NGPCatalog(CatalogBase): else: self.rotation = 0. # Update output .cat header with catalog specific output columns self._add_output_columns_header() self._get_healpix_list() self._load() def _add_output_columns_header(self): self.add_hdr = " model_tag teff logg feh" self.add_fmt = " %10s %8.4f %8.4f %8.4f" self.chip_output.update_ouptut_header(additional_column_names=self.add_hdr) def _get_healpix_list(self): self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax Loading Loading @@ -133,10 +139,8 @@ class NGPCatalog(CatalogBase): param['dec'] = dec_arr[igals] param['ra_orig'] = gals['ra_true'][igals] param['dec_orig'] = gals['dec_true'][igals] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] if param['mag_use_normal'] >= 24: if param['mag_use_normal'] >= 26.5: continue param['z'] = gals['redshift_true'][igals] param['model_tag'] = 'None' Loading Loading @@ -173,15 +177,22 @@ class NGPCatalog(CatalogBase): param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im param['star'] = 2 # Quasar # NOTE: this cut cannot be put before the SED type has been assigned if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue self.ids += 1 # param['id'] = self.ids param['id'] = gals['galaxyID'][igals] if param['star'] == 0: obj = Galaxy(param, self.rotation, logger=self.logger) self.objs.append(obj) if param['star'] == 2: obj = Quasar(param, logger=self.logger) # Need to deal with additional output columns obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.) self.objs.append(obj) def _load_stars(self, stars, pix_id=None): Loading Loading @@ -234,7 +245,7 @@ class NGPCatalog(CatalogBase): if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = stars['app_sdss_g'][istars] if param['mag_use_normal'] >= 26: if param['mag_use_normal'] >= 26.5: continue self.ids += 1 # param['id'] = self.ids Loading @@ -247,6 +258,10 @@ class NGPCatalog(CatalogBase): param['z'] = 0.0 param['star'] = 1 # Star obj = Star(param, logger=self.logger) # Append additional output columns to the .cat file obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh']) self.objs.append(obj) def _load(self, **kwargs): Loading Loading @@ -302,7 +317,8 @@ class NGPCatalog(CatalogBase): raise ValueError("Object type not known") speci = interpolate.interp1d(wave, flux) # lamb = np.arange(2500, 10001 + 0.5, 0.5) lamb = np.arange(2400, 11001 + 0.5, 0.5) # lamb = np.arange(2400, 11001 + 0.5, 0.5) lamb = np.arange(2000, 18001 + 0.5, 0.5) y = speci(lamb) # erg/s/cm2/A --> photo/s/m2/A all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 Loading
ObservationSim/Config/ChipOutput.py +20 −55 Original line number Diff line number Diff line Loading @@ -11,6 +11,7 @@ class ChipOutput(object): self.imgKey2 = imgKey2 self.exptime = exptime self.mjdTime = mjdTime self.pointing_type = pointing_type if (ra_cen is not None) and (dec_cen is not None): self.ra_cen = ra_cen self.dec_cen = dec_cen Loading @@ -23,95 +24,59 @@ class ChipOutput(object): # self.cat_name = 'MSC_' + config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat" self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat" self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat" self.subdir = subdir # Setup logger for each chip logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(7, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log" logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log" self.logger = logging.getLogger() fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8') fh.setLevel(logging.DEBUG) self.logger.setLevel(logging.DEBUG) formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') logging.getLogger('numba').setLevel(logging.WARNING) # formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s') formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s') fh.setFormatter(formatter) self.logger.addHandler(fh) hdr1 = "# obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type " hdr2 = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge g1 g2 " hdr3 = "sed_type av redden " hdr4 = "pm_ra pm_dec RV parallax\n" hdr4 = "pm_ra pm_dec RV parallax" fmt1 = "%10d %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s " fmt2 = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f " fmt3 = "%2d %8.4f %8.4f " fmt4 = "%15.8f %15.8f %15.8f %15.8f\n" fmt4 = "%15.8f %15.8f %15.8f %15.8f" self.hdr = hdr1 + hdr2 + hdr3 + hdr4 self.fmt = fmt1 + fmt2 + fmt3 + fmt4 self.logger.info("pointing_type = %s\n"%(pointing_type)) if pointing_type == 'MS': def update_ouptut_header(self, additional_column_names=""): self.hdr += additional_column_names def create_output_file(self): if self.pointing_type == 'MS': self.cat = open(os.path.join(self.subdir, self.cat_name), "w") self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name))) if not self.hdr.endswith("\n"): self.hdr += "\n" self.cat.write(self.hdr) # def updateHDR(self, hdr): # hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise, "comment":"read noise in e-/pixel"}, # {"name":"DARK", "value":self.chip.dark_noise, "comment":"Dark noise (e-/pixel/s)"}, # {"name":"EXPTIME", "value":self.exptime, "comment":"exposure time in second"}, # {"name":"GAIN", "value":self.chip.gain, "comment":"CCD gain in e-/ADU"}, # {"name":"SATURATE","value":65535.0, "comment":"saturation level"}, # {"name":"CCDCHIP", "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"}, # {"name":"FILTER", "value":self.filt.filter_type, "comment":"filter name"}, # {"name":"MJD-OBS", "value":self.mjdTime, "comment":"Modified Julian Date (MJD) of observation"}, # {"name":"DATE-OBS","value":self.imgKey1, "comment":"Date of observation"}, # {"name":"EQUINOX", "value":2000.0}, # {"name":"RADECSYS","value":"ICRS"}, # {"name":"RA", "value":self.ra_cen, "comment":"telescope pointing center"}, # {"name":"DEC", "value":self.dec_cen, "comment":"telescope pointing center"}, # {"name":"OBJECT", "value":"CSS-OS"}, # {"name":"WCSDIM", "value":2.0, "comment":"WCS Dimensionality"}, # {"name":"EXTNAME", "value":"IM1", "comment":"Extension name"}, # {"name":"BSCALE", "value":1.0}, # {"name":"BZERO", "value":0.0}, # {"name":"OBSID", "value":self.imgKey0, "comment":"Observation ID"}, # {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"}, # {"name":"RSPEED", "value":10.0, "comment":"Read speed"}, # {"name":"CHIPTEMP","value":-100.0, "comment":"Chip temperature"}, # {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"}, # {"name":"CCDSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, # {"name":"NSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, # {"name":"AUTHOR", "value":"CSST-Sim Group"}, # {"name":"GROUP", "value":"Weak Lensing Working Group for CSST"}] # for item in hdrNew: # hdr.add_record(item) # return hdr # def cat_add_obj(self, obj, pos_img, snr, pos_shear, g1, g2): def cat_add_obj(self, obj, pos_img, pos_shear): ximg = pos_img.x - self.chip.bound.xmin + 1.0 yimg = pos_img.y - self.chip.bound.ymin + 1.0 # if obj.type == 'galaxy': # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, # obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0) # elif obj.type == "quasar": # line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, # obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0) # else: # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh']) # print( # obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, # obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2, # obj.sed_type, obj.av, obj.redden, # obj.pmra, obj.pmdec, obj.rv, obj.parallax) line = self.fmt%( obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2, obj.sed_type, obj.av, obj.redden, obj.pmra, obj.pmdec, obj.rv, obj.parallax) line += obj.additional_output_str if not line.endswith("\n"): line += "\n" self.cat.write(line)
ObservationSim/Instrument/Chip/Chip.py +16 −14 Original line number Diff line number Diff line Loading @@ -33,6 +33,7 @@ class Chip(FocalPlane): self.dark_exptime = float(config["ins_effects"]['dark_exptime']) self.flat_exptime = float(config["ins_effects"]['flat_exptime']) self.readout_time = float(config["ins_effects"]['readout_time']) self.full_well = int(config["ins_effects"]["full_well"]) self.logger = logger self.direct_Img_sls = direct_Img_sls Loading Loading @@ -96,16 +97,17 @@ class Chip(FocalPlane): if filter_type in ['g', 'r', 'GV']: filename = 'Astro_MB.txt' if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt' # Mirror efficiency: if filter_type == 'nuv': mirror_eff = 0.54 if filter_type == 'u': mirror_eff = 0.68 if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8 if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right # if filter_type == 'nuv': mirror_eff = 0.54 # if filter_type == 'u': mirror_eff = 0.68 # if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8 # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right # path = os.path.join(self.ccdEffCurve_dir, filename) # table = Table.read(path, format='ascii') with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path: table = Table.read(ccd_path, format='ascii') throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear') # throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear') throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear') bandpass = galsim.Bandpass(throughput, wave_type='nm') return bandpass Loading Loading @@ -303,7 +305,6 @@ class Chip(FocalPlane): hdu1.writeto(fname, output_verify='ignore', overwrite=True) def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None, logger=None): SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"]) SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"]) SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"]) Loading Loading @@ -340,7 +341,7 @@ class Chip(FocalPlane): img.array[img.array > fullwell] = fullwell img = img - sky_map - self.dark_noise * (self.exptime + 0.5 * self.readout_time) return img return img/self.exptime if config["ins_effects"]["add_hotpixels"] == True: Loading Loading @@ -533,7 +534,7 @@ class Chip(FocalPlane): # Output images for calibration pointing ###################################################################################### # Bias output if config["output_setting"]["bias_output"] == True and pointing_type=='CAL': if config["ins_effects"]["add_bias"] == True and config["output_setting"]["bias_output"] == True and pointing_type=='CAL': if self.logger is not None: self.logger.info(" Output N frame Bias files") else: Loading @@ -546,6 +547,7 @@ class Chip(FocalPlane): ncombine=1, read_noise=self.read_noise, gain=1, seed=SeedBiasNonuni+self.chipID, logger=self.logger) # Readout noise for Biases is not generated with random seeds. So readout noise for bias images can't be reproduced. if config["ins_effects"]["cosmic_ray"] == True: if config["ins_effects"]["cray_differ"] == True: cr_map, cr_event_num = effects.produceCR_Map( Loading Loading @@ -601,7 +603,7 @@ class Chip(FocalPlane): del BiasCombImg # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file if config["output_setting"]["flat_output"] == True and pointing_type=='CAL': if config["ins_effects"]["flat_fielding"] == True and config["output_setting"]["flat_output"] == True and pointing_type=='CAL': if self.logger is not None: self.logger.info(" Output N frame Flat-Field files") else: Loading Loading @@ -673,7 +675,7 @@ class Chip(FocalPlane): # Add Read-out Noise if config["ins_effects"]["add_readout"] == True: seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 3 rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) FlatCombImg.addNoise(readout_noise) Loading Loading @@ -710,7 +712,7 @@ class Chip(FocalPlane): del flat_img # Export Dark current images if config["output_setting"]["dark_output"] == True and pointing_type=='CAL': if config["ins_effects"]["add_dark"] == True and config["output_setting"]["dark_output"] == True and pointing_type=='CAL': if self.logger is not None: self.logger.info(" Output N frame Dark Current files") else: Loading Loading @@ -795,7 +797,7 @@ class Chip(FocalPlane): # Add Read-out Noise if config["ins_effects"]["add_readout"] == True: seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID + 2 rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) DarkCombImg.addNoise(readout_noise) Loading
ObservationSim/Instrument/Chip/Effects.py +1 −1 Original line number Diff line number Diff line Loading @@ -197,7 +197,7 @@ def MakeFlatSmooth(GSBounds, seed): Flty, Fltx = np.mgrid[ymin:(ymax+1), xmin:(xmax+1)] rg = Generator(PCG64(int(seed))) p1,p2,bg=rg.poisson(1000, 3) Fltz = 1e-6*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20 Fltz = 0.6*1e-7*(a1 * (Fltx-p1) ** 2 + a2 * (Flty-p2) ** 2 - a3*Fltx - a4*Flty) + bg*20 FlatImg = galsim.ImageF(Fltz) return FlatImg Loading
ObservationSim/Instrument/Filter.py +26 −6 Original line number Diff line number Diff line Loading @@ -3,7 +3,7 @@ import pylab as pl import os import numpy as np from ObservationSim.Instrument._util import photonEnergy from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag from ObservationSim.Instrument.FilterParam import FilterParam try: Loading @@ -13,7 +13,7 @@ except ImportError: import importlib_resources as pkg_resources class Filter(object): def __init__(self, filter_id, filter_type, filter_param, ccd_bandpass): def __init__(self, filter_id, filter_type, filter_param, ccd_bandpass=None): self.filter_id = filter_id self.filter_type = filter_type self.ccd_bandpass = ccd_bandpass Loading @@ -36,14 +36,15 @@ class Filter(object): self.efficiency = filter_param.param[filter_type][4] self.sky_background = filter_param.param[filter_type][5] self.mag_saturation = filter_param.param[filter_type][6] self.mag_dim = filter_param.param[filter_type][7] self.mag_limiting = filter_param.param[filter_type][7] # self.filter_dir = filter_param.filter_dir def is_too_bright(self, mag): return mag <= self.mag_saturation - 1.0 return mag <= self.mag_saturation - 2.5 # return mag <= 14.0 def is_too_dim(self, mag): return mag >= self.mag_dim + 1.0 return mag >= self.mag_limiting + 1.0 def _get_bandpasses(self, filter_dir=None, unit='A'): if self.filter_id < 7: # Photometric Loading @@ -51,8 +52,10 @@ class Filter(object): # filter_file = os.path.join(filter_dir, self.filter_type+".dat") # bandpass_full = galsim.Bandpass(filter_file, wave_type=unit) with pkg_resources.path('ObservationSim.Instrument.data.filters', self.filter_type.lower() + '.txt') as filter_file: self.filter_bandpass = galsim.Bandpass(str(filter_file), wave_type=unit) with pkg_resources.path('ObservationSim.Instrument.data.throughputs', self.filter_type.lower() + '_throughput.txt') as filter_file: bandpass_full = galsim.Bandpass(str(filter_file), wave_type=unit) bandpass_full = bandpass_full * self.ccd_bandpass # bandpass_full = bandpass_full * self.ccd_bandpass # Get sub-bandpasses bandpass_sub_list = [] Loading Loading @@ -95,3 +98,20 @@ class Filter(object): def getSkyNoise(self, exptime, gain=1.): return self.sky_background * exptime / gain def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None): if self.filter_type in ["GI", "GV", "GU"]: return if chip is not None: pix_scale = chip.pix_scale read_noise = chip.read_noise dark_noise = chip.dark_noise full_well = chip.full_well else: pix_scale = 0.074 read_noise = 5.0 dark_noise = 0.02 full_well = 90000 throughput_file = self.filter_type.lower() + '_throughput.txt' self.mag_limiting, self.mag_saturation = calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well)
