Merge remote-tracking branch 'origin/develop' (8a2f7687) · Commits · csst-sims / csst_msc_sim

catalog/C10_Catalog.py

+318 −180

File changed.

Preview size limit exceeded, changes collapsed.

catalog/C9_Catalog.py

+392 −191

File changed.

Preview size limit exceeded, changes collapsed.

catalog/Catalog_example.py

+21 −19

Original line number	Diff line number	Diff line
		@@ -4,7 +4,7 @@ import astropy.constants as cons
		from astropy.table import Table
		from scipy import interpolate

		from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar
		from observation_sim.mock_objects import CatalogBase, Star


		class Catalog(CatalogBase):
		@@ -54,14 +54,17 @@ class Catalog(CatalogBase):

		super().__init__()
		self.cat_dir = os.path.join(
		config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
		config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"]
		)
		self.chip = chip
		if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"]:
		if (
		"star_cat" in config["catalog_options"]["input_path"]
		and config["catalog_options"]["input_path"]["star_cat"]
		):
		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.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
		# NOTE: must call _load() method here to read in all objects
		self.objs = []
		self._load()
		@@ -142,17 +145,17 @@ class Catalog(CatalogBase):
		"""

		stars = Table.read(self.star_path)
		nstars = stars['sourceID'].size
		nstars = stars["sourceID"].size
		for istars in range(nstars):
		param = self.initialize_param()
		param['id'] = istars + 1
		param['ra'] = stars['RA'][istars]
		param['dec'] = stars['Dec'][istars]
		param['sed_type'] = stars['sourceID'][istars]
		param['model_tag'] = stars['model_tag'][istars]
		param['z'] = 0.0
		param['star'] = 1 # Star
		param['mag_use_normal'] = stars['app_sdss_g'][istars]
		param["id"] = istars + 1
		param["ra"] = stars["RA"][istars]
		param["dec"] = stars["Dec"][istars]
		param["sed_type"] = stars["sourceID"][istars]
		param["model_tag"] = stars["model_tag"][istars]
		param["z"] = 0.0
		param["star"] = 1 # Star
		param["mag_use_normal"] = stars["app_sdss_g"][istars]
		obj = Star(param)
		self.objs.append(obj)

		@@ -176,11 +179,10 @@ class Catalog(CatalogBase):
		NOTE: the range of wavelengthes must at least cover [2450 - 11000] Angstorms
		"""

		if obj.type == 'star':
		wave = Table.read(self.star_SED_path,
		path=f"/SED/wave_{obj.model_tag}")
		if obj.type == "star":
		wave = Table.read(self.star_SED_path, path=f"/SED/wave_{obj.model_tag}")
		flux = Table.read(self.star_SED_path, path=f"/SED/{obj.sed_type}")
		wave, flux = wave['col0'].data, flux['col0'].data
		wave, flux = wave["col0"].data, flux["col0"].data
		else:
		raise ValueError("Object type not known")
		speci = interpolate.interp1d(wave, flux)
		@@ -188,7 +190,7 @@ class Catalog(CatalogBase):
		y = speci(lamb)
		# erg/s/cm^2/A --> photons/s/m^2/A
		all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
		sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
		sed = Table(np.array([lamb, all_sed]).T, names=("WAVELENGTH", "FLUX"))
		return sed

		def load_norm_filt(self, obj):

catalog/init.py

0 → 100644

+0 −0

Empty file added.

catalog/testCat_fits.py

+68 −52

