import os import numpy as np import random import galsim import h5py as h5 import healpy as hp from .Star import Star from .Galaxy import Galaxy from .Quasar import Quasar from ._util import seds, sed_assign, extAv from astropy.table import Table from astropy.coordinates import spherical_to_cartesian 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"]: 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() 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) # print(ngals) 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']) # print(nstars) 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 = Table.read(self.galaxy_path) # stars = Table.read(self.star_path) # self.ngals = gals['galaxyID'].size # self.nstars = stars['sourceID'].size 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 # for igals in range(self.ngals): # param = {} # param['id'] = igals + 1 # 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'] = avGal[int(self.ud()*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 # param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] # if param['star'] == 0: # obj = Galaxy(param, self.rotation) # self.objs.append(obj) # if param['star'] == 2: # obj = Quasar(param) # self.objs.append(obj) # for istars in range(self.nstars): # param = {} # param['id'] = self.ngals + istars + 1 # param['ra'] = stars['RA'][istars] # param['dec'] = stars['Dec'][istars] # param['sed_type'] = stars['sourceID'][istars] # param['model_tag'] = stars['model_tag'][istars] # param['z'] = 0.0 # param['star'] = 1 # Star # param['mag_use_normal'] = stars['app_sdss_g'][istars] # obj = Star(param) # self.objs.append(obj) # del gals # del stars # del avGal