diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index 685e12d6c4f684553d10e1bafa6f5191207af808..4004d2ac647108157b7ed075724642776d0d60b3 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -14,13 +14,13 @@ class ChipOutput(object): self.ra_cen = ra_cen self.dec_cen = dec_cen else: - self.ra_cen = config["ra_center"] - self.dec_cen = config["dec_center"] + self.ra_cen = config["obs_setting"]["ra_center"] + self.dec_cen = config["obs_setting"]["dec_center"] exp_name = imgKey0 + "_%s_%s.fits" self.chipLabel = focal_plane.getChipLabel(chip.chipID) self.img_name = prefix + exp_name%(self.chipLabel, filt.filter_type) # self.cat_name = self.img_name[:-5] + ".cat" - self.cat_name = 'MSC_' + config["date_obs"] + config["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat" + 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.subdir = subdir # hdr1 = "#ID ID_chip filter xImage yImage ra dec z mag flag SNR " diff --git a/ObservationSim/Config/Config.py b/ObservationSim/Config/Config.py index 140f143b2432256c61ca4dafa218bac787749895..a17f0a0d9460b4ac2750b256728dd7756d7e4701 100755 --- a/ObservationSim/Config/Config.py +++ b/ObservationSim/Config/Config.py @@ -2,7 +2,7 @@ import galsim import os from astropy.time import Time as asTime -def ConfigDir(cat_dir=None, work_dir=None, data_dir=None, config_file_path=None): +def ConfigDir(config, work_dir=None, data_dir=None): path_dict = {} # Working directory if work_dir == None: @@ -11,48 +11,35 @@ def ConfigDir(cat_dir=None, work_dir=None, data_dir=None, config_file_path=None) else: path_dict["work_dir"] = work_dir - # Configuration file - if config_file_path is not None: - path_dict["config_file"] = config_file_path - else: - path_dict["config_file"] = os.path.join(path_dict["work_dir"], "ObservationSim.cfg") - - # Output directories - # path_dict["output_fig_dir"] = os.path.join(path_dict["work_dir"], "figure/") - # if not os.path.exists(path_dict["output_fig_dir"]): - # os.system("mkdir %s"%path_dict["output_fig_dir"]) - # path_dict["output_cat_dir"] = os.path.join(path_dict["work_dir"], "scat/") - # if not os.path.exists(path_dict["output_cat_dir"]): - # os.system("mkdir %s"%path_dict["output_cat_dir"]) - # path_dict["output_img_dir"] = os.path.join(path_dict["work_dir"], "simg/") - # if not os.path.exists(path_dict["output_img_dir"]): - # os.system("mkdir %s"%path_dict["output_img_dir"]) - # Data directory if data_dir == None: + # Assume all input datasets are in the work directory path_dict["data_dir"] =os.path.join(path_dict["work_dir"], "data/") else: path_dict["data_dir"] = data_dir # Data sub-catalogs # Object catalog direcotry # path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], "catalog_points_7degree2/", cat_dir) - path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], "Catalog_20210126") + path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], config["input_path"]["cat_dir"]) # PSF data directory - path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], "csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90_proc") + path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_dir"]) # SED catalog directory + # TODO: SED_dir is deprecated path_dict["SED_dir"] = os.path.join(path_dict["data_dir"], "imageSims/Catalog/SEDObject") - path_dict["template_dir"] = path_dict["data_dir"] + "Templates/" + # path_dict["template_dir"] = path_dict["data_dir"] + "Templates/" + # path_dict["template_dir"] = os.path.join(path_dict["data_dir"], config["SED_templates_path"]["galaxy_SED"]) # Directories/files for instrument parameters, e.g. efficiency curves. - path_dict["filter_dir"] = os.path.join(path_dict["data_dir"], "Filters") - path_dict["ccd_dir"] = os.path.join(path_dict["data_dir"], "Filter_CCD_Mirror/ccd") - path_dict["mirror_file"] = os.path.join(path_dict["data_dir"], "Filter_CCD_Mirror/mirror_ccdnote.txt") + path_dict["filter_dir"] = os.path.join(path_dict["data_dir"], config["Efficiency_curve_path"]["filter_eff"]) + path_dict["ccd_dir"] = os.path.join(path_dict["data_dir"], config["Efficiency_curve_path"]["ccd_eff"]) + path_dict["mirror_file"] = os.path.join(path_dict["data_dir"], config["Efficiency_curve_path"]["mirror_eff"]) # Cosmic-ray data directory: - path_dict["CRdata_dir"] = os.path.join(path_dict["data_dir"], "CRdata") + path_dict["CRdata_dir"] = os.path.join(path_dict["data_dir"], config["CR_data_path"]) + path_dict["sky_file"] = os.path.join(path_dict["data_dir"], config["sky_data_path"]) # Slitless spectroscopy realted - path_dict["sls_dir"] = os.path.join(path_dict["data_dir"], "CONF/") - path_dict["normalize_dir"] = os.path.join(path_dict["data_dir"], "normalize_filter/") + path_dict["sls_dir"] = os.path.join(path_dict["data_dir"], config["SLS_path"]["SLS_conf"]) + path_dict["normalize_dir"] = os.path.join(path_dict["data_dir"], config["SLS_path"]["SLS_norm"]) return path_dict @@ -118,4 +105,4 @@ def ParseConfig(config): config["seed_Av"] = int(config["seed_Av"]) config["bias_level"] = int(config["bias_level"]) config["df_strength"] = float(config["df_strength"]) - return config \ No newline at end of file + return config diff --git a/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc b/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc index a4d3370b5d08194750f3cc5f26b9f3621d478ced..5fa40c48881f6052fa33c5dcaef8bd21e388b296 100644 Binary files a/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc and b/ObservationSim/Config/Header/__pycache__/ImageHeader.cpython-38.pyc differ diff --git a/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc index 95498b910bc3cbf213678df11cd50137784c98ce..1341eda0d8a2615b2a576d9fc08910af4b3b5ddd 100644 Binary files a/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/Config/Header/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc b/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc index 0af438058be76f3f7d83400ab223b8ab5e2f2c7e..f155567189d29c105786f72576465fb0589c6466 100644 Binary files a/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc and b/ObservationSim/Config/__pycache__/ChipOutput.cpython-38.pyc differ diff --git a/ObservationSim/Config/__pycache__/Config.cpython-38.pyc b/ObservationSim/Config/__pycache__/Config.cpython-38.pyc index bbaea4d7f2b9e0b83a4d3e16dfe00d3780980e93..ab25d1fc4ed1628753d8f0d812fd18ce6da7509b 100644 Binary files a/ObservationSim/Config/__pycache__/Config.cpython-38.pyc and b/ObservationSim/Config/__pycache__/Config.cpython-38.pyc differ diff --git a/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc index dcce485ab8de6ed97359abf898f9c603d8b1fb91..c334df4b112cf42398995ec284caf5859e8f31ff 100644 Binary files a/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/Config/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 911b29f483d77ab10b57c901ba792dca99a88632..a4348570429be8f062adfd07c7bc90efced3dead 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -30,13 +30,13 @@ class Chip(FocalPlane): self.read_noise = 5.0 # e/pix self.dark_noise = 0.02 # e/pix/s self.pix_scale = 0.074 # pixel scale - self.gain = float(config["gain"]) - self.bias_level = float(config["bias_level"]) + self.gain = float(config["ins_effects"]["gain"]) + self.bias_level = float(config["ins_effects"]["bias_level"]) self.overscan = 1000 - self.exptime = float(config["exp_time"]) # second - self.dark_exptime = float(config['dark_exptime']) - self.flat_exptime = float(config['flat_exptime']) - self.readout_time = float(config['readout_time']) + self.exptime = float(config["obs_setting"]["exp_time"]) # second + 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']) # A chip ID must be assigned self.chipID = int(chipID) @@ -65,8 +65,8 @@ class Chip(FocalPlane): self._getCRdata() # Define the sensor - if config["bright_fatter"].lower() == "y": - self.sensor = galsim.SiliconSensor(strength=config["df_strength"], treering_func=treering_func) + if config["ins_effects"]["bright_fatter"] == True: + self.sensor = galsim.SiliconSensor(strength=config["ins_effects"]["df_strength"], treering_func=treering_func) else: self.sensor = galsim.Sensor() @@ -295,18 +295,18 @@ 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): - SeedGainNonuni=int(config["seed_gainNonUniform"]) - SeedBiasNonuni=int(config["seed_biasNonUniform"]) - SeedRnNonuni = int(config["seed_rnNonUniform"]) - SeedBadColumns = int(config["seed_badcolumns"]) - SeedDefective = int(config["seed_defective"]) - SeedCosmicRay = int(config["seed_CR"]) - fullwell = int(config["full_well"]) - if config["add_hotpixels"].lower() == "y": + SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"]) + SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"]) + SeedRnNonuni = int(config["random_seeds"]["seed_rnNonUniform"]) + SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"]) + SeedDefective = int(config["random_seeds"]["seed_defective"]) + SeedCosmicRay = int(config["random_seeds"]["seed_CR"]) + fullwell = int(config["ins_effects"]["full_well"]) + if config["ins_effects"]["add_hotpixels"] == True: BoolHotPix = True else: BoolHotPix = False - if config["add_deadpixels"].lower() == "y": + if config["ins_effects"]["add_deadpixels"]== True: BoolDeadPix = True else: BoolDeadPix = False @@ -318,47 +318,47 @@ class Chip(FocalPlane): raise ValueError("The shape img and sky_map must be equal.") elif tel is not None: # If sky_map is given in flux sky_map = sky_map * tel.pupil_area * self.exptime - if config["add_back"].lower() == "y": + if config["ins_effects"]["add_back"] == True: img += sky_map del sky_map # Apply flat-field large scale structure for one chip - if config["flat_fielding"].lower() == "y": + if config["ins_effects"]["flat_fielding"] == True: print(" Creating and applying Flat-Fielding", flush=True) print(img.bounds, flush=True) flat_img = effects.MakeFlatSmooth( img.bounds, - int(config["seed_flat"])) + int(config["random_seeds"]["seed_flat"])) flat_normal = flat_img / np.mean(flat_img.array) img *= flat_normal del flat_normal - if config["flat_output"].lower() == "n": + if config["output_setting"]["flat_output"] == False: del flat_img # Apply Shutter-effect for one chip - if config["shutter_effect"].lower() == "y": + if config["ins_effects"]["shutter_effect"] == True: print(" Apply shutter effect", flush=True) shuttimg = effects.ShutterEffectArr(img, t_shutter=1.3, dist_bearing=735, dt=1E-3) # shutter effect normalized image for this chip img *= shuttimg - if config["shutter_output"].lower() == "y": # output 16-bit shutter effect image with pixel value <=65535 + if config["output_setting"]["shutter_output"] == True: # output 16-bit shutter effect image with pixel value <=65535 shutt_gsimg = galsim.ImageUS(shuttimg*6E4) shutt_gsimg.write("%s/ShutterEffect_%s_1.fits" % (chip_output.subdir, self.chipID)) del shutt_gsimg del shuttimg # Add Poisson noise - seed = int(config["seed_poisson"]) + pointing_ID*30 + self.chipID + seed = int(config["random_seeds"]["seed_poisson"]) + pointing_ID*30 + self.chipID rng_poisson = galsim.BaseDeviate(seed) poisson_noise = galsim.PoissonNoise(rng_poisson, sky_level=0.) img.addNoise(poisson_noise) # Add cosmic-rays - if config["cosmic_ray"].lower() == "y" and pointing_type=='MS': + if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS': print(" Adding Cosmic-Ray", flush=True) - cr_map = effects.produceCR_Map( + cr_map, cr_event_num = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=self.exptime+0.5*self.readout_time, - cr_pixelRatio=0.003*(1+0.5*self.readout_time/self.exptime), + cr_pixelRatio=0.003*(self.exptime+0.5*self.readout_time)/150., gain=self.gain, attachedSizes=self.attachedSizes, seed=SeedCosmicRay+pointing_ID*30+self.chipID) # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3; @@ -382,25 +382,25 @@ class Chip(FocalPlane): date_obs=date_obs, time_obs=time_obs, output_dir=chip_output.subdir, - exptime=150.) + exptime=self.exptime) del crmap_gsimg # Apply PRNU effect and output PRNU flat file: - if config["prnu_effect"].lower() == "y": + if config["ins_effects"]["prnu_effect"] == True: print(" Applying PRNU effect", flush=True) prnu_img = effects.PRNU_Img( xsize=self.npix_x, ysize=self.npix_y, sigma=0.01, - seed=int(config["seed_prnu"]+self.chipID)) + seed=int(config["random_seeds"]["seed_prnu"]+self.chipID)) img *= prnu_img - if config["prnu_output"].lower() == "y": + if config["output_setting"]["prnu_output"] == True: prnu_img.write("%s/FlatImg_PRNU_%s.fits" % (chip_output.subdir,self.chipID)) - if config["flat_output"].lower() == "n": + if config["output_setting"]["flat_output"] == False: del prnu_img # Add dark current - if config["add_dark"].lower() == "y": + if config["ins_effects"]["add_dark"] == True: dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(self.exptime+0.5*self.readout_time))) img.addNoise(dark_noise) @@ -410,35 +410,35 @@ class Chip(FocalPlane): img = effects.DefectivePixels(img, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0) # Apply Bad lines - if config["add_badcolumns"].lower() == "y": + if config["ins_effects"]["add_badcolumns"] == True: img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID) # Add Bias level - if config["add_bias"].lower() == "y": + if config["ins_effects"]["add_bias"] == True: print(" Adding Bias level and 16-channel non-uniformity") img = effects.AddBiasNonUniform16(img, - bias_level=float(config["bias_level"]), + bias_level=float(config["ins_effects"]["bias_level"]), nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID) # Apply Nonlinearity on the chip image - if config["non_linear"].lower() == "y": + if config["ins_effects"]["non_linear"] == True: print(" Applying Non-Linearity on the chip image", flush=True) img = effects.NonLinearity(GSImage=img, beta1=5.e-7, beta2=0) # Apply CCD Saturation & Blooming - if config["saturbloom"].lower() == "y": + if config["ins_effects"]["saturbloom"] == True: print(" Applying CCD Saturation & Blooming") img = effects.SaturBloom(GSImage=img, nsect_x=1, nsect_y=1, fullwell=fullwell) # Apply CTE Effect - if config["cte_trail"].lower() == "y": + if config["ins_effects"]["cte_trail"] == True: print(" Apply CTE Effect") img = effects.CTE_Effect(GSImage=img, threshold=27) # Add Read-out Noise - if config["add_readout"].lower() == "y": - seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID + if config["ins_effects"]["add_readout"] == True: + seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) img.addNoise(readout_noise) @@ -458,21 +458,21 @@ class Chip(FocalPlane): # Output images for calibration pointing ###################################################################################### # Bias output - if config["bias_output"].lower() == "y" and pointing_type=='CAL': + if config["output_setting"]["bias_output"] == True and pointing_type=='CAL': print(" Output N frame Bias files", flush=True) - NBias = int(config["NBias"]) + NBias = int(config["ins_effects"]["NBias"]) for i in range(NBias): BiasCombImg, BiasTag = effects.MakeBiasNcomb( self.npix_x, self.npix_y, - bias_level=float(config["bias_level"]), + bias_level=float(config["ins_effects"]["bias_level"]), ncombine=1, read_noise=self.read_noise, gain=1, seed=SeedBiasNonuni+self.chipID) - if config["cosmic_ray"].lower() == "y": - if config["cray_differ"].lower() == "y": - cr_map = effects.produceCR_Map( + if config["ins_effects"]["cosmic_ray"] == True: + if config["ins_effects"]["cray_differ"] == True: + cr_map, cr_event_num = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=0.01, - cr_pixelRatio=0.003*(0.01+0.5*self.readout_time)/(self.exptime+0.5*self.readout_time), + cr_pixelRatio=0.003*(0.01+0.5*self.readout_time)/150., gain=self.gain, attachedSizes=self.attachedSizes, seed=SeedCosmicRay+pointing_ID*30+self.chipID+1) @@ -481,20 +481,20 @@ class Chip(FocalPlane): del cr_map # Non-Linearity for Bias - if config["non_linear"].lower() == "y": + if config["ins_effects"]["non_linear"] == True: print(" Applying Non-Linearity on the Bias image", flush=True) BiasCombImg = effects.NonLinearity(GSImage=BiasCombImg, beta1=5.e-7, beta2=0) # Apply Bad lines - if config["add_badcolumns"].lower() == "y": - BiasCombImg = effects.BadColumns(BiasCombImg-float(config["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID) + float(config["bias_level"])-5 + if config["ins_effects"]["add_badcolumns"] == True: + BiasCombImg = effects.BadColumns(BiasCombImg-float(config["ins_effects"]["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID) + float(config["ins_effects"]["bias_level"])-5 BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain, nsecy = 2, nsecx=8, seed=SeedGainNonuni+self.chipID) # BiasCombImg = effects.AddOverscan( # BiasCombImg, - # overscan=float(config["bias_level"])-2, gain=self.gain, + # overscan=float(config["ins_effects"]["bias_level"])-2, gain=self.gain, # widthl=27, widthr=27, widtht=8, widthb=8) BiasCombImg.replaceNegative(replace_value=0) BiasCombImg.quantize() @@ -518,13 +518,13 @@ class Chip(FocalPlane): del BiasCombImg # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file - if config["flat_output"].lower() == "y" and pointing_type=='CAL': + if config["output_setting"]["flat_output"] == True and pointing_type=='CAL': print(" Output N frame Flat-Field files", flush=True) - NFlat = int(config["NFlat"]) - if config["add_bias"].lower() == "y": + NFlat = int(config["ins_effects"]["NFlat"]) + if config["ins_effects"]["add_bias"] == True: biaslevel = self.bias_level overscan = biaslevel-2 - elif config["add_bias"].lower() == "n": + elif config["ins_effects"]["add_bias"] == True: biaslevel = 0 overscan = 0 darklevel = self.dark_noise*(self.flat_exptime+0.5*self.readout_time) @@ -539,12 +539,12 @@ class Chip(FocalPlane): biaslevel=0, seed_bias=SeedDefective+self.chipID ) - if config["cosmic_ray"].lower() == "y": - if config["cray_differ"].lower() == "y": - cr_map = effects.produceCR_Map( + if config["ins_effects"]["cosmic_ray"] == True: + if config["ins_effects"]["cray_differ"] == True: + cr_map, cr_event_num = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=self.flat_exptime+0.5*self.readout_time, - cr_pixelRatio=0.003*(self.flat_exptime+0.5*self.readout_time)/(self.exptime+0.5*self.readout_time), + cr_pixelRatio=0.003*(self.flat_exptime+0.5*self.readout_time)/150., gain=self.gain, attachedSizes=self.attachedSizes, seed=SeedCosmicRay+pointing_ID*30+self.chipID+3) @@ -552,11 +552,11 @@ class Chip(FocalPlane): FlatCombImg += cr_map del cr_map - if config["non_linear"].lower() == "y": + if config["ins_effects"]["non_linear"] == True: print(" Applying Non-Linearity on the Flat image", flush=True) FlatCombImg = effects.NonLinearity(GSImage=FlatCombImg, beta1=5.e-7, beta2=0) - if config["cte_trail"].lower() == "y": + if config["ins_effects"]["cte_trail"] == True: FlatCombImg = effects.CTE_Effect(GSImage=FlatCombImg, threshold=3) # Add Hot Pixels or/and Dead Pixels @@ -565,21 +565,21 @@ class Chip(FocalPlane): FlatCombImg = effects.DefectivePixels(FlatCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0) # Apply Bad lines - if config["add_badcolumns"].lower() == "y": + if config["ins_effects"]["add_badcolumns"] == True: FlatCombImg = effects.BadColumns(FlatCombImg, seed=SeedBadColumns, chipid=self.chipID) # Add Bias level - if config["add_bias"].lower() == "y": + if config["ins_effects"]["add_bias"] == True: print(" Adding Bias level and 16-channel non-uniformity") - # img += float(config["bias_level"]) + # img += float(config["ins_effects"]["bias_level"]) FlatCombImg = effects.AddBiasNonUniform16(FlatCombImg, bias_level=biaslevel, nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID) # Add Read-out Noise - if config["add_readout"].lower() == "y": - seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID + if config["ins_effects"]["add_readout"] == True: + seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) FlatCombImg.addNoise(readout_noise) @@ -615,13 +615,13 @@ class Chip(FocalPlane): del flat_img # Export Dark current images - if config["dark_output"].lower() == "y" and pointing_type=='CAL': + if config["output_setting"]["dark_output"] == True and pointing_type=='CAL': print(" Output N frame Dark Current files", flush=True) - NDark = int(config["NDark"]) - if config["add_bias"].lower() == "y": + NDark = int(config["ins_effects"]["NDark"]) + if config["ins_effects"]["add_bias"] == True: biaslevel = self.bias_level overscan = biaslevel-2 - elif config["add_bias"].lower() == "n": + elif config["ins_effects"]["add_bias"] == True: biaslevel = 0 overscan = 0 for i in range(NDark): @@ -630,12 +630,12 @@ class Chip(FocalPlane): overscan=overscan, bias_level=0, darkpsec=0.02, exptime=self.dark_exptime+0.5*self.readout_time, ncombine=1, read_noise=self.read_noise, gain=1, seed_bias=SeedBiasNonuni+self.chipID) - if config["cosmic_ray"].lower() == "y": - if config["cray_differ"].lower() == "y": - cr_map = effects.produceCR_Map( + if config["ins_effects"]["cosmic_ray"] == True: + if config["ins_effects"]["cray_differ"] == True: + cr_map, cr_event_num = effects.produceCR_Map( xLen=self.npix_x, yLen=self.npix_y, exTime=self.dark_exptime+0.5*self.readout_time, - cr_pixelRatio=0.003*(self.dark_exptime+0.5*self.readout_time)/(self.exptime+0.5*self.readout_time), + cr_pixelRatio=0.003*(self.dark_exptime+0.5*self.readout_time)/150., gain=self.gain, attachedSizes=self.attachedSizes, seed=SeedCosmicRay+pointing_ID*30+self.chipID+2) @@ -662,11 +662,11 @@ class Chip(FocalPlane): del crmap_gsimg # Non-Linearity for Dark - if config["non_linear"].lower() == "y": + if config["ins_effects"]["non_linear"] == True: print(" Applying Non-Linearity on the Dark image", flush=True) DarkCombImg = effects.NonLinearity(GSImage=DarkCombImg, beta1=5.e-7, beta2=0) - if config["cte_trail"].lower() == "y": + if config["ins_effects"]["cte_trail"] == True: DarkCombImg = effects.CTE_Effect(GSImage=DarkCombImg, threshold=3) # Add Hot Pixels or/and Dead Pixels @@ -675,21 +675,21 @@ class Chip(FocalPlane): DarkCombImg = effects.DefectivePixels(DarkCombImg, IfHotPix=BoolHotPix, IfDeadPix=BoolDeadPix, fraction=badfraction, seed=SeedDefective+self.chipID, biaslevel=0) # Apply Bad lines - if config["add_badcolumns"].lower() == "y": + if config["ins_effects"]["add_badcolumns"] == True: DarkCombImg = effects.BadColumns(DarkCombImg, seed=SeedBadColumns, chipid=self.chipID) # Add Bias level - if config["add_bias"].lower() == "y": + if config["ins_effects"]["add_bias"] == True: print(" Adding Bias level and 16-channel non-uniformity") - # img += float(config["bias_level"]) + # img += float(config["ins_effects"]["bias_level"]) DarkCombImg = effects.AddBiasNonUniform16(DarkCombImg, bias_level=biaslevel, nsecy = 2, nsecx=8, seed=SeedBiasNonuni+self.chipID) # Add Read-out Noise - if config["add_readout"].lower() == "y": - seed = int(config["seed_readout"]) + pointing_ID*30 + self.chipID + if config["ins_effects"]["add_readout"] == True: + seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID rng_readout = galsim.BaseDeviate(seed) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) DarkCombImg.addNoise(readout_noise) @@ -725,7 +725,7 @@ class Chip(FocalPlane): # img = galsim.ImageUS(img) # # 16 output channel, with each a single image file - # if config["readout16"].lower() == "y": + # if config["ins_effects"]["readout16"] == True: # print(" 16 Output Channel simulation") # for coli in [0, 1]: # for rowi in range(8): diff --git a/ObservationSim/Instrument/Chip/Effects.py b/ObservationSim/Instrument/Chip/Effects.py index 441359ef9fa82e90bcb0f6995fb3a02c70d7e21f..49496ee24e7ec04bd53b811a9e4927f264d21e59 100644 --- a/ObservationSim/Instrument/Chip/Effects.py +++ b/ObservationSim/Instrument/Chip/Effects.py @@ -669,7 +669,7 @@ def produceCR_Map(xLen, yLen, exTime, cr_pixelRatio, gain, attachedSizes, seed=2 CRmap[sly, slx] += crMatrix_n[oky,:][:,okx] - return CRmap.astype(np.int32); + return CRmap.astype(np.int32), cr_event_size def ShutterEffectArr(GSImage, t_shutter=1.3, dist_bearing=735, dt=1E-3): diff --git a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc index 5db77ca069ccb0db010ad81ff1dc9acc6472fdea..494be8e7feb822165cd70af27a4429697328180a 100644 Binary files a/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc and b/ObservationSim/Instrument/Chip/__pycache__/Chip.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc index 303968f15726e04a2c7b9af955df831995dc922e..5e2fc0e5926b1ffd19e65dd2b8aebc6b4f963435 100644 Binary files a/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc and b/ObservationSim/Instrument/Chip/__pycache__/Effects.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc index 277c47b11ccb2236e498dfb7ae94d0e2b5c6e16e..08322b16d722afbc1f14d8449b5c2b3cb2df92ff 100644 Binary files a/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/Instrument/Chip/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc index 2b37ab185f518e1771a98f36a77a51e6d17f926b..3a3a9a6609e559f8145ea9a1a0883548e44037d0 100644 Binary files a/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/Filter.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc index bb08b43d2cb7a5e622d265d9963bf473660d28d6..881c1cf95e5e6c1a90ac9b8e38cc20e4e57f83e8 100644 Binary files a/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/FilterParam.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc index b43df585778cb90566db7f9e8b1333f5abff2b1e..7997be3e03e48f50d9ffaf6c74af880e7f256906 100644 Binary files a/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/FocalPlane.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc index 6bf29c20e588c6038023de8ce4567086ed8e14be..dae2816701485d0eebfc28213839a1e2a6ac0086 100644 Binary files a/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/Telescope.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc index 76ee1b645269ab0afe006079e97121d5fe1137ca..92ea9cd16903a0ace81a3b6a93432ec42993d90e 100644 Binary files a/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc b/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc index 048e4c873f60d594ec5b3d266f22fc7814e98b2c..59c3215c5df8d78daeb96409c9ab3b54c20e1fbd 100644 Binary files a/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc and b/ObservationSim/Instrument/__pycache__/_util.cpython-38.pyc differ diff --git a/ObservationSim/Instrument/_util.py b/ObservationSim/Instrument/_util.py index 9442cb4ba1358605317435507fa20e8a775dc826..f93f385639d58ba4896e3dac45f535e26e45f51b 100755 --- a/ObservationSim/Instrument/_util.py +++ b/ObservationSim/Instrument/_util.py @@ -1,9 +1,9 @@ import numpy as np import os -vc_A = 2.99792458e+18 # speed of light: A/s -vc_m = 2.99792458e+8 # speed of light: m/s -h_Plank = 6.626196e-27 # Plank constant: erg s +VC_A = 2.99792458e+18 # speed of light: A/s +VC_M = 2.99792458e+8 # speed of light: m/s +H_PLANK = 6.626196e-27 # Plank constant: erg s def photonEnergy(lambd): """ The energy of photon at a given wavelength @@ -13,6 +13,6 @@ def photonEnergy(lambd): Return: eph: energy of photon in unit of erg """ - nu = vc_A / lambd - eph = h_Plank * nu + nu = VC_A / lambd + eph = H_PLANK * nu return eph \ No newline at end of file diff --git a/ObservationSim/MockObject/Catalog.py b/ObservationSim/MockObject/Catalog.py index 6da49ec8785ad4bfd6f07aedc3d382f887442d55..2ebbcdd697d759150566f054e25f243dd3ccc64c 100644 --- a/ObservationSim/MockObject/Catalog.py +++ b/ObservationSim/MockObject/Catalog.py @@ -11,32 +11,43 @@ from ._util import seds, sed_assign, extAv from astropy.table import Table from astropy.coordinates import spherical_to_cartesian -nside = 128 +NSIDE = 128 class Catalog(object): - def __init__(self, config, chip, cat_dir, sed_dir, pRa, pDec, rotation=None, template_dir=None): - self.cat_dir = cat_dir + def __init__(self, config, chip, cat_dir=None, pRa=None, pDec=None, sed_dir=None, rotation=None): + if cat_dir is not None: + self.cat_dir = cat_dir + else: + self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"]) self.sed_dir = sed_dir self.chip = chip - self.template_dir = template_dir - self.seed_Av = config["seed_Av"] + self.seed_Av = config["random_seeds"]["seed_Av"] - self.pRa = float('{:8.4f}'.format(pRa)) - self.pDec = float('{:8.4f}'.format(pDec)) + if pRa is not None: + self.pRa = float('{:8.4f}'.format(pRa)) + if pDec is not None: + self.pDec = float('{:8.4f}'.format(pDec)) # star_file = 'stars_ccd' + chip.getChipLabel(chip.chipID) + '_p_RA'+str(self.pRa) + '_DE' + str(self.pDec) + '.hdf5' # galaxy_file = 'galaxies_ccd' + chip.getChipLabel(chip.chipID) + '_p_RA'+str(self.pRa) + '_DE' + str(self.pDec) + '.hdf5' - # star_file = 'MMW_Cluster_D20_10_healpix.hdf5' - star_file = 'MMW_Cluster_D20_healpix_21.hdf5' - galaxy_file = 'galaxyCats_r_10.0_healpix.hdf5' - star_SED_file = 'SpecLib.hdf5' - self.star_path = os.path.join(cat_dir, star_file) - self.galaxy_path = os.path.join(cat_dir, galaxy_file) - self.star_SED_path = os.path.join(cat_dir, star_SED_file) - self.rotation = rotation - self._load_SED_lib_gals() - self._load_SED_lib_star() + + if "star_cat" in config["input_path"]: + star_file = config["input_path"]["star_cat"] + star_SED_file = config["SED_templates_path"]["star_SED"] + self.star_path = os.path.join(self.cat_dir, star_file) + self.star_SED_path = os.path.join(config["data_dir"], star_SED_file) + self._load_SED_lib_star() + if "galaxy_cat" in config["input_path"]: + + galaxy_file = config["input_path"]["galaxy_cat"] + self.galaxy_path = os.path.join(self.cat_dir, galaxy_file) + self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"]) + self._load_SED_lib_gals() + if "rotateEll" in config["shear_setting"]: + self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.)) + else: + self.rotation = 0. # self._load_SED_info() self._get_healpix_list() self._load() @@ -48,20 +59,20 @@ class Catalog(object): ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min])) dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min])) # phi, theta = ra, np.pi/2. - dec - # vertices_pix = np.unique(hp.ang2pix(nside, theta, phi)) + # vertices_pix = np.unique(hp.ang2pix(NSIDE, theta, phi)) # print(vertices_pix) vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon(nside, np.array(vertices).T, inclusive=True) + self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True) print("HEALPix List: ", self.pix_list) def _load_SED_lib_star(self): self.tempSED_star = h5.File(self.star_SED_path,'r') def _load_SED_lib_gals(self): - tempdir_gal = os.path.join(self.template_dir, "Galaxy/") + # tempdir_gal = os.path.join(self.template_dir, "Galaxy/") # tempdir_star = os.path.join(self.template_dir, "PicklesStars/") - self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=tempdir_gal) # self.tempSed_star, self.tempRed_star = seds("star.list", seddir=tempdir_star) + self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=self.galaxy_SED_path) def _load_SED_info(self): sed_info_file = os.path.join(self.sed_dir, "sed.info") diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 4c75c3fb7acf1700bfdb75d84433b51d2f8c7e89..6d89c15efc64c564dce3967bea349d412b1f87a9 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -27,7 +27,7 @@ class Galaxy(MockObject): if rotation is not None: self.rotateEllipticity(rotation) - def load_SED(self, survey_type, sed_path, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None): + def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None): if survey_type == "photometric": norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 if sed_templates is None: diff --git a/ObservationSim/MockObject/Quasar.py b/ObservationSim/MockObject/Quasar.py index 128318569a7c3346854d11c970efe9172ab4f5d3..d88659e064e371aa5cddca36818e375d9025259a 100755 --- a/ObservationSim/MockObject/Quasar.py +++ b/ObservationSim/MockObject/Quasar.py @@ -12,7 +12,7 @@ class Quasar(MockObject): def __init__(self, param): super().__init__(param) - def load_SED(self, survey_type, sed_path, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None): + def load_SED(self, survey_type, sed_path=None, cosids=None, objtypes=None, sed_templates=None, normFilter=None, target_filt=None): if survey_type == "photometric": norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 if sed_templates is None: diff --git a/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc b/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc index b19854db741e875348d2f21151daef7d417a8f52..4b3396d630ade0fecffa135e8888c359181ca56f 100644 Binary files a/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc and b/ObservationSim/MockObject/SpecDisperser/__pycache__/SpecDisperser.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc b/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc index bb9ae4ff8db2343afa4cf8f97e893b84a08de064..13c4d80fd308e0ffa85ad98ed98efc58c7b55605 100644 Binary files a/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/MockObject/SpecDisperser/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc b/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc index b97e9a26fd9c012c40ea2433061431b57d551d33..60a53d5729681648bf79e00b9e0a2b4cad6ace70 100644 Binary files a/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/MockObject/SpecDisperser/disperse_c/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc index d34d179ec40c26b378c0cdcea67b8fb0410f5073..2aec03940114bfd44f9bbf11f84bca836deb02e1 100644 Binary files a/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/Catalog.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc index 4eaa9589bf1fc54e04217d2675c89b5651ace89c..a663442ad59b1faa9182ac376dddb0907d44c9b6 100644 Binary files a/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/Galaxy.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc index e2afe2b63b9bccbc017e6eb977bd759859edb614..4fe9639e24557f0c661e2900a958163b47b1359b 100644 Binary files a/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/MockObject.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc index 7b75a46094c6b6b94ba3f5e576621bd2aea6a2ae..e3f189462efdcd7f72c2e6525c890b62a279f505 100644 Binary files a/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/Quasar.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc index 84666cd5e170eb86c2fb25334012611289c2ef66..f1c7e805f9c62cd7bfabda603c47a82abf5b13dd 100644 Binary files a/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/SkybackgroundMap.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc index 3f06002d807ac881f578e31fe9d0addd0e881571..512a0c3e954127c7ac3736af3e1e13775937b48a 100644 Binary files a/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/Star.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc index 4aeb75f96f2fb0bf2f27fd4bf6183b956a7926d1..556670e94d9879c41d266edc520fa4d6e51685b4 100644 Binary files a/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc b/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc index a3c2bdaba0cafe57c6dd7d1ec6e1cce576bc459f..d1e72047f5abebeb3fa7f57bcda249428a676637 100644 Binary files a/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc and b/ObservationSim/MockObject/__pycache__/_util.cpython-38.pyc differ diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index d49957f877f1efe3eb0ef99abb2378bf8a7737e9..1824cf76b2defbde11ac889736f61ddd52d6b358 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -1,9 +1,9 @@ -from Config import ConfigDir, ReadConfig, ChipOutput +from Config import ConfigDir, ChipOutput from Config.Header import generatePrimaryHeader, generateExtensionHeader from Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip from MockObject import Catalog, MockObject, Star, Galaxy, Quasar, calculateSkyMap_split_g -from PSF import PSFGauss, PSFInterp, FieldDistortion -from _util import makeSubDir, getShearFiled, makeSubDir_PointingList +from PSF import PSFGauss, FieldDistortion, PSFInterp +from _util import getShearFiled, makeSubDir_PointingList from astropy.io import fits from datetime import datetime import numpy as np @@ -14,17 +14,18 @@ import logging import psutil class Observation(object): - def __init__(self, input_cat_dir=None, work_dir=None, data_dir=None): - self.path_dict = ConfigDir(input_cat_dir, work_dir, data_dir) - self.config = ReadConfig(self.path_dict["config_file"]) + def __init__(self, config, work_dir=None, data_dir=None): + self.path_dict = ConfigDir(config=config, work_dir=work_dir, data_dir=data_dir) + # self.config = ReadConfig(self.path_dict["config_file"]) + self.config = config self.tel = Telescope(optEffCurve_path=self.path_dict["mirror_file"]) # Currently the default values are hard coded in - self.focal_plane = FocalPlane(survey_type=self.config["survey_type"]) # Currently the default values are hard coded in + self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) # Currently the default values are hard coded in self.filter_param = FilterParam(filter_dir=self.path_dict["filter_dir"]) # Currently the default values are hard coded in self.chip_list = [] self.filter_list = [] # if we want to apply field distortion? - if self.config["field_dist"].lower() == "y": + if self.config["ins_effects"]["field_dist"] == True: self.fd_model = FieldDistortion() else: self.fd_model = None @@ -47,17 +48,17 @@ class Observation(object): self.g1_field, self.g2_field, self.nshear = getShearFiled(config=self.config) - def runOneChip(self, chip, filt, chip_output, wcs_fp=None, psf_model=None, pointing_ID=0, ra_cen=None, dec_cen=None, img_rot=None, exptime=150., timestamp_obs=1621915200, pointing_type='MS', input_cat_name=None, shear_cat_file=None, cat_dir=None, sed_dir=None): + def runOneChip(self, chip, filt, chip_output, wcs_fp=None, psf_model=None, pointing_ID=0, ra_cen=None, dec_cen=None, img_rot=None, exptime=150., timestamp_obs=1621915200, pointing_type="MS", shear_cat_file=None, cat_dir=None, sed_dir=None): if (ra_cen is None) or (dec_cen is None): - ra_cen = self.config["ra_center"] - dec_cen = self.config["dec_center"] + ra_cen = self.config["obs_setting"]["ra_center"] + dec_cen = self.config["obs_setting"]["dec_center"] if img_rot is None: - img_rot = self.config["image_rot"] + img_rot = self.config["obs_setting"]["image_rot"] - if self.config["psf_model"] == "Gauss": + if self.config["psf_setting"]["psf_model"] == "Gauss": psf_model = PSFGauss(chip=chip) - elif self.config["psf_model"] == "Interp": + elif self.config["psf_setting"]["psf_model"] == "Interp": psf_model = PSFInterp(chip=chip) else: print("unrecognized PSF model type!!", flush=True) @@ -77,12 +78,11 @@ class Observation(object): if chip.survey_type == "photometric": sky_map = None elif chip.survey_type == "spectroscopic": - skyfile = os.path.join(self.path_dict["data_dir"], 'skybackground/sky_emiss_hubble_50_50_A.dat') - sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=skyfile, conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0) + sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0) if pointing_type == 'MS': # Load catalogues and templates - self.cat = Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir, pRa=ra_cen, pDec=dec_cen, rotation=img_rot, template_dir=self.path_dict["template_dir"]) + self.cat = Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir) self.nobj = len(self.cat.objs) # Loop over objects @@ -97,7 +97,8 @@ class Observation(object): # Load SED if obj.type == 'star': normF = chip.normF_star - #continue + if self.config["galaxy_only"]: + continue try: obj.load_SED( survey_type=chip.survey_type, @@ -109,18 +110,18 @@ class Observation(object): continue elif obj.type == 'galaxy': # or obj.type == quasar normF = chip.normF_galaxy - # continue + if self.config["star_only"]: + continue obj.load_SED( - sed_path=sed_dir, survey_type=chip.survey_type, sed_templates=self.cat.tempSed_gal, normFilter=normF, target_filt=filt) elif obj.type == 'quasar': normF = chip.normF_galaxy - # continue + if self.config["star_only"]: + continue obj.load_SED( - sed_path=sed_dir, survey_type=chip.survey_type, sed_templates=self.cat.tempSed_gal, normFilter=normF, @@ -140,14 +141,18 @@ class Observation(object): # print(obj.getMagFilter(filt)) continue - if self.config["shear_method"] == "constant": + if self.config["shear_setting"]["shear_type"] == "constant": if obj.type == 'star': g1, g2 = 0, 0 else: g1, g2 = self.g1_field, self.g2_field - elif self.config["shear_method"] == "extra": - # TODO: every object with individual shear from input catalog(s) - g1, g2 = self.g1_field[j], self.g2_field[j] + elif self.config["shear_setting"]["shear_type"] == "extra": + try: + # TODO: every object with individual shear from input catalog(s) + g1, g2 = self.g1_field[j], self.g2_field[j] + except: + print("failed to load external shear.") + pass pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False) if pos_img.x == -1 or pos_img.y == -1: @@ -159,7 +164,9 @@ class Observation(object): # Draw object & update output catalog try: - if chip.survey_type == "photometric": + if self.config["out_cat_only"]: + isUpdated = True + if chip.survey_type == "photometric" and not self.config["out_cat_only"]: isUpdated, pos_shear = obj.drawObj_multiband( tel=self.tel, pos_img=pos_img, @@ -170,7 +177,7 @@ class Observation(object): g1=g1, g2=g2, exptime=exptime) - elif chip.survey_type == "spectroscopic": + elif chip.survey_type == "spectroscopic" and not self.config["out_cat_only"]: isUpdated, pos_shear = obj.drawObj_slitless( tel=self.tel, pos_img=pos_img, @@ -203,8 +210,6 @@ class Observation(object): # Detector Effects # =========================================================== - # chip.img = chip.addNoise(config=self.config, tel=self.tel, filt=filt, img=chip.img, sky_map=sky_map) - # chip.img = chip.addEffects(config=self.config, img=chip.img, chip_output=chip_output, filt=filt, pointing_ID=pointing_ID) # whether to output zero, dark, flat calibration images. chip.img = chip.addEffects( @@ -253,6 +258,7 @@ class Observation(object): col_num=chip.colID, extName='raw') chip.img = galsim.Image(chip.img.array, dtype=np.uint16) + # chip.img = galsim.Image(chip.img.array, dtype=np.uint32) hdu1 = fits.PrimaryHDU(header=h_prim) hdu2 = fits.ImageHDU(chip.img.array, header=h_ext) hdu1 = fits.HDUList([hdu1, hdu2]) @@ -266,13 +272,13 @@ class Observation(object): print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing_ID, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True) - def runExposure(self, ra_cen=None, dec_cen=None, timestamp_obs=1621915200, pointing_ID=0, pointing_type='MS', img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None, oneChip=None): + def runExposure(self, ra_cen=None, dec_cen=None, timestamp_obs=1621915200, pointing_ID=0, pointing_type='MS', img_rot=None, exptime=150., shear_cat_file=None, chips=None): if (ra_cen == None) or (dec_cen == None): - ra_cen = self.config["ra_center"] - dec_cen = self.config["dec_center"] + ra_cen = self.config["obs_setting"]["ra_center"] + dec_cen = self.config["obs_setting"]["dec_center"] if img_rot == None: - img_rot = self.config["image_rot"] + img_rot = self.config["obs_setting"]["image_rot"] sub_img_dir, prefix = makeSubDir_PointingList(path_dict=self.path_dict, config=self.config, pointing_ID=pointing_ID) @@ -281,9 +287,9 @@ class Observation(object): chip = self.chip_list[i] filt = self.filter_list[i] - # Just run one chip - if oneChip is not None: - if chip.chipID != oneChip: + # Only run a particular set of chips + if chips is not None: + if chip.chipID not in chips: continue # Prepare output files @@ -309,16 +315,30 @@ class Observation(object): exptime=exptime, timestamp_obs=timestamp_obs, pointing_type=pointing_type, - cat_dir=self.path_dict["cat_dir"], - sed_dir=self.path_dict["SED_dir"]) + cat_dir=self.path_dict["cat_dir"]) print("finished running chip#%d..."%(chip.chipID), flush=True) - def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=np.array([1621915200]), pointing_type=np.array(['MS']), img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None): - comm = MPI.COMM_WORLD - ind_thread = comm.Get_rank() - num_thread = comm.Get_size() + def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=np.array([1621915200]), pointing_type=np.array(['MS']), img_rot=None, exptime=150., shear_cat_file=None, chips=None, use_mpi=False): + if use_mpi: + comm = MPI.COMM_WORLD + ind_thread = comm.Get_rank() + num_thread = comm.Get_size() - nchips_per_fp = len(self.chip_list) + 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) # TEMP if len(timestamp_obs) == 1: @@ -343,15 +363,18 @@ class Observation(object): pointing_ID = pStart + ipoint else: pointing_ID = pRange[ipoint] - if i % num_thread != ind_thread: - continue + if use_mpi: + if i % 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 = self.chip_list[ichip] - filt = self.filter_list[ichip] + # chip = self.chip_list[ichip] + # filt = self.filter_list[ichip] + chip = run_chips[ichip] + filt = run_filts[ichip] print("running pointing#%d, chip#%d, at PID#%d..."%(pointing_ID, chip.chipID, pid), flush=True) chip_output = ChipOutput( config=self.config, @@ -374,6 +397,5 @@ class Observation(object): exptime=exptime, timestamp_obs=timestamp_obs[ipoint], pointing_type=pointing_type[ipoint], - cat_dir=self.path_dict["cat_dir"], - sed_dir=self.path_dict["SED_dir"]) + cat_dir=self.path_dict["cat_dir"]) print("finished running chip#%d..."%(chip.chipID), flush=True) diff --git a/ObservationSim/PSF/PSFInterp/PSFConfig.py b/ObservationSim/PSF/PSFInterp/PSFConfig.py index 3f2fef8a97c71f7d60015ca896c19812c9956185..cdbc7ead2b60e4eefd9187bdded1b71fbdbab752 100644 --- a/ObservationSim/PSF/PSFInterp/PSFConfig.py +++ b/ObservationSim/PSF/PSFInterp/PSFConfig.py @@ -3,7 +3,8 @@ import numpy as np import scipy.io from itertools import islice -from PSF.PSFInterp.PSFUtil import * +# from PSF.PSFInterp.PSFUtil import * +# from PSFUtil import findNeighbors, findNeighbors_hoclist @@ -141,6 +142,8 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=Tru npsf = PSFMat[:, :, :].shape[0] psfWeight = np.zeros([npsf]) + from .PSFUtil import findNeighbors, findNeighbors_hoclist + if OnlyNeighbors == True: if hoc is None: neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=4, OnlyDistance=False) diff --git a/ObservationSim/PSF/PSFInterp/PSFInterp.py b/ObservationSim/PSF/PSFInterp/PSFInterp.py index 92d79ca7efa517f25d641fdfe0833a85d3c53be8..dc3ba25916506cb02bef4dffa1a8d498475b1033 100644 --- a/ObservationSim/PSF/PSFInterp/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp/PSFInterp.py @@ -10,9 +10,10 @@ import numpy as np import os import time import copy -import PSF.PSFInterp.PSFConfig as myConfig -import PSF.PSFInterp.PSFUtil as myUtil -from PSF.PSFModel import PSFModel + +from . import PSFConfig as myConfig +from . import PSFUtil as myUtil +from ..PSFModel import PSFModel LOG_DEBUG = False #***# NPSF = 900 #***# 30*30 @@ -29,6 +30,7 @@ class PSFInterp(PSFModel): params: Other parameters? PSF_data_file: The file for PSF data matrix (optional). """ + from . import PSFUtil as myUtil if LOG_DEBUG: print('===================================================') print('DEBUG: psf module for csstSim ' \ @@ -55,6 +57,7 @@ class PSFInterp(PSFModel): self.cen_row= np.zeros([self.nwave, NPSF], dtype=np.float32) self.hoc =[] self.hoclist=[] + for twave in range(self.nwave): for tpsf in range(NPSF): self.psfMat[twave, tpsf, :, :] = self.PSF_data[twave][tpsf]['psfMat'] @@ -100,6 +103,7 @@ class PSFInterp(PSFModel): Returns: psfSet: The matrix of the csst-psf """ + from . import PSFConfig as myConfig psfSet = [] for ii in range(self.nwave): iwave = ii+1 @@ -144,7 +148,7 @@ class PSFInterp(PSFModel): Returns: PSF: A 'galsim.GSObject'. """ - + from . import PSFConfig as myConfig pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize #***# assert self.iccd == chip.chipID, 'ERROR: self.iccd != chip.chipID' @@ -200,8 +204,6 @@ class PSFInterp(PSFModel): return self.psf.shear(PSFshear), PSFshear - - if __name__ == '__main__': if True: psfPath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90' @@ -256,6 +258,8 @@ if __name__ == '__main__': from scipy import ndimage y, x = ndimage.center_of_mass(img) plt.plot([x-dng],[y-dng], 'rx', ms = 10, mew=4) + + from . import PSFUtil as myUtil x, y = myUtil.findMaxPix(img) plt.plot([x-dng],[y-dng], 'b+', ms = 10, mew=4) plt.savefig('test/figs/test2.jpg') diff --git a/ObservationSim/PSF/PSFInterp/PSFProcess.py b/ObservationSim/PSF/PSFInterp/PSFProcess.py index 6bf54466a88071ff1e0bfe3f5946793d3a46dd74..39addcea688d3e2554cf482b021aa3ac64d1cf18 100644 --- a/ObservationSim/PSF/PSFInterp/PSFProcess.py +++ b/ObservationSim/PSF/PSFInterp/PSFProcess.py @@ -6,8 +6,10 @@ import mpi4py.MPI as MPI #sys.path.append("/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_20201222") -import PSFConfig as myConfig -import PSFUtil as myUtil +# import PSFConfig as myConfig +# import PSFUtil as myUtil +import PSFConfig +import PSFUtil def mkdir(path): isExists = os.path.exists(path) @@ -55,8 +57,7 @@ for iccd in range(1, 31): # continue ipsfOutput = psfMatPath + '/psf_{:}_centroidWgt.mat'.format(ipsf) - - psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True) + psfInfo = PSFConfig.LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True) ipsfMat = psfInfo['psfMat'] npix_y, npix_x = ipsfMat.shape @@ -67,7 +68,7 @@ for iccd in range(1, 31): ncut = int(npix_y*0.90) ncut = ncut + ncut%2 #even pixs - img, cx, cy = myUtil.centroidWgt(ipsfMat, nt=ncut) + img, cx, cy = PSFUtil.centroidWgt(ipsfMat, nt=ncut) dcx = cx - npix_x/2 #pixel coords -> global coords dcy = cy - npix_y/2 #pixel coords -> global coords dcx*= psfInfo['pixsize']*1e-3 #5e-3 #pixels -> mm diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc index 76f505d0428a825e0ccdf624a54c0f9f248d0ee2..468f0294e213a98c00d3ca280a3712c028f70af0 100644 Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc and b/ObservationSim/PSF/PSFInterp/__pycache__/PSFConfig.cpython-38.pyc differ diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc index f5be4d2ac0cdf119f4f1b41592eff70f742b2f88..324cb82e4d50e76b23b6a794ffb6f267f8534f8a 100644 Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc and b/ObservationSim/PSF/PSFInterp/__pycache__/PSFInterp.cpython-38.pyc differ diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc index 91ce15d443edbf546e8556af513f5f41117f52c8..14727bc555a3d3e9ba3958ceafb9cf9c05d5b2aa 100644 Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc and b/ObservationSim/PSF/PSFInterp/__pycache__/PSFUtil.cpython-38.pyc differ diff --git a/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc b/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc index 0261ab3e66d147ed69ec746d72e76afec0120863..8fac220f537665f2e6018313955a178434ff72fa 100644 Binary files a/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/PSF/PSFInterp/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/PSF/__init__.py b/ObservationSim/PSF/__init__.py index fc26fbfa0339aaed44805f2461fb3d2c34d479cb..e8823072b0f52b720d37229253a3924c4c7a5e81 100755 --- a/ObservationSim/PSF/__init__.py +++ b/ObservationSim/PSF/__init__.py @@ -1,4 +1,4 @@ from .PSFModel import PSFModel from .PSFGauss import PSFGauss -from .PSFInterp import PSFInterp +from .PSFInterp.PSFInterp import PSFInterp from .FieldDistortion import FieldDistortion \ No newline at end of file diff --git a/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc b/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc index eb4cc3241920e91255bd16ed47010b40ba4f6144..574f7adb3acb80c5dfcddcb139e383443c9a990a 100644 Binary files a/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc and b/ObservationSim/PSF/__pycache__/FieldDistortion.cpython-38.pyc differ diff --git a/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc b/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc index 743b59e2ca6bc80bb165530155391c95f89d8f17..526253816d6c36de939653768ad8c5a8825cc110 100644 Binary files a/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc and b/ObservationSim/PSF/__pycache__/PSFGauss.cpython-38.pyc differ diff --git a/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc b/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc index 69f51a786878cd27bf55ed0eaedfe8f1cf39545b..c7adf737f573d401e0b954effab184402bfb5686 100644 Binary files a/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc and b/ObservationSim/PSF/__pycache__/PSFModel.cpython-38.pyc differ diff --git a/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc b/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc index 130eaafe650cbc3e12701b4e63a5f8f07d61d290..953fc965bdbc4b947834ddd9befdfa08787b2090 100644 Binary files a/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc and b/ObservationSim/PSF/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc b/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc index 758e8ec39517bc7aba7295398f3cda415a577705..4c27078d55d9396bd26e6cd82fa92d62f2c05050 100644 Binary files a/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc and b/ObservationSim/__pycache__/ObservationSim.cpython-38.pyc differ diff --git a/ObservationSim/__pycache__/Pointing.cpython-38.pyc b/ObservationSim/__pycache__/Pointing.cpython-38.pyc index 2dccd0258a9dc7d83260f981da6f71d290870cbe..49945ad93f669fe1b8b86ad00e905644ced0c89d 100644 Binary files a/ObservationSim/__pycache__/Pointing.cpython-38.pyc and b/ObservationSim/__pycache__/Pointing.cpython-38.pyc differ diff --git a/ObservationSim/__pycache__/__init__.cpython-38.pyc b/ObservationSim/__pycache__/__init__.cpython-38.pyc new file mode 100644 index 0000000000000000000000000000000000000000..611fe47a01b30a6df8913fdc8beaf5b0bf349d99 Binary files /dev/null and b/ObservationSim/__pycache__/__init__.cpython-38.pyc differ diff --git a/ObservationSim/__pycache__/_util.cpython-38.pyc b/ObservationSim/__pycache__/_util.cpython-38.pyc index 0df2540b77bb84486a2b4eb24a8201af6902ec97..3f3b03b1491efd33880b8588ae9f3f65d17b535c 100644 Binary files a/ObservationSim/__pycache__/_util.cpython-38.pyc and b/ObservationSim/__pycache__/_util.cpython-38.pyc differ diff --git a/ObservationSim/_util.py b/ObservationSim/_util.py index b2159cbb98224dd6fc023e0397f5a1d156ad76bf..3f839ef29f92096e5a3b8b4df448360068f00d4f 100755 --- a/ObservationSim/_util.py +++ b/ObservationSim/_util.py @@ -1,67 +1,89 @@ import numpy as np import os from datetime import datetime +import argparse + +def parse_args(): + ''' + Parse command line arguments. Many of the following + can be set in the .yaml config file as well. + ''' + parser = argparse.ArgumentParser() + parser.add_argument('config_file', help='.yaml config file for simulation settings.') + parser.add_argument('-c', '--config_dir', help='Directory that houses the ,yaml config file.') + parser.add_argument('-d', '--data_dir', help='Directory that houses the input data.') + parser.add_argument('-w', '--work_dir', help='The path for output.') + return parser.parse_args() 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 + ut = datetime.utcnow() + eye, emo, eda, eho, emi, ese = str(ut.year), str(ut.month), str(ut.day), str(ut.hour), str(ut.minute), str(ut.second) + emse = str(ut.microsecond) + if int(emo)<10: emo = "0%s"%emo + if int(eda)<10: eda = "0%s"%eda + if int(eho)<10: eho = "0%s"%eho + if int(emi)<10: emi = "0%s"%emi + if int(ese)<10: ese = "0%s"%ese - if tt==0: - namekey = "CSST%s%s%sT%s%s%s"%(eye,emo,eda,eho,emi,ese) - elif tt==1: - namekey = "%s-%s-%sT%s:%s:%s.%s"%(eye,emo,eda,eho,emi,ese,emse) - elif tt==2: - namekey = "%s%s%s%s%s%s"%(eye,emo,eda,eho,emi,ese) - else: - raise ValueError("!!! Give a right 'tt' value.") + if tt==0: + namekey = "CSST%s%s%sT%s%s%s"%(eye,emo,eda,eho,emi,ese) + elif tt==1: + namekey = "%s-%s-%sT%s:%s:%s.%s"%(eye,emo,eda,eho,emi,ese,emse) + elif tt==2: + namekey = "%s%s%s%s%s%s"%(eye,emo,eda,eho,emi,ese) + else: + raise ValueError("!!! Give a right 'tt' value.") - return namekey + return namekey -def makeSubDir(path_dict, config): - subImgdir = path_dict["output_img_dir"] + config["mockImgDir"] + '/' - if not os.path.exists(subImgdir): - os.system("mkdir %s"%subImgdir) - imgKey0 = imgName(tt=0) - imgKey1 = imgName(tt=1) - imgKey2 = imgName(tt=2) - subImgdir = subImgdir + imgKey0 + config["reIndex"] + config["reShear"] + "/" - if not os.path.exists(subImgdir): - os.system("mkdir %s"%subImgdir) - return subImgdir, imgKey0, imgKey1, imgKey2 +# def makeSubDir(path_dict, config): +# subImgdir = path_dict["output_img_dir"] + config["mockImgDir"] + '/' +# if not os.path.exists(subImgdir): +# os.system("mkdir %s"%subImgdir) +# imgKey0 = imgName(tt=0) +# imgKey1 = imgName(tt=1) +# imgKey2 = imgName(tt=2) +# subImgdir = subImgdir + imgKey0 + config["reIndex"] + config["reShear"] + "/" +# if not os.path.exists(subImgdir): +# os.system("mkdir %s"%subImgdir) +# return subImgdir, imgKey0, imgKey1, imgKey2 def makeSubDir_PointingList(path_dict, config, pointing_ID=0): - imgDir = path_dict["work_dir"] + config["mockImgDir"] + '/' - if not os.path.exists(imgDir): - os.system("mkdir %s"%imgDir) - # prefix = "MSC_" + config["date_obs"] + config["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') - prefix = "MSC_" + str(pointing_ID).rjust(7, '0') - subImgdir = os.path.join(imgDir, prefix) - if not os.path.exists(subImgdir): - os.system("mkdir %s"%subImgdir) - return subImgdir, prefix + 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') + subImgdir = os.path.join(imgDir, prefix) + if not os.path.exists(subImgdir): + try: + os.makedirs(subImgdir, exist_ok=True) + except OSError: + pass + return subImgdir, prefix def getShearFiled(config, shear_cat_file=None): - if not config["shear_method"] in ["constant", "extra"]: - raise ValueError("Please set a right 'shear_method' parameter.") + if not config["shear_setting"]["shear_type"] in ["constant", "extra"]: + raise ValueError("Please set a right 'shear_method' parameter.") - if config["shear_method"] == "constant": - g1 = config["reduced_g1"] - g2 = config["reduced_g2"] - reduced_shear = np.sqrt(g1**2 + g2**2) - nshear = 1 - # TODO logging - else: - # TODO logging - # if not os.path.exists(shear_cat_file): - # raise ValueError("Cannot find shear catalog file.") - # shearCat = np.loadtxt(shear_cat_file) - # nshear = shearCat.shape[0] - # g1, g2 = shearCat[:, 0], shearCat[:, 1] - pass - return g1, g2, nshear \ No newline at end of file + if config["shear_setting"]["shear_type"] == "constant": + g1 = config["shear_setting"]["reduced_g1"] + g2 = config["shear_setting"]["reduced_g2"] + reduced_shear = np.sqrt(g1**2 + g2**2) + nshear = 1 + # TODO logging + else: + # TODO logging + if not os.path.exists(shear_cat_file): + raise ValueError("Cannot find shear catalog file.") + try: + shearCat = np.loadtxt(shear_cat_file) + nshear = shearCat.shape[0] + g1, g2 = shearCat[:, 0], shearCat[:, 1] + except: + print("Failed to the shear catalog file.") + print("Setting to no shear.") + g1, g2 = 0., 0. + return g1, g2, nshear \ No newline at end of file diff --git a/ObservationSim/run.pbs b/ObservationSim/run.pbs index 34274e334c5a5211e9247eda8472f25e1ea05d00..e199d491aafcf8c9cc65506582f32ad98578fce7 100755 --- a/ObservationSim/run.pbs +++ b/ObservationSim/run.pbs @@ -12,11 +12,12 @@ #PBS -q batch #PBS -u fangyuedong -NP=54 +NP=40 date echo $NP # mpirun -np $NP --oversubscribe -H comput101 python /public/home/fangyuedong/CSST/ObservationSim/runExposure.py +# python /public/home/fangyuedong/sim_code_release/CSSTSim/ObservationSim/preprocess.py +# mpirun -np $NP python /public/home/fangyuedong/sim_code_release/CSSTSim/ObservationSim/runExposure.py -python /public/home/fangyuedong/test/CSST/ObservationSim/preprocess.py -mpirun -np $NP python /public/home/fangyuedong/test/CSST/ObservationSim/runExposure.py +mpirun -np $NP python /public/home/fangyuedong/sim_code_release/CSSTSim/ObservationSim/runExposure.py config_sim.yaml -c /public/home/fangyuedong/sim_code_release/CSSTSim/config diff --git a/ObservationSim/runExposure.py b/ObservationSim/runExposure.py index a90d30e0cd79e736c4be0c01b7c8929dde640ee5..8ae357c8f396ee7de17c221aea2ea0a2909e0e91 100755 --- a/ObservationSim/runExposure.py +++ b/ObservationSim/runExposure.py @@ -1,12 +1,135 @@ from ObservationSim import Observation -from Pointing import * +from _util import parse_args +from datetime import datetime import os import numpy as np import galsim +import yaml import gc gc.enable() +def Pointing(config, pointing_filename, data_dir): + pointing_file = os.path.join(data_dir, pointing_filename) + f = open(pointing_file, 'r') + for _ in range(1): + header = f.readline() + iline = 0 + pRA = [] + pDEC = [] + for line in f: + line = line.strip() + columns = line.split() + pRA.append(float(columns[0])) + pDEC.append(float(columns[1])) + f.close() + pRA = np.array(pRA) + pDEC = np.array(pDEC) + + # Create calibration pointings + # NOTE: temporary implementation + ncal = config['obs_setting']['np_cal'] + pointing_type = ['MS']*len(pRA) + pRA = np.append([pRA[0]]*ncal, pRA) + pDEC = np.append([pDEC[0]]*ncal, pDEC) + pointing_type = ['CAL']*ncal + pointing_type + + # Calculate starting time(s) + # NOTE: temporary implementation + t0 = datetime(2021, 5, 25, 12, 0, 0) + t = datetime.timestamp(t0) + timestamp_obs = [] + delta_t = 10 # Time elapsed between exposures (minutes) + for i in range(len(pointing_type)): + timestamp_obs.append(t) + if pointing_type[i] == 'CAL': + t += 3 * delta_t * 60 # 3 calibration exposure + elif pointing_type[i] == 'MS': + t += delta_t * 60 + timestamp_obs = np.array(timestamp_obs) + pointing_type = np.array(pointing_type) + return pRA, pDEC, timestamp_obs, pointing_type + +def MakeDirectories(work_dir, run_name, nPointings, pRange=None): + if not os.path.exists(work_dir): + try: + os.makedirs(work_dir, exist_ok=True) + except OSError: + pass + imgDir = os.path.join(work_dir, run_name) + if not os.path.exists(imgDir): + try: + os.makedirs(imgDir, exist_ok=True) + except OSError: + pass + prefix = "MSC_" + for pointing_ID in range(nPointings): + if pRange is not None: + if pointing_ID not in pRange: + continue + 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 + +def runSim(): + """ + Main simulation call. + """ + args = parse_args() + if args.config_dir is None: + args.config_dir = '' + args.config_dir = os.path.abspath(args.config_dir) + args.config_file = os.path.join(args.config_dir, args.config_file) + with open(args.config_file, "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) + config["obs_setting"]["image_rot"] = float(config["obs_setting"]["image_rot"])*galsim.degrees + + if args.data_dir is not None: + config['data_dir'] = args.data_dir + if args.work_dir is not None: + config['work_dir'] = args.work_dir + + pRA, pDEC, timestamp_obs, pointing_type = Pointing(config=config, pointing_filename=config['pointing_file'], data_dir=config['data_dir']) + + MakeDirectories(work_dir=config['work_dir'], run_name=config['run_name'], nPointings=len(pRA), pRange=config['obs_setting']['run_pointings']) + if "run_chips" in config["obs_setting"]: + run_chips = config["obs_setting"]["run_chips"] + else: + run_chips = None + + # from Config import ConfigDir + # path_dict = ConfigDir(config=config, work_dir=config["work_dir"], data_dir=config["data_dir"]) + # for key, value in path_dict.items(): + # print (key + " : " + str(value)) + + obs = Observation(config=config, work_dir=config['work_dir'], data_dir=config['data_dir']) + if config["pointing_file"] is None: + obs.runExposure(chips=run_chips) + else: + obs.runExposure_MPI_PointingList( + ra_cen=pRA, + dec_cen=pDEC, + pRange=config["obs_setting"]["run_pointings"], + timestamp_obs=timestamp_obs, + pointing_type=pointing_type, + exptime=config["obs_setting"]["exp_time"], + use_mpi=config["run_option"]["use_mpi"], + chips=run_chips + ) + print("run finished") + +if __name__=='__main__': + runSim() + ############################################# # Testing run one exposure (NOT using MPI) # ipoint = 2 @@ -18,6 +141,6 @@ gc.enable() ############################################# # Testing run pointing list (using MPI) -obs = Observation(work_dir=work_dir, data_dir=data_dir) -# obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type) -obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type, exptime=exp_time) +# obs = Observation(work_dir=work_dir, data_dir=data_dir) +# # obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type) +# obs.runExposure_MPI_PointingList(ra_cen=pRA, dec_cen=pDEC, pRange=pRange, timestamp_obs=timestamp_obs, pointing_type=pointing_type, exptime=exp_time) diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/config/config_sim.yaml b/config/config_sim.yaml new file mode 100644 index 0000000000000000000000000000000000000000..f8e51f45f3d7e5d01db81e78fe4053295896fbb0 --- /dev/null +++ b/config/config_sim.yaml @@ -0,0 +1,209 @@ +--- +############################################### +# +# Configuration file for CSST simulation +# CSST-Sim Group, 2021/10/07, version 0.3 +# +############################################### + +# 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: "/public/home/fangyuedong/sim_code_release/CSST/test/" +work_dir: "/public/home/fangyuedong/sim_code_release/CSSTSim/workplace/" +data_dir: "/data/simudata/CSSOSDataProductsSims/data/" +run_name: "TEST" + +# (Optional) a file of point list +pointing_file: "pointing10_20210202.dat" + +# Whether to use MPI +run_option: + use_mpi: YES + n_cores: 40 + +# Output catalog only? +# If yes, no imaging simulation will run +out_cat_only: NO + +# Only simulate stars? +star_only: NO + +# Only simulate galaxies? +galaxy_only: NO + +############################################### +# Observation setting +############################################### +obs_setting: + + # Options for survey types: + # "Photometric": simulate photometric chips only + # "Spectroscoplic": simulate slitless spectroscopic chips only + # "All": simulate full focal plane + survey_type: "Photometric" + + # Exposure time [seconds] + exp_time: 150. + + # Observation starting date & time + # (Subject to change) + date_obs: "210525" # [yymmdd] + time_obs: "120000" # [hhmmss] + + # Default Pointing [degrees] + # Note: NOT valid when a pointing list file is specified + ra_center: 60. + dec_center: -40. + + # Image rotation [degree] + image_rot: -113.4333 + + # Number of calibration pointings + # Note: only valid when a pointing list is specified + np_cal: 1 + + # (Optional) only run specific pointing(s). + # Note: only valid when a pointing list is specified + run_pointings: [0, 1, 2, 3, 4] + + # (Optional) only run specific chip(s) + # Note: for all pointings + run_chips: [6, 25] + +############################################### +# Input path setting +############################################### + +# Default path settings for WIDE survey simulation +input_path: + cat_dir: "Catalog_20210126/" + star_cat: "MMW_Cluster_D20_healpix_21.hdf5" + galaxy_cat: "galaxyCats_r_10.0_healpix.hdf5" + +SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Templates/Galaxy/" + +# TODO: should the following path settings be made hidden from user? +Efficiency_curve_path: + filter_eff: "Filters/" + ccd_eff: "Filter_CCD_Mirror/ccd/" + mirror_eff: "Filter_CCD_Mirror/mirror_ccdnote.txt" + +CR_data_path: "CRdata/" +sky_data_path: "skybackground/sky_emiss_hubble_50_50_A.dat" + +SLS_path: + SLS_conf: "CONF" + SLS_norm: "normalize_filter" + + +############################################### +# PSF setting +############################################### +psf_setting: + + # Which PSF model to use: + # "Gauss": simple gaussian profile + # "Interp": Interpolated PSF from sampled ray-tracing data + psf_model: "Interp" + + # PSF size [arcseconds] + # radius of 80% energy encircled + # NOTE: only valid for "Gauss" PSF + psf_rcont: 0.15 + + # path to PSF data + # NOTE: only valid for "Interp" PSF + psf_dir: "csstPSFdata/CSSOS_psf_20210108/CSST_psf_ciomp_2p5um_cycle3_ccr90_proc" + + # sigma_spin: 0.0 # psf spin? + +############################################### +# Shear setting +############################################### + +shear_setting: + # Options to generate mock shear field: + # "constant": all galaxies are assigned a constant reduced shear + # "catalog": from catalog (not available yet) + # "extra": from seprate file + shear_type: "constant" + + # For constant shear filed + reduced_g1: 0.026 + reduced_g2: 0.015 + + # Representation of the shear vector? + reShear: "E" + + # rotate galaxy ellipticity + rotateEll: 0. # [degree] + + # Extra shear catalog + # (currently not used) + # shear_cat: "mockShear.cat" + +############################################### +# Instrumental effects setting +############################################### +ins_effects: + # switches + field_dist: ON # Whether to add field distortions + add_back: ON # Whether to add sky background + 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 + shutter_effect: ON # Whether to add shutter effect + flat_fielding: ON # Whether to add flat-fielding effect + prnu_effect: ON # Whether to add PRNU effect + non_linear: OFF # 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 + 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 + + # values + dark_exptime: 300 # Exposure time for dark current frames [seconds] + flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] + readout_time: 40 # The read-out time for each channel [seconds] + df_strength: 2.3 # Sillicon sensor diffusion strength + bias_level: 500 # bias level [e-/pixel] + gain: 1.1 # Gain + full_well: 90000 # Full well depth [e-] + NBias: 1 # Number of bias frames to be exported for each exposure + NDark: 1 # Number of dark frames to be exported for each exposure + NFlat: 1 # Number of flat frames to be exported for each exposure + +############################################### +# Output options +############################################### +output_setting: + readout16: OFF # Whether to export as 16 channels (subimages) with pre- and over-scan + shutter_output: OFF # Whether to export shutter effect 16-bit image + bias_output: ON # Whether to export bias frames + dark_output: ON # Whether to export the combined dark current files + flat_output: ON # Whether to export the combined flat-fielding files + prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files + +############################################### +# Random seeds +############################################### +random_seeds: + seed_Av: 121212 # Seed for generating random intrinsic extinction + seed_poisson: 20210601 # Seed for Poisson noise + seed_CR: 20210317 # Seed for generating random cosmic ray maps + seed_flat: 20210101 # Seed for generating random flat fields + seed_prnu: 20210102 # Seed for photo-response non-uniformity + seed_gainNonUniform: 20210202 # Seed for gain nonuniformity + seed_biasNonUniform: 20210203 # Seed for bias nonuniformity + seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity + 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 diff --git a/config/test_yaml.py b/config/test_yaml.py new file mode 100644 index 0000000000000000000000000000000000000000..e44a5ce19870788adb78214236a2370cab51ef31 --- /dev/null +++ b/config/test_yaml.py @@ -0,0 +1,18 @@ +import yaml + +if __name__ == '__main__': + + stream = open("config_sim.yaml", 'r') + # dictionary = yaml.load(stream) + dictionary = yaml.load(stream, Loader=yaml.SafeLoader) + # print(dictionary) + for key, value in dictionary.items(): + print (key + " : " + str(value)) + + with open("config_sim.yaml", "r") as stream: + try: + dictionary = yaml.safe_load(stream) + for key, value in dictionary.items(): + print (key + " : " + str(value)) + except yaml.YAMLError as exc: + print(exc) \ No newline at end of file diff --git a/setup.py b/setup.py new file mode 100644 index 0000000000000000000000000000000000000000..e750979c8d4a970aa310e42b99d5b0102991c8fb --- /dev/null +++ b/setup.py @@ -0,0 +1,19 @@ +from setuptools import setup, find_packages + +setup(name='CSSTSim', + version='0.3', + packages=find_packages(), + # install_requires=[ + # 'numpy>=1.18.5', + # 'galsim>=2.2.4', + # 'yaml>=5.3.1', + # 'astropy>=4.0.1', + # 'scipy>=1.5.0', + # 'mpi4py>=3.0.3', + # 'sep>=1.0.3', + # 'healpy>=1.14.0', + # 'h5py>=2.10.0', + # 'Cython>=0.29.21' + # 'numba>=0.50.1' + # ] +) \ No newline at end of file diff --git a/tests/.coverage b/tests/.coverage new file mode 100644 index 0000000000000000000000000000000000000000..27384e3d018368c4995104be55393efb96727a2d --- /dev/null +++ b/tests/.coverage @@ -0,0 +1 @@ +}q(U collectorqUcoverage v3.6b3qUlinesq}u. \ No newline at end of file diff --git a/tests/PSFMatsEll_coverage.py b/tests/PSFMatsEll_coverage.py new file mode 100644 index 0000000000000000000000000000000000000000..de2d924445cecd5fccff14d613d5e4484e893a31 --- /dev/null +++ b/tests/PSFMatsEll_coverage.py @@ -0,0 +1,209 @@ +import unittest + +import sys,os,math +from itertools import islice + +#import mpi4py.MPI as MPI + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +mpl.use('Agg') + +import scipy.io +from scipy import ndimage + +#sys.path.append("/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_v4_20210326") +sys.path.append("../") +from ObservationSim.PSF.PSFInterp import PSFConfig as myConfig +#import PSFUtil as myUtil + +NPSF = 400 +############################## + +###计算PSF椭率### +def psfSecondMoments(psfMat, cenX, cenY, pixSize=1): + apr = 0.5 #arcsec, 0.5角秒内测量 + fl = 28. #meters + pxs = 2.5 #microns + apr = np.deg2rad(apr/3600.)*fl*1e6 + apr = apr/pxs + apr = np.int(np.ceil(apr)) + + I = psfMat + ncol = I.shape[1] + nrow = I.shape[0] + w = 0.0 + w11 = 0.0 + w12 = 0.0 + w22 = 0.0 + for icol in range(ncol): + for jrow in range(nrow): + x = icol*pixSize - cenX + y = jrow*pixSize - cenY + rr = np.sqrt(x*x + y*y) + wgt= 0.0 + if rr <= apr: + wgt = 1.0 + w += I[jrow, icol]*wgt + w11 += x*x*I[jrow, icol]*wgt + w12 += x*y*I[jrow, icol]*wgt + w22 += y*y*I[jrow, icol]*wgt + w11 /= w + w12 /= w + w22 /= w + sz = w11 + w22 + e1 = (w11 - w22)/sz + e2 = 2.0*w12/sz + + return sz, e1, e2 + +############################## +############################## +''' +def assignTasks(npsf, NTasks, ThisTask): + npsfPerTasks = int(npsf/NTasks) + iStart= 0 + npsfPerTasks*ThisTask + iEnd = npsfPerTasks + npsfPerTasks*ThisTask + if ThisTask == NTasks: + iEnd = npsf + return iStart, iEnd +''' + +#def test_psfEll(iccd, iwave, psfPath, ThisTask, NTasks): +def test_psfEll(iccd, iwave, psfPath): + nccd = 30 + npsf = NPSF + + #iStart, iEnd = assignTasks(npsf, NTasks, ThisTask) + + imx = np.zeros(npsf) + imy = np.zeros(npsf) + psf_e1 = np.zeros(npsf) + psf_e2 = np.zeros(npsf) + psf_sz = np.zeros(npsf) + + #for ipsf in range(iStart+1, iEnd+1): + for ipsf in range(1, 401): + if ipsf != 1: + continue + print('ipsf-{:}'.format(ipsf), end='\r') + psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True, PSFCentroidWgt=True, VPSF=False) + imx[ipsf-1] = psfInfo['image_x']+psfInfo['centroid_x'] + imy[ipsf-1] = psfInfo['image_y']+psfInfo['centroid_y'] + psfMat = psfInfo['psfMat'] + cenX = 256 + cenY = 256 + sz, e1, e2 = psfSecondMoments(psfMat, cenX, cenY, pixSize=1) + psf_e1[ipsf-1] = e1 + psf_e2[ipsf-1] = e2 + psf_sz[ipsf-1] = sz + print('test:' ,sz, e1, e2) + + ####### + #comm.barrier() + #imx = comm.allreduce(imx, op=MPI.SUM) + #imy = comm.allreduce(imy, op=MPI.SUM) + + #psf_e1 = comm.allreduce(psf_e1, op=MPI.SUM) + #psf_e2 = comm.allreduce(psf_e2, op=MPI.SUM) + #psf_sz = comm.allreduce(psf_sz, op=MPI.SUM) + + #comm.barrier() + #if ThisTask == 0: + # arr = [imx, imy, psf_e1, psf_e2, psf_sz] + # np.save('/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_v4_20210326/test4report/data/psfEll20_{:}_{:}'.format(iccd, iwave), arr) + + +def test_psfEllPlot(OVERPLOT=False): + #if ThisTask == 0: + if True: + prefix = 'psfEll30' + iccd = 1 + iwave= 1 + data = np.load('data/'+prefix+'_1_1.npy') + imx= data[0] + imy= data[1] + psf_e1 = data[2] + psf_e2 = data[3] + print(np.shape(imx)) + npsf = np.shape(imx)[0] + + plt.cla() + plt.close("all") + fig = plt.figure(figsize=(12,12)) + plt.plot(imx, imy, 'r.') + plt.savefig('figs/psfPos.pdf') + + ####### + + fig = plt.figure(figsize=(12, 12)) + plt.subplots_adjust(wspace=0.1, hspace=0.1) + ax = plt.subplot(1, 1, 1) + for ipsf in range(npsf): + plt.plot(imx[ipsf], imy[ipsf], 'r.') + + ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 + ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) + ell *= 15 + lcos = ell*np.cos(ang) + lsin = ell*np.sin(ang) + plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'r', lw=2) + ########### + ang = 0. + ell = 0.05 + ell*= 15 + lcos = ell*np.cos(ang) + lsin = ell*np.sin(ang) + plt.plot([imx[898]-lcos, imx[898]+lcos],[imy[898]+5.-lsin, imy[898]+5.+lsin],'k', lw=2) + plt.annotate('{:}'.format(ell/15), (imx[898]-2., imy[898]+6.), xycoords='data', fontsize=10) + plt.xlabel('CCD X (mm)') + plt.ylabel('CCD Y (mm)') + + if OVERPLOT == True: + prefix = 'psfEll20' + data = np.load('data/'+prefix+'_1_1.npy') + imx= data[0] + imy= data[1] + psf_e1 = data[2] + psf_e2 = data[3] + + npsf = np.shape(imx)[0] + for ipsf in range(npsf): + plt.plot(imx[ipsf], imy[ipsf], 'b.') + + ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 + ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) + ell *= 15 + lcos = ell*np.cos(ang) + lsin = ell*np.sin(ang) + plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'b', lw=2) + + plt.gca().set_aspect(1) + if OVERPLOT == True: + prefix = 'psfEllOP' + plt.savefig('figs/'+prefix+'_iccd{:}.pdf'.format(iccd)) + + +class PSFMatsEll_coverage(unittest.TestCase): + def test_psfEll_(self): + #comm = MPI.COMM_WORLD + #ThisTask = comm.Get_rank() + #NTasks = comm.Get_size() + print('#####haha#####') + + iccd = 1 + iwave= 1 + psfPath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_20210326/CSST_psf_ciomp_20x20field' + #test_psfEll(iccd, iwave, psfPath, ThisTask, NTasks) + test_psfEll(iccd, iwave, psfPath) + + test_psfEllPlot(OVERPLOT=True) + + +############################## +############################## +############################## +if __name__=='__main__': + unittest.main() + diff --git a/tests/PSFMatsIDW_coverage.py b/tests/PSFMatsIDW_coverage.py new file mode 100644 index 0000000000000000000000000000000000000000..c63ce096dd3fa7a13cdcc43350b056f5d3e2d7f8 --- /dev/null +++ b/tests/PSFMatsIDW_coverage.py @@ -0,0 +1,468 @@ +import unittest + +import sys,os,math +from itertools import islice + +import mpi4py.MPI as MPI + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +mpl.use('Agg') + +import scipy.io +from scipy import ndimage + +#sys.path.append("/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_v4_20210326") +# sys.path.append("../") +from ObservationSim.PSF.PSFInterp import PSFConfig as myConfig +from ObservationSim.PSF.PSFInterp import PSFUtil as myUtil + +NPSF = 400 +############################## + +###计算PSF椭率### +def psfSecondMoments(psfMat, cenX, cenY, pixSize=1): + apr = 0.5 #arcsec, 0.5角秒内测量 + fl = 28. #meters + pxs = 2.5 #microns + apr = np.deg2rad(apr/3600.)*fl*1e6 + apr = apr/pxs + apr = np.int(np.ceil(apr)) + + I = psfMat + ncol = I.shape[1] + nrow = I.shape[0] + w = 0.0 + w11 = 0.0 + w12 = 0.0 + w22 = 0.0 + for icol in range(ncol): + for jrow in range(nrow): + x = icol*pixSize - cenX + y = jrow*pixSize - cenY + rr = np.sqrt(x*x + y*y) + wgt= 0.0 + if rr <= apr: + wgt = 1.0 + w += I[jrow, icol]*wgt + w11 += x*x*I[jrow, icol]*wgt + w12 += x*y*I[jrow, icol]*wgt + w22 += y*y*I[jrow, icol]*wgt + w11 /= w + w12 /= w + w22 /= w + sz = w11 + w22 + e1 = (w11 - w22)/sz + e2 = 2.0*w12/sz + + return sz, e1, e2 + +############################## +############################## +''' +def assignTasks(npsf, NTasks, ThisTask): + npsfPerTasks = int(npsf/NTasks) + iStart= 0 + npsfPerTasks*ThisTask + iEnd = npsfPerTasks + npsfPerTasks*ThisTask + if ThisTask == NTasks: + iEnd = npsf + return iStart, iEnd +''' + +#def test_psfIDW(iccd, iwave, psfPath, ThisTask, NTasks): +def test_psfIDW(iccd, iwave, psfPath): + nccd = 30 + npsfA = 900 + npsfB = 400 + + #iStart, iEnd = assignTasks(400, NTasks, ThisTask) + + psfPathA = psfPath+'_30x30field' + psfPathB = psfPath+'_20x20field' + + imxA = np.zeros(npsfA) + imyA = np.zeros(npsfA) + psfA = np.zeros([npsfA, 512, 512]) + + imxB = np.zeros(npsfB) + imyB = np.zeros(npsfB) + psfB = np.zeros([npsfB, 512, 512]) + + #for ipsf in range(iStart+1, iEnd+1): + for ipsf in range(1, npsfA+1): + print('ipsfA:', ipsf, end='\r', flush=True) + psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPathA, InputMaxPixelPos=True, PSFCentroidWgt=True, VPSF=False) + imxA[ipsf-1] = psfInfo['image_x'] + psfInfo['centroid_x'] + imyA[ipsf-1] = psfInfo['image_y'] + psfInfo['centroid_y'] + psfA[ipsf-1, :, :] = psfInfo['psfMat'] + + for ipsf in range(1, npsfB+1): + print('ipsfB:', ipsf, end='\r', flush=True) + psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPathB, InputMaxPixelPos=True, PSFCentroidWgt=True, VPSF=False) + imxB[ipsf-1] = psfInfo['image_x'] + psfInfo['centroid_x'] + imyB[ipsf-1] = psfInfo['image_y'] + psfInfo['centroid_y'] + psfB[ipsf-1, :, :] = psfInfo['psfMat'] + + #myConfig.psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False) + for ipsf in range(npsfB): + print('ipsf:', ipsf, end='\r', flush=True) + px = imxB[ipsf] + py = imyB[ipsf] + cen_col = imxA + cen_row = imyA + PSFMat = psfA + psfIDW = myConfig.psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False) + np.save('figs/psfIDW/psfIDW_{:}_{:}_{:}'.format(iccd, iwave, ipsf+1), psfIDW) + + +#def test_psfEll(iccd, iwave, ThisTask, NTasks): +def test_psfEll(iccd, iwave): + nccd = 30 + npsf = 400 + + #iStart, iEnd = assignTasks(npsf, NTasks, ThisTask) + + #imx = np.zeros(npsf) + #imy = np.zeros(npsf) + psf_e1 = np.zeros(npsf) + psf_e2 = np.zeros(npsf) + psf_sz = np.zeros(npsf) + + #for ipsf in range(iStart+1, iEnd+1): + for ipsf in range(1, 401): + if ipsf > 1: + continue + print('ipsf-{:}'.format(ipsf), end='\r') + #psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPath, InputMaxPixelPos=True, PSFCentroidWgt=True, VPSF=False) + #psfMat = psfInfo['psfMat'] + #imx[ipsf-1] = psfInfo['image_x']+psfInfo['centroid_x'] + #imy[ipsf-1] = psfInfo['image_y']+psfInfo['centroid_y'] + psfMat = np.load('figs/psfIDW/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf)) #psfInfo['psfMat'] + cenX = 256 + cenY = 256 + sz, e1, e2 = psfSecondMoments(psfMat, cenX, cenY, pixSize=1) + psf_e1[ipsf-1] = e1 + psf_e2[ipsf-1] = e2 + psf_sz[ipsf-1] = sz + + ####### + #comm.barrier() + #imx = comm.allreduce(imx, op=MPI.SUM) + #imy = comm.allreduce(imy, op=MPI.SUM) + + #psf_e1 = comm.allreduce(psf_e1, op=MPI.SUM) + #psf_e2 = comm.allreduce(psf_e2, op=MPI.SUM) + #psf_sz = comm.allreduce(psf_sz, op=MPI.SUM) + + #comm.barrier() + #if ThisTask == 0: + # arr = [psf_e1, psf_e2, psf_sz] + # np.save('/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_v4_20210326/test4report/data/psfEll20IDW_{:}_{:}'.format(iccd, iwave), arr) + +''' +def test_psfResidualCalc(iccd, iwave, ipsf, psfPath): + psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPath+'_20x20field', InputMaxPixelPos=True, PSFCentroidWgt=True, VPSF=False) + psfMatORG = psfInfo['psfMat'] + psfMatIDW = np.load('/public/home/weichengliang/lnData/CSST_new_framwork/csstPSF_v4_20210326/test4report/figs/psfIDW/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf)) + + _,tREE80 = myUtil.psfEncircle(psfMatORG, fraction=0.8, psfSampleSizeInMicrons=2.5, focalLengthInMeters=28, cenPix=None) + tREE80_pix = np.int32(np.ceil(tREE80/(0.074/2/2))[0]) + #print(tREE80, np.ceil(tREE80/(0.074/2/2)), tREE80_pix) + + timg0 = psfMatORG[256-tREE80_pix:256+tREE80_pix, 256-tREE80_pix:256+tREE80_pix] + timg1 = psfMatIDW[256-tREE80_pix:256+tREE80_pix, 256-tREE80_pix:256+tREE80_pix] + #print("residual::", np.max((timg1-timg0)/timg0), np.min((timg1-timg0)/timg0), np.mean((timg1-timg0)/timg0)) + + return np.mean((timg1-timg0)/timg0) +''' + +def test_psfResidualPlot(iccd, iwave, ipsf, psfPath): + psfInfo = myConfig.LoadPSF(iccd, iwave, ipsf, psfPath+'_20x20field', InputMaxPixelPos=True, PSFCentroidWgt=True, VPSF=False) + psfMatORG = psfInfo['psfMat'] + psfMatIDW = np.load('figs/psfIDW/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf)) + + npix = psfMatORG.shape[0] + pixCutEdge= int(npix/2-15) + + img0 = psfMatORG[pixCutEdge:npix-pixCutEdge, pixCutEdge:npix-pixCutEdge] + img1 = psfMatIDW[pixCutEdge:npix-pixCutEdge, pixCutEdge:npix-pixCutEdge] + imgX = (img1 - img0)/img0 + + img0 = np.log10(img0) + img1 = np.log10(img1) + imgX = np.log10(np.abs(imgX)) + + fig = plt.figure(figsize=(18,4)) + ax = plt.subplot(1,3,1) + plt.imshow(img0, origin='lower', vmin=-7, vmax=-1.3) + plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--') + plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--') + plt.annotate('ORG', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15) + cticks=[-7, -6, -5, -4, -3, -2, -1] + cbar = plt.colorbar(ticks=cticks) + cbar.ax.set_yticklabels(['$10^{-7}$', '$10^{-6}$', '$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$']) + print(img0.min(), img0.max()) + + + ax = plt.subplot(1,3,2) + plt.imshow(img1, origin='lower', vmin=-7, vmax=-1.3) + plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--') + plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--') + plt.annotate('IDW', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15) + cticks=[-7, -6, -5, -4, -3, -2, -1] + cbar = plt.colorbar(ticks=cticks) + cbar.ax.set_yticklabels(['$10^{-7}$', '$10^{-6}$', '$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$']) + print(img1.min(), img1.max()) + + + ax = plt.subplot(1,3,3) + plt.imshow(imgX, origin='lower', vmin =-3, vmax =np.log10(3e-1)) + plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--') + plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--') + #plt.annotate('(IDW-ORG)/ORG', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15) + cticks=[-5, -4, -3, -2, -1] + cbar = plt.colorbar(ticks=cticks) + cbar.ax.set_yticklabels(['$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$']) + + print(np.max((psfMatORG-psfMatIDW))) + plt.savefig('figs/psfResidual_iccd{:}.pdf'.format(iccd)) + + +def test_psfEllPlot(OVERPLOT=False): + #if ThisTask == 0: + if True: + prefix = 'psfEll20' + iccd = 1 + iwave= 1 + data = np.load('data/'+prefix+'_1_1.npy') + imx= data[0] + imy= data[1] + psf_e1 = data[2] + psf_e2 = data[3] + print(np.shape(imx)) + npsf = np.shape(imx)[0] + + plt.cla() + plt.close("all") + fig = plt.figure(figsize=(12,12)) + plt.plot(imx, imy, 'r.') + plt.savefig('figs/psfPos.pdf') + + ####### + + fig = plt.figure(figsize=(12, 12)) + plt.subplots_adjust(wspace=0.1, hspace=0.1) + ax = plt.subplot(1, 1, 1) + for ipsf in range(npsf): + plt.plot(imx[ipsf], imy[ipsf], 'b.') + + ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 + ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) + ell *= 15 + lcos = ell*np.cos(ang) + lsin = ell*np.sin(ang) + plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'b', lw=2) + ########### + ang = 0. + ell = 0.05 + ell*= 15 + lcos = ell*np.cos(ang) + lsin = ell*np.sin(ang) + #plt.plot([imx[898]-lcos, imx[898]+lcos],[imy[898]+5.-lsin, imy[898]+5.+lsin],'k', lw=2) + #plt.annotate('{:}'.format(ell/15), (imx[898]-2., imy[898]+6.), xycoords='data', fontsize=10) + plt.xlabel('CCD X (mm)') + plt.ylabel('CCD Y (mm)') + + if OVERPLOT == True: + prefix = 'psfEll20IDW' + data = np.load('data/'+prefix+'_1_1.npy') + #imx= data[0] + #imy= data[1] + psf_e1 = data[0] + psf_e2 = data[1] + + npsf = np.shape(imx)[0] + for ipsf in range(npsf): + #plt.plot(imx[ipsf], imy[ipsf], 'r.') + + ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 + ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) + ell *= 15 + lcos = ell*np.cos(ang) + lsin = ell*np.sin(ang) + plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'r', lw=1) + + plt.gca().set_aspect(1) + if OVERPLOT == True: + prefix = 'psfEllOPIDW' + plt.savefig('figs/'+prefix+'_iccd{:}.pdf'.format(iccd)) + +''' +def test_psfdEllPlot(): + if ThisTask == 0: + prefix = 'psfEll20' + iccd = 1 + iwave= 1 + data = np.load('data/'+prefix+'_1_1.npy') + imx= data[0] + imy= data[1] + psf_e1 = data[2] + psf_e2 = data[3] + psf_sz = data[4] + print(np.shape(imx)) + npsf = np.shape(imx)[0] + + ellX = np.sqrt(psf_e1**2 + psf_e2**2) + angX = np.arctan2(psf_e2, psf_e1)/2 + angX = np.rad2deg(angX) + szX = psf_sz + + ############################## + prefix = 'psfEll20IDW' + data = np.load('data/'+prefix+'_1_1.npy') + #imx= data[0] + #imy= data[1] + psf_e1 = data[0] + psf_e2 = data[1] + psf_sz = data[2] + + ellY = np.sqrt(psf_e1**2 + psf_e2**2) + angY = np.arctan2(psf_e2, psf_e1)/2 + angY = np.rad2deg(angY) + szY = psf_sz + + ############################## + fig=plt.figure(figsize=(15, 4)) + ax = plt.subplot(1,3,1) + plt.hist(ellX, bins=20, color='b', alpha=0.5) + plt.hist(ellY, bins=20, color='r', alpha=0.5) + plt.xlabel('$\epsilon$') + plt.ylabel('PDF') + + ax = plt.subplot(1,3,2) + plt.hist((ellY-ellX)/ellX, bins=20, color='r', alpha=0.5) + plt.xlabel('$(\epsilon_{\\rm IDW}-\epsilon_{\\rm ORG})/\epsilon_{\\rm ORG}$') + plt.ylabel('PDF') + + + ax = plt.subplot(1,3,3) + plt.hist((angY-angX)/angX, bins=20, color='r', alpha=0.5, range=[-0.1, 0.1]) + plt.xlabel('$(\\alpha_{\\rm IDW}-\\alpha_{\\rm ORG})/\\alpha_{\\rm ORG}$') + plt.ylabel('PDF') + plt.savefig('figs/psfEllOPIDWPDF.pdf') + + fig=plt.figure(figsize=(4, 4)) + plt.hist((szY-szX)/szX, bins=20, color='r', alpha=0.5) + plt.xlabel('$(R_{\\rm IDW}-R_{\\rm ORG})/R_{\\rm ORG}$') + plt.ylabel('PDF') + plt.savefig('figs/psfEllOPIDWPDF_dsz.pdf') +''' + +def test_psfdEllabsPlot(iccd): + #if ThisTask == 0: + if True: + prefix = 'psfEll20' + #iccd = 1 + #iwave= 1 + data = np.load('data/'+prefix+'_{:}_1.npy'.format(iccd)) + imx= data[0] + imy= data[1] + psf_e1 = data[2] + psf_e2 = data[3] + psf_sz = data[4] + print(np.shape(imx)) + npsf = np.shape(imx)[0] + + ellX = np.sqrt(psf_e1**2 + psf_e2**2) + angX = np.arctan2(psf_e2, psf_e1)/2 + angX = np.rad2deg(angX) + szX = psf_sz + + ############################## + prefix = 'psfEll20IDW' + data = np.load('data/'+prefix+'_{:}_1.npy'.format(iccd)) + #imx= data[0] + #imy= data[1] + psf_e1 = data[0] + psf_e2 = data[1] + psf_sz = data[2] + + ellY = np.sqrt(psf_e1**2 + psf_e2**2) + angY = np.arctan2(psf_e2, psf_e1)/2 + angY = np.rad2deg(angY) + szY = psf_sz + + ############################## + fig=plt.figure(figsize=(6, 5)) + grid = plt.GridSpec(3,1,left=0.15, right=0.95, wspace=None, hspace=0.02) + #plt.subplots_adjust(left=None,bottom=None,right=None,top=None,wspace=None,hspace=0.02) + ax = plt.subplot(grid[0:2,0]) + plt.plot([0.01,0.1],[0.01,0.1], 'k--', lw=1. ) + plt.scatter(ellX, ellY, color='b', alpha=1., s=3., edgecolors='None') + plt.xlim([0.015, 0.085]) + plt.ylim([0.015, 0.085]) + plt.ylabel('$\epsilon_{\\rm IDW}$') + plt.gca().axes.get_xaxis().set_visible(False) + + + ax = plt.subplot(grid[2,0]) + plt.plot([0.015,0.085],[0.,0.], 'k--', lw=1. ) + plt.scatter(ellX, (ellY-ellX), color='b', s=3., edgecolors='None') + plt.xlim([0.015, 0.085]) + plt.ylim([-0.0018, 0.0018]) + plt.xlabel('$\epsilon_{\\rm ORG}$') + plt.ylabel('$\Delta$') + plt.savefig('figs/psfEllOPIDWPDF_{:}.pdf'.format(iccd)) + + fig=plt.figure(figsize=(4, 4)) + plt.hist((szY-szX)/szX, bins=20, color='r', alpha=0.5) + plt.xlabel('$(R_{\\rm IDW}-R_{\\rm ORG})/R_{\\rm ORG}$') + plt.ylabel('PDF') + plt.savefig('figs/psfEllOPIDWPDF_dsz_{:}.pdf'.format(iccd)) + + + +class PSFMatsIDW_coverage(unittest.TestCase): + def test_psfIDW_(self): + #comm = MPI.COMM_WORLD + #ThisTask = comm.Get_rank() + #NTasks = comm.Get_size() + + iccd = 1 + iwave= 1 + ipsf = 400 + psfPath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/CSSOS_psf_20210326/CSST_psf_ciomp' + #test_psfIDW(iccd, iwave, psfPath, ThisTask, NTasks) + test_psfIDW(iccd, iwave, psfPath) + + ''' + for iccd in range(7, 10): + res = np.zeros(400) + for ipsf in range(1,401): + print(ipsf, end="\r") + res[ipsf-1] = test_psfResidualCalc(iccd, iwave, ipsf, psfPath) + #fig = plt.figure(figsize=(6,6)) + #plt.hist(np.abs(res), bins=50) + #plt.xlim([0,1]) + #plt.savefig('figs/psfResidualREE80PDF.pdf') + print("{:}:".format(iccd), res[res<=0.01].size/400*100) + ''' + + test_psfResidualPlot(iccd, iwave, ipsf, psfPath) + + #test_psfEll(iccd, iwave, ThisTask, NTasks) + test_psfEll(iccd, iwave) + + test_psfEllPlot(OVERPLOT=True) + #test_psfdEllPlot() + + test_psfdEllabsPlot(iccd) + + + +############################## +############################## +############################## +if __name__=='__main__': + unittest.main() diff --git a/tests/TestSpecDisperse.py b/tests/TestSpecDisperse.py new file mode 100644 index 0000000000000000000000000000000000000000..8e0f9cbe39168d3fae09cdd1d7bba446bb38a463 --- /dev/null +++ b/tests/TestSpecDisperse.py @@ -0,0 +1,435 @@ +import unittest +from ObservationSim.MockObject.SpecDisperser import rotate90, SpecDisperser + +from ObservationSim.Config import ConfigDir, ReadConfig, ChipOutput +from ObservationSim.Instrument import Telescope, Chip, FilterParam, Filter, FocalPlane +from ObservationSim.MockObject import MockObject, Star +from ObservationSim.PSF import PSFGauss + +import numpy as np +import galsim +from astropy.table import Table +from scipy import interpolate + +import matplotlib.pyplot as plt + +from lmfit.models import LinearModel, GaussianModel + + +def fit_SingleGauss(xX, yX, contmX, iHa0): + background = LinearModel(prefix='line_') + pars = background.make_params(intercept=yX.max(), slope=0) + pars = background.guess(yX, x=xX) + + gauss = GaussianModel(prefix='g_') + pars.update(gauss.make_params()) + pars['g_center'].set(iHa0, min=iHa0 - 3, max=iHa0 + 3) + pars['g_amplitude'].set(50, min=0) + pars['g_sigma'].set(12, min=0.01) + + mod = gauss + background + init = mod.eval(pars, x=xX) + outX = mod.fit(yX, pars, x=xX) + compsX = outX.eval_components(x=xX) + # print(outX.fit_report(min_correl=0.25)) + # print outX.params['g_center'] + outX.fit_report(min_correl=0.25) + # print(outX.fit_report(min_correl=0.25)) + line_slopeX = float(outX.fit_report(min_correl=0.25).split('line_slope:')[1].split('+/-')[0]) * contmX + err_line_slopeX = float( + outX.fit_report(min_correl=0.25).split('line_slope:')[1].split('+/-')[1].split('(')[0]) * contmX + + line_interceptX = float(outX.fit_report(min_correl=0.25).split('line_intercept:')[1].split('+/-')[0]) * contmX + err_line_interceptX = float( + outX.fit_report(min_correl=0.25).split('line_intercept:')[1].split('+/-')[1].split('(')[0]) * contmX + + sigmaX = float(outX.fit_report(min_correl=0.25).split('g_sigma:')[1].split('+/-')[0]) + err_sigmaX = float(outX.fit_report(min_correl=0.25).split('g_sigma:')[1].split('+/-')[1].split('(')[0]) + + fwhmX = float(outX.fit_report(min_correl=0.25).split('g_fwhm:')[1].split('+/-')[0]) + err_fwhmX = float(outX.fit_report(min_correl=0.25).split('g_fwhm:')[1].split('+/-')[1].split('(')[0]) + + centerX = float(outX.fit_report(min_correl=0.25).split('g_center:')[1].split('+/-')[0]) + err_centerX = float(outX.fit_report(min_correl=0.25).split('g_center:')[1].split('+/-')[1].split('(')[0]) + + return sigmaX, err_sigmaX, fwhmX, err_fwhmX, centerX, err_centerX + +def produceObj(x,y,chip): + pos_img = galsim.PositionD(chip.bound.xmin + x, chip.bound.ymin + y) + + param = {} + param["star"] = 1 + param["id"] = 1 + param["ra"] = chip.img.wcs.posToWorld(pos_img).ra.deg + param["dec"] = chip.img.wcs.posToWorld(pos_img).dec.deg + param["z"] = 0 + param["sed_type"] = 1 + param["model_tag"] = 1 + param["mag_use_normal"] = 10 + + obj = Star(param) + + pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, chip=chip) + wave = np.arange(2500, 11000.5, 0.5) + # sedLen = wave.shape[0] + flux = pow(wave, -2) * 1e8 + flux[200] = flux[200] * 10 + flux[800] = flux[800] * 30 + flux[2000] = flux[2000] * 5 + + obj.sed = Table(np.array([wave, flux]).T, + names=('WAVELENGTH', 'FLUX')) + return obj, pos_img + + +class TestSpecDisperse(unittest.TestCase): + + def __init__(self, methodName='runTest',conff = '', throughputf = ''): + super(TestSpecDisperse,self).__init__(methodName) + self.conff = conff + self.throughputf = throughputf + + def test_rotate901(self): + m = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]]) + m1 = np.array([[21,16,11,6,1],[22,17,12,7,2],[23,18,13,8,3],[24,19,14,9,4],[25,20,15,10,5]]) + m2 = np.array([[5,10,15,20,25],[4,9,14,19,24],[3,8,13,18,23],[2,7,12,17,22],[1,6,11,16,21]]) + xc = 2 + yc = 2 + isClockwise = 0 + m1, xc1, yc1 = rotate90(array_orig=m, xc=xc, yc=yc, isClockwise=isClockwise) + self.assertTrue(xc1-xc == 0) + self.assertTrue(yc1-yc == 0) + self.assertTrue(np.sum(m-m1) == 0) + + def test_rotate902(self): + m = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]]) + m1 = np.array([[21,16,11,6,1],[22,17,12,7,2],[23,18,13,8,3],[24,19,14,9,4],[25,20,15,10,5]]) + m2 = np.array([[5,10,15,20,25],[4,9,14,19,24],[3,8,13,18,23],[2,7,12,17,22],[1,6,11,16,21]]) + xc = 2 + yc = 2 + isClockwise =1 + m1, xc1, yc1 = rotate90(array_orig=m, xc=xc, yc=yc, isClockwise=isClockwise) + self.assertTrue(xc1-xc == 0) + self.assertTrue(yc1-yc == 0) + self.assertTrue(np.sum(m-m2) == 0) + + + def test_Specdistperse1(self): + star = galsim.Gaussian(fwhm=0.39) + g_img = galsim.Image(100, 100, scale=0.074) + starImg = star.drawImage(image=g_img) + + wave = np.arange(6200, 10000.5, 0.5) + # sedLen = wave.shape[0] + flux = pow(wave, -2) * 1e8 + # flux[200] = flux[200] * 10 + # flux[800] = flux[800] * 30 + # flux[2000] = flux[2000] * 5 + + sed = Table(np.array([wave, flux]).T, + names=('WAVELENGTH', 'FLUX')) + conff = self.conff + throughput_f = self.throughputf + thp = Table.read(throughput_f) + thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY']) + sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200, + band_end=10000, isAlongY=0, conf=conff, gid=0) + spec = sdp.compute_spec_orders() + Aimg = spec['A'][0] + wave_pix = spec['A'][5] + wave_pos = spec['A'][3] + sens = spec['A'][6] + sh = Aimg.shape + spec_pix = np.zeros(sh[1]) + for i in range(sh[1]): + spec_pix[i] = sum(Aimg[:, i]) + # figure() + # imshow(Aimg) + + wave_flux = np.zeros(wave_pix.shape[0]) + for i in np.arange(1, wave_pix.shape[0] - 1): + w = wave_pix[i] + thp_w = thp_i(w) + deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2 + f = spec_pix[wave_pos[0] - 1 + i] + + if 6200 <= w <= 10000: + f = f / thp_w + else: + f = 0 + wave_flux[i] = f / deltW + sed_i = interpolate.interp1d(sed['WAVELENGTH'], sed['FLUX']) + ids = wave_pix < 9700 + ids1 = wave_pix[ids] > 6500 + print('Spec disperse flux test') + self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004) + plt.figure() + plt.plot(wave_pix, wave_flux) + plt.plot(sed['WAVELENGTH'], sed['FLUX']) + plt.xlim(6200, 10000) + plt.ylim(1, 3) + plt.yscale('log') + plt.xlabel('$\lambda$') + plt.ylabel('$F\lambda$') + plt.legend(['extracted', 'input']) + plt.show() + + def test_Specdistperse2(self): + + psf_fwhm = 0.39 + pix_scale = 0.074 + star = galsim.Gaussian(fwhm=psf_fwhm) + g_img = galsim.Image(100, 100, scale=pix_scale) + starImg = star.drawImage(image=g_img) + + wave = np.arange(6200, 10000.5, 0.5) + # sedLen = wave.shape[0] + flux = pow(wave, -2) * 1e8 + flux[200] = flux[200] * 10 + flux[800] = flux[800] * 30 + flux[2000] = flux[2000] * 5 + + sed = Table(np.array([wave, flux]).T, + names=('WAVELENGTH', 'FLUX')) + conff = self.conff + throughput_f = self.throughputf + thp = Table.read(throughput_f) + thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY']) + sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200, + band_end=10000, isAlongY=0, conf=conff, gid=0) + spec = sdp.compute_spec_orders() + Aimg = spec['A'][0] + wave_pix = spec['A'][5] + wave_pos = spec['A'][3] + sens = spec['A'][6] + sh = Aimg.shape + spec_pix = np.zeros(sh[1]) + for i in range(sh[1]): + spec_pix[i] = sum(Aimg[:, i]) + # figure() + # imshow(Aimg) + + wave_flux = np.zeros(wave_pix.shape[0]) + for i in np.arange(1, wave_pix.shape[0] - 1): + w = wave_pix[i] + thp_w = thp_i(w) + deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2 + f = spec_pix[wave_pos[0] - 1 + i] + + if 6200 <= w <= 10000: + f = f / thp_w + else: + f = 0 + wave_flux[i] = f / deltW + sed_i = interpolate.interp1d(sed['WAVELENGTH'], sed['FLUX']) + input_em_lam = 6600 + ids = wave_pix < input_em_lam+200 + ids1 = wave_pix[ids] > input_em_lam-200 + deltLamda_pix = (max(wave_pix[ids][ids1]) - min(wave_pix[ids][ids1])) / (wave_pix[ids][ids1].shape[0] - 1) + _, _, fwhmx, fwhmx_err, center, center_err = fit_SingleGauss(wave_pix[ids][ids1], wave_flux[ids][ids1], 1.0, 6600) + + print('Emission line position and shape test') + + self.assertTrue(input_em_lam-center < deltLamda_pix) + + self.assertTrue(fwhmx/deltLamda_pix*pix_scale - psf_fwhm < np.abs(0.01)) + # print('error is ',np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))) + # self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004) + plt.figure() + plt.plot(wave_pix, wave_flux) + plt.plot(sed['WAVELENGTH'], sed['FLUX']) + plt.xlim(6200, 10000) + plt.ylim(1, 75) + plt.yscale('log') + plt.xlabel('$\lambda$') + plt.ylabel('$F\lambda$') + plt.legend(['extracted', 'input']) + plt.show() + + def test_Specdistperse3(self): + + psf_fwhm = 0.39 + pix_scale = 0.074 + star = galsim.Gaussian(fwhm=psf_fwhm) + g_img = galsim.Image(100, 100, scale=pix_scale) + starImg = star.drawImage(image=g_img) + + wave = np.arange(6200, 10000.5, 0.5) + # sedLen = wave.shape[0] + flux = pow(wave, -2) * 1e8 + flux[200] = flux[200] * 10 + flux[800] = flux[800] * 30 + flux[2000] = flux[2000] * 5 + + sed = Table(np.array([wave, flux]).T, + names=('WAVELENGTH', 'FLUX')) + conff = self.conff + throughput_f = self.throughputf + thp = Table.read(throughput_f) + thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY']) + sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200, + band_end=8000, isAlongY=0, conf=conff, gid=0) + sdp1 = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=8000, + band_end=10000, isAlongY=0, conf=conff, gid=0) + spec = sdp.compute_spec_orders() + spec1 = sdp1.compute_spec_orders() + Aimg = spec['A'][0] + spec1['A'][0] + wave_pix = spec['A'][5] + wave_pos = spec['A'][3] + sens = spec['A'][6] + sh = Aimg.shape + spec_pix = np.zeros(sh[1]) + for i in range(sh[1]): + spec_pix[i] = sum(Aimg[:, i]) + + + wave_flux = np.zeros(wave_pix.shape[0]) + for i in np.arange(1, wave_pix.shape[0] - 1): + w = wave_pix[i] + thp_w = thp_i(w) + deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2 + f = spec_pix[wave_pos[0] - 1 + i] + + if 6200 <= w <= 10000: + f = f / thp_w + else: + f = 0 + wave_flux[i] = f / deltW + + sdp2 = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200, + band_end=10000, isAlongY=0, conf=conff, gid=0) + + spec2 = sdp2.compute_spec_orders() + Aimg2 = spec2['A'][0] + + spec_pix2 = np.zeros(sh[1]) + for i in range(sh[1]): + spec_pix2[i] = sum(Aimg2[:, i]) + + wave_flux2 = np.zeros(wave_pix.shape[0]) + for i in np.arange(1, wave_pix.shape[0] - 1): + w = wave_pix[i] + thp_w = thp_i(w) + deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2 + f = spec_pix2[wave_pos[0] - 1 + i] + + if 6200 <= w <= 10000: + f = f / thp_w + else: + f = 0 + wave_flux2[i] = f / deltW + + r1_i = interpolate.interp1d(wave_pix, wave_flux2) + r2_i = interpolate.interp1d(wave_pix, wave_flux) + + print('Spec Splicing test') + self.assertTrue(r1_i(8000)-r2_i(8000) < np.abs(0.0001)) + + # self.assertTrue(fwhmx/deltLamda_pix*pix_scale - psf_fwhm < np.abs(0.01)) + # print('error is ',np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))) + # self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004) + plt.figure() + plt.plot(wave_pix, wave_flux2) + plt.plot(wave_pix, wave_flux) + # plt.plot(sed['WAVELENGTH'], sed['FLUX']) + plt.xlim(6200, 10000) + plt.ylim(1, 4) + plt.yscale('log') + plt.xlabel('$\lambda$') + plt.ylabel('$F\lambda$') + plt.legend(['one spec', 'split in 8000 A']) + plt.show() + + + + def test_double_disperse(self): + work_dir = "/public/home/fangyuedong/CSST_unittest/CSST/test/" + # data_dir = "/Volumes/Extreme SSD/SimData/" + data_dir = "/data/simudata/CSSOSDataProductsSims/data/" + path_dict = ConfigDir(work_dir=work_dir, data_dir=data_dir) + config = ReadConfig(path_dict["config_file"]) + + filter_param = FilterParam(filter_dir=path_dict["filter_dir"]) + focal_plane = FocalPlane(survey_type=config["survey_type"]) + chip = Chip(1, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"], + normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict["sls_dir"], + config=config) + filter_id, filter_type = chip.getChipFilter() + filt = Filter(filter_id=filter_id, filter_type=filter_type, filter_param=filter_param, + ccd_bandpass=chip.effCurve) + tel = Telescope(optEffCurve_path=path_dict["mirror_file"]) + + psf_model = PSFGauss(chip=chip) + + + wcs_fp = focal_plane.getTanWCS(config["ra_center"], config["dec_center"], config["image_rot"], chip.pix_scale) + chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + chip.img.wcs = wcs_fp + + obj, pos_img = produceObj(2000,4500, chip) + obj.drawObj_slitless( + tel=tel, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=0, + g2=0, + exptime=150, + normFilter=chip.normF_star) + + obj, pos_img = produceObj(3685, 6500, chip) + obj.drawObj_slitless( + tel=tel, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=0, + g2=0, + exptime=150, + normFilter=chip.normF_star) + + obj, pos_img = produceObj(5000, 2500, chip) + obj.drawObj_slitless( + tel=tel, + pos_img=pos_img, + psf_model=psf_model, + bandpass_list=filt.bandpass_sub_list, + filt=filt, + chip=chip, + g1=0, + g2=0, + exptime=150, + normFilter=chip.normF_star) + + print('Spec double disperse test') + from astropy.io import fits + fits.writeto('test.fits',chip.img.array, overwrite = True) + + # plt.figure() + # plt.imshow(chip.img.array) + # plt.show() + + + +if __name__ == '__main__': + conff='/data/simudata/CSSOSDataProductsSims/data/CONF/CSST_GI2.conf' + throughputf='/data/simudata/CSSOSDataProductsSims/data/CONF/GI.Throughput.1st.fits' + + suit = unittest.TestSuite() + case1 = TestSpecDisperse('test_Specdistperse1',conff,throughputf) + suit.addTest(case1) + case2 = TestSpecDisperse('test_Specdistperse2', conff, throughputf) + suit.addTest(case2) + case3 = TestSpecDisperse('test_Specdistperse3', conff, throughputf) + suit.addTest(case3) + case4 = TestSpecDisperse('test_double_disperse', conff, throughputf) + suit.addTest(case4) + + unittest.TextTestRunner(verbosity=2).run(suit) + # runner = unittest.TextTestRunner() + # runner.run(suit) \ No newline at end of file diff --git a/tests/test.fits b/tests/test.fits new file mode 100644 index 0000000000000000000000000000000000000000..f6b8b81706219c422b18bd4237d88fe26254c726 Binary files /dev/null and b/tests/test.fits differ