Original line number	Diff line number	Diff line
		@@ -5,23 +5,21 @@ 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 observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar, Stamp
		from observation_sim.mock_objects._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
		from observation_sim.mock_objects import CatalogBase, Stamp
		from observation_sim.mock_objects._util import (
		getObservedSED,
		)
		from observation_sim.astrometry.Astrometry_util import on_orbit_obs_position

		import astropy.io.fits as fitsio
		from observation_sim.mock_objects._util import seds, sed_assign, extAv
		from observation_sim.mock_objects._util import seds, sed_assign

		# (TEST)
		from astropy.cosmology import FlatLambdaCDM
		from astropy import constants
		from astropy import units as U

		try:
		import importlib.resources as pkg_resources
		@@ -45,24 +43,30 @@ class Catalog(CatalogBase):
		self.filt = filt
		self.logger = chip_output.logger

		with pkg_resources.path('catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
		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:
		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.
		self.max_size = 0.0

		if "stamp_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["stamp_cat"] and config["catalog_options"]["stamp_yes"]:
		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.stamp_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["stamp_SED"]) ###shoule be stamp-SED
		# self._load_SED_lib_stamps() ###shoule be stamp-SED
		self.tempSed_gal, self.tempRed_gal = seds(
		"galaxy.list", seddir="/public/home/chengliang/CSSOSDataProductsSims/testCats/Templates/Galaxy/") # only for test
		"galaxy.list",
		seddir="/public/home/chengliang/CSSOSDataProductsSims/testCats/Templates/Galaxy/",
		) # only for test

		self._add_output_columns_header()
		self._get_healpix_list()
		@@ -70,24 +74,31 @@ class Catalog(CatalogBase):

		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_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)
		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
		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
		NSIDE, hp.ang2vec(np.radians(90.0) - dec, ra), inclusive=True
		)
		# self.pix_list = hp.query_polygon(NSIDE, np.array(vertices).T, inclusive=True)
		if self.logger is not None:
		@@ -104,7 +115,7 @@ class Catalog(CatalogBase):

		def _load_stamps(self, stamps, pix_id=None):
		print("debug:: load_stamps")
		nstamps = len(stamps['filename'])
		nstamps = len(stamps["filename"])
		self.rng_sedGal = random.Random()
		# Use healpix index as the random seed
		self.rng_sedGal.seed(float(pix_id))
		@@ -113,31 +124,32 @@ class Catalog(CatalogBase):
		for istamp in range(nstamps):
		print("debug::", istamp)
		fitsfile = os.path.join(
		self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8'))
		self.cat_dir, "stampCats/" + stamps["filename"][istamp].decode("utf-8")
		)
		print("debug::", istamp, fitsfile)
		hdu = fitsio.open(fitsfile)

		param = self.initialize_param()
		param['id'] = hdu[0].header['index'] # istamp
		param['star'] = 3 # Stamp type in .cat file
		param['ra'] = hdu[0].header['ra']
		param['dec'] = hdu[0].header['dec']
		param['pixScale'] = hdu[0].header['pixScale']
		param["id"] = hdu[0].header["index"] # istamp
		param["star"] = 3 # Stamp type in .cat file
		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']
		# gals['mag_true_g_lsst']
		param['mag_use_normal'] = hdu[0].header['mag_g']
		param["mag_use_normal"] = hdu[0].header["mag_g"]

		# Apply astrometric modeling
		# in C3 case only aberration
		param['ra_orig'] = param['ra']
		param['dec_orig'] = param['dec']
		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()
		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()
		@@ -161,25 +173,26 @@ class Catalog(CatalogBase):
		input_vz=self.pointing.sat_vz,
		input_epoch="J2015.5",
		input_date_str=date_str,
		input_time_str=time_str
		input_time_str=time_str,
		)
		param['ra'] = ra_arr[0]
		param['dec'] = dec_arr[0]
		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['mu'] = 1
		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["mu"] = 1

		# 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']))
		param["image"] = hdu[0].data
		param["image"] = param["image"] / (np.sum(param["image"]))
		obj = Stamp(param)
		self.objs.append(obj)

		@@ -187,8 +200,12 @@ class Catalog(CatalogBase):
		self.objs = []
		self.ids = 0

		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']
		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"]
		print("debug::", stamps_cat.keys())

		for pix in self.pix_list:
		@@ -203,18 +220,17 @@ class Catalog(CatalogBase):

		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)))
		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 == 'stamp':
		if obj.type == "stamp":
		sed_data = getObservedSED(
		sedCat=self.tempSed_gal[obj.sed_type],
		redshift=obj.z,
		av=obj.param["av"],
		redden=obj.param["redden"]
		redden=obj.param["redden"],
		)
		wave, flux = sed_data[0], sed_data[1]
		else:
		@@ -224,7 +240,7 @@ class Catalog(CatalogBase):
		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'))
		sed = Table(np.array([lamb, all_sed]).T, names=("WAVELENGTH", "FLUX"))

		del wave
		del flux