diff --git a/Catalog/C3Catalog.py b/Catalog/C3Catalog.py index b36b7223596a99a85354bf05050b3cb64c55d215..08d9ecf352821c27074bdf824c755c7c0b223b49 100644 --- a/Catalog/C3Catalog.py +++ b/Catalog/C3Catalog.py @@ -84,10 +84,46 @@ class C3Catalog(CatalogBase): self.rng_sedGal = random.Random() self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed self.ud = galsim.UniformDeviate(pix_id) + + # Apply astrometric modeling + # in C3 case only aberration + ra_arr = gals['ra_true'][:] + dec_arr = gals['dec_true'][:] + if self.config["obs_setting"]["enable_astrometric_model"]: + ra_list = ra_arr.tolist() + dec_list = dec_arr.tolist() + pmra_list = np.zeros(ngals).tolist() + pmdec_list = np.zeros(ngals).tolist() + rv_list = np.zeros(ngals).tolist() + parallax_list = [1e-9] * ngals + dt = datetime.fromtimestamp(self.pointing.timestamp) + date_str = dt.date().isoformat() + time_str = dt.time().isoformat() + ra_arr, dec_arr = on_orbit_obs_position( + input_ra_list=ra_list, + input_dec_list=dec_list, + input_pmra_list=pmra_list, + input_pmdec_list=pmdec_list, + input_rv_list=rv_list, + input_parallax_list=parallax_list, + input_nstars=ngals, + input_x=self.pointing.sat_x, + input_y=self.pointing.sat_y, + input_z=self.pointing.sat_z, + input_vx=self.pointing.sat_vx, + input_vy=self.pointing.sat_vy, + input_vz=self.pointing.sat_vz, + input_epoch="J2015.5", + input_date_str=date_str, + input_time_str=time_str + ) + for igals in range(ngals): param = self.initialize_param() - param['ra'] = gals['ra_true'][igals] - param['dec'] = gals['dec_true'][igals] + param['ra'] = ra_arr[igals] + param['dec'] = dec_arr[igals] + param['ra_orig'] = gals['ra_true'][igals] + param['dec_orig'] = gals['dec_true'][igals] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] @@ -129,7 +165,8 @@ class C3Catalog(CatalogBase): param['star'] = 2 # Quasar self.ids += 1 - param['id'] = self.ids + # param['id'] = self.ids + param['id'] = gals['galaxyID'][igals] if param['star'] == 0: obj = Galaxy(param, self.rotation) @@ -141,9 +178,9 @@ class C3Catalog(CatalogBase): def _load_stars(self, stars, pix_id=None): nstars = len(stars['sourceID']) # Apply astrometric modeling + # in C3 case only aberration ra_arr = stars["RA"][:] dec_arr = stars["Dec"][:] - # if "astrometric_lib" in self.config["obs_setting"] and self.config["obs_setting"]["enable_astrometric_model"]: if self.config["obs_setting"]["enable_astrometric_model"]: ra_list = ra_arr.tolist() dec_list = dec_arr.tolist() @@ -170,20 +207,22 @@ class C3Catalog(CatalogBase): input_vz=self.pointing.sat_vz, input_epoch="J2015.5", input_date_str=date_str, - input_time_str=time_str, - # lib_path=lib_path + input_time_str=time_str ) for istars in range(nstars): param = self.initialize_param() param['ra'] = ra_arr[istars] param['dec'] = dec_arr[istars] + param['ra_orig'] = stars["RA"][istars] + param['dec_orig'] = stars["Dec"][istars] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue param['mag_use_normal'] = stars['app_sdss_g'][istars] if param['mag_use_normal'] >= 26.5: continue self.ids += 1 - param['id'] = self.ids + # param['id'] = self.ids + param['id'] = stars['sourceID'][istars] param['sed_type'] = stars['sourceID'][istars] param['model_tag'] = stars['model_tag'][istars] param['teff'] = stars['teff'][istars] diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index dc7e465cd4d6915ebedf14f78cf5dc66d35c51a4..a2285f36688da23c99062852fdb1e59cba2fd4b6 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -23,73 +23,81 @@ class ChipOutput(object): self.cat_name = 'MSC_' + config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat" self.subdir = subdir - # hdr1 = "#ID ID_chip filter xImage yImage ra dec z mag flag SNR " - hdr1 = "#ID ID_chip filter xImage yImage ra dec z mag flag " - hdr2 = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge " - hdr3 = "e1PSF e2PSF e1 e2 g1 g2 e1OBS e2OBS " - hdr4 = "sed_type av redden " - hdr5 = "star_model teff logg feh\n" - # fmt1 = "%10d %4d %5s %10.3f %10.3f %15.6f %15.6f %7.4f %8.4f %2d %9.2f " - fmt1 = "%10d %4d %5s %10.3f %10.3f %15.6f %15.6f %7.4f %8.4f %2d " - fmt2 = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f " - fmt3 = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f " - fmt4 = "%2d %8.4f %8.4f " - fmt5 = "%10s %8.4f %8.4f %8.4f\n" - self.hdr = hdr1 + hdr2 + hdr3 + hdr4 + hdr5 - self.fmt = fmt1 + fmt2 + fmt3 + fmt4 + fmt5 + hdr1 = "obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type " + hdr2 = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge g1 g2 " + hdr3 = "sed_type av redden " + hdr4 = "pm_ra pm_dec RV parallax\n" + + fmt1 = "%10d %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s " + fmt2 = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f " + fmt3 = "%2d %8.4f %8.4f " + fmt4 = "%15.8f %15.8f %15.8f %15.8f\n" + + self.hdr = hdr1 + hdr2 + hdr3 + hdr4 + self.fmt = fmt1 + fmt2 + fmt3 + fmt4 + print("pointing_type = %s\n"%(pointing_type)) if pointing_type == 'MS': self.cat = open(os.path.join(self.subdir, self.cat_name), "w") print("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name))) self.cat.write(self.hdr) - def updateHDR(self, hdr): - hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise, "comment":"read noise in e-/pixel"}, - {"name":"DARK", "value":self.chip.dark_noise, "comment":"Dark noise (e-/pixel/s)"}, - {"name":"EXPTIME", "value":self.exptime, "comment":"exposure time in second"}, - {"name":"GAIN", "value":self.chip.gain, "comment":"CCD gain in e-/ADU"}, - {"name":"SATURATE","value":65535.0, "comment":"saturation level"}, - {"name":"CCDCHIP", "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"}, - {"name":"FILTER", "value":self.filt.filter_type, "comment":"filter name"}, - {"name":"MJD-OBS", "value":self.mjdTime, "comment":"Modified Julian Date (MJD) of observation"}, - {"name":"DATE-OBS","value":self.imgKey1, "comment":"Date of observation"}, - {"name":"EQUINOX", "value":2000.0}, - {"name":"RADECSYS","value":"ICRS"}, - {"name":"RA", "value":self.ra_cen, "comment":"telescope pointing center"}, - {"name":"DEC", "value":self.dec_cen, "comment":"telescope pointing center"}, - {"name":"OBJECT", "value":"CSS-OS"}, - {"name":"WCSDIM", "value":2.0, "comment":"WCS Dimensionality"}, - {"name":"EXTNAME", "value":"IM1", "comment":"Extension name"}, - {"name":"BSCALE", "value":1.0}, - {"name":"BZERO", "value":0.0}, - {"name":"OBSID", "value":self.imgKey0, "comment":"Observation ID"}, - {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"}, - {"name":"RSPEED", "value":10.0, "comment":"Read speed"}, - {"name":"CHIPTEMP","value":-100.0, "comment":"Chip temperature"}, - {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"}, - {"name":"CCDSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, - {"name":"NSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, - {"name":"AUTHOR", "value":"CSST-Sim Group"}, - {"name":"GROUP", "value":"Weak Lensing Working Group for CSST"}] - for item in hdrNew: - hdr.add_record(item) - return hdr + # def updateHDR(self, hdr): + # hdrNew = [{"name":"RDNOISE", "value":self.chip.read_noise, "comment":"read noise in e-/pixel"}, + # {"name":"DARK", "value":self.chip.dark_noise, "comment":"Dark noise (e-/pixel/s)"}, + # {"name":"EXPTIME", "value":self.exptime, "comment":"exposure time in second"}, + # {"name":"GAIN", "value":self.chip.gain, "comment":"CCD gain in e-/ADU"}, + # {"name":"SATURATE","value":65535.0, "comment":"saturation level"}, + # {"name":"CCDCHIP", "value":int(self.chipLabel), "comment":"chip ID in the CCD mosaic"}, + # {"name":"FILTER", "value":self.filt.filter_type, "comment":"filter name"}, + # {"name":"MJD-OBS", "value":self.mjdTime, "comment":"Modified Julian Date (MJD) of observation"}, + # {"name":"DATE-OBS","value":self.imgKey1, "comment":"Date of observation"}, + # {"name":"EQUINOX", "value":2000.0}, + # {"name":"RADECSYS","value":"ICRS"}, + # {"name":"RA", "value":self.ra_cen, "comment":"telescope pointing center"}, + # {"name":"DEC", "value":self.dec_cen, "comment":"telescope pointing center"}, + # {"name":"OBJECT", "value":"CSS-OS"}, + # {"name":"WCSDIM", "value":2.0, "comment":"WCS Dimensionality"}, + # {"name":"EXTNAME", "value":"IM1", "comment":"Extension name"}, + # {"name":"BSCALE", "value":1.0}, + # {"name":"BZERO", "value":0.0}, + # {"name":"OBSID", "value":self.imgKey0, "comment":"Observation ID"}, + # {"name":"CCDNAME", "value":"ccd"+self.chipLabel,"comment":"CCD name"}, + # {"name":"RSPEED", "value":10.0, "comment":"Read speed"}, + # {"name":"CHIPTEMP","value":-100.0, "comment":"Chip temperature"}, + # {"name":"DATASEC", "value":"1:%d,1:%d"%(self.chip.npix_x,self.chip.npix_y), "comment":"Data section"}, + # {"name":"CCDSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, + # {"name":"NSUM", "value":self.chip.npix_x*self.chip.npix_y, "comment":"CCD pixel summing"}, + # {"name":"AUTHOR", "value":"CSST-Sim Group"}, + # {"name":"GROUP", "value":"Weak Lensing Working Group for CSST"}] + # for item in hdrNew: + # hdr.add_record(item) + # return hdr # def cat_add_obj(self, obj, pos_img, snr, pos_shear, g1, g2): - def cat_add_obj(self, obj, pos_img, pos_shear, g1, g2): + def cat_add_obj(self, obj, pos_img, pos_shear): ximg = pos_img.x - self.chip.bound.xmin + 1.0 yimg = pos_img.y - self.chip.bound.ymin + 1.0 - e1, e2, g1, g2, e1OBS, e2OBS = obj.getObservedEll(g1, g2) - if obj.type == 'galaxy': - line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, - obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, - pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0) - elif obj.type == "quasar": - line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, - obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0) - else: - line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, - pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh']) - # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS) + # if obj.type == 'galaxy': + # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, + # obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, + # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0, 0, 0) + # elif obj.type == "quasar": + # line = self.fmt % (obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, + # obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, obj.sed_type, obj.param['av'], obj.param['redden'], 'n', 0.0, 0.0, 0.0) + # else: + # line = self.fmt%(obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.z, obj.getMagFilter(self.filt), obj.param["star"], 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + # pos_shear.g1, pos_shear.g2, e1, e2, g1, g2, e1OBS, e2OBS, 0, 0.0, 0.0, obj.param['model_tag'], obj.param['teff'], obj.param['logg'],obj.param['feh']) + # print( + # obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, + # obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2, + # obj.sed_type, obj.av, obj.redden, + # obj.pmra, obj.pmdec, obj.rv, obj.parallax) + + line = self.fmt%( + obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type, + obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2, + obj.sed_type, obj.av, obj.redden, + obj.pmra, obj.pmdec, obj.rv, obj.parallax) self.cat.write(line) diff --git a/ObservationSim/MockObject/Catalog.py b/ObservationSim/MockObject/Catalog.py deleted file mode 100644 index a86d36fe2f0c115b1e38e2b37e03b07479772b7b..0000000000000000000000000000000000000000 --- a/ObservationSim/MockObject/Catalog.py +++ /dev/null @@ -1,179 +0,0 @@ -import os -import numpy as np -import random -import galsim -import h5py as h5 -import healpy as hp -from astropy.table import Table -from astropy.coordinates import spherical_to_cartesian - -from ObservationSim.MockObject.Star import Star -from ObservationSim.MockObject.Galaxy import Galaxy -from ObservationSim.MockObject.Quasar import Quasar -from ObservationSim.MockObject._util import seds, sed_assign, extAv - -NSIDE = 128 - -class Catalog(object): - 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.seed_Av = config["random_seeds"]["seed_Av"] - - 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' - - if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]: - 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"] and config["input_path"]["galaxy_cat"]: - 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() - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - # self.sky_coverage_enlarged = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.5) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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)) - # print(vertices_pix) - vertices = spherical_to_cartesian(1., dec, ra) - 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_star = os.path.join(self.template_dir, "PicklesStars/") - # 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") - sed_info = Table.read(sed_info_file, format="ascii") - self.cosids = sed_info["IDCosmoDC2"] - self.objtypes = sed_info["objType"] - - def _load_gals(self, gals, pix_id=None): - ngals = len(gals['galaxyID']) - self.rng_sedGal = random.Random() - self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed - self.ud = galsim.UniformDeviate(pix_id) - for igals in range(ngals): - param = {} - param['ra'] = gals['ra_true'][igals] - param['dec'] = gals['dec_true'][igals] - param['z'] = gals['redshift_true'][igals] - param['model_tag'] = 'None' - param['gamma1'] = 0 - param['gamma2'] = 0 - param['kappa'] = 0 - param['delta_ra'] = 0 - param['delta_dec'] = 0 - sersicB = gals['sersic_bulge'][igals] - hlrMajB = gals['size_bulge_true'][igals] - hlrMinB = gals['size_minor_bulge_true'][igals] - sersicD = gals['sersic_disk'][igals] - hlrMajD = gals['size_disk_true'][igals] - hlrMinD = gals['size_minor_disk_true'][igals] - aGal = gals['size_true'][igals] - bGal = gals['size_minor_true'][igals] - param['bfrac'] = gals['bulge_to_total_ratio_i'][igals] - param['theta'] = gals['position_angle_true'][igals] - param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB) - param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD) - param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB) - param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD) - param['ell_tot'] = (aGal - bGal) / (aGal + bGal) - - # Assign each galaxy a template SED - param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal) - param['redden'] = self.tempRed_gal[param['sed_type']] - param['av'] = self.avGal[int(self.ud()*self.nav)] - if param['sed_type'] <= 5: - param['av'] = 0.0 - param['redden'] = 0 - param['star'] = 0 # Galaxy - if param['sed_type'] >= 29: - param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im - param['star'] = 2 # Quasar - - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] - if param['mag_use_normal'] >= 26.5: - continue - self.ids += 1 - param['id'] = self.ids - - if param['star'] == 0: - obj = Galaxy(param, self.rotation) - self.objs.append(obj) - if param['star'] == 2: - obj = Quasar(param) - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - for istars in range(nstars): - param = {} - param['ra'] = stars['RA'][istars] - param['dec'] = stars['Dec'][istars] - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - param['mag_use_normal'] = stars['app_sdss_g'][istars] - if param['mag_use_normal'] >= 26.5: - continue - self.ids += 1 - param['id'] = self.ids - param['sed_type'] = stars['sourceID'][istars] - param['model_tag'] = stars['model_tag'][istars] - param['teff'] = stars['teff'][istars] - param['logg'] = stars['grav'][istars] - param['feh'] = stars['feh'][istars] - param['z'] = 0.0 - param['star'] = 1 # Star - obj = Star(param) - self.objs.append(obj) - - def _load(self): - self.nav = 15005 - self.avGal = extAv(self.nav, seed=self.seed_Av) - gals_cat = h5.File(self.galaxy_path, 'r')['galaxies'] - star_cat = h5.File(self.star_path, 'r')['catalog'] - self.objs = [] - self.ids = 0 - for pix in self.pix_list: - gals = gals_cat[str(pix)] - stars = star_cat[str(pix)] - self._load_gals(gals, pix_id=pix) - self._load_stars(stars, pix_id=pix) - print("number of objects in catalog: ", len(self.objs)) - del self.avGal \ No newline at end of file diff --git a/ObservationSim/MockObject/CatalogBase.py b/ObservationSim/MockObject/CatalogBase.py index 08e82ca110ad837bb105849cced1dbc73f918b96..442c325dc237fde7536fa4d0b28dd7694bf1ede1 100644 --- a/ObservationSim/MockObject/CatalogBase.py +++ b/ObservationSim/MockObject/CatalogBase.py @@ -29,29 +29,33 @@ class CatalogBase(metaclass=ABCMeta): "star":-1, "id":-1, "ra":0, - "dec":0, - "z":0, + "dec":0., + "ra_orig":0., + "dec_orig":0., + "z":0., "sed_type":-1, "model_tag":"unknown", - "mag_use_normal":100, - "theta":0, - "kappa":0, - "gamma1":0, - "gamma2":0, - "bfrac":0, - "hlr_bulge":0, - "hlr_disk":0, - "ell_bulge":0, - "ell_disk":0, - "ell_tot":0, - "teff":0, - "logg":0, - "feh":0, - "g1":0, - "g2":0, - "pmra":0, - "pmdec":0, - "rv":0, + "mag_use_normal":100., + "theta":0., + "kappa":0., + "gamma1":0., + "gamma2":0., + "bfrac":0., + "av":0., + "redden":0., + "hlr_bulge":0., + "hlr_disk":0., + "ell_bulge":0., + "ell_disk":0., + "ell_tot":0., + "teff":0., + "logg":0., + "feh":0., + "g1":0., + "g2":0., + "pmra":0., + "pmdec":0., + "rv":0., "parallax":1e-9 } return param diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 5386f490bf8ce5a0d49c5dc1a1899de7644ab09e..d0cdcef2eb11aa4bcf1e8a1701a1ba2782064c98 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -131,12 +131,12 @@ class Galaxy(MockObject): bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) bulge = bulge.shear(bulge_shape) - # gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk + gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk - # (TEST) Random knots - knots = galsim.RandomKnots(npoints=100, profile=disk) - kfrac = np.random.random()*(1.0 - self.bfrac) - gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots + # # (TEST) Random knots + # knots = galsim.RandomKnots(npoints=100, profile=disk) + # kfrac = np.random.random()*(1.0 - self.bfrac) + # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots gal = gal.withFlux(nphotons) gal_shear = galsim.Shear(g1=g1, g2=g2) @@ -220,12 +220,12 @@ class Galaxy(MockObject): bulge_shape = galsim.Shear(g1=self.e1_bulge, g2=self.e2_bulge) bulge = bulge.shear(bulge_shape) - # gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk + gal = self.bfrac * bulge + (1.0 - self.bfrac) * disk # (TEST) Random knots - knots = galsim.RandomKnots(npoints=100, profile=disk) - kfrac = np.random.random()*(1.0 - self.bfrac) - gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots + # knots = galsim.RandomKnots(npoints=100, profile=disk) + # kfrac = np.random.random()*(1.0 - self.bfrac) + # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots gal = gal.withFlux(tel.pupil_area * exptime) gal_shear = galsim.Shear(g1=g1, g2=g2) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index a0ebdc4b58d973bef8df8285971da88d55dd16fb..055bab0aecbe5bf9d9f7abb6ee7b6f6d3922e9dc 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -17,15 +17,34 @@ class MockObject(object): self.type = "star" elif self.param["star"] == 2: self.type = "quasar" + self.id = self.param["id"] self.ra = self.param["ra"] self.dec = self.param["dec"] + self.ra_orig = self.param["ra_orig"] + self.dec_orig = self.param["dec_orig"] self.z = self.param["z"] self.sed_type = self.param["sed_type"] self.model_tag = self.param["model_tag"] self.mag_use_normal = self.param["mag_use_normal"] self.sed = None + # Place holder for outputs + self.av = self.param["av"] + self.redden = self.param["redden"] + self.pmra = self.param["pmra"] + self.pmdec = self.param["pmdec"] + self.rv = self.param["rv"] + self.parallax = self.param["parallax"] + self.g1 = self.param["g1"] + self.g2 = self.param["g2"] + self.thetaR = self.param["theta"] + self.bfrac = self.param["bfrac"] + self.hlr_disk = self.param["hlr_disk"] + self.hlr_bulge = self.param["hlr_bulge"] + self.e1_disk, self.e2_disk = 0., 0. + self.e1_bulge, self.e2_bulge = 0., 0. + def getMagFilter(self, filt): if filt.filter_type in ["GI", "GV", "GU"]: return self.param["mag_use_normal"] @@ -304,5 +323,5 @@ class MockObject(object): snr_obj = img_flux / sig_obj return snr_obj - def getObservedEll(self, g1=0, g2=0): - return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + # def getObservedEll(self, g1=0, g2=0): + # return 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 diff --git a/ObservationSim/MockObject/__init__.py b/ObservationSim/MockObject/__init__.py index 4868239200c55c30ee405cad6d0b341a647a1859..63734c6e4faef8672477eeca8f3d64b03c9aa207 100755 --- a/ObservationSim/MockObject/__init__.py +++ b/ObservationSim/MockObject/__init__.py @@ -3,6 +3,5 @@ from .Galaxy import Galaxy from .CatalogBase import CatalogBase from .Quasar import Quasar from .Star import Star -from .Catalog import Catalog from .SkybackgroundMap import * # from .CosmicRay import CosmicRay \ No newline at end of file diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 137107ca1b38663c31a6ae3beceecb7e547413f6..5485463e075722d8aa39729350d1e2d252ad62c4 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -19,10 +19,8 @@ class Observation(object): def __init__(self, config, Catalog, work_dir=None, data_dir=None): self.path_dict = config_dir(config=config, work_dir=work_dir, data_dir=data_dir) self.config = config - # self.tel = Telescope(optEffCurve_path=self.path_dict["mirror_file"]) self.tel = Telescope() self.focal_plane = FocalPlane(survey_type=self.config["obs_setting"]["survey_type"]) - # self.filter_param = FilterParam(filter_dir=self.path_dict["filter_dir"]) self.filter_param = FilterParam() self.chip_list = [] self.filter_list = [] @@ -44,10 +42,6 @@ class Observation(object): # Make Chip & Filter lists chip = Chip( chipID=chipID, - # ccdEffCurve_dir=self.path_dict["ccd_dir"], - # CRdata_dir=self.path_dict["CRdata_dir"], - # normalize_dir=self.path_dict["normalize_dir"], - # sls_dir=self.path_dict["sls_dir"], config=self.config) filter_id, filter_type = chip.getChipFilter() filt = Filter(filter_id=filter_id, @@ -114,7 +108,6 @@ class Observation(object): skyMap=flat_normal, blueLimit=filt.blue_limit, redLimit=filt.red_limit, - # skyfn=self.path_dict["sky_file"], conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0) @@ -130,7 +123,7 @@ class Observation(object): bright_obj = 0 dim_obj = 0 for j in range(self.nobj): - # if j >= 100: + # if j >= 10: # break obj = self.cat.objs[j] if obj.type == 'star' and self.config["run_option"]["galaxy_only"]: @@ -170,13 +163,13 @@ class Observation(object): if self.config["shear_setting"]["shear_type"] == "constant": if obj.type == 'star': - obj.param["g1"], obj.param["g2"] = 0, 0 + obj.g1, obj.g2 = 0., 0. else: - obj.param["g1"], obj.param["g2"] = self.g1_field, self.g2_field + obj.g1, obj.g2 = self.g1_field, self.g2_field elif self.config["shear_setting"]["shear_type"] == "extra": try: # TODO: every object with individual shear from input catalog(s) - obj.param["g1"], obj.param["g2"] = self.g1_field[j], self.g2_field[j] + obj.g1, obj.g2 = self.g1_field[j], self.g2_field[j] except: print("failed to load external shear.") pass @@ -205,8 +198,8 @@ class Observation(object): bandpass_list=filt.bandpass_sub_list, filt=filt, chip=chip, - g1=obj.param["g1"], - g2=obj.param["g2"], + g1=obj.g1, + g2=obj.g2, exptime=pointing.exp_time ) elif chip.survey_type == "spectroscopic" and not self.config["run_option"]["out_cat_only"]: @@ -217,14 +210,15 @@ class Observation(object): bandpass_list=filt.bandpass_sub_list, filt=filt, chip=chip, - g1=obj.param["g1"], - g2=obj.param["g2"], + g1=obj.g1, + g2=obj.g2, exptime=pointing.exp_time, normFilter=norm_filt, ) if isUpdated: # TODO: add up stats - chip_output.cat_add_obj(obj, pos_img, pos_shear, obj.param["g1"], obj.param["g2"]) + # print("updating output catalog...") + chip_output.cat_add_obj(obj, pos_img, pos_shear) pass else: # print("object omitted", flush=True) diff --git a/config/config_C3.yaml b/config/config_C3.yaml index 1e084a698c19f82de2788eddf830940869d0b65a..4726ae30df8b45dbc72435788c4f3dd369578bc9 100644 --- a/config/config_C3.yaml +++ b/config/config_C3.yaml @@ -67,7 +67,6 @@ obs_setting: image_rot: -113.4333 # Number of calibration pointings - # Note: only valid when a pointing list is specified np_cal: 0 # Run specific pointing(s): @@ -100,7 +99,6 @@ SED_templates_path: star_SED: "Catalog_20210126/SpecLib.hdf5" galaxy_SED: "Templates/Galaxy/" - ############################################### # PSF setting ############################################### diff --git a/config/config_NGP.yaml b/config/config_NGP.yaml new file mode 100644 index 0000000000000000000000000000000000000000..e0e78d5b68459951ded1064ea2c7f7e8bf195fc4 --- /dev/null +++ b/config/config_NGP.yaml @@ -0,0 +1,213 @@ +--- +############################################### +# +# 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/temp/CSST/workplace/" +data_dir: "/data/simudata/CSSOSDataProductsSims/data/" +run_name: "NGP_20220327" + +# (Optional) a file of point list +# if you just want to run default pointing: +# - pointing_dir: null +# - pointing_file: null +pointing_dir: null +pointing_file: null + +# Whether to use MPI +run_option: + use_mpi: YES + # NOTE: "n_threads" paramters is currently not used in the backend + # simulation codes. It should be implemented later in the web frontend + # in order to config the number of threads to request from NAOC cluster + n_threads: 80 + + # 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 + # "Spectroscopic": simulate slitless spectroscopic chips only + # "All": simulate full focal plane + survey_type: "All" + + # 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: 192.8595 + dec_center: 27.1283 + + # Image rotation [degree] + image_rot: -113.4333 + + # Number of calibration pointings + np_cal: 0 + + # Run specific pointing(s): + # - give a list of indexes of pointings: [ip_1, ip_2...] + # - run all pointings: null + # Note: only valid when a pointing list is specified + run_pointings: null + + # Run specific chip(s): + # - give a list of indexes of chips: [ip_1, ip_2...] + # - run all chips: null + # Note: for all pointings + run_chips: null + + # Whether to enable astrometric modeling + # astrometric_lib: "libshao.so" + enable_astrometric_model: True + +############################################### +# Input path setting +############################################### + +# Default path settings for WIDE survey simulation +input_path: + cat_dir: "OnOrbitCalibration/CTargets20211231" + star_cat: "CT-NGP_V2.2.hdf5" + galaxy_cat: "galaxyCats_r_10.0_healpix_shift_192.859500_27.128300.hdf5" + +SED_templates_path: + star_SED: "Catalog_20210126/SpecLib.hdf5" + galaxy_SED: "Templates/Galaxy/" + +############################################### +# 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: "/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube" + + # path to field-distortion model + # Note: only valid when ins_effects: field_dist is "ON" + fd_path: "FieldDistModelGlobal_v1.0.pickle" + + # 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: ON # Whether to add non-linearity + cosmic_ray: ON # Whether to add cosmic-ray + cray_differ: ON # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: OFF # Whether to simulate CTE trails + saturbloom: ON # Whether to simulate Saturation & Blooming + add_badcolumns: ON # Whether to add bad columns + add_hotpixels: ON # Whether to add hot pixels + add_deadpixels: ON # Whether to add dead(dark) pixels + bright_fatter: ON # Whether to simulate Brighter-Fatter (also diffusion) effect + + # 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/run_NGP.pbs b/run_NGP.pbs new file mode 100755 index 0000000000000000000000000000000000000000..8d7ed893611944d1445c9c581b063555619015a6 --- /dev/null +++ b/run_NGP.pbs @@ -0,0 +1,19 @@ +#!/bin/bash + +#PBS -N SIMS +##PBS -l walltime=70:00:00 + +##mpdallexit +##mpdboot -n 10 -f ./mpd.hosts + +##PBS -j oe +#PBS -l nodes=comput110:ppn=80 +#####PBS -q longq +#PBS -q batch +#PBS -u fangyuedong + +NP=30 +date +echo $NP + +mpirun -np $NP python /public/home/fangyuedong/temp/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/temp/CSST/config diff --git a/tests/PSFInterpTest/PSFInterpModule_coverage.py b/tests/PSFInterpTest/PSFInterpModule_coverage.py new file mode 100644 index 0000000000000000000000000000000000000000..079623803dbb62efd0870aaa3dd78da00f554881 --- /dev/null +++ b/tests/PSFInterpTest/PSFInterpModule_coverage.py @@ -0,0 +1,402 @@ +import unittest + +import sys,os,math +from itertools import islice + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib as mpl +mpl.use('Agg') + +import yaml +from ObservationSim.Config import Config +from ObservationSim.Config.Config import config_dir +from ObservationSim.Instrument import Chip + +from ObservationSim.PSF.PSFInterp import PSFInterp + + +def defineCCD(iccd, config_file): + with open(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) + path_dict = config_dir(config=config, work_dir=config['work_dir'], data_dir=config['data_dir']) + chip = Chip(chipID=iccd, config=config) + #chip = Chip(chipID=iccd, 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) + return chip + +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 test_psfEll(iccd, iwave, psfMat): + psfMat_iwave = psfMat.psfMat[iwave-1, :,:,:] + npsf = np.shape(psfMat_iwave)[0] + + 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(1, npsf+1): + print('ipsf-{:}'.format(ipsf), end='\r') + imx[ipsf-1] = psfMat.cen_col[iwave-1, ipsf-1] + imy[ipsf-1] = psfMat.cen_row[iwave-1, ipsf-1] + psfMat_iwave_ipsf = psfMat_iwave[ipsf-1, :, :] + cenX = 256 + cenY = 256 + sz, e1, e2 = psfSecondMoments(psfMat_iwave_ipsf, cenX, cenY, pixSize=1) + psf_e1[ipsf-1] = e1 + psf_e2[ipsf-1] = e2 + psf_sz[ipsf-1] = sz + #print('ell======', ipsf, np.sqrt(e1**2 + e2**2)) + ####### + arr = [imx, imy, psf_e1, psf_e2, psf_sz] + np.save('data/psfEll{:}_{:}_{:}'.format(int(np.sqrt(npsf)),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.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)) + + +def test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB): + bandpass = iwave-1 + class pos_img(): + def __init__(self,x, y): + self.x = x*1e3/10. #in unit of pixels + self.y = y*1e3/10. + + psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:] + npsf = np.shape(psfMat_iwave)[0] + psf_e1 = np.zeros(npsf) + psf_e2 = np.zeros(npsf) + psf_sz = np.zeros(npsf) + + for ipsf in range(1, npsf+1): + print('ipsf:', ipsf, end='\r', flush=True) + tpos_img = pos_img(psfMatA.cen_col[iwave-1, ipsf-1], psfMatA.cen_row[iwave-1, ipsf-1]) + psfIDW = psfMatB.get_PSF(chip, tpos_img, bandpass, galsimGSObject=False, findNeighMode='treeFind') + np.save('figs/psfIDW/psfIDW_{:}_{:}_{:}'.format(iccd, iwave, ipsf), psfIDW) + + cenX = 256 + cenY = 256 + sz, e1, e2 = psfSecondMoments(psfIDW, cenX, cenY, pixSize=1) + psf_e1[ipsf-1] = e1 + psf_e2[ipsf-1] = e2 + psf_sz[ipsf-1] = sz + arr = [psf_e1, psf_e2, psf_sz] + np.save('data/psfEll20IDW_{:}_{:}'.format(iccd, iwave), arr) + + +def test_psfResidualPlot(iccd, iwave, ipsf, psfMatA): + psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:] + + psfMatORG = psfMat_iwave[ipsf-1, :, :] + 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_psfEllIDWPlot(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.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_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=(6, 6)) + 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 PSFInterpModule_coverage(unittest.TestCase): + def test_psfEll_(self): + iccd = 1 + iwave= 1 + + config_file = "/public/home/weichengliang/CSST_git/newVersion/CSST/config/config_C3.yaml" + chip = defineCCD(iccd, config_file) + print(chip.chipID) + print(chip.cen_pix_x, chip.cen_pix_y) + + ipath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCubeTest' + psfMatA = PSFInterp(chip, npsf=400, PSF_data_file=ipath, PSF_data_prefix="S20x20_") + psfMatB = PSFInterp(chip, npsf=900, PSF_data_file=ipath, PSF_data_prefix="S30x30_") + + test_psfEll(iccd, iwave, psfMatA) + test_psfEll(iccd, iwave, psfMatB) + test_psfEllPlot(OVERPLOT=True) + + test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB) + ipsf = 1 + test_psfResidualPlot(iccd, iwave, ipsf, psfMatA) + test_psfEllIDWPlot(OVERPLOT=True) + test_psfdEllabsPlot(iccd) + + +if __name__ == '__main__': + unittest.main() + print('#####haha#####') + + + diff --git a/tests/PSFInterpTest/loadPSFSet.py b/tests/PSFInterpTest/loadPSFSet.py new file mode 100644 index 0000000000000000000000000000000000000000..5ebef6f265d9d30707b263da3ec1e843364e754a --- /dev/null +++ b/tests/PSFInterpTest/loadPSFSet.py @@ -0,0 +1,63 @@ +import unittest + +import sys,os,math +from itertools import islice + +import numpy as np + +import yaml +from ObservationSim.Config import Config +from ObservationSim.Config.Config import config_dir +from ObservationSim.Instrument import Chip + +from ObservationSim.PSF.PSFInterp import PSFInterp + + +def defineCCD(iccd, config_file): + with open(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) + path_dict = config_dir(config=config, work_dir=config['work_dir'], data_dir=config['data_dir']) + chip = Chip(chipID=iccd, config=config) + #chip = Chip(chipID=iccd, 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) + return chip + +def loadPSFSet(iccd): + config_file = "/public/home/weichengliang/CSST_git/newVersion/CSST/config/config_C3.yaml" + chip = defineCCD(iccd, config_file) + print(chip.chipID) + print(chip.cen_pix_x, chip.cen_pix_y) + + ipath = '/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube' + psfMat= PSFInterp(chip, npsf=900, PSF_data_file=ipath, PSF_data_prefix="") + psfSet= psfMat._loadPSF(iccd, ipath, PSF_data_prefix="") + + twave = 0 #[0...3] + tpsf = 0 #[0...899] + field_x= psfSet[twave][tpsf]['field_x'] + field_y= psfSet[twave][tpsf]['field_y'] + image_x= psfSet[twave][tpsf]['image_x'] + image_y= psfSet[twave][tpsf]['image_y'] + centroid_x= psfSet[twave][tpsf]['centroid_x'] + centroid_y= psfSet[twave][tpsf]['centroid_y'] + print("pos_info:", field_x, field_y, image_x, image_y, centroid_x, centroid_y) + return psfSet + + + +class PSFInterpModule_coverage(unittest.TestCase): + def test_psfEll_(self): + iccd = 1 #[1...30] + psfSet = loadPSFSet(iccd) + + +if __name__ == '__main__': + unittest.main() + print('#####haha#####') + + +