Catalog.py

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')