diff --git a/Catalog/C6_50sqdeg.py b/Catalog/C6_50sqdeg.py deleted file mode 100644 index 4ce9f3ce95089cd7377957d856bcc6efa5fa73c9..0000000000000000000000000000000000000000 --- a/Catalog/C6_50sqdeg.py +++ /dev/null @@ -1,476 +0,0 @@ -import os -import galsim -import random -import copy -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U -from astropy.coordinates import SkyCoord -from astropy.io import fits - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5','galaxies_C6_bundle000241.h5','galaxies_C6_bundle000242.h5','galaxies_C6_bundle000287.h5','galaxies_C6_bundle000288.h5','galaxies_C6_bundle000714.h5','galaxies_C6_bundle000715.h5','galaxies_C6_bundle000778.h5','galaxies_C6_bundle000779.h5','galaxies_C6_bundle000842.h5','galaxies_C6_bundle000843.h5','galaxies_C6_bundle002046.h5','galaxies_C6_bundle002110.h5','galaxies_C6_bundle002111.h5','galaxies_C6_bundle002173.h5','galaxies_C6_bundle002174.h5','galaxies_C6_bundle002238.h5','galaxies_C6_bundle002596.h5','galaxies_C6_bundle002597.h5','galaxies_C6_bundle002656.h5','galaxies_C6_bundle002657.h5','galaxies_C6_bundle002711.h5','galaxies_C6_bundle002712.h5','galaxies_C6_bundle002844.h5','galaxies_C6_bundle002845.h5','galaxies_C6_bundle002884.h5','galaxies_C6_bundle002885.h5','galaxies_C6_bundle002921.h5','galaxies_C6_bundle002922.h5'] - -qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits'] - -star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5'] -star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)] - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -def get_agnsed_file(bundle_file_name): - return qsosed_file_list[bundle_file_list.index(bundle_file_name)] - -def get_star_cat(ra_pointing, dec_pointing): - pointing_c = SkyCoord(ra=ra_pointing*U.deg, dec=dec_pointing*U.deg) - max_dist = 10 - return_star_path = None - for star_file, center in zip(star_file_list, star_center_list): - center_c = SkyCoord(ra=center[0]*U.deg, dec=center[1]*U.deg) - dist = pointing_c.separation(center_c).to(U.deg).value - if dist < max_dist: - return_star_path = star_file - max_dist = dist - return return_star_path - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = config["catalog_options"]["input_path"]["cat_dir"] - - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - # Get the cloest star catalog file - star_file_name = get_star_cat(ra_pointing=self.pointing.ra, dec_pointing=self.pointing.dec) - star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name) - self.star_path = os.path.join(self.cat_dir, star_path) - self.star_SED_path = config["catalog_options"]["SED_templates_path"]["star_SED"] - self._load_SED_lib_star() - - if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = config["catalog_options"]["SED_templates_path"]["galaxy_SED"] - self._load_SED_lib_gals() - self.agn_seds = {} - - if "AGN_SED" in config["catalog_options"]["SED_templates_path"] and not config["catalog_options"]["star_only"]: - self.AGN_SED_path = config["catalog_options"]["SED_templates_path"]["AGN_SED"] - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_gals(self, gals, pix_id=None, cat_id=0, agnsed_file=""): - ngals = len(gals['ra']) - - # Apply astrometric modeling - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - if self.filt.filter_type == 'NUV': - param['mag_use_normal'] = gals['mag_csst_nuv'][igals] - else: - param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - param['g1'] = gals['shear'][igals][0] - param['g2'] = gals['shear'][igals][1] - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity( - e1=gals['ellipticity_true'][igals][0], - e2=gals['ellipticity_true'][igals][1], - rotation=self.rotation, - unit='radians') - # param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - # phi_e = cmath.phase(complex(param['e1'], param['e2'])) - # param['e1'] = param['ell_total'] * np.cos(phi_e + 2*self.rotation) - # param['e2'] = param['ell_total'] * np.sin(phi_e + 2*self.rotation) - - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = 1. - param['bulge_sersic_idx'] = 4. - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - # TEMP - self.ids += 1 - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - # Is this an Quasar? - param['qsoindex'] = gals['qsoindex'][igals] - if param['qsoindex'] == -1: - param['star'] = 0 # Galaxy - param['agnsed_file'] = "" - obj = Galaxy(param, logger=self.logger) - else: - param_qso = copy.deepcopy(param) - param_qso['star'] = 2 # Quasar - param_qso['agnsed_file'] = agnsed_file - # First add QSO model - obj = Quasar(param_qso, logger=self.logger) - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0, 0.) - self.objs.append(obj) - # Then add host galaxy model - param['star'] = 0 # Galaxy - param['agnsed_file'] = "" - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns for (host) galaxy - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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] - self.ids += 1 - param['id'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - bundle_file = "galaxies_C6_bundle{:06}.h5".format(bundleID) - file_path = os.path.join(self.galaxy_path, bundle_file) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - - # Get corresponding AGN SED file - agnsed_file = get_agnsed_file(bundle_file) - agnsed_path = os.path.join(self.AGN_SED_path, agnsed_file) - self.agn_seds[agnsed_file] = fits.open(agnsed_path)[0].data - - self._load_gals(gals, pix_id=pix, cat_id=bundleID, agnsed_file=agnsed_file) - - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.agn_seds[obj.agnsed_file][int(obj.qsoindex)] * 1e-17 - flux[flux < 0] = 0. - wave = self.lamb_gal * (1.0 + obj.z) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/C6_Catalog.py b/Catalog/C6_Catalog.py deleted file mode 100644 index 481494504d0964dc106c76960bd4d1deee3e6e6a..0000000000000000000000000000000000000000 --- a/Catalog/C6_Catalog.py +++ /dev/null @@ -1,514 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: - AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] - self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) - self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - param['g1'] = gals['shear'][igals][0] - param['g2'] = gals['shear'][igals][1] - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = 1. - param['bulge_sersic_idx'] = 4. - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: - try: - self._load_AGNs() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/C6_fits_Catalog.py b/Catalog/C6_fits_Catalog.py deleted file mode 100644 index d8eea0c6b53a834dddfb3ab226e1718e5940e727..0000000000000000000000000000000000000000 --- a/Catalog/C6_fits_Catalog.py +++ /dev/null @@ -1,636 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar, Stamp -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -import astropy.io.fits as fitsio -from ObservationSim.MockObject._util import seds, sed_assign, extAv - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - # (TEST) - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path: - self.normF_galaxy = Table.read(str(filter_path)) - - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: - AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] - self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) - self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - ###mock_stamp_START - if "stamp_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["stamp_cat"] and config["catalog_options"]["stamp_yes"]: - stamp_file = config["catalog_options"]["input_path"]["stamp_cat"] - self.stamp_path = os.path.join(self.cat_dir, stamp_file) - self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir="/share/simudata/CSSOSDataProductsSims/data/Templates/Galaxy/") #only for test - ###mock_stamp_END - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - # self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - ###mock_stamp_START - elif obj.type == "stamp": - return self.normF_galaxy ###normalize_filter for stamp - #return None - ###mock_stamp_END - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - if igals > 2000: - break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - param['g1'] = gals['shear'][igals][0] - param['g2'] = gals['shear'][igals][1] - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - if istars > 100: - break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - if iAGNs > 100: - break - - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - ###mock_stamp_START - def _load_stamps(self, stamps, pix_id=None): - nstamps = len(stamps['filename']) - 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 istamp in range(nstamps): - print('DEBUG:::istamp=', istamp) - - fitsfile = os.path.join(self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8')) - hdu=fitsio.open(fitsfile) - - param = self.initialize_param() - param['id'] = hdu[0].header['index'] #istamp - param['star'] = 3 # Stamp type in .cat file - ###param['lensGalaxyID'] = hdu[0].header['lensGID'] - param['ra'] = hdu[0].header['ra'] - param['dec']= hdu[0].header['dec'] - param['pixScale']= hdu[0].header['pixScale'] - #param['srcGalaxyID'] = hdu[0].header['srcGID'] - #param['mu']= hdu[0].header['mu'] - #param['PA']= hdu[0].header['PA'] - #param['bfrac']= hdu[0].header['bfrac'] - #param['z']= hdu[0].header['z'] - param['mag_use_normal'] = 20 #hdu[0].header['m_normal'] #gals['mag_true_g_lsst'] - - ###assert(stamps['lensGID'][istamp] == param['lensGalaxyID']) - - # Apply astrometric modeling - # in C3 case only aberration - param['ra_orig'] = param['ra'] - param['dec_orig']= param['dec'] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = [param['ra']] #ra_arr.tolist() - dec_list= [param['dec']] #dec_arr.tolist() - pmra_list = np.zeros(1).tolist() - pmdec_list = np.zeros(1).tolist() - rv_list = np.zeros(1).tolist() - parallax_list = [1e-9] * 1 - 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=1, - 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 - ) - param['ra'] = ra_arr[0] - param['dec']= dec_arr[0] - - # 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'] = 0.0 - param['redden'] = 0 - - #param["CSSTmag"]= True - #param["mag_r"] = 20. - #param[''] - ###more keywords for stamp### - param['image'] = hdu[0].data - param['image'] = param['image']/(np.sum(param['image'])) - obj = Stamp(param) - self.objs.append(obj) - ###mock_stamp_END - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: - try: - self._load_AGNs() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - ###mock_stamp_START - if "stamp_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["stamp_cat"] and self.config["catalog_options"]["stamp_yes"]: - stamps_cat = h5.File(self.stamp_path, 'r')['Stamps'] - for pix in self.pix_list: - try: - stamps = stamps_cat[str(pix)] - self._load_stamps(stamps, pix_id=pix) - del stamps - except Exception as e: - self.logger.error(str(e)) - print(e) - ###mock_stamp_END - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - ###mock_stamp_START - elif obj.type == 'stamp': - sed_data = getObservedSED( - sedCat=self.tempSed_gal[obj.sed_type], - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - ###mock_stamp_END - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/C6_50sqdeg_ns.py b/Catalog/C9_Catalog.py similarity index 100% rename from Catalog/C6_50sqdeg_ns.py rename to Catalog/C9_Catalog.py diff --git a/Catalog/Calib_Catalog.py b/Catalog/Calib_Catalog.py deleted file mode 100644 index 308a9d5fd08146b1e266ba41f692cb760f9395b4..0000000000000000000000000000000000000000 --- a/Catalog/Calib_Catalog.py +++ /dev/null @@ -1,673 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: - AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] - self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) - self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - if "CALIB_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"][ - "CALIB_cat"] and not config["catalog_options"]["star_only"]: - self.CALIB_cat_path = os.path.join(config["data_dir"], - config["catalog_options"]["input_path"]["CALIB_cat"]) - self.CALIB_SED_path = os.path.join(config["data_dir"], - config["catalog_options"]["SED_templates_path"]["CALIB_SED"]) - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - param['g1'] = gals['shear'][igals][0] - param['g2'] = gals['shear'][igals][1] - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = 1. - param['bulge_sersic_idx'] = 4. - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - def _load_calibObj(self): - - data = Table.read(self.CALIB_cat_path) - ra_arr = data['RA'] - dec_arr = data['DEC'] - - ngals = len(data) - - # Apply astrometric modeling - # in C3 case only aberration - # ra_arr = gals['ra'][:] - # dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = data['RA'][igals] - param['dec_orig'] = data['DEC'][igals] - param['mag_use_normal'] = data['MAG_g'][igals] - # if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - # continue - - param['z'] = -99 - param['model_tag'] = 'None' - param['g1'] = 0 - param['g2'] = 0 - param['kappa'] = 0 - param['e1'] = 0 - param['e2'] = 0 - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1'] ** 2 + param['e2'] ** 2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = 0 - param['e2_disk'] = 0 - param['e1_bulge'] = 0 - param['e2_bulge'] = 0 - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - # param['bulgemass'] = gals['bulgemass'][igals] - # param['diskmass'] = gals['diskmass'][igals] - - # param['size'] = gals['size'][igals] - # if param['size'] > self.max_size: - # self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = data['SERSIC_N'][igals] - param['bulge_sersic_idx'] = 1. - param['hlr_bulge'] = data['RE'][igals] - param['hlr_disk'] = data['RE'][igals] - param['bfrac'] = 0 - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 4 - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = data['SPEC_FN'][igals][0:-5] - - if param['star'] == 4: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - # obj.additional_output_str = self.add_fmt % ("n", 0., 0., 0., - # param['bulgemass'], param['diskmass'], param['detA'], - # param['e1'], param['e2'], param['kappa'], param['g1'], - # param['g2'], param['size'], - # param['galType'], param['veldisp']) - - self.objs.append(obj) - - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: - try: - self._load_AGNs() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if "CALIB_cat" in self.config["catalog_options"]["input_path"] and \ - self.config["catalog_options"]["input_path"][ - "CALIB_cat"] and not self.config["catalog_options"]["star_only"]: - try: - self._load_calibObj() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - elif obj.type == 'calib': - data = Table.read(os.path.join(self.CALIB_SED_path,obj.id+'.fits')) - obj_w = data['WAVELENGTH'] - obj_f = data['FLUX'] - input_delt_w = np.min(obj_w[1:]-obj_w[0:-1]) - if input_delt_w > 0.5: - lamb = np.arange(2000, 11000 + 0.5, 0.5) - speci = interpolate.interp1d(obj_w, obj_f) - y1 = speci(lamb) - else: - lamb = obj_w - y1 = obj_f - # erg/s/cm2/A --> photon/s/m2/A - y1_phot = y1 * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, y1_phot]).T, - names=('WAVELENGTH', 'FLUX')) - - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), - interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - return sed - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/FGS_Catalog.py b/Catalog/FGS_Catalog.py deleted file mode 100644 index bdab18cd497a1b5b72a719961c83557cfe2db225..0000000000000000000000000000000000000000 --- a/Catalog/FGS_Catalog.py +++ /dev/null @@ -1,516 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - # (TEST) - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: - AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] - self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) - self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - # if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - # continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - param['g1'] = gals['shear'][igals][0] - param['g2'] = gals['shear'][igals][1] - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = 1. - param['bulge_sersic_idx'] = 4. - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: - try: - self._load_AGNs() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - # if obj.type == 'quasar': - if obj.type == 'galaxy' or obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/NGPCatalog.py b/Catalog/NGPCatalog.py deleted file mode 100644 index cf6ae97a2758d98ab827f5e94f06222c5a9cdd0b..0000000000000000000000000000000000000000 --- a/Catalog/NGPCatalog.py +++ /dev/null @@ -1,331 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import seds, sed_assign, extAv, tag_sed, getObservedSED -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - self.chip_output = chip_output - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path: - self.normF_galaxy = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_file = config["catalog_options"]["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["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - return self.normF_galaxy - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=self.galaxy_SED_path) - - 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) - - # 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.utcfromtimestamp(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'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra_true'][igals] - param['dec_orig'] = gals['dec_true'][igals] - param['mag_use_normal'] = gals['mag_true_g_lsst'][igals] - # if param['mag_use_normal'] >= 26.5: - # continue - param['z'] = gals['redshift_true'][igals] - param['model_tag'] = 'None' - param['g1'] = 0 - param['g2'] = 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)] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - 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 - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - self.ids += 1 - # param['id'] = self.ids - param['id'] = gals['galaxyID'][igals] - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh']) - - self.objs.append(obj) - - def _load(self, **kwargs): - self.nav = 15005 - self.avGal = extAv(self.nav, seed=self.seed_Av) - self.objs = [] - self.ids = 0 - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - gals_cat = h5.File(self.galaxy_path, 'r')['galaxies'] - for pix in self.pix_list: - try: - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix) - del gals - except Exception as e: - self.logger.error(str(e)) - print(e) - if self.logger is not None: - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - del self.avGal - - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - sed_data = getObservedSED( - sedCat=self.tempSed_gal[obj.sed_type], - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 18001 + 0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photo/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - del wave - del flux - return sed diff --git a/Catalog/fd_test_C6.py b/Catalog/fd_test_C6.py deleted file mode 100644 index cbb41e88bcdeed54074a4ef2bca09f053ed1e31c..0000000000000000000000000000000000000000 --- a/Catalog/fd_test_C6.py +++ /dev/null @@ -1,524 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: - AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] - self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) - self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - # for igals in range(ngals): - for igals in range(0, ngals, 5): - # # (TEST) - # if igals > 1000: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - # param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - param['mag_use_normal'] = 20. - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - # param['g1'] = gals['shear'][igals][0] - # param['g2'] = gals['shear'][igals][1] - param['g1'] = 0. - param['g2'] = 0. - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - # param['e1_disk'] = param['e1'] - # param['e2_disk'] = param['e2'] - # param['e1_bulge'] = param['e1'] - # param['e2_bulge'] = param['e2'] - param['e1_disk'] = 0. - param['e2_disk'] = 0. - param['e1_bulge'] = 0. - param['e2_bulge'] = 0. - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - # param['size'] = gals['size'][igals] - param['size'] = 1. - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = 1. - param['bulge_sersic_idx'] = 4. - - # Sizes - # param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - param['bfrac'] = 0. - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - # if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - # star_cat = h5.File(self.star_path, 'r')['catalog'] - # for pix in self.pix_list: - # try: - # stars = star_cat[str(pix)] - # self._load_stars(stars, pix_id=pix) - # del stars - # except Exception as e: - # self.logger.error(str(e)) - # print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - # if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: - # try: - # self._load_AGNs() - # except Exception as e: - # traceback.print_exc() - # self.logger.error(str(e)) - # print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/testCat_galaxy.py b/Catalog/testCat_galaxy.py deleted file mode 100644 index fd9a5a41c686ebee91ab5633e760dc6902325551..0000000000000000000000000000000000000000 --- a/Catalog/testCat_galaxy.py +++ /dev/null @@ -1,447 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar, Stamp -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -import astropy.io.fits as fitsio -from ObservationSim.MockObject._util import seds, sed_assign, extAv - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - # (TEST) - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and config["catalog_options"]["galaxy_yes"]: - self.galaxy_path = os.path.join(self.cat_dir, config["catalog_options"]["input_path"]["galaxy_cat"]) - - self.galaxy_SED_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and config["catalog_options"]["galaxy_yes"]: - self.AGN_path = os.path.join(self.cat_dir, config["catalog_options"]["input_path"]["AGN_cat"]) - - self.AGN_SED_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - # self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - # if igals < 3447: - # continue - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - # if param['mag_use_normal'] >= 26.5: - # continue - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - param['g1'] = gals['shear'][igals][0] - param['g2'] = gals['shear'][igals][1] - param['kappa'] = gals['kappa'][igals] - param['e1'] = gals['ellipticity_true'][igals][0] - param['e2'] = gals['ellipticity_true'][igals][1] - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and self.config["catalog_options"]["galaxy_yes"]: - # for i in range(76): - # file_path = os.path.join(self.galaxy_path, "galaxies%04d_C6.h5"%(i)) - # gals_cat = h5.File(file_path, 'r')['galaxies'] - # for pix in self.pix_list: - # try: - # gals = gals_cat[str(pix)] - # self._load_gals(gals, pix_id=pix, cat_id=i) - # del gals - # except Exception as e: - # traceback.print_exc() - # self.logger.error(str(e)) - # print(e) - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - # # (TEST C6): - # if pix != 35421 or bundleID != 523: - # continue - gals_cat = h5.File(self.galaxy_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - #if "AGN_cat" in self.config["input_path"] and self.config["input_path"]["AGN_cat"] and not self.config["run_option"]["star_only"]: - if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and self.config["catalog_options"]["galaxy_yes"]: - try: - self._load_AGNs() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - # (TEST) - # def convert_mags(self, obj): - # spec = np.matmul(obj.coeff, self.pcs.T) - # lamb = self.lamb_gal * U.angstrom - # unit_sed = ((lamb * lamb)/constants.c*U.erg/U.second/U.cm**2/U.angstrom).to(U.jansky) - # unit_sed = unit_sed.value - # lamb = lamb.value - # lamb *= (1 + obj.z) - # spec *= (1 + obj.z) - # mags = -2.5 * np.log10(unit_sed*spec) + 8.9 - # mags += self.cosmo.distmod(obj.z).value - # return lamb, mags - - - def load_sed(self, obj, **kwargs): - if obj.type == 'galaxy' or obj.type == 'quasar': - # dist_L_pc = (1 + obj.z) * comoving_dist(z=obj.z)[0] - # factor = (10 / dist_L_pc)**2 - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - # flux = self.SED_AGN[int(obj.id)] * factor * 1e-17 - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/Catalog/testCat_star.py b/Catalog/testCat_star.py deleted file mode 100644 index 458e41cd9756c974204633bf4c1ea454726603c5..0000000000000000000000000000000000000000 --- a/Catalog/testCat_star.py +++ /dev/null @@ -1,378 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar, Stamp -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -import astropy.io.fits as fitsio -from ObservationSim.MockObject._util import seds, sed_assign, extAv - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -class Catalog(CatalogBase): - """An user customizable class for reading in catalog(s) of objects and SEDs. - - NOTE: must inherit the "CatalogBase" abstract class - - ... - - Attributes - ---------- - cat_dir : str - a directory that contains the catalog file(s) - star_path : str - path to the star catalog file - star_SED_path : str - path to the star SED data - objs : list - a list of ObservationSim.MockObject (Star, Galaxy, or Quasar) - NOTE: must have "obj" list when implement your own Catalog - - Methods - ---------- - load_sed(obj, **kwargs): - load the corresponding SED data for one object - load_norm_filt(obj): - load the filter throughput for the input catalog's photometric system. - """ - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - """Constructor method. - - Parameters - ---------- - config : dict - configuration dictionary which is parsed from the input YAML file - chip: ObservationSim.Instrument.Chip - an ObservationSim.Instrument.Chip instance, can be used to identify the band etc. - pointing: ObservationSim.Config.Pointing - an ObservationSim.Config.Pointing instance, can be used to configure the astrometry module - chip_output: ObservationSim.Config.ChipOutput - an ObservationSim.Config.ChipOutput instance, can be used to setup the output format - filt: ObservationSim.Instrument.Filter - an ObservationSim.Instrument.Filter instance, can be used to identify the filter type - **kwargs : dict - other needed input parameters (in key-value pairs), please modify corresponding - initialization call in "ObservationSim.py" as you need. - - Returns - ---------- - None - """ - - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - #self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and config["catalog_options"]["star_yes"]: - self.star_path = os.path.join(self.cat_dir, config["catalog_options"]["input_path"]["star_cat"]) - - self.star_SED_path = os.path.join(self.cat_dir, config["catalog_options"]["SED_templates_path"]["star_SED"]) - self._load_SED_lib_star() - - - if "rotateEll" in config["catalog_options"]: - self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - # self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - """Load the corresponding thourghput for the input magnitude "param["mag_use_normal"]". - - NOTE: if the input magnitude is already in CSST magnitude, simply return None - - Parameters - ---------- - obj : ObservationSim.MockObject - the object to get thourghput data for - - Returns - ---------- - norm_filt : Astropy.Table - the throughput Table with two columns (namely, "WAVELENGTH", "SENSITIVITY"): - norm_filt["WAVELENGTH"] : wavelengthes in Angstroms - norm_filt["SENSITIVITY"] : efficiencies - """ - if obj.type == "star": - return self.normF_star - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load(self, **kwargs): - """Read in all objects in from the catalog file(s). - - This is a must implemented method which is used to read in all objects, and - then convert them to ObservationSim.MockObject (Star, Galaxy, or Quasar). - - Currently, - the model of ObservationSim.MockObject.Star class requires: - param["star"] : int - specify the object type: 0: galaxy, 1: star, 2: quasar - param["id"] : int - ID number of the object - param["ra"] : float - Right ascension (in degrees) - param["dec"] : float - Declination (in degrees) - param["mag_use_normal"]: float - the absolute magnitude in a particular filter - NOTE: if that filter is not the corresponding CSST filter, the - load_norm_filt(obj) function must be implemented to load the filter - throughput of that particular photometric system - - the model of ObservationSim.MockObject.Galaxy class requires: - param["star"] : int - specify the object type: 0: galaxy, 1: star, 2: quasar - param["id"] : int - ID number of the object - param["ra"] : float - Right ascension (in degrees) - param["dec"] : float - Declination (in degrees) - param["mag_use_normal"]: float - the absolute magnitude in a particular filter - NOTE: if that filter is not the corresponding CSST filter, the - load_norm_filt(obj) function must be implemented to load the filter - throughput of that particular photometric system - param["bfrac"] : float - the bulge fraction - param["hlr_bulge"] : float - the half-light-radius of the bulge - param["hlr_disk"] : float - the half-light-radius of the disk - param["e1_bulge"], param["e2_bulge"] : float - the ellipticity of the bulge components - param["e1_disk"], param["e2_disk"] : float - the ellipticity of the disk components - (Optional parameters): - param['disk_sersic_idx']: float - Sersic index for galaxy disk component - param['bulge_sersic_idx']: float - Sersic index for galaxy bulge component - param['g1'], param['g2']: float - Reduced weak lensing shear components (valid for shear type: catalog) - the model of ObservationSim.MockObject.Galaxy class requires: - Currently a Quasar is modeled as a point source, just like a Star. - - NOTE: To construct an object, according to its type, just call: - Star(param), Galaxy(param), or Quasar(param) - - NOTE: All constructed objects should be appened to "self.objs". - - NOTE: Any other parameters can also be set within "param" dict: - Used to calculate required quantities and/or SEDs etc. - - Parameters - ---------- - **kwargs : dict - other needed input parameters (in key-value pairs), please modify corresponding - initialization call in "ObservationSim.py" as you need. - - Returns - ---------- - None - """ - self.objs = [] - self.ids = 0 - - #if "star_cat" in self.config["input_path"] and self.config["input_path"]["star_cat"] and not self.config["run_option"]["galaxy_only"]: - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and self.config["catalog_options"]["star_yes"]: - star_cat = h5.File(self.star_path, 'r')['stars'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - """Load the corresponding SED data for a particular obj. - - Parameters - ---------- - obj : ObservationSim.MockObject - the object to get SED data for - **kwargs : dict - other needed input parameters (in key-value pairs), please modify corresponding - initialization call in "ObservationSim.py" as you need. - - Returns - ---------- - sed : Astropy.Table - the SED Table with two columns (namely, "WAVELENGTH", "FLUX"): - sed["WAVELENGTH"] : wavelength in Angstroms - sed["FLUX"] : fluxes in photons/s/m^2/A - NOTE: the range of wavelengthes must at least cover [2450 - 11000] Angstorms - """ - - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - del wave - del flux - return sed diff --git a/Catalog/wcs_test_C6.py b/Catalog/wcs_test_C6.py deleted file mode 100644 index 5623ac23e71c4795dea569959058aa9957a2b042..0000000000000000000000000000000000000000 --- a/Catalog/wcs_test_C6.py +++ /dev/null @@ -1,520 +0,0 @@ -import os -import galsim -import random -import numpy as np -import h5py as h5 -import healpy as hp -import astropy.constants as cons -import traceback -from astropy.coordinates import spherical_to_cartesian -from astropy.table import Table -from scipy import interpolate -from datetime import datetime - -from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar -from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist -from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position - -# (TEST) -from astropy.cosmology import FlatLambdaCDM -from astropy import constants -from astropy import units as U - -try: - import importlib.resources as pkg_resources -except ImportError: - # Try backported to PY<37 'importlib_resources' - import importlib_resources as pkg_resources - -NSIDE = 128 - -def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7): - assert NSIDE == 2**healpixOrder - shift = healpixOrder - bundleOrder - shift = 2*shift - - nside_bundle = 2**bundleOrder - nside_healpix= 2**healpixOrder - - healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring) - bundleID_nest = (healpixID_nest >> shift) - bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest) - - return bundleID_ring - -class Catalog(CatalogBase): - def __init__(self, config, chip, pointing, chip_output, filt, **kwargs): - super().__init__() - self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]) - self.seed_Av = config["catalog_options"]["seed_Av"] - - self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111) - - self.chip_output = chip_output - self.filt = filt - self.logger = chip_output.logger - - with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: - self.normF_star = Table.read(str(filter_path)) - - self.config = config - self.chip = chip - self.pointing = pointing - - self.max_size = 0. - - if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]: - star_file = config["catalog_options"]["input_path"]["star_cat"] - star_SED_file = config["catalog_options"]["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["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]: - galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"] - self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir) - self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"]) - self._load_SED_lib_gals() - - if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]: - AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"] - self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"]) - self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"]) - self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"]) - self._load_SED_lib_AGN() - - if "rotateEll" in config["catalog_options"]: - self.rotation = np.radians(float(config["catalog_options"]["rotateEll"])) - else: - self.rotation = 0. - - # Update output .cat header with catalog specific output columns - self._add_output_columns_header() - - self._get_healpix_list() - self._load() - - def _add_output_columns_header(self): - self.add_hdr = " model_tag teff logg feh" - self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp " - - self.add_fmt = " %10s %8.4f %8.4f %8.4f" - self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f " - self.chip_output.update_output_header(additional_column_names=self.add_hdr) - - def _get_healpix_list(self): - self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2) - ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax - 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])) - # vertices = spherical_to_cartesian(1., dec, ra) - self.pix_list = hp.query_polygon( - NSIDE, - hp.ang2vec(np.radians(90.) - dec, ra), - inclusive=True - ) - if self.logger is not None: - msg = str(("HEALPix List: ", self.pix_list)) - self.logger.info(msg) - else: - print("HEALPix List: ", self.pix_list) - - def load_norm_filt(self, obj): - if obj.type == "star": - return self.normF_star - elif obj.type == "galaxy" or obj.type == "quasar": - # return self.normF_galaxy - return None - else: - return None - - def _load_SED_lib_star(self): - self.tempSED_star = h5.File(self.star_SED_path,'r') - - def _load_SED_lib_gals(self): - pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r") - lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r") - self.lamb_gal = lamb['lamb'][()] - self.pcs = pcs['pcs'][()] - - def _load_SED_lib_AGN(self): - from astropy.io import fits - self.SED_AGN = fits.open(self.AGN_SED_path)[0].data - self.lamb_AGN = np.load(self.AGN_SED_wave_path) - - - def _load_gals(self, gals, pix_id=None, cat_id=0): - ngals = len(gals['ra']) - - # Apply astrometric modeling - # in C3 case only aberration - ra_arr = gals['ra'][:] - dec_arr = gals['dec'][:] - 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.utcfromtimestamp(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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - - for igals in range(ngals): - # # (TEST) - # if igals > 100: - # break - - param = self.initialize_param() - param['ra'] = ra_arr[igals] - param['dec'] = dec_arr[igals] - param['ra_orig'] = gals['ra'][igals] - param['dec_orig'] = gals['dec'][igals] - param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals] - # param['mag_use_normal'] = 20. - if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]): - continue - - param['z'] = gals['redshift'][igals] - param['model_tag'] = 'None' - # param['g1'] = gals['shear'][igals][0] - # param['g2'] = gals['shear'][igals][1] - param['g1'] = 0. - param['g2'] = 0. - param['kappa'] = gals['kappa'][igals] - # param['e1'] = gals['ellipticity_true'][igals][0] - # param['e2'] = gals['ellipticity_true'][igals][1] - param['e1'] = 0. - param['e2'] = 0. - - # For shape calculation - - param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2) - if param['ell_total'] > 0.9: - continue - param['e1_disk'] = param['e1'] - param['e2_disk'] = param['e2'] - param['e1_bulge'] = param['e1'] - param['e2_bulge'] = param['e2'] - - - param['delta_ra'] = 0 - param['delta_dec'] = 0 - - # Masses - param['bulgemass'] = gals['bulgemass'][igals] - param['diskmass'] = gals['diskmass'][igals] - - param['size'] = gals['size'][igals] - if param['size'] > self.max_size: - self.max_size = param['size'] - - # Sersic index - param['disk_sersic_idx'] = 1. - param['bulge_sersic_idx'] = 4. - - # Sizes - param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass']) - if param['bfrac'] >= 0.6: - param['hlr_bulge'] = param['size'] - param['hlr_disk'] = param['size'] * (1. - param['bfrac']) - else: - param['hlr_disk'] = param['size'] - param['hlr_bulge'] = param['size'] * param['bfrac'] - - # SED coefficients - param['coeff'] = gals['coeff'][igals] - param['detA'] = gals['detA'][igals] - - # Others - param['galType'] = gals['type'][igals] - param['veldisp'] = gals['veldisp'][igals] - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 0 # Galaxy - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEMP - self.ids += 1 - # param['id'] = self.ids - param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) - - if param['star'] == 0: - obj = Galaxy(param, logger=self.logger) - - # Need to deal with additional output columns - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - param['bulgemass'], param['diskmass'], param['detA'], - param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'], - param['galType'], param['veldisp']) - - self.objs.append(obj) - - def _load_stars(self, stars, pix_id=None): - nstars = len(stars['sourceID']) - # Apply astrometric modeling - ra_arr = stars["RA"][:] - dec_arr = stars["Dec"][:] - pmra_arr = stars['pmra'][:] - pmdec_arr = stars['pmdec'][:] - rv_arr = stars['RV'][:] - parallax_arr = stars['parallax'][:] - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = pmra_arr.tolist() - pmdec_list = pmdec_arr.tolist() - rv_list = rv_arr.tolist() - parallax_list = parallax_arr.tolist() - dt = datetime.utcfromtimestamp(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=nstars, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for istars in range(nstars): - # # (TEST) - # if istars > 100: - # break - - 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] - param['pmra'] = pmra_arr[istars] - param['pmdec'] = pmdec_arr[istars] - param['rv'] = rv_arr[istars] - param['parallax'] = parallax_arr[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] - param['mag_use_normal'] = 20. - # if param['mag_use_normal'] >= 26.5: - # continue - self.ids += 1 - param['id'] = stars['sourceID'][istars] - 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, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'], - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - - self.objs.append(obj) - - def _load_AGNs(self): - data = Table.read(self.AGN_path) - ra_arr = data['ra'] - dec_arr = data['dec'] - nAGNs = len(data) - if self.config["obs_setting"]["enable_astrometric_model"]: - ra_list = ra_arr.tolist() - dec_list = dec_arr.tolist() - pmra_list = np.zeros(nAGNs).tolist() - pmdec_list = np.zeros(nAGNs).tolist() - rv_list = np.zeros(nAGNs).tolist() - parallax_list = [1e-9] * nAGNs - dt = datetime.utcfromtimestamp(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=nAGNs, - 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="J2000", - input_date_str=date_str, - input_time_str=time_str - ) - for iAGNs in range(nAGNs): - param = self.initialize_param() - param['ra'] = ra_arr[iAGNs] - param['dec'] = dec_arr[iAGNs] - param['ra_orig'] = data['ra'][iAGNs] - param['dec_orig'] = data['dec'][iAGNs] - param['z'] = data['z'][iAGNs] - param['appMag'] = data['appMag'][iAGNs] - param['absMag'] = data['absMag'][iAGNs] - - # NOTE: this cut cannot be put before the SED type has been assigned - if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): - continue - - # TEST no redening and no extinction - param['av'] = 0.0 - param['redden'] = 0 - - param['star'] = 2 # Quasar - param['id'] = data['igmlos'][iAGNs] - - if param['star'] == 2: - obj = Quasar(param, logger=self.logger) - - # Append additional output columns to the .cat file - obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.) - self.objs.append(obj) - - def _load(self, **kwargs): - self.objs = [] - self.ids = 0 - - if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]: - star_cat = h5.File(self.star_path, 'r')['catalog'] - for pix in self.pix_list: - try: - stars = star_cat[str(pix)] - self._load_stars(stars, pix_id=pix) - del stars - except Exception as e: - self.logger.error(str(e)) - print(e) - - if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]: - for pix in self.pix_list: - try: - bundleID = get_bundleIndex(pix) - file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID)) - gals_cat = h5.File(file_path, 'r')['galaxies'] - gals = gals_cat[str(pix)] - self._load_gals(gals, pix_id=pix, cat_id=bundleID) - del gals - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]: - try: - self._load_AGNs() - except Exception as e: - traceback.print_exc() - self.logger.error(str(e)) - print(e) - - if self.logger is not None: - self.logger.info("maximum galaxy size: %.4f"%(self.max_size)) - self.logger.info("number of objects in catalog: %d"%(len(self.objs))) - else: - print("number of objects in catalog: ", len(self.objs)) - - def load_sed(self, obj, **kwargs): - if obj.type == 'star': - _, wave, flux = tag_sed( - h5file=self.tempSED_star, - model_tag=obj.param['model_tag'], - teff=obj.param['teff'], - logg=obj.param['logg'], - feh=obj.param['feh'] - ) - elif obj.type == 'galaxy' or obj.type == 'quasar': - factor = 10**(-.4 * self.cosmo.distmod(obj.z).value) - if obj.type == 'galaxy': - flux = np.matmul(self.pcs, obj.coeff) * factor - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - sedcat = np.vstack((self.lamb_gal, flux)).T - sed_data = getObservedSED( - sedCat=sedcat, - redshift=obj.z, - av=obj.param["av"], - redden=obj.param["redden"] - ) - wave, flux = sed_data[0], sed_data[1] - elif obj.type == 'quasar': - flux = self.SED_AGN[int(obj.id)] * 1e-17 - # if np.any(flux < 0): - # raise ValueError("Glaxy %s: negative SED fluxes"%obj.id) - flux[flux < 0] = 0. - # sedcat = np.vstack((self.lamb_AGN, flux)).T - wave = self.lamb_AGN - # print("sed (erg/s/cm2/A) = ", sed_data) - # np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat) - else: - raise ValueError("Object type not known") - speci = interpolate.interp1d(wave, flux) - lamb = np.arange(2000, 11001+0.5, 0.5) - y = speci(lamb) - # erg/s/cm2/A --> photon/s/m2/A - all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13 - sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX')) - - if obj.type == 'quasar': - # integrate to get the magnitudes - sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T - sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest') - sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False) - interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full) - obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full) - # if obj.param['mag_use_normal'] >= 30: - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]])) - # print("obj ID = %d"%obj.id) - # print("mag_use_normal = %.3f"%obj.param['mag_use_normal']) - # print("integrated flux = %.7f"%(interFlux)) - # print("app mag = %.3f"%obj.param['appMag']) - # print("abs mag = %.3f"%obj.param['absMag']) - # mag = getABMAG(interFlux, self.filt.bandpass_full) - # print("mag diff = %.3f"%(mag - obj.param['mag_use_normal'])) - del wave - del flux - return sed diff --git a/config/config_50sqdeg.yaml b/config/config_50sqdeg.yaml deleted file mode 100644 index f28a4cc01bfdcc891b036a8be3fcde11a86d20cc..0000000000000000000000000000000000000000 --- a/config/config_50sqdeg.yaml +++ /dev/null @@ -1,222 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# CSST-Sim Group, 2023/04/25 -# -############################################### - -# 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: "/share/home/weichengliang/CSST_git/test_new_sim/outputs/" -data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "testRun0" - -# Whether to use MPI -run_option: - # 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 - -############################################### -# Catalog setting -############################################### -# Configure your catalog: options to be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "" - # star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" - star_cat: "starcat/" - galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" - # AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - - SED_templates_path: - star_SED: "Catalog_20210126/SpecLib.hdf5" - galaxy_SED: "Catalog_C6_20221212/sedlibs/" - AGN_SED: "qsocat/qsosed/" - # AGN_SED_WAVE: "wave_ross13.npy" - - # Only simulate stars? - star_only: YES - - # Only simulate galaxies? - galaxy_only: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - - # Options for survey types: - # "Photometric": simulate photometric chips only - # "Spectroscopic": simulate slitless spectroscopic chips only - # "FGS": simulate FGS chips only (31-42) - # "All": simulate full focal plane - survey_type: "Photometric" - - # Exposure time [seconds] - exp_time: 150. - - # Observation starting date & time - 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 - - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" - pointing_file: "pointing_50_combined.dat" - - # 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: [0] - - # 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: [8] - - # Whether to enable astrometric modeling - enable_astrometric_model: True - - # Whether to enable straylight model - enable_straylight_model: True - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - mag_sat_margin: -2.5 - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": from catalog - shear_type: "catalog" - - # For constant shear filed - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Instrumental effects setting -############################################### -ins_effects: - # switches - # Note: bias_16channel, gain_16channel, and shutter_effect - # is currently not applicable to "FGS" observations - 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 - add_prescan: OFF - bias_16channel: ON # Whether to add different biases for 16 channels - gain_16channel: ON # Whether to make different gains for 16 channels - 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 - format_output: OFF - - # Values: - # default values have been defined individually for each chip in: - # ObservationSim/Instrument/data/ccd/chip_definition.json - # Set them here will override the default 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-] - -############################################### -# Output options (for calibration pointings only) -############################################### -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 - 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 - -############################################### -# Random seeds -############################################### -random_seeds: - 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 -... diff --git a/config/config_C6.yaml b/config/config_C6.yaml deleted file mode 100644 index 5c7ebe6ede17d86d41d3366c3dfdd5627698844d..0000000000000000000000000000000000000000 --- a/config/config_C6.yaml +++ /dev/null @@ -1,223 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# CSST-Sim Group, 2023/04/25 -# -############################################### - -# 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: "/share/home/fangyuedong/new_sim/workplace/" -# work_dir: "/share/C6_new_sim_2sq" -data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "C6_new_sim_2sq_run1" -project_cycle: 6 -run_counter: 1 - -# Whether to use MPI -run_option: - # 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 - -############################################### -# Catalog setting -############################################### -# Configure your catalog: options to be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "Catalog_C6_20221212" - star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" - galaxy_cat: "cat2CSSTSim_bundle/" - AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - - SED_templates_path: - star_SED: "Catalog_20210126/SpecLib.hdf5" - galaxy_SED: "Catalog_C6_20221212/sedlibs/" - AGN_SED: "quickspeclib_ross13.fits" - AGN_SED_WAVE: "wave_ross13.npy" - - # Only simulate stars? - star_only: YES - - # Only simulate galaxies? - galaxy_only: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - - # Options for survey types: - # "Photometric": simulate photometric chips only - # "Spectroscopic": simulate slitless spectroscopic chips only - # "FGS": simulate FGS chips only (31-42) - # "All": simulate full focal plane - survey_type: "Photometric" - - # Exposure time [seconds] - exp_time: 150. - - # Observation starting date & time - 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 - - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" - pointing_file: "pointing_radec_246.5_40.dat" - - # 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: [0] - - # 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: [8] - - # Whether to enable astrometric modeling - enable_astrometric_model: True - - # Whether to enable straylight model - enable_straylight_model: True - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - # mag_sat_margin: -2.5 - mag_sat_margin: -15. - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": from catalog - shear_type: "catalog" - - # For constant shear filed - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Instrumental effects setting -############################################### -ins_effects: - # switches - # Note: bias_16channel, gain_16channel, and shutter_effect - # is currently not applicable to "FGS" observations - field_dist: YES # Whether to add field distortions - add_back: YES # Whether to add sky background - add_dark: YES # Whether to add dark noise - add_readout: YES # Whether to add read-out (Gaussian) noise - add_bias: YES # Whether to add bias-level to images - bias_16channel: YES # Whether to add different biases for 16 channels - gain_16channel: YES # Whether to make different gains for 16 channels - shutter_effect: YES # Whether to add shutter effect - flat_fielding: YES # Whether to add flat-fielding effect - prnu_effect: YES # Whether to add PRNU effect - non_linear: YES # Whether to add non-linearity - cosmic_ray: YES # Whether to add cosmic-ray - cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: YES # Whether to simulate CTE trails - saturbloom: YES # Whether to simulate Saturation & Blooming - add_badcolumns: YES # Whether to add bad columns - add_hotpixels: YES # Whether to add hot pixels - add_deadpixels: YES # Whether to add dead(dark) pixels - bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect - - # Values: - # default values have been defined individually for each chip in: - # ObservationSim/Instrument/data/ccd/chip_definition.json - # Set them here will override the default 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-] - -############################################### -# Output options (for calibration pointings only) -############################################### -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 - 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 - -############################################### -# Random seeds -############################################### -random_seeds: - 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 -... diff --git a/config/config_C6_dev.yaml b/config/config_C6_dev.yaml deleted file mode 100644 index dd619d58f6f4607253e399cccccaa61167147d14..0000000000000000000000000000000000000000 --- a/config/config_C6_dev.yaml +++ /dev/null @@ -1,234 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# CSST-Sim Group, 2023/04/25 -# -############################################### - -# 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: "/share/home/fangyuedong/20231211/workplace/" -data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "test_spec" -project_cycle: 6 -run_counter: 1 - -# Whether to use MPI -run_option: - # 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 - -############################################### -# Catalog setting -############################################### -# Configure your catalog: options to be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "Catalog_C6_20221212" - star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" - galaxy_cat: "cat2CSSTSim_bundle/" - AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - - SED_templates_path: - star_SED: "Catalog_20210126/SpecLib.hdf5" - galaxy_SED: "Catalog_C6_20221212/sedlibs/" - AGN_SED: "quickspeclib_ross13.fits" - AGN_SED_WAVE: "wave_ross13.npy" - - # Only simulate stars? - star_only: NO - - # Only simulate galaxies? - galaxy_only: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - - # Options for survey types: - # "Photometric": simulate photometric chips only - # "Spectroscopic": simulate slitless spectroscopic chips only - # "FGS": simulate FGS chips only (31-42) - # "All": simulate full focal plane - # "CALIBRATION": falt, bias, dark with or without postflash - survey_type: "Spectroscopic" - #"LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null - #'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', - #'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050','LED13': '340', 'LED14': '365' - LED_TYPE: ['LED5'] - LED_TIME: [1.] - # unit e- ,flat level - FLAT_LEVEL: 20000 - FLAT_LEVEL_FIL: 'g' - - # Exposure time [seconds] - exp_time: 150. - - # Observation starting date & time - 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 - - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" - pointing_file: "pointing_radec_246.5_40.dat" - - # 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: [0] - - # 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: [1] - - # Whether to enable astrometric modeling - enable_astrometric_model: True - - # Whether to enable straylight model - enable_straylight_model: True - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - # mag_sat_margin: -2.5 - mag_sat_margin: -15. - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": from catalog - shear_type: "catalog" - - # For constant shear filed - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Instrumental effects setting -############################################### -ins_effects: - # switches - # Note: bias_16channel, gain_16channel, and shutter_effect - # is currently not applicable to "FGS" observations - field_dist: YES # Whether to add field distortions - add_back: YES # Whether to add sky background - add_dark: YES # Whether to add dark noise - add_readout: YES # Whether to add read-out (Gaussian) noise - add_bias: YES # Whether to add bias-level to images - add_prescan: YES - bias_16channel: YES # Whether to add different biases for 16 channels - gain_16channel: YES # Whether to make different gains for 16 channels - shutter_effect: YES # Whether to add shutter effect - flat_fielding: YES # Whether to add flat-fielding effect - prnu_effect: YES # Whether to add PRNU effect - non_linear: YES # Whether to add non-linearity - cosmic_ray: YES # Whether to add cosmic-ray - cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: YES # Whether to simulate CTE trails, CTI_lgl_v0.3.tar.gz - saturbloom: YES # Whether to simulate Saturation & Blooming - add_badcolumns: YES # Whether to add bad columns - add_hotpixels: YES # Whether to add hot pixels - add_deadpixels: YES # Whether to add dead(dark) pixels - bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect - add_prescan: YES # Whether to add pre/over-scan - format_output: YES ##1*16 output - - # Values: - # default values have been defined individually for each chip in: - # ObservationSim/Instrument/data/ccd/chip_definition.json - # Set them here will override the default 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-] - -############################################### -# Output options (for calibration pointings only) -############################################### -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 - 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 - -############################################### -# Random seeds -############################################### -random_seeds: - 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 -... diff --git a/config/config_C6_fits.yaml b/config/config_C6_fits.yaml deleted file mode 100644 index 2c91057b4b51c376a6905c11f2daabad6bc92ebc..0000000000000000000000000000000000000000 --- a/config/config_C6_fits.yaml +++ /dev/null @@ -1,241 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# CSST-Sim Group, 2023/04/25 -# -############################################### - -# 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: "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C8/" -data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "testRun_FGS" -project_cycle: 6 -run_counter: 1 - -# Whether to use MPI -run_option: - # 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 - - #output PSF - out_psf: YES - -############################################### -# Catalog setting -############################################### -# Configure your catalog: options to be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "Catalog_C6_20221212" - star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" - galaxy_cat: "cat2CSSTSim_bundle/" - AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - stamp_cat: "stampCatsIndex.hdf5" - - SED_templates_path: - star_SED: "Catalog_20210126/SpecLib.hdf5" - galaxy_SED: "Catalog_C6_20221212/sedlibs/" - AGN_SED: "quickspeclib_ross13.fits" - AGN_SED_WAVE: "wave_ross13.npy" - #stamp_SED: - - # Only simulate stars? - star_only: YES - - # Only simulate galaxies? - galaxy_only: NO - - # With stamp? - stamp_yes: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - - # Options for survey types: - # "Photometric": simulate photometric chips only - # "Spectroscopic": simulate slitless spectroscopic chips only - # "FGS": simulate FGS chips only (31-42) - # "All": simulate full focal plane - # "CALIBRATION": falt, bias, dark with or without postflash - survey_type: "All" - # "LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null - # 'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', - # 'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050','LED13': '340', 'LED14': '365' - # LED_TYPE: ['LED5'] - # LED_TIME: [1.] - # # unit e- ,flat level - # FLAT_LEVEL: 20000 - # FLAT_LEVEL_FIL: 'g' - - # Exposure time [seconds] - exp_time: 150. - - # Observation starting date & time - 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 - - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" - pointing_file: "pointing_radec_246.5_40.dat" - - # 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: [0] - - # 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: [39,40,41,42] - - # Whether to enable astrometric modeling - enable_astrometric_model: True - - # Whether to enable straylight model - enable_straylight_model: True - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - # mag_sat_margin: -2.5 - mag_sat_margin: -15. - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": from catalog - shear_type: "catalog" - - # For constant shear filed - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Instrumental effects setting -############################################### -ins_effects: - # switches - # Note: bias_16channel, gain_16channel, and shutter_effect - # is currently not applicable to "FGS" observations - field_dist: YES # Whether to add field distortions - add_back: YES # Whether to add sky background - add_dark: YES # Whether to add dark noise - add_readout: YES # Whether to add read-out (Gaussian) noise - add_bias: YES # Whether to add bias-level to images - bias_16channel: YES # Whether to add different biases for 16 channels - gain_16channel: YES # Whether to make different gains for 16 channels - shutter_effect: NO # Whether to add shutter effect - flat_fielding: YES # Whether to add flat-fielding effect - prnu_effect: YES # Whether to add PRNU effect - non_linear: YES # Whether to add non-linearity - cosmic_ray: YES # Whether to add cosmic-ray - cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: NO # Whether to simulate CTE trails, CTI_lgl_v0.3.tar.gz - saturbloom: YES # Whether to simulate Saturation & Blooming - add_badcolumns: YES # Whether to add bad columns - add_hotpixels: YES # Whether to add hot pixels - add_deadpixels: YES # Whether to add dead(dark) pixels - bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect - add_prescan: NO # Whether to add pre/over-scan - format_output: NO ##1*16 output - - # Values: - # default values have been defined individually for each chip in: - # ObservationSim/Instrument/data/ccd/chip_definition.json - # Set them here will override the default values - # dark_exptime: 300 # Exposure time for dark current frames [seconds] - # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] - # readout_time: 0.01 # The read-out time for each channel [seconds] - # df_strength: 2.3 # Sillicon sensor diffusion strength - # bias_level: 2000 # bias level [e-/pixel] - # gain: 1. # Gain - # full_well: 90000 # Full well depth [e-] - -############################################### -# Output options (for calibration pointings only) -############################################### -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 - 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 - -############################################### -# Random seeds -############################################### -random_seeds: - 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 -... diff --git a/config/config_ooc_c6.yaml b/config/config_ooc_c6.yaml deleted file mode 100644 index 54d6ef5a0809c08f3a37f949e3e1839aac0e0626..0000000000000000000000000000000000000000 --- a/config/config_ooc_c6.yaml +++ /dev/null @@ -1,233 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# CSST-Sim Group, 2023/04/25 -# -############################################### - -# 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: "/share/home/fangyuedong/20231211/workplace/" -data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "test20231218_c6_onlyCat" -project_cycle: 6 -run_counter: 1 - -# Whether to use MPI -run_option: - # 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: YES - -############################################### -# Catalog setting -############################################### -# Configure your catalog: options to be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "Catalog_C6_20221212" - star_cat: "C6_MMW_GGC_Astrometry_healpix.hdf5" - galaxy_cat: "cat2CSSTSim_bundle/" - AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - - SED_templates_path: - star_SED: "Catalog_20210126/SpecLib.hdf5" - galaxy_SED: "Catalog_C6_20221212/sedlibs/" - AGN_SED: "quickspeclib_ross13.fits" - AGN_SED_WAVE: "wave_ross13.npy" - - # Only simulate stars? - star_only: YES - - # Only simulate galaxies? - galaxy_only: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - - # Options for survey types: - # "Photometric": simulate photometric chips only - # "Spectroscopic": simulate slitless spectroscopic chips only - # "FGS": simulate FGS chips only (31-42) - # "All": simulate full focal plane - # "CALIBRATION": falt, bias, dark with or without postflash - survey_type: "Photometric" - #"LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null - #'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', - #'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050','LED13': '340', 'LED14': '365' - LED_TYPE: ['LED5'] - LED_TIME: [1.] - # unit e- ,flat level - FLAT_LEVEL: 20000 - FLAT_LEVEL_FIL: 'g' - - # Exposure time [seconds] - exp_time: 150. - - # Observation starting date & time - 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 - - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" - pointing_file: "pointing_radec_246.5_40.dat" - - # 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: [0] - - # 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: [8] - - # Whether to enable astrometric modeling - enable_astrometric_model: False - - # Whether to enable straylight model - enable_straylight_model: True - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - mag_sat_margin: -2.5 - # mag_sat_margin: -15. - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": from catalog - shear_type: "constant" - - # For constant shear filed - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Instrumental effects setting -############################################### -ins_effects: - # switches - # Note: bias_16channel, gain_16channel, and shutter_effect - # is currently not applicable to "FGS" observations - field_dist: YES # Whether to add field distortions - add_back: YES # Whether to add sky background - add_dark: YES # Whether to add dark noise - add_readout: YES # Whether to add read-out (Gaussian) noise - add_bias: YES # Whether to add bias-level to images - bias_16channel: YES # Whether to add different biases for 16 channels - gain_16channel: YES # Whether to make different gains for 16 channels - shutter_effect: YES # Whether to add shutter effect - flat_fielding: YES # Whether to add flat-fielding effect - prnu_effect: YES # Whether to add PRNU effect - non_linear: YES # Whether to add non-linearity - cosmic_ray: YES # Whether to add cosmic-ray - cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: YES # Whether to simulate CTE trails, CTI_lgl_v0.3.tar.gz - saturbloom: YES # Whether to simulate Saturation & Blooming - add_badcolumns: YES # Whether to add bad columns - add_hotpixels: YES # Whether to add hot pixels - add_deadpixels: YES # Whether to add dead(dark) pixels - bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect - add_prescan: YES # Whether to add pre/over-scan - format_output: YES ##1*16 output - - # Values: - # default values have been defined individually for each chip in: - # ObservationSim/Instrument/data/ccd/chip_definition.json - # Set them here will override the default 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-] - -############################################### -# Output options (for calibration pointings only) -############################################### -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 - 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 - -############################################### -# Random seeds -############################################### -random_seeds: - 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 -... diff --git a/config/config_overall.yaml b/config/config_overall.yaml index 85b6e933a83fe2a7a7abb0af644f993bf78a0530..70985703552d19c40615e7c148b5bbd09192fe42 100644 --- a/config/config_overall.yaml +++ b/config/config_overall.yaml @@ -32,12 +32,10 @@ run_option: catalog_options: input_path: cat_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/" - # star_cat: "starcat/" star_cat: "starcat_C9/" galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" SED_templates_path: - # star_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/SpecLib.hdf5" star_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/starcat_C9/" galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/" AGN_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/qsocat/qsosed/" @@ -59,10 +57,9 @@ obs_setting: # if you just want to run default pointing: # - pointing_dir: null # - pointing_file: null - # pointing_file: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing50_C9/pointing_50_1_n.dat" - pointing_file: "/public/home/fangyuedong/project/unit_test_data/pointing_for_test.dat" + pointing_file: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing50_C9/pointing_50_1_n.dat" - obs_config_file: "/public/home/fangyuedong/project/csst_msc_sim/config/obs_config_SCI_WIDE_phot.yaml" + obs_config_file: "/public/home/fangyuedong/project/csst_msc_sim/config/obs_config_SCI.yaml" # Run specific pointing(s): # - give a list of indexes of pointings: [ip_1, ip_2...] diff --git a/config/config_overall_newStar.yaml b/config/config_overall_newStar.yaml deleted file mode 100644 index fb6fbbb26e38c1060e38514287916cba8a70d375..0000000000000000000000000000000000000000 --- a/config/config_overall_newStar.yaml +++ /dev/null @@ -1,141 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# Overall settings -# CSST-Sim Group, 2024/01/08 -# -############################################### - -# 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: "/home/zhangxin/CSST_SIM/CSST_sim_scheduler/" -run_name: "QSO_50sqdeg_test" - -# Project cycle and run counter are used to name the outputs -project_cycle: 9 -run_counter: 1 - -# Run options -run_option: - # Output catalog only? - # If yes, no imaging simulation will be run. Only the catalogs - # of corresponding footprints will be generated. - out_cat_only: NO - -############################################### -# Catalog setting -############################################### -# Configure the input catalog: options should be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "/nfsdata/share/CSSOSDataProductsSims/data/Catalog_C6_20221212/" - star_cat: "star_catalog/" - # galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" - - SED_templates_path: - star_SED: "/nfsdata/share/CSSOSDataProductsSims/data/Catalog_C6_20221212/star_catalog/" - # galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/" - # AGN_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/qsocat/qsosed/" - - # Only simulate stars? - star_only: YES - - # Only simulate galaxies? - galaxy_only: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - -############################################### -# Observation setting -############################################### -obs_setting: - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_file: "/nfsdata/share/CSSOSDataProductsSims/data/pointing_50_1_n.dat" - - obs_config_file: "/home/zhangxin/CSST_SIM/CSST_sim_scheduler/csst-simulation/config/obs_config_SCI_WIDE_phot.yaml" - - # 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: [1] - - # Whether to enable astrometric modeling - enable_astrometric_model: NO - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - mag_sat_margin: -2.5 - # mag_sat_margin: -15. - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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 models for photometry survey simulation - psf_pho_dir: "/nfsdata/share/CSSOSDataProductsSims/data/psfCube1" - # PSF models for slitless spectrum survey simulation - psf_sls_dir: "/nfsdata/share/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" - -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": get shear values from catalog - shear_type: "constant" - - # For constant shear field - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Output options -############################################### -output_setting: - output_format: "channels" # Whether to export as 16 channels (subimages) with pre- and over-scan ("image"/"channels") - shutter_output: NO # Whether to export shutter effect 16-bit image - prnu_output: NO # Whether to export the PRNU (pixel-to-pixel flat-fielding) files - -############################################### -# Random seeds -############################################### -random_seeds: - 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/config_overall_newStar_fits.yaml b/config/config_overall_newStar_fits.yaml index 91312e052facd509ae78d0b68653e1674e286e8a..cac80ea6c35676f921789e14cc4bb944cfab7c31 100644 --- a/config/config_overall_newStar_fits.yaml +++ b/config/config_overall_newStar_fits.yaml @@ -64,7 +64,7 @@ obs_setting: # - pointing_file: null pointing_file: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing50_C9/pointing_50_1_n.dat" - obs_config_file: "/public/home/chengliang/CSSOSDataProductsSims/csst_msc_sim/config/obs_config_SCI_WIDE_phot.yaml" + obs_config_file: "/public/home/chengliang/CSSOSDataProductsSims/csst_msc_sim/config/obs_config_SCI.yaml" # Run specific pointing(s): # - give a list of indexes of pointings: [ip_1, ip_2...] diff --git a/config/config_testCASE.yaml b/config/config_testCASE.yaml deleted file mode 100644 index f44cf5c972c1d0d53a5772a1b97c94e306917653..0000000000000000000000000000000000000000 --- a/config/config_testCASE.yaml +++ /dev/null @@ -1,221 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# CSST-Sim Group, 2023/04/30 -# -############################################### - -# 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: "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C6/" -data_dir: "/share/home/weichengliang/CSST_git/csst-simulation/tools/" -run_name: "C6_fits_testRun" - -# Whether to use MPI -run_option: - # 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 - -############################################### -# Catalog setting -############################################### -# Configure your catalog: options to be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - cat_dir: "Catalog_test" - star_cat: "testStarCat.h5" - galaxy_cat: "testGalaxyCat.h5" - AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - stamp_cat: "stampCatsIndex.hdf5" - - SED_templates_path: - star_SED: "SpecLib.hdf5" - galaxy_SED: "sedlibs/" - AGN_SED: "quickspeclib_ross13.fits" - AGN_SED_WAVE: "wave_ross13.npy" - #stamp_SED: - - # simulate stars? - star_yes: NO - - # simulate galaxies? - galaxy_yes: NO - - # With stamp? - stamp_yes: YES - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - - # Options for survey types: - # "Photometric": simulate photometric chips only - # "Spectroscopic": simulate slitless spectroscopic chips only - # "FGS": simulate FGS chips only (31-42) - # "All": simulate full focal plane - survey_type: "Photometric" - - # Exposure time [seconds] - exp_time: 150. - - # Observation starting date & time - date_obs: "210525" # [yymmdd] - time_obs: "120000" # [hhmmss] - - # Default Pointing [degrees] - # Note: NOT valid when a pointing list file is specified - ra_center: 60.11828178499743 #244.972773 - dec_center: -39.75897455697294 #39.895901 - # Image rotation [degree] - image_rot: 112.23685624012192 #109.59 - - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: null #"/share/simudata/CSSOSDataProductsSims/data/" - pointing_file: null #"pointing_radec_246.5_40.dat" - - # 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: [0] - - # 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: [7, 8, 9] - - # Whether to enable astrometric modeling - enable_astrometric_model: False - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - mag_sat_margin: -2.5 - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": from catalog - shear_type: "catalog" - - # For constant shear filed - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Instrumental effects setting -############################################### -ins_effects: - # switches - # Note: bias_16channel, gain_16channel, and shutter_effect - # is currently not applicable to "FGS" observations - field_dist: ON # Whether to add field distortions - add_back: OFF # Whether to add sky background - add_dark: OFF # Whether to add dark noise - add_readout: OFF # Whether to add read-out (Gaussian) noise - add_bias: OFF # Whether to add bias-level to images - bias_16channel: OFF # Whether to add different biases for 16 channels - gain_16channel: OFF # Whether to make different gains for 16 channels - shutter_effect: OFF # Whether to add shutter effect - flat_fielding: OFF # Whether to add flat-fielding effect - prnu_effect: OFF # Whether to add PRNU effect - non_linear: OFF # Whether to add non-linearity - cosmic_ray: OFF # Whether to add cosmic-ray - cray_differ: OFF # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: OFF # Whether to simulate CTE trails - saturbloom: OFF # Whether to simulate Saturation & Blooming - add_badcolumns: OFF # Whether to add bad columns - add_hotpixels: OFF # Whether to add hot pixels - add_deadpixels: OFF # Whether to add dead(dark) pixels - bright_fatter: OFF # Whether to simulate Brighter-Fatter (also diffusion) effect - - # Values: - # default values have been defined individually for each chip in: - # ObservationSim/Instrument/data/ccd/chip_definition.json - # Set them here will override the default values - # dark_exptime: 300 # Exposure time for dark current frames [seconds] - # flat_exptime: 150 # Exposure time for flat-fielding frames [seconds] - # readout_time: 0.01 # The read-out time for each channel [seconds] - # df_strength: 2.3 # Sillicon sensor diffusion strength - # bias_level: 2000 # bias level [e-/pixel] - # gain: 1. # Gain - # full_well: 90000 # Full well depth [e-] - -############################################### -# Output options (for calibration pointings only) -############################################### -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 - 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 - -############################################### -# Random seeds -############################################### -random_seeds: - 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/obs_config_SCI_WIDE_phot.yaml b/config/obs_config_SCI.yaml similarity index 100% rename from config/obs_config_SCI_WIDE_phot.yaml rename to config/obs_config_SCI.yaml diff --git a/example.sh b/example.sh deleted file mode 100755 index 5e93491ffc5cf673dd25f3b49f681f1eecd9f47e..0000000000000000000000000000000000000000 --- a/example.sh +++ /dev/null @@ -1,8 +0,0 @@ -#!/bin/bash - -date - -python3 /share/home/fangyuedong/fgs_sim/csst-simulation/run_sim.py \ - --config_file config_example.yaml \ - --catalog Catalog_example \ - -c /share/home/fangyuedong/fgs_sim/csst-simulation/config diff --git a/mpi_example.pbs b/mpi_example.pbs deleted file mode 100755 index 2dbcc53e7ca19aca1437b42afb5d6fe97f0f6446..0000000000000000000000000000000000000000 --- a/mpi_example.pbs +++ /dev/null @@ -1,22 +0,0 @@ -#!/bin/bash - -#PBS -N SIMS -##PBS -l walltime=70:00:00 - -##mpdallexit -##mpdboot -n 10 -f ./mpd.hosts - -##PBS -j oe -#PBS -l nodes=comput108:ppn=80 -#####PBS -q longq -#PBS -q batch -#PBS -u fangyuedong - -NP=40 -date -echo $NP - -mpirun -np $NP python /share/home/fangyuedong/fgs_sim/csst-simulation/run_sim.py \ - --config_file config_example.yaml \ - --catalog Catalog_example \ - -c /share/home/fangyuedong/fgs_sim/csst-simulation/config \ No newline at end of file diff --git a/profile_C6.sh b/profile_C6.sh deleted file mode 100755 index c7d03048dc0c4b4f85238f79d3c8e8d1f88e4710..0000000000000000000000000000000000000000 --- a/profile_C6.sh +++ /dev/null @@ -1,25 +0,0 @@ -#!/bin/bash - -date - -python -m cProfile -o C6_profiler_test.pstats /share/home/fangyuedong/new_sim/csst-simulation/run_sim.py \ - --config_file config_C6.yaml \ - --catalog C6_Catalog \ - -c /share/home/fangyuedong/new_sim/csst-simulation/config - - # --config_file config_test_new_sim.yaml \ - # --catalog New_sim_Test \ - # -c /share/home/fangyuedong/new_sim/csst-simulation/config - - # --config_file config_50sqdeg.yaml \ - # --catalog C6_50sqdeg \ - # -c /share/home/fangyuedong/new_sim/csst-simulation/config - - # --config_file config_fgs.yaml \ - # --catalog FGS_Catalog \ - # -c /share/home/fangyuedong/csst-simulation/config - - # --config_file test_fd_C6.yaml \ - # --catalog fd_test_C6 \ - # --config_file config_C6_test_wcs.yaml \ - # --catalog wcs_test_C6 \ diff --git a/profile_run.sh b/profile_run.sh new file mode 100755 index 0000000000000000000000000000000000000000..bde7199c25f44227bf106083a6b252afe7a765bb --- /dev/null +++ b/profile_run.sh @@ -0,0 +1,8 @@ +#!/bin/bash + +date + +python -m cProfile -o profiler_test.pstats /share/home/fangyuedong/new_sim/csst-simulation/run_sim.py \ + --config_file config_overall.yaml \ + --catalog C9_Catalog \ + -c /share/home/fangyuedong/new_sim/csst-simulation/config diff --git a/run_C6.pbs b/submit_jobs.pbs similarity index 85% rename from run_C6.pbs rename to submit_jobs.pbs index 32cfa73df9123ce6714f5d91c12b6e7cfdcff32b..34867ae2a24730ab5f485dd69725642343272687 100755 --- a/run_C6.pbs +++ b/submit_jobs.pbs @@ -12,4 +12,4 @@ date #限定单节点任务数 srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes -mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 6 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py --config_file config_overall.yaml --catalog C6_50sqdeg_ns -c /public/home/fangyuedong/project/csst_msc_sim/config \ No newline at end of file +mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 6 --map-by node python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py --config_file config_overall.yaml --catalog C9_Catalog -c /public/home/fangyuedong/project/csst_msc_sim/config \ No newline at end of file diff --git a/test_C6.sh b/test_C6.sh deleted file mode 100755 index 1ce5be63ec599042c43f680178ef7a0ce877721f..0000000000000000000000000000000000000000 --- a/test_C6.sh +++ /dev/null @@ -1,19 +0,0 @@ -#!/bin/bash - -date - -python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py \ - --config_file config_overall.yaml \ - -c /public/home/fangyuedong/project/csst_msc_sim/config \ - --catalog C6_50sqdeg - -# python3 /share/home/fangyuedong/20231211/csst-simulation/run_sim.py \ -# --config_file config_C6_dev.yaml \ -# --catalog C6_Catalog \ -# -c /share/home/fangyuedong/20231211/csst-simulation/config - -# python3 /share/home/fangyuedong/20231211/csst-simulation/run_sim.py \ -# --config_file config_ooc_c6.yaml \ -# --catalog C6_Catalog \ -# -c /share/home/fangyuedong/20231211/csst-simulation/config - diff --git a/test_run.sh b/test_run.sh new file mode 100755 index 0000000000000000000000000000000000000000..436bfaa90c61cce627534e4923440b06fe19c46e --- /dev/null +++ b/test_run.sh @@ -0,0 +1,9 @@ +#!/bin/bash + +date + +python3 /public/home/fangyuedong/project/csst_msc_sim/run_sim.py \ + --config_file config_overall.yaml \ + -c /public/home/fangyuedong/project/csst_msc_sim/config \ + --catalog C9_Catalog +