Skip to content
GitLab
Commit fd6c3108 authored by Fang Yuedong's avatar Fang Yuedong Committed by Zhang Xin
Browse files

Release v2.0

parent ca7fb5db
Hide whitespace changes
Inline Side-by-side
Catalog/C3Catalog.py → Catalog/C6_Catalog.py
View file @ fd6c3108
......@@ -5,15 +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 ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
from ObservationSim.MockObject._util import seds, sed_assign, extAv, tag_sed, getObservedSED
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:
......@@ -22,52 +28,90 @@ except ImportError:
NSIDE = 128
class C3Catalog(CatalogBase):
def __init__(self, config, chip, pointing, **kwargs):
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["input_path"]["cat_dir"])
self.seed_Av = config["random_seeds"]["seed_Av"]
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.seed_Av = config["catalog_options"]["seed_Av"]
if "logger" in kwargs:
self.logger = kwargs["logger"]
else:
self.logger = None
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.normF_star = Table.read(str(filter_path))
self.config = config
self.chip = chip
self.pointing = pointing
if "star_cat" in config["input_path"] and config["input_path"]["star_cat"] and not config["run_option"]["galaxy_only"]:
star_file = config["input_path"]["star_cat"]
star_SED_file = config["SED_templates_path"]["star_SED"]
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["input_path"] and config["input_path"]["galaxy_cat"] and not config["run_option"]["star_only"]:
galaxy_file = config["input_path"]["galaxy_cat"]
self.galaxy_path = os.path.join(self.cat_dir, galaxy_file)
self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"])
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 "rotateEll" in config["shear_setting"]:
self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.))
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 = 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_ouptut_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)
# 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)
......@@ -78,7 +122,8 @@ class C3Catalog(CatalogBase):
if obj.type == "star":
return self.normF_star
elif obj.type == "galaxy" or obj.type == "quasar":
return self.normF_galaxy
# return self.normF_galaxy
return None
else:
return None
......@@ -86,18 +131,24 @@ class C3Catalog(CatalogBase):
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)
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):
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)
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_true'][:]
dec_arr = gals['dec_true'][:]
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()
......@@ -105,7 +156,7 @@ class C3Catalog(CatalogBase):
pmdec_list = np.zeros(ngals).tolist()
rv_list = np.zeros(ngals).tolist()
parallax_list = [1e-9] * ngals
dt = datetime.fromtimestamp(self.pointing.timestamp)
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
......@@ -122,84 +173,119 @@ class C3Catalog(CatalogBase):
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2015.5",
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_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:
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_true'][igals]
param['z'] = gals['redshift'][igals]
param['model_tag'] = 'None'
param['gamma1'] = 0
param['gamma2'] = 0
param['kappa'] = 0
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
# sersicB = gals['sersic_bulge'][igals]
hlrMajB = gals['size_bulge_true'][igals]
hlrMinB = gals['size_minor_bulge_true'][igals]
# sersicD = gals['sersic_disk'][igals]
hlrMajD = gals['size_disk_true'][igals]
hlrMinD = gals['size_minor_disk_true'][igals]
aGal = gals['size_true'][igals]
bGal = gals['size_minor_true'][igals]
param['bfrac'] = gals['bulge_to_total_ratio_i'][igals]
param['theta'] = gals['position_angle_true'][igals]
param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB)
param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD)
param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB)
param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD)
param['ell_tot'] = (aGal - bGal) / (aGal + bGal)
# Assign each galaxy a template SED
param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
param['redden'] = self.tempRed_gal[param['sed_type']]
param['av'] = self.avGal[int(self.ud()*self.nav)]
if param['sed_type'] <= 5:
param['av'] = 0.0
param['redden'] = 0
param['star'] = 0 # Galaxy
if param['sed_type'] >= 29:
param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
param['star'] = 2 # Quasar
# 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'] = gals['galaxyID'][igals]
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger)
self.objs.append(obj)
if param['star'] == 2:
obj = Quasar(param, logger=self.logger)
self.objs.append(obj)
# 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
# in C3 case only aberration
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 = np.zeros(nstars).tolist()
pmdec_list = np.zeros(nstars).tolist()
rv_list = np.zeros(nstars).tolist()
parallax_list = [1e-9] * nstars
dt = datetime.fromtimestamp(self.pointing.timestamp)
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(
......@@ -216,23 +302,30 @@ class C3Catalog(CatalogBase):
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2015.5",
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
# 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]
......@@ -242,31 +335,117 @@ class C3Catalog(CatalogBase):
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.nav = 15005
self.avGal = extAv(self.nav, seed=self.seed_Av)
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 not self.config["catalog_options"]["galaxy_only"]:
star_cat = h5.File(self.star_path, 'r')['catalog']
for pix in self.pix_list:
stars = star_cat[str(pix)]
self._load_stars(stars, pix_id=pix)
del stars
if "galaxy_cat" in self.config["input_path"] and self.config["input_path"]["galaxy_cat"] and not self.config["run_option"]["star_only"]:
gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
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:
gals = gals_cat[str(pix)]
self._load_gals(gals, pix_id=pix)
del gals
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))
del self.avGal
def load_sed(self, obj, **kwargs):
if obj.type == 'star':
......@@ -278,22 +457,58 @@ class C3Catalog(CatalogBase):
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]
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(2500, 10001 + 0.5, 0.5)
lamb = np.arange(2400, 11001 + 0.5, 0.5)
lamb = np.arange(2000, 11001+0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
# 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
Catalog/C6_fits_Catalog.py 0 → 100644
View file @ fd6c3108
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 "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 = 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_ouptut_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 > 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']
# 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, self.rotation, 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)
###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):
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'] = 22 #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
Catalog/Catalog_example.py
View file @ fd6c3108
......@@ -6,7 +6,7 @@ from scipy import interpolate
from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
class Catalog_example(CatalogBase):
class Catalog(CatalogBase):
"""An user customizable class for reading in catalog(s) of objects and SEDs.
NOTE: must inherit the "CatalogBase" abstract class
......@@ -32,8 +32,7 @@ class Catalog_example(CatalogBase):
load_norm_filt(obj):
load the filter throughput for the input catalog's photometric system.
"""
def __init__(self, config, chip, pointing, **kwargs):
def __init__(self, config, chip, **kwargs):
"""Constructor method.
Parameters
......@@ -41,7 +40,7 @@ class Catalog_example(CatalogBase):
config : dict
configuration dictionary which is parsed from the input YAML file
chip: ObservationSim.Instrument.Chip
a ObservationSim.Instrument.Chip instance, can be used to identify the band etc.
an ObservationSim.Instrument.Chip instance, can be used to identify the band etc.
**kwargs : dict
other needed input parameters (in key-value pairs), please modify corresponding
initialization call in "ObservationSim.py" as you need.
......@@ -52,12 +51,11 @@ class Catalog_example(CatalogBase):
"""
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.chip = chip
self.pointing = pointing
if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]:
star_file = config["input_path"]["star_cat"]
star_SED_file = config["SED_templates_path"]["star_SED"]
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)
# NOTE: must call _load() method here to read in all objects
......@@ -100,28 +98,34 @@ class Catalog_example(CatalogBase):
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["theta"] : float
the position angle (in degrees)
param["bfrac"] : float
the bulge fraction
param["hlr_bulge] : float
param["hlr_bulge"] : float
the half-light-radius of the bulge
param["hlr_disk] : float
param["hlr_disk"] : float
the half-light-radius of the disk
param["ell_bulge] : float
the ellipticity of the bulge
param["ell_disk] : float
the ellipticity of the disk
param["ell_tot] : float
the total ellipticity
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)
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
......
Catalog/FGS_Catalog.py 0 → 100644
View file @ fd6c3108
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 = 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_ouptut_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, self.rotation, 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
Catalog/NGPCatalog.py
View file @ fd6c3108
......@@ -22,11 +22,11 @@ except ImportError:
NSIDE = 128
class NGPCatalog(CatalogBase):
class Catalog(CatalogBase):
def __init__(self, config, chip, pointing, chip_output, **kwargs):
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
self.seed_Av = config["random_seeds"]["seed_Av"]
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
......@@ -40,19 +40,19 @@ class NGPCatalog(CatalogBase):
self.chip = chip
self.pointing = pointing
if "star_cat" in config["input_path"] and config["input_path"]["star_cat"] and not config["run_option"]["galaxy_only"]:
star_file = config["input_path"]["star_cat"]
star_SED_file = config["SED_templates_path"]["star_SED"]
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["input_path"] and config["input_path"]["galaxy_cat"] and not config["run_option"]["star_only"]:
galaxy_file = config["input_path"]["galaxy_cat"]
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["SED_templates_path"]["galaxy_SED"])
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["shear_setting"]:
self.rotation = float(int(config["shear_setting"]["rotateEll"]/45.))
if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
else:
self.rotation = 0.
......@@ -111,7 +111,7 @@ class NGPCatalog(CatalogBase):
pmdec_list = np.zeros(ngals).tolist()
rv_list = np.zeros(ngals).tolist()
parallax_list = [1e-9] * ngals
dt = datetime.fromtimestamp(self.pointing.timestamp)
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
......@@ -140,12 +140,12 @@ class NGPCatalog(CatalogBase):
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
# if param['mag_use_normal'] >= 26.5:
# continue
param['z'] = gals['redshift_true'][igals]
param['model_tag'] = 'None'
param['gamma1'] = 0
param['gamma2'] = 0
param['g1'] = 0
param['g2'] = 0
param['kappa'] = 0
param['delta_ra'] = 0
param['delta_dec'] = 0
......@@ -167,8 +167,13 @@ class NGPCatalog(CatalogBase):
# 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)]
# 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
......@@ -211,7 +216,7 @@ class NGPCatalog(CatalogBase):
pmdec_list = pmdec_arr.tolist()
rv_list = rv_arr.tolist()
parallax_list = parallax_arr.tolist()
dt = datetime.fromtimestamp(self.pointing.timestamp)
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
......@@ -245,8 +250,8 @@ class NGPCatalog(CatalogBase):
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
# if param['mag_use_normal'] >= 26.5:
# continue
self.ids += 1
# param['id'] = self.ids
param['id'] = stars['sourceID'][istars]
......@@ -269,7 +274,7 @@ class NGPCatalog(CatalogBase):
self.avGal = extAv(self.nav, seed=self.seed_Av)
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 not self.config["catalog_options"]["galaxy_only"]:
star_cat = h5.File(self.star_path, 'r')['catalog']
for pix in self.pix_list:
try:
......@@ -279,7 +284,7 @@ class NGPCatalog(CatalogBase):
except Exception as e:
self.logger.error(str(e))
print(e)
if "galaxy_cat" in self.config["input_path"] and self.config["input_path"]["galaxy_cat"] and not self.config["run_option"]["star_only"]:
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:
......@@ -316,8 +321,6 @@ class NGPCatalog(CatalogBase):
else:
raise ValueError("Object type not known")
speci = interpolate.interp1d(wave, flux)
# lamb = np.arange(2500, 10001 + 0.5, 0.5)
# lamb = np.arange(2400, 11001 + 0.5, 0.5)
lamb = np.arange(2000, 18001 + 0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photo/s/m2/A
......
Catalog/testCat_fits.py 0 → 100644
View file @ fd6c3108
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):
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 "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="/share/simudata/CSSOSDataProductsSims/data/Templates/Galaxy/") #only for test
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_ouptut_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 == "stamp":
#return self.normF_galaxy ###normalize_filter for stamp
return None
else:
return None
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):
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['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'] = hdu[0].header['mag_g'] #gals['mag_true_g_lsst']
# 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)
def _load(self, **kwargs):
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']
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)
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 == '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]
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
Catalog/testCat_galaxy.py 0 → 100644
View file @ fd6c3108
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 = 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_ouptut_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, self.rotation, 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
Catalog/testCat_star.py 0 → 100644
View file @ fd6c3108
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_ouptut_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
ObservationSim/Config/ChipOutput.py
View file @ fd6c3108
......@@ -18,43 +18,39 @@ class ChipOutput(object):
else:
self.ra_cen = config["obs_setting"]["ra_center"]
self.dec_cen = config["obs_setting"]["dec_center"]
exp_name = imgKey0 + "_%s_%s.fits"
self.chipLabel = focal_plane.getChipLabel(chip.chipID)
self.img_name = prefix + exp_name%(self.chipLabel, filt.filter_type)
# self.cat_name = 'MSC_' + config["obs_setting"]["date_obs"] + config["obs_setting"]["time_obs"] + "_" + str(pointing_ID).rjust(7, '0') + "_" + self.chipLabel.rjust(2,'0') + ".cat"
self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".cat"
self.chipLabel = focal_plane.getChipLabel(chip.chipID)
self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat"
self.subdir = subdir
# Setup logger for each chip
logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), focal_plane.getChipLabel(chip.chipID), filt.filter_type) + ".log"
logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".log"
self.logger = logging.getLogger()
fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
fh.setLevel(logging.DEBUG)
self.logger.setLevel(logging.DEBUG)
logging.getLogger('numba').setLevel(logging.WARNING)
# formatter = logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s')
formatter = logging.Formatter('%(asctime)s - %(msecs)d - %(levelname)-8s - [%(filename)s:%(lineno)d] - %(message)s')
fh.setFormatter(formatter)
self.logger.addHandler(fh)
hdr1 = "# obj_ID ID_chip filter xImage yImage ra dec ra_orig dec_orig z mag obj_type "
hdr2 = "thetaR bfrac hlr_disk hlr_bulge e1_disk e2_disk e1_bulge e2_bulge g1 g2 "
hdr3 = "sed_type av redden "
hdr4 = "pm_ra pm_dec RV parallax"
fmt1 = "%10d %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s "
fmt2 = "%8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f "
fmt3 = "%2d %8.4f %8.4f "
fmt4 = "%15.8f %15.8f %15.8f %15.8f"
self.hdr = hdr1 + hdr2 + hdr3 + hdr4
self.fmt = fmt1 + fmt2 + fmt3 + fmt4
hdr2 = "pm_ra pm_dec RV parallax"
fmt1 = "%20s %4d %5s %10.3f %10.3f %15.8f %15.8f %15.8f %15.8f %7.4f %8.4f %15s "
fmt2 = "%15.8f %15.8f %15.8f %15.8f"
self.hdr = hdr1 + hdr2
self.fmt = fmt1 + fmt2
self.logger.info("pointing_type = %s\n"%(pointing_type))
def Log_info(self, message):
print(message)
self.logger.info(message)
def Log_error(self, message):
print(message)
self.logger.error(message)
def update_ouptut_header(self, additional_column_names=""):
self.hdr += additional_column_names
......@@ -70,16 +66,14 @@ class ChipOutput(object):
def cat_add_obj(self, obj, pos_img, pos_shear):
# ximg = pos_img.x - self.chip.bound.xmin + 1.0
# yimg = pos_img.y - self.chip.bound.ymin + 1.0
# print('-------------debug-----------------')
# print('from global',ximg, yimg)
# self.logger.info('-------------debug-----------------')
# self.logger.info('from global',ximg, yimg)
ximg = obj.real_pos.x + 1.0
yimg = obj.real_pos.y + 1.0
# print('from loacl',ximg, yimg)
# self.logger.info('from loacl',ximg, yimg)
line = self.fmt%(
obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type,
obj.thetaR, obj.bfrac, obj.hlr_disk, obj.hlr_bulge, obj.e1_disk, obj.e2_disk, obj.e1_bulge, obj.e2_bulge, obj.g1, obj.g2,
obj.sed_type, obj.av, obj.redden,
obj.pmra, obj.pmdec, obj.rv, obj.parallax)
line += obj.additional_output_str
if not line.endswith("\n"):
......
ObservationSim/Config/Config.py
View file @ fd6c3108
......@@ -17,14 +17,9 @@ def config_dir(config, work_dir=None, data_dir=None):
path_dict["data_dir"] =os.path.join(path_dict["work_dir"], "data/")
else:
path_dict["data_dir"] = data_dir
# Data sub-catalogs
# Object catalog direcotry
path_dict["cat_dir"] = os.path.join(path_dict["data_dir"], config["input_path"]["cat_dir"])
# PSF data directory
if config["psf_setting"]["psf_model"] == "Interp":
path_dict["psf_dir"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["psf_dir"])
if config["ins_effects"]["field_dist"] == True:
path_dict["fd_path"] = os.path.join(path_dict["data_dir"], config["psf_setting"]["fd_path"])
return path_dict
......
ObservationSim/Config/Header/ImageHeader.py
View file @ fd6c3108
......@@ -13,7 +13,10 @@ import random
import os
import sys
import astropy.coordinates as coord
from astropy.coordinates import SkyCoord
from astropy.wcs.utils import fit_wcs_from_points
from astropy.time import Time
from astropy import wcs
def chara2digit(char):
""" Function to judge and convert characters to digitals
......@@ -164,7 +167,8 @@ def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232
##9232 9216 898 534 1309 60 -40 -23.4333
def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'):
def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra_ref = 60, dec_ref = -40, pa = -23.433, pixel_scale = 0.074, pixel_size=1e-2,
rotate_chip=0., filter = 'GI', row_num = None, col_num = None, xcen = None, ycen = None):
""" Creat a wcs frame for CCST with multiple extensions
......@@ -172,141 +176,172 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
----------
"""
r_dat = OrderedDict()
r_dat['EQUINOX'] = 2000.0
r_dat['WCSDIM'] = 2.0
r_dat['CTYPE1'] = 'RA---TAN'
r_dat['CTYPE2'] = 'DEC--TAN'
r_dat['CRVAL1'] = ra_ref
r_dat['CRVAL2'] = dec_ref
flag_x = [0, 1, -1, 1, -1]
flag_y = [0, 1, 1, -1, -1]
flag_ext_x = [0,-1,1,-1,1]
flag_ext_y = [0,-1,-1,1,1]
x_num = 6
y_num = 5
detector_num = x_num*y_num
detector_size_x = xlen
detector_size_y = ylen
gap_y = gapy
gap_x = [gapx1,gapx2]
ra_ref = ra
dec_ref = dec
pa_aper = pa
pixel_size = psize
gap_x1_num = 3
gap_x2_num = 2
y_center = (detector_size_y*y_num+gap_y*(y_num-1))/2
if (row_num is not None) and (col_num is not None):
x_center = (detector_size_x*x_num+gap_x[0]*gap_x1_num+gap_x[1]*gap_x2_num)/2
x_num = 6
y_num = 5
detector_num = x_num*y_num
gap_x_map = np.array([[0,0,0,0,0],[gap_x[0],gap_x[1],gap_x[1],gap_x[1],gap_x[1]],[gap_x[1],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[1]],[gap_x[1],gap_x[1],gap_x[1],gap_x[1],gap_x[0]]])
detector_size_x = xlen
detector_size_y = ylen
gap_y = gapy
gap_x = [gapx1,gapx2]
gap_x1_num = 3
gap_x2_num = 2
# frame_array = np.empty((5,6),dtype=np.float64)
# print(x_center,y_center)
j = row_num
i = col_num
# ccdnum = str((j-1)*5+i)
x_ref, y_ref = detector_size_x*i + sum(gap_x_map[0:i,j-1]) - detector_size_x/2. , (detector_size_y+gap_y)*j-gap_y-detector_size_y/2
y_center = (detector_size_y*y_num+gap_y*(y_num-1))/2
# print(i,j,x_ref,y_ref,ra_ref,dec_ref)
x_center = (detector_size_x*x_num+gap_x[0]*gap_x1_num+gap_x[1]*gap_x2_num)/2
r_dat = OrderedDict()
gap_x_map = np.array([[0,0,0,0,0],[gap_x[0],gap_x[1],gap_x[1],gap_x[1],gap_x[1]],[gap_x[1],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[0]],[gap_x[0],gap_x[0],gap_x[0],gap_x[0],gap_x[1]],[gap_x[1],gap_x[1],gap_x[1],gap_x[1],gap_x[0]]])
j = row_num
i = col_num
# ccdnum = str((j-1)*5+i)
# name = []
# value = []
# description = []
for k in range(1,2):
x_ref, y_ref = detector_size_x*i + sum(gap_x_map[0:i,j-1]) - detector_size_x/2. , (detector_size_y+gap_y)*j-gap_y-detector_size_y/2
cd = np.array([[ pixel_size, 0], [0, pixel_size]])/3600.*flag_x[k]
for k in range(1,2):
cd = np.array([[ pixel_scale, 0], [0, pixel_scale]])/3600.*flag_x[k]
cd_rot = rotate_CD_matrix(cd, pa_aper)
# f = open("CCD"+ccdnum.rjust(2,'0')+"_extension"+str(k)+"_wcs.param","w")
r_dat['CRPIX1'] = flag_ext_x[k]*((x_ref+flag_ext_x[k]*detector_size_x/2)-x_center)
r_dat['CRPIX2'] = flag_ext_y[k]*((y_ref+flag_ext_y[k]*detector_size_y/2)-y_center)
r_dat['CD1_1'] = cd_rot[0,0]
r_dat['CD1_2'] = cd_rot[0,1]
r_dat['CD2_1'] = cd_rot[1,0]
r_dat['CD2_2'] = cd_rot[1,1]
if filter in ['GU', 'GV', 'GI']:
w = wcs.WCS(naxis=2)
w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
w.wcs.cd = cd_rot
w.wcs.crval = [ra_ref, dec_ref]
w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']]
# test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
sls_rot = 1
if i > 2:
sls_rot = -sls_rot
sn_x = 30
sn_y = 30
x_pixs = np.zeros(sn_y * sn_x)
y_pixs = np.zeros(sn_y * sn_x)
xpixs_line = np.linspace(1, xlen, sn_x)
ypixs_line = np.linspace(1, ylen, sn_y)
sky_coors = []
for n1, y in enumerate(ypixs_line):
for n2, x in enumerate(xpixs_line):
i_pix = n1 * sn_x + n2
x_pixs[i_pix] = x
y_pixs[i_pix] = y
pix_coor = np.array([x, y])
sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, xlen=xlen, ylen=ylen, w=w)
# print(sc1[0,0],sc1[0,1])
sky_coors.append((sc1[0, 0], sc1[0, 1]))
wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]),
world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN')
# print(wcs_new)
# test_center = wcs_new.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
#
# print(test_center - test_center_o)
r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0]
r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1]
r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0]
r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1]
r_dat['CRPIX1'] = wcs_new.wcs.crpix[0]
r_dat['CRPIX2'] = wcs_new.wcs.crpix[1]
r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
r_dat['CRVAL2'] = wcs_new.wcs.crval[1]
elif (xcen is not None) and (ycen is not None):
xcen, ycen = xcen/pixel_size, ycen/pixel_size
x1, y1 = xcen - xlen/2., ycen - ylen/2.
r_dat['CRPIX1'] = -x1
r_dat['CRPIX2'] = -y1
cd = np.array([[ pixel_scale, 0], [0, pixel_scale]])/3600.*flag_x[1]
cd_rot = rotate_CD_matrix(cd, pa_aper)
# f = open("CCD"+ccdnum.rjust(2,'0')+"_extension"+str(k)+"_wcs.param","w")
r_dat['EQUINOX'] = 2000.0
r_dat['WCSDIM'] = 2.0
r_dat['CTYPE1'] = 'RA---TAN'
r_dat['CTYPE2'] = 'DEC--TAN'
r_dat['CRVAL1'] = ra_ref
r_dat['CRVAL2'] = dec_ref
r_dat['CRPIX1'] = flag_ext_x[k]*((x_ref+flag_ext_x[k]*detector_size_x/2)-x_center)
r_dat['CRPIX2'] = flag_ext_y[k]*((y_ref+flag_ext_y[k]*detector_size_y/2)-y_center)
r_dat['CD1_1'] = cd_rot[0,0]
r_dat['CD1_2'] = cd_rot[0,1]
r_dat['CD2_1'] = cd_rot[1,0]
r_dat['CD2_2'] = cd_rot[1,1]
if filter in ['GU', 'GV', 'GI']:
from astropy import wcs
w = wcs.WCS(naxis=2)
w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
w.wcs.cd = cd_rot
w.wcs.crval = [ra_ref, dec_ref]
w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']]
# test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
sls_rot = 1
if i > 2:
sls_rot = -sls_rot
sn_x = 30
sn_y = 30
x_pixs = np.zeros(sn_y * sn_x)
y_pixs = np.zeros(sn_y * sn_x)
xpixs_line = np.linspace(1, xlen, sn_x)
ypixs_line = np.linspace(1, ylen, sn_y)
sky_coors = []
for n1, y in enumerate(ypixs_line):
for n2, x in enumerate(xpixs_line):
i_pix = n1 * sn_x + n2
x_pixs[i_pix] = x
y_pixs[i_pix] = y
pix_coor = np.array([x, y])
sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, w=w)
# print(sc1[0,0],sc1[0,1])
sky_coors.append((sc1[0, 0], sc1[0, 1]))
from astropy.coordinates import SkyCoord
from astropy.wcs.utils import fit_wcs_from_points
wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]),
world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN')
# print(wcs_new)
# test_center = wcs_new.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
#
# print(test_center - test_center_o)
r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0]
r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1]
r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0]
r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1]
r_dat['CRPIX1'] = wcs_new.wcs.crpix[0]
r_dat['CRPIX2'] = wcs_new.wcs.crpix[1]
r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
r_dat['CRVAL2'] = wcs_new.wcs.crval[1]
return r_dat
w = wcs.WCS(naxis=2)
w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
w.wcs.cd = cd_rot
w.wcs.crval = [ra_ref, dec_ref]
w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']]
sn_x = 30
sn_y = 30
x_pixs = np.zeros(sn_y * sn_x)
y_pixs = np.zeros(sn_y * sn_x)
xpixs_line = np.linspace(1, xlen, sn_x)
ypixs_line = np.linspace(1, ylen, sn_y)
sky_coors = []
for n1, y in enumerate(ypixs_line):
for n2, x in enumerate(xpixs_line):
i_pix = n1 * sn_x + n2
x_pixs[i_pix] = x
y_pixs[i_pix] = y
pix_coor = np.array([x, y])
sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=rotate_chip, xlen=xlen, ylen=ylen, w=w)
sky_coors.append((sc1[0, 0], sc1[0, 1]))
wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]),
world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"),
projection='TAN')
r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0]
r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1]
r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0]
r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1]
r_dat['CRPIX1'] = wcs_new.wcs.crpix[0]
r_dat['CRPIX2'] = wcs_new.wcs.crpix[1]
r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
r_dat['CRVAL2'] = wcs_new.wcs.crval[1]
else:
raise ValueError('In function WCS_def(): Either (row_num, col_num) or (xcen, ycen, pixel_size) should be given')
return r_dat
def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, psize = 0.074, row_num = 1, col_num = 1, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.]):
#TODO project_cycle is temporary， is not in header defined, delete in future
def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, chip_name="01"):
# array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0
k = (row_num-1)*6+col_num
# k = (row_num-1)*6+col_num
# ccdnum = str(k)
g_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/global_header.header'
......@@ -325,48 +360,64 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
h_prim = fits.Header()
h_prim = fits.Header.fromfile(g_header_fn)
h_prim['PIXSIZE1'] = xlen
h_prim['PIXSIZE2'] = ylen
h_prim['DATE'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6]
h_prim['TIME'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
h_prim['DATE-OBS'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6]
h_prim['TIME-OBS'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
h_prim['DETECTOR'] = 'CCD'+CCDID[k-1].rjust(2,'0')
h_prim['RA_OBJ'] = ra
h_prim['DEC_OBJ'] = dec
h_prim['OBJECT'] = '1'+ pointNum.rjust(8,'0')
h_prim['OBSID'] = '1'+ pointNum.rjust(8,'0')
h_prim['TELFOCUS'] = 'f/14'
# h_prim['PIXSIZE1'] = xlen
# h_prim['PIXSIZE2'] = ylen
h_prim['DATE'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6] + 'T' + time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
# h_prim['TIME'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
h_prim['DATE-OBS'] = '20'+date[0:2]+'-' + date[2:4]+'-'+date[4:6] + 'T' + time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
# h_prim['TIME-OBS'] = time_obs[0:2]+':'+time_obs[2:4]+':'+time_obs[4:6]
# h_prim['DETECTOR'] = 'CHIP'+CCDID[k-1].rjust(2,'0')
h_prim['OBJ_RA'] = ra
h_prim['OBJ_DEC'] = dec
obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}
OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')
h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7,'0')
h_prim['OBSID'] = OBS_id
# h_prim['TELFOCUS'] = 'f/14'
h_prim['EXPTIME'] = exptime
# Define file types
file_type = {'SCI':'sci', 'BIAS':'zero', 'DARK':'dark', 'FLAT':'flat', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray'}
file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray','CALS':'CALS','CALF':'CALF'}
h_prim['FILETYPE'] = file_type[im_type]
co = coord.SkyCoord(ra, dec, unit='deg')
ra_hms = format(co.ra.hms.h, '02.0f') + ':' + format(co.ra.hms.m, '02.0f') + ':' + format(co.ra.hms.s, '05.2f')
dec_hms = format(co.dec.dms.d, '02.0f') + ':' + format(abs(co.dec.dms.m), '02.0f') + ':' + format(abs(co.dec.dms.s),
'05.2f')
ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m, '02.0f') + format(co.ra.hms.s, '02.0f')
dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m), '02.0f') + format(abs(co.dec.dms.s), '02.0f')
h_prim['TARGET'] = ra_hms + '+' + dec_hms
#
# h_prim['RA_NOM'] = ra_hms
# h_prim['DEC_NOM'] = dec_hms
h_prim['RA_NOM'] = ra_hms
h_prim['DEC_NOM'] = dec_hms
h_prim['PIXSCAL1'] = psize
h_prim['PIXSCAL2'] = psize
h_prim['RA_PNT0'] = ra
h_prim['DEC_PNT0'] = dec
h_prim['RA_PNT1'] = ra
h_prim['DEC_PNT1'] = dec
ttt = h_prim['DATE'] + 'T' + h_prim['TIME']
# h_prim['PIXSCAL1'] = pixel_scale
# h_prim['PIXSCAL2'] = pixel_scale
ttt = h_prim['DATE']
tstart = Time(ttt)
h_prim['EXPSTART'] = round(tstart.mjd, 5)
h_prim['CABSTART'] = h_prim['EXPSTART']
# tend = Time(tstart.cxcsec + h_prim['EXPTIME'], format="cxcsec")
tend = Time(tstart.mjd + h_prim['EXPTIME']/86400., format="mjd")
h_prim['EXPEND'] = round(tend.mjd, 5)
h_prim['CABEND'] = h_prim['EXPEND']
file_start_time = '20' + date[0:6] + time_obs[0:6]
end_time_str = str(tend.datetime)
file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_1' + pointNum.rjust(8, '0') + '_' + CCDID[
k - 1].rjust(2, '0') + '_L0_1'
# h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
# k - 1].rjust(2, '0') + '_L0_V01'
h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
h_prim['POSI0_X'] = sat_pos[0]
......@@ -376,17 +427,18 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
h_prim['VELO0_X'] = sat_vel[0]
h_prim['VELO0_Y'] = sat_vel[1]
h_prim['VELO0_Z'] = sat_vel[2]
h_prim['RA_PNT0'] = ra_hms
h_prim['DEC_PNT0'] = dec_hms
# h_prim['RA_PNT0'] = ra_hms
# h_prim['DEC_PNT0'] = dec_hms
# Get version of CSSTSim Package
from pkg_resources import get_distribution
h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
# h_prim['SIM_VER'] = (get_distribution("CSSTSim").version, "Version of CSST MSC simulation software")
h_prim['FITSCREA'] = get_distribution("CSSTSim").version
return h_prim
def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1, extName='SCI'):
def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, pixel_scale = 0.074, pixel_size=1e-2,
extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None):
e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header'
f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')
......@@ -401,21 +453,58 @@ def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -2
s = s.strip("\n")
CCDID = s.split()
k = (row_num - 1) * 6 + col_num
# k = (row_num - 1) * 6 + col_num
h_ext = fits.Header.fromfile(e_header_fn)
h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
h_ext['FILTER'] = filters[k-1]
# h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
# h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
# h_ext['FILTER'] = filters[k-1]
h_ext['CCDCHIP'] = chip.chip_name
h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0')
h_ext['FILTER'] = chip.filter_type
h_ext['NAXIS1'] = xlen
h_ext['NAXIS2'] = ylen
h_ext['EXTNAME'] = extName
h_ext['GAIN1'] = gain
h_ext['RDNOISE1'] = readout
h_ext['CCDCHIP'] = 'ccd' + CCDID[k-1].rjust(2,'0')
h_ext['POS_ANG'] = pa
header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize,
row_num=row_num, col_num=col_num, filter = h_ext['FILTER'])
h_ext['GAIN2'] = gain
h_ext['GAIN3'] = gain
h_ext['GAIN4'] = gain
h_ext['GAIN5'] = gain
h_ext['GAIN6'] = gain
h_ext['GAIN7'] = gain
h_ext['GAIN8'] = gain
h_ext['GAIN9'] = gain
h_ext['GAIN10'] = gain
h_ext['GAIN11'] = gain
h_ext['GAIN12'] = gain
h_ext['GAIN13'] = gain
h_ext['GAIN14'] = gain
h_ext['GAIN15'] = gain
h_ext['GAIN16'] = gain
h_ext['RDNOIS1'] = readout
h_ext['RDNOIS2'] = readout
h_ext['RDNOIS3'] = readout
h_ext['RDNOIS4'] = readout
h_ext['RDNOIS5'] = readout
h_ext['RDNOIS6'] = readout
h_ext['RDNOIS7'] = readout
h_ext['RDNOIS8'] = readout
h_ext['RDNOIS9'] = readout
h_ext['RDNOIS10'] = readout
h_ext['RDNOIS11'] = readout
h_ext['RDNOIS12'] = readout
h_ext['RDNOIS13'] = readout
h_ext['RDNOIS14'] = readout
h_ext['RDNOIS15'] = readout
h_ext['RDNOIS16'] = readout
h_ext['PIXSCAL1'] = pixel_scale
h_ext['PIXSCAL2'] = pixel_scale
# h_ext['POS_ANG'] = pa
header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra_ref=ra, dec_ref=dec, pa=pa, pixel_scale=pixel_scale, pixel_size=pixel_size,
rotate_chip=chip.rotate_angle, filter = h_ext['FILTER'], row_num=row_num, col_num=col_num, xcen = xcen, ycen = ycen)
h_ext['CRPIX1'] = header_wcs['CRPIX1']
h_ext['CRPIX2'] = header_wcs['CRPIX2']
......
ObservationSim/Config/Header/extension_header.header
View file @ fd6c3108
XTENSION= 'IMAGE ' / Image extension BITPIX = 16 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 9216 / length of first array axis NAXIS2 = 9232 / length of second array axis PCOUNT = 0 GCOUNT = 1 EXTNAME = 'SCI ' EXTVER = 1 BSCALE = 1 BZERO = 32768 BUNIT = 'ADU ' COMMENT ================================================================== COMMENT CCD chip information COMMENT ================================================================== NCHAN = 16 / Number of readout channels NCHAN1 = 8 / Number of horizintal channels NCHAN2 = 2 / Number of vertical channels DETSIZE = '[1:9216:9232]' / detector size BIASSEC = '[9216:9217,9232:9234]' / bias section CCDCHIP = 6 / CCD chip ID CCDLABEL= 'y-1 ' / CCD chip label PSCAN1 = 27 / horizontal prescan width PSCAN2 = 8 / vertical prescan height OSCAN1 = 16 / horizontal overscan width OSCAN2 = 8 / vertical overscan height FILTER = 'y ' / filter name COMMENT ================================================================== COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS COMMENT ================================================================== WCSDIM = 2 / Number of World Coordinate System axes EQUINOX = 2000.0 / Epoch (year) CRPIX1 = -10017.0 / Coordinate reference pixel of x CRPIX2 = 24876.0 / Coordinate reference pixel of y CRVAL1 = 62.228226 / Coordinate reference value of x CRVAL2 = -42.316932 / Coordinate reference value of y CTYPE1 = 'RA---TAN' / the coordinate type CTYPE2 = 'DEC--TAN' / the coordinate type CD1_1 = 1.88602083707394E-05 / partial of first axis coordinate of x CD2_1 = -8.1745583617600E-06 / partial of first axis coordinate of y CD1_2 = 8.17455836176000E-06 / partial of second axis coordinate of x CD2_2 = 1.88602083707394E-05 / partial of second axis coordinate of y POS_ANG = 30.0 / Positon angle of detector array COMMENT ================================================================== COMMENT Readout information COMMENT ================================================================== GAIN1 = 1.1 / CCD gain RDNOISE1= 5.0 / read noise READTIME= 40.0 / read time (sec) RDSPEED = 10.0 / read speed (in MHz) CHIPTEMP= -100.0 / chip temperature (in K) COMMENT ================================================================== COMMENT Other information COMMENT ================================================================== CHECKSUM= '''abcde''' / SHA256 checksum of global headers END
XTENSION= 'IMAGE ' / extension type BITPIX = 16 / bits per data value NAXIS = 2 / number of data axes NAXIS1 = 9216 / length of first array axis NAXIS2 = 9232 / length of second array axis PCOUNT = 0 GCOUNT = 1 EXTNAME = 'SCI ' EXTVER = 1 BSCALE = 1 BZERO = 32768 BUNIT = 'ADU ' / physical unit of array values COMMENT ================================================================== COMMENT Detector information COMMENT ================================================================== DETECTOR= 'CCD' / detector name DETTEMP1= 173 / detector temperature at EXPSTART DETTEMP2= 173 / detector temperature at EXPEND DETTEMP3= 173 / detector temperature when readout is finished DETSIZE = '9560x9264' / detector size DATASEC = '9216x9232' / data size PIXSCAL1= 0.074 / pixel scale for axis 1 PIXSCAL2= 0.074 / pixel scale for axis 2 PIXSIZE1= 10 / pixel size for axis 1 PIXSIZE2= 10 / pixel size for axis 1 COMMENT ================================================================== COMMENT CCD chip information COMMENT ================================================================== CCDCHIP = 8 / CCD chip ID CCDLABEL= 'y-1' / CCD chip label FILTER = 'y' / filter name NCHAN = 16 / number of readout channels NCHAN1 = 8 / number of horizintal channels NCHAN2 = 2 / number of vertical channels PSCAN1 = 27 / horizontal prescan width, per readout channel PSCAN2 = 8 / horizontal prescan width, per readout channel OSCAN1 = 16 / horizontal overscan width,per readout channel OSCAN2 = 16 / horizontal overscan width,per readout channel COMMENT ================================================================== COMMENT WORLD COORDINATE SYSTEM AND RELATED PARAMETERS COMMENT ================================================================== WCSAXES = 2 / number of World Coordinate System axes CRPIX1 = -10017.0 / x-coordinate of reference pixel CRPIX2 = 24876.0 / y-coordinate of reference pixel CRVAL1 = 62.228226 / first axis value at reference pixel CRVAL2 = -42.316932 / second axis value at reference pixel CTYPE1 = 'RA---TAN' / the coordinate type for the first axis CTYPE2 = 'DEC--TAN' / the coordinate type for the second axis CD1_1 = 1.88602083707394E-05 / partial of first axis coordinate w.r.t. x CD2_1 = -8.1745583617600E-06 / partial of first axis coordinate w.r.t. x CD1_2 = 8.17455836176000E-06 / partial of second axis coordinate w.r.t. x CD2_2 = 1.88602083707394E-05 / partial of second axis coordinate w.r.t. x OTHERS = '' / COMMENT ================================================================== COMMENT Readout information COMMENT ================================================================== GAIN1 = 1.1 / CCD gain (channel 1) GAIN2 = 1.1 / CCD gain (channel 2) GAIN3 = 1.1 / CCD gain (channel 3) GAIN4 = 1.1 / CCD gain (channel 4) GAIN5 = 1.1 / CCD gain (channel 5) GAIN6 = 1.1 / CCD gain (channel 6) GAIN7 = 1.1 / CCD gain (channel 7) GAIN8 = 1.1 / CCD gain (channel 8) GAIN9 = 1.1 / CCD gain (channel 9) GAIN10 = 1.1 / CCD gain (channel 10) GAIN11 = 1.1 / CCD gain (channel 11) GAIN12 = 1.1 / CCD gain (channel 12) GAIN13 = 1.1 / CCD gain (channel 13) GAIN14 = 1.1 / CCD gain (channel 14) GAIN15 = 1.1 / CCD gain (channel 15) GAIN16 = 1.1 / CCD gain (channel 16) RDNOIS1 = 5.0 / read noise (channel 1) RDNOIS2 = 5.0 / read noise (channel 2) RDNOIS3 = 5.0 / read noise (channel 3) RDNOIS4 = 5.0 / read noise (channel 4) RDNOIS5 = 5.0 / read noise (channel 5) RDNOIS6 = 5.0 / read noise (channel 6) RDNOIS7 = 5.0 / read noise (channel 7) RDNOIS8 = 5.0 / read noise (channel 8) RDNOIS9 = 5.0 / read noise (channel 9) RDNOIS10= 5.0 / read noise (channel 10) RDNOIS11= 5.0 / read noise (channel 11) RDNOIS12= 5.0 / read noise (channel 12) RDNOIS13= 5.0 / read noise (channel 13) RDNOIS14= 5.0 / read noise (channel 14) RDNOIS15= 5.0 / read noise (channel 15) RDNOIS16= 5.0 / read noise (channel 16) RDSPEED = 10.0 / read speed (in MHz) COMMENT ================================================================== COMMENT Shutter information COMMENT ================================================================== SHUTSTAT= T / shutter status SHTOPEN0= 0.0 / shutter open time (begin) SHTOPEN1= 0.0 / shutter open time (end) SHTCLOS0= 0.0 / shutter close time (begin) SHTCLOS1= 0.0 / shutter close time (end) COMMENT ================================================================== COMMENT LED information COMMENT ================================================================== LEDSTAT1= '0000000000000000000000000000' / main LED status LEDEXP1 = 0.0 / main LED flash time (s) LEDSTAT2= '0000000000000000000000000000' / backup LED status LEDEXP2 = 0.0 / backup LED flash time (s) LEDTEMPA= 173.0 / LED temperature (main LED in K) LEDTEMPB= 173.0 / LED temperature (backup LED in K) COMMENT ================================================================== COMMENT Other information COMMENT ================================================================== CHECKSUM= '''abcde''' / SHA256 checksum of global headers END
ObservationSim/Config/Header/global_header.header
View file @ fd6c3108
SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T NEXTEND = 1 GROUPS = F DATE = '2021-03-04' / Date this file was written TIME = '09:30:00' / Time this file was written FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / Name of file FILETYPE= 'sci ' / type of data found in data data file TELESCOP= 'CSST ' / telescope used to acquire data INSTRUME= 'MSC ' / identifier for instrument used to acquire data RADECSYS= 'ICRS ' / frame of reference of coordinates EQUINOX = 2000.0 / Default equinox COMMENT ================================================================== COMMENT Object information COMMENT ================================================================== OBJECT = '00000000' / name of the object observed TARGET = '+000000000000' / Name of the target (RA & Dec of the object) OBSID = '00000000' / ID of this observation RA_OBJ = 62.228226 / R.A. of the object (degrees) DEC_OBJ = -42.316932 / declination of the object (degrees) COMMENT ================================================================== COMMENT Exposure information COMMENT ================================================================== SUNANGLE= 30.0 / angle between sun and direction perpendicular MOONANGL= 30.0 / angle between moon and direction perpendicular SUN_ALT = 30.0 / Sun altitude REFFRAME= 'CSSTGSC-1.0' / guide star catalog version DATE-OBS= '2021-03-04' / UTC date of start of observation (yyyy-mm-dd) TIME-OBS= '09:30:00' / UTC time of start of observation (hh:mm:ss) EXPSTART= 58849.0 / exposure start time (MJD) EXPEND = 58849.0 / exposure end time (MJD) EXPTIME = 150.0 / exposure duration EPOCH = 2021.0 / coordinate epoch COMMENT ================================================================== COMMENT Detector information COMMENT ================================================================== DETECTOR= 'CCD06 ' / deterctor name PIXSCAL1= 0.074 / Pixel scale for axis 1 PIXSCAL2= 0.074 / Pixel scale for axis 2 PIXSIZE1= 9216.0 / Pixel size for axis 1 PIXSIZE2= 9232.0 / Pixel size for aixs 2 CCDTEMP = -100.0 / CCD chip / dewar temperature (in K) COMMENT ================================================================== COMMENT Telescope information COMMENT ================================================================== TELSTAT = T / Tracking status TELTEMP1= -100.0 / Telescope temperature (in K) POSI0_X = 0.0 / The orbital position in X (shuttr open) POSI0_Y = 0.0 / The orbital position in Y (shuttr open) POSI0_Z = 0.0 / The orbital position in Z (shuttr open) VELO0_X = 0.0 / The orbital velocity in X (shuttr open) VELO0_Y = 0.0 / The orbital velocity in Y (shuttr open) VELO0_Z = 0.0 / The orbital velocity in Z (shuttr open) RA_PNT0 = '00:00:00.00' / R.A. at shutter open time DEC_PNT0= '00:00:00.00' / declination at shutter open time POSI1_X = 0.0 / The orbital position in X (shuttr close) POSI1_Y = 0.0 / The orbital position in Y (shuttr close) POSI1_Z = 0.0 / The orbital position in Z (shuttr close) VELO1_X = 0.0 / The orbital velocity in X (shuttr close) VELO1_Y = 0.0 / The orbital velocity in Y (shuttr close) VELO1_Z = 0.0 / The orbital velocity in Z (shuttr close) RA_PNT1 = 0.0 / R.A. at shutter close time DEC_PNT1= 0.0 / declination at shutter close time RA_NOM = '00:00:00.00' / nominal proposed R.A. of the pointing DEC_NOM = '00:00:00.00' / nominal proposed DEC of the pointing TELFOCUS= 'abcdefg ' / telescope focus AOSTAT = T / Active opticas status HOODSTAT= T / Lens hood status COMMENT ================================================================== COMMENT Shutter information COMMENT ================================================================== SHUTHWV = 'shutter-h-1.0' / shutter hardware version SHUTSWV = 'shutter-s-1.0' / shutter software version SHUTSTAT= T / shutter status SHTOPEN0= 0.0 / shutter open time (begin) SHTOPEN1= 0.0 / shutter open time (end) SHTCLOS0= 0.0 / shutter close time (begin) SHTCLOS1= 0.0 / shutter close time (end) COMMENT ================================================================== COMMENT Guide information COMMENT ================================================================== STARNO1 = 1.0 / guide star number STARPOX1= 0.0 / guide star average position STARPOY1= 0.0 / guide star average position SECMOM1 = 0.0 / guide star second momentum GAFLAG1 = T / guide flag indicate whether adjust position STARNO2 = 1.0 / guide star number STARPOX2= 0.0 / guide star average position STARPOY2= 0.0 / guide star average position SECMOM2 = 0.0 / guide star second momentum GAFLAG2 = T / guide flag indicate whether adjust position STARNO3 = 1.0 / guide star number STARPOX3= 0.0 / guide star average position STARPOY3= 0.0 / guide star average position SECMOM3 = 0.0 / guide star second momentum GAFLAG3 = T / guide flag indicate whether adjust position STARNO4 = 1.0 / guide star number STARPOX4= 0.0 / guide star average position STARPOY4= 0.0 / guide star average position SECMOM4 = 0.0 / guide star second momentum GAFLAG4 = T / guide flag indicate whether adjust position COMMENT ================================================================== COMMENT LED information COMMENT ================================================================== LEDSTAT1= '0000000000000000000000000000' / LED status LEDEXP1 = 0.0 / LED flash time (second) LEDSTAT2= '0000000000000000000000000000' / LED status LEDEXP2 = 0.0 / LED flash time (second) LEDTEMPA= -100.0 / LED temperature (main LED in K) LEDTEMPB= -100.0 / LED temperature (backup LED in K) COMMENT ================================================================== COMMENT Other information COMMENT ================================================================== CHECKSUM= 'abcdefg ' / SHA256 checksum of global headers END
SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T NEXTEND = 1 GROUPS = F DATE = '2021-03-04T09:30:00'/ date this file was written FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name FILETYPE= 'SCI ' / observation type TELESCOP= 'CSST ' / telescope used to acquire data INSTRUME= 'MSC ' / instrument used to acquire data RADECSYS= 'ICRS ' / frame of reference of coordinates EQUINOX = 2000.0 / FITSCREA= 'C6' / FITS create software version COMMENT ================================================================== COMMENT Object information COMMENT ================================================================== OBJECT = '00000000' / object name TARGET = '+000000000000' / target name (hhmmss+ddmmss) OBSID = '00000000' / observation ID OBJ_RA = 62.228226 / R.A. of the object (degrees) OBJ_DEC = -42.316932 / declination of the object (degrees) COMMENT ================================================================== COMMENT Telescope information COMMENT ================================================================== REFFRAME= 'CSSTGSC-1.0' / guide star catalog version DATE-OBS= '2021-03-04T09:30:00'/ date of the observation (yyyy-mm-dd hh:mm:ss) SATESWV = '0001' / software version in the satellite EXPSTART= 59130.5 / exposure start time (MJD) CABSTART= 59130.5 / (MJD) SUNANGL0= 50.0 / angle between sun and optical axis at CABSTART MOONANG0= 30.0 / angle moon and opt axis at CABSTART at CABST TEL_ALT0= 20.0 / angle opt axis and the ground-piston at CABST POS_ANG0= 20.0 / angle bwt y axis and the North Pole at CABST POSI0_X = 0.0 / the orbital position in X at CABSTART POSI0_Y = 0.0 / the orbital position in Y at CABSTART POSI0_Z = 0.0 / the orbital position in Z at CABSTART VELO0_X = 0.0 / the orbital velocity in X at CABSTART VELO0_Y = 0.0 / the orbital velocity in Y at CABSTART VELO0_Z = 0.0 / the orbital velocity in Z at CABSTART EULER0_1= 0.0 / Euler angle 1 at CABSTART EULER0_2= 0.0 / Euler angle 2 at CABSTART EULER0_3= 0.0 / Euler angle 3 at CABSTART RA_PNT0 = 0.0 / RA of the pointing (degrees) at CABSTART DEC_PNT0= 0.0 / DEC of the pointing (degrees) at CABSTART EXPEND = 0.0 / exposure end time (MJD) CABEND = 0.0 / (MJD) SUNANGL1= 50.0 / angle between sun and optical axis at CABEND MOONANG1= 30.0 / angle btw moon and optical axis at CABEND TEL_ALT1= 20.0 / angle opt axis and the ground-piston at CAEND POS_ANG1= 20.0 / angle bwt y axis and the North Pole at CABEND POSI1_X = 0.0 / the orbital position in X at CABEND POSI1_Y = 0.0 / the orbital position in Y at CABEND POSI1_Z = 0.0 / the orbital position in Z at CABEND VELO1_X = 0.0 / The orbital velocity in X at CABEND VELO1_Y = 0.0 / The orbital velocity in Y at CABEND VELO1_Z = 0.0 / The orbital velocity in Z at CABEND EULER1_1= 0.0 / Euler angle 1 at CABEND EULER1_2= 0.0 / Euler angle 2 at CABEND EULER1_3= 0.0 / Euler angle 3 at CABEND RA_PNT1 = 0.0 / RA of the pointing (degrees) at CABEND DEC_PNT1= 0.0 / DEC of the pointing (degrees) at CABEND EXPTIME = 150.0 / exposure duration EPOCH = 2000 / coordinate epoch COMMENT Other information COMMENT ================================================================== CHECKSUM= 'abcdefg ' / SHA256 checksum of global headers END
ObservationSim/Instrument/Chip/Chip.py
View file @ fd6c3108
import galsim
import os
import numpy as np
import pickle
import json
from astropy.table import Table
from numpy.random import Generator, PCG64
from astropy.io import fits
......@@ -9,6 +11,7 @@ from datetime import datetime
from ObservationSim.Instrument.Chip import Effects as effects
from ObservationSim.Instrument.FocalPlane import FocalPlane
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Instrument._util import rotate_conterclockwise
try:
import importlib.resources as pkg_resources
......@@ -21,30 +24,65 @@ class Chip(FocalPlane):
# Get focal plane (instance of paraent class) info
# TODO: use chipID to config individual chip?
super().__init__()
self.npix_x = 9216
self.npix_y = 9232
self.read_noise = 5.0 # e/pix
self.dark_noise = 0.02 # e/pix/s
self.pix_scale = 0.074 # pixel scale
self.gain = float(config["ins_effects"]["gain"])
self.bias_level = float(config["ins_effects"]["bias_level"])
self.overscan = 1000
self.exptime = float(config["obs_setting"]["exp_time"]) # second
self.dark_exptime = float(config["ins_effects"]['dark_exptime'])
self.flat_exptime = float(config["ins_effects"]['flat_exptime'])
self.readout_time = float(config["ins_effects"]['readout_time'])
self.full_well = int(config["ins_effects"]["full_well"])
# self.npix_x = 9216
# self.npix_y = 9232
# self.pix_scale = 0.074 # pixel scale
self._set_attributes_from_config(config)
self.logger = logger
# A chip ID must be assigned
self.chipID = int(chipID)
self._getChipRowCol()
self.chip_name = str(chipID).rjust(2, '0')
# Get corresponding filter info
self.filter_id, self.filter_type = self.getChipFilter()
self.survey_type = self._getSurveyType()
# [TODO]
if self.filter_type != "FGS":
self._getChipRowCol()
# Set the relavent specs for FGS detectors
# [TODO]
try:
with pkg_resources.files('ObservationSim.Instrument.data.ccd').joinpath("chip_definition.json") as chip_definition:
with open(chip_definition, "r") as f:
chip_dict = json.load(f)[str(self.chipID)]
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.ccd', "chip_definition.json") as chip_definition:
with open(chip_definition, "r") as f:
chip_dict = json.load(f)[str(self.chipID)]
for key in chip_dict:
setattr(self, key, chip_dict[key])
if self.filter_type == "FGS":
if config["ins_effects"]["field_dist"] == False:
self.fdModel = None
else:
fgs_name = self.chip_name[0:4]
try:
with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
with open(field_distortion, "rb") as f:
self.fdModel = pickle.load(f)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_pr4_%s.pickle"%(fgs_name.lower())) as field_distortion:
with open(field_distortion, "rb") as f:
self.fdModel = pickle.load(f)
else:
# Get the corresponding field distortion model
if config["ins_effects"]["field_dist"] == False:
self.fdModel = None
else:
try:
with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
with open(field_distortion, "rb") as f:
self.fdModel = pickle.load(f)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.field_distortion', "FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
with open(field_distortion, "rb") as f:
self.fdModel = pickle.load(f)
# Get boundary (in pix)
self.bound = self.getChipLim()
self.ccdEffCurve_dir = ccdEffCurve_dir
......@@ -76,17 +114,25 @@ class Chip(FocalPlane):
self.effCurve = self._getChipEffCurve(self.filter_type)
self._getCRdata()
# Define the sensor
# Define the sensor model
if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
self.sensor = galsim.SiliconSensor(strength=config["ins_effects"]["df_strength"], treering_func=treering_func)
self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
else:
self.sensor = galsim.Sensor()
self.flat_cube = None # for spectroscopic flat field cube simulation
# def _getChipRowCol(self):
# self.rowID = (self.chipID - 1) // 5 + 1
# self.colID = (self.chipID - 1) % 5 + 1
def _set_attributes_from_config(self, config):
# Default setting
self.read_noise = 5.0 # e/pix
self.dark_noise = 0.02 # e/pix/s
self.rotate_angle = 0.
self.overscan = 1000
# Override default values
for key in ["gain", "bias_level, dark_exptime", "flat_exptime", "readout_time", "full_well", "read_noise", "dark_noise", "overscan"]:
if key in config["ins_effects"]:
setattr(self, key, config["ins_effects"][key])
def _getChipRowCol(self):
self.rowID, self.colID = self.getChipRowCol(self.chipID)
......@@ -98,20 +144,21 @@ class Chip(FocalPlane):
def _getSurveyType(self):
if self.filter_type in ["GI", "GV", "GU"]:
return "spectroscopic"
else:
elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
return "photometric"
# elif self.filter_type in ["FGS"]:
# return "FGS"
def _getChipEffCurve(self, filter_type):
# CCD efficiency curves
if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt'
if filter_type in ['g', 'r', 'GV']: filename = 'Astro_MB.txt'
if filter_type in ['g', 'r', 'GV', 'FGS']: filename = 'Astro_MB.txt' # TODO, need to switch to the right efficiency curvey for FGS CMOS
if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
# Mirror efficiency:
# if filter_type == 'nuv': mirror_eff = 0.54
# if filter_type == 'u': mirror_eff = 0.68
# if filter_type in ['g', 'r', 'i', 'z', 'y']: mirror_eff = 0.8
# if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
# path = os.path.join(self.ccdEffCurve_dir, filename)
# table = Table.read(path, format='ascii')
try:
......@@ -120,14 +167,11 @@ class Chip(FocalPlane):
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.ccd', filename) as ccd_path:
table = Table.read(ccd_path, format='ascii')
# throughput = galsim.LookupTable(x=table['col1'], f=table['col2']*mirror_eff, interpolant='linear')
throughput = galsim.LookupTable(x=table['col1'], f=table['col2'], interpolant='linear')
bandpass = galsim.Bandpass(throughput, wave_type='nm')
return bandpass
def _getCRdata(self):
# path = os.path.join(self.CRdata_dir, 'wfc-cr-attachpixel.dat')
# self.attachedSizes = np.loadtxt(path)
try:
with pkg_resources.files('ObservationSim.Instrument.data').joinpath("wfc-cr-attachpixel.dat") as cr_path:
self.attachedSizes = np.loadtxt(cr_path)
......@@ -138,97 +182,82 @@ class Chip(FocalPlane):
def getChipFilter(self, chipID=None, filter_layout=None):
"""Return the filter index and type for a given chip #(chipID)
"""
filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI"]
# TODO: maybe a more elegent way other than hard coded?
# e.g. use something like a nested dict:
filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
if filter_layout is not None:
return filter_layout[chipID][0], filter_layout[chipID][1]
if chipID == None:
chipID = self.chipID
# updated configurations
# if chipID>30 or chipID<1: raise ValueError("!!! Chip ID: [1,30]")
# if chipID in [10, 15, 16, 21]: filter_type = 'y'
# if chipID in [11, 20]: filter_type = "z"
# if chipID in [9, 22]: filter_type = "i"
# if chipID in [12, 19]: filter_type = "u"
# if chipID in [7, 24]: filter_type = "r"
# if chipID in [14, 13, 18, 17]: filter_type = "nuv"
# if chipID in [8, 23]: filter_type = "g"
# if chipID in [6, 5, 25, 26]: filter_type = "GI"
# if chipID in [27, 30, 1, 4]: filter_type = "GV"
# if chipID in [28, 29, 2, 3]: filter_type = "GU"
if chipID in [6, 15, 16, 25]: filter_type = "y"
if chipID>42 or chipID<1: raise ValueError("!!! Chip ID: [1,42]")
if chipID in [6, 15, 16, 25]: filter_type = "y"
if chipID in [11, 20]: filter_type = "z"
if chipID in [7, 24]: filter_type = "i"
if chipID in [7, 24]: filter_type = "i"
if chipID in [14, 17]: filter_type = "u"
if chipID in [9, 22]: filter_type = "r"
if chipID in [12, 13, 18, 19]: filter_type = "nuv"
if chipID in [8, 23]: filter_type = "g"
if chipID in [1, 10, 21, 30]: filter_type = "GI"
if chipID in [2, 5, 26, 29]: filter_type = "GV"
if chipID in [3, 4, 27, 28]: filter_type = "GU"
if chipID in [9, 22]: filter_type = "r"
if chipID in [12, 13, 18, 19]: filter_type = "nuv"
if chipID in [8, 23]: filter_type = "g"
if chipID in [1, 10, 21, 30]: filter_type = "GI"
if chipID in [2, 5, 26, 29]: filter_type = "GV"
if chipID in [3, 4, 27, 28]: filter_type = "GU"
if chipID in range(31, 43): filter_type = 'FGS'
filter_id = filter_type_list.index(filter_type)
return filter_id, filter_type
def getChipLim(self, chipID=None):
"""Calculate the edges in pixel for a given CCD chip on the focal plane
NOTE: There are 5*4 CCD chips in the focus plane for photometric observation.
NOTE: There are 5*4 CCD chips in the focus plane for photometric / spectroscopic observation.
Parameters:
chipID: int
the index of the chip
Returns:
A galsim BoundsD object
"""
# if chipID == None:
# chipID = self.chipID
# gx = self.npix_gap_x
# gy1, gy2 = self.npix_gap_y
# # xlim of a given ccd chip
# xrem = (chipID-1)%self.nchip_x - self.nchip_x // 2
# xcen = (self.npix_x + gx) * xrem
# nx0 = xcen - self.npix_x//2 + 1
# nx1 = xcen + self.npix_x//2
# # ylim of a given ccd chip
# yrem = 2*((chipID-1)//self.nchip_x) - (self.nchip_y-1)
# ycen = (self.npix_y//2 + gy1//2) * yrem
# if chipID <= 6: ycen = (self.npix_y//2 + gy1//2) * yrem - (gy2-gy1)
# if chipID >= 25: ycen = (self.npix_y//2 + gy1//2) * yrem + (gy2-gy1)
# ny0 = ycen - self.npix_y//2 + 1
# ny1 = ycen + self.npix_y//2
if chipID == None:
chipID = self.chipID
rowID, colID = self.rowID, self.colID
if ((chipID is not None) and (int(chipID) <= 30)) or (self.chipID <= 30):
# [TODO]
if chipID == None:
chipID = self.chipID
rowID, colID = self.rowID, self.colID
else:
rowID, colID = self.getChipRowCol(chipID)
gx1, gx2 = self.npix_gap_x
gy = self.npix_gap_y
# xlim of a given CCD chip
xrem = 2*(colID - 1) - (self.nchip_x - 1)
xcen = (self.npix_x//2 + gx1//2) * xrem
if chipID >= 26 or chipID == 21:
xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1)
if chipID <= 5 or chipID == 10:
xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1)
nx0 = xcen - self.npix_x//2 + 1
nx1 = xcen + self.npix_x//2
# ylim of a given CCD chip
yrem = (rowID - 1) - self.nchip_y // 2
ycen = (self.npix_y + gy) * yrem
ny0 = ycen - self.npix_y//2 + 1
ny1 = ycen + self.npix_y//2
return galsim.BoundsD(nx0-1, nx1-1, ny0-1, ny1-1)
else:
rowID, colID = self.getChipRowCol(chipID)
gx1, gx2 = self.npix_gap_x
gy = self.npix_gap_y
# xlim of a given CCD chip
xrem = 2*(colID - 1) - (self.nchip_x - 1)
xcen = (self.npix_x//2 + gx1//2) * xrem
if chipID >= 26 or chipID == 21:
xcen = (self.npix_x//2 + gx1//2) * xrem - (gx2-gx1)
if chipID <= 5 or chipID == 10:
xcen = (self.npix_x//2 + gx1//2) * xrem + (gx2-gx1)
nx0 = xcen - self.npix_x//2 + 1
nx1 = xcen + self.npix_x//2
xmin, xmax, ymin, ymax = 1e10, -1e10, 1e10, -1e10
xcen = self.x_cen / self.pix_size
ycen = self.y_cen / self.pix_size
sign_x = [-1., 1., -1., 1.]
sign_y = [-1., -1., 1., 1.]
for i in range(4):
x = xcen + sign_x[i] * self.npix_x / 2.
y = ycen + sign_y[i] * self.npix_y / 2.
x, y = rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
xmin, xmax = min(xmin, x), max(xmax, x)
ymin, ymax = min(ymin, y), max(ymax, y)
return galsim.BoundsD(xmin, xmax, ymin, ymax)
# ylim of a given CCD chip
yrem = (rowID - 1) - self.nchip_y // 2
ycen = (self.npix_y + gy) * yrem
ny0 = ycen - self.npix_y//2 + 1
ny1 = ycen + self.npix_y//2
return galsim.BoundsD(nx0-1, nx1-1, ny0-1, ny1-1)
def getSkyCoverage(self, wcs):
# print("In getSkyCoverage: xmin = %.3f, xmax = %.3f, ymin = %.3f, ymax = %.3f"%(self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax))
return super().getSkyCoverage(wcs, self.bound.xmin, self.bound.xmax, self.bound.ymin, self.bound.ymax)
......@@ -239,14 +268,19 @@ class Chip(FocalPlane):
bound = self.getSkyCoverage(wcs)
return galsim.BoundsD(bound.xmin - margin, bound.xmax + margin, bound.ymin - margin, bound.ymax + margin)
def isContainObj(self, ra_obj, dec_obj, wcs=None, margin=1):
def isContainObj(self, ra_obj=None, dec_obj=None, x_image=None, y_image=None, wcs=None, margin=1):
# magin in number of pix
if wcs is None:
wcs = self.img.wcs
pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
if (ra_obj is not None) and (dec_obj is not None):
if wcs is None:
wcs = self.img.wcs
pos_obj = wcs.toImage(galsim.CelestialCoord(ra=ra_obj*galsim.degrees, dec=dec_obj*galsim.degrees))
x_image, y_image = pos_obj.x, pos_obj.y
elif (x_image is None) or (y_image is None):
raise ValueError("Either (ra_obj, dec_obj) or (x_image, y_image) should be given")
xmin, xmax = self.bound.xmin - margin, self.bound.xmax + margin
ymin, ymax = self.bound.ymin - margin, self.bound.ymax + margin
if (pos_obj.x - xmin) * (pos_obj.x - xmax) > 0.0 or (pos_obj.y - ymin) * (pos_obj.y - ymax) > 0.0:
if (x_image - xmin) * (x_image - xmax) > 0.0 or (y_image - ymin) * (y_image - ymax) > 0.0:
return False
return True
......@@ -277,27 +311,28 @@ class Chip(FocalPlane):
pointNum = str(pointing_ID),
ra=ra_cen,
dec=dec_cen,
psize=self.pix_scale,
row_num=self.rowID,
col_num=self.colID,
pixel_scale=self.pix_scale,
date=date_obs,
time_obs=time_obs,
im_type = im_type,
exptime=exptime
exptime=exptime,
chip_name=str(self.chipID).rjust(2, '0')
)
h_ext = generateExtensionHeader(
chip=self,
xlen=self.npix_x,
ylen=self.npix_y,
ra=ra_cen,
dec=dec_cen,
dec=dec_cen,
pa=img_rot.deg,
gain=self.gain,
readout=self.read_noise,
dark=self.dark_noise,
saturation=90000,
psize=self.pix_scale,
row_num=self.rowID,
col_num=self.colID,
pixel_scale=self.pix_scale,
pixel_size=self.pix_size,
xcen=self.x_cen,
ycen=self.y_cen,
extName='raw')
return h_prim, h_ext
......@@ -312,7 +347,9 @@ class Chip(FocalPlane):
time_obs=time_obs,
exptime=exptime)
hdu1 = fits.PrimaryHDU(header=h_prim)
hdu1.add_checksum()
hdu2 = fits.ImageHDU(img.array, header=h_ext)
hdu2.add_checksum()
hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
......@@ -324,7 +361,7 @@ class Chip(FocalPlane):
SeedBadColumns = int(config["random_seeds"]["seed_badcolumns"])
SeedDefective = int(config["random_seeds"]["seed_defective"])
SeedCosmicRay = int(config["random_seeds"]["seed_CR"])
fullwell = int(config["ins_effects"]["full_well"])
fullwell = int(self.full_well)
if config["ins_effects"]["add_hotpixels"] == True:
BoolHotPix = True
else:
......@@ -342,7 +379,7 @@ class Chip(FocalPlane):
# Add sky background
if sky_map is None:
sky_map = filt.getSkyNoise(exptime=self.exptime)
sky_map = filt.getSkyNoise(exptime=exptime)
sky_map = sky_map * np.ones_like(img.array)
sky_map = galsim.Image(array=sky_map)
# Apply Poisson noise to the sky map
......@@ -354,7 +391,7 @@ class Chip(FocalPlane):
elif img.array.shape != sky_map.shape:
raise ValueError("The shape img and sky_map must be equal.")
elif tel is not None: # If sky_map is given in flux
sky_map = sky_map * tel.pupil_area * self.exptime
sky_map = sky_map * tel.pupil_area * exptime
if config["ins_effects"]["add_back"] == True:
img += sky_map
del sky_map
......@@ -407,8 +444,8 @@ class Chip(FocalPlane):
print(" Adding Cosmic-Ray", flush=True)
cr_map, cr_event_num = effects.produceCR_Map(
xLen=self.npix_x, yLen=self.npix_y,
exTime=self.exptime+0.5*self.readout_time,
cr_pixelRatio=0.003*(self.exptime+0.5*self.readout_time)/600.,
exTime=exptime+0.5*self.readout_time,
cr_pixelRatio=0.003*(exptime+0.5*self.readout_time)/600.,
gain=self.gain,
attachedSizes=self.attachedSizes,
seed=SeedCosmicRay+pointing_ID*30+self.chipID) # seed: obj-imaging:+0; bias:+1; dark:+2; flat:+3;
......@@ -419,7 +456,7 @@ class Chip(FocalPlane):
del cr_map
# crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID))
# crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
datetime_obs = datetime.fromtimestamp(timestamp_obs)
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal(
......@@ -432,7 +469,7 @@ class Chip(FocalPlane):
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir,
exptime=self.exptime)
exptime=exptime)
del crmap_gsimg
# Apply PRNU effect and output PRNU flat file:
......@@ -454,7 +491,7 @@ class Chip(FocalPlane):
# Add dark current
if config["ins_effects"]["add_dark"] == True:
dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(self.exptime+0.5*self.readout_time)))
dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(rng_poisson, self.dark_noise*(exptime+0.5*self.readout_time)))
img.addNoise(dark_noise)
# Add Hot Pixels or/and Dead Pixels
......@@ -474,7 +511,7 @@ class Chip(FocalPlane):
print(" Adding Bias level and 16-channel non-uniformity")
if config["ins_effects"]["bias_16channel"] == True:
img = effects.AddBiasNonUniform16(img,
bias_level=float(config["ins_effects"]["bias_level"]),
bias_level=float(self.bias_level),
nsecy = 2, nsecx=8,
seed=SeedBiasNonuni+self.chipID,
logger=self.logger)
......@@ -539,11 +576,11 @@ class Chip(FocalPlane):
self.logger.info(" Output N frame Bias files")
else:
print(" Output N frame Bias files", flush=True)
NBias = int(config["ins_effects"]["NBias"])
NBias = int(config["output_setting"]["NBias"])
for i in range(NBias):
BiasCombImg, BiasTag = effects.MakeBiasNcomb(
self.npix_x, self.npix_y,
bias_level=float(config["ins_effects"]["bias_level"]),
bias_level=float(self.bias_level),
ncombine=1, read_noise=self.read_noise, gain=1,
seed=SeedBiasNonuni+self.chipID,
logger=self.logger)
......@@ -571,7 +608,7 @@ class Chip(FocalPlane):
# Apply Bad lines
if config["ins_effects"]["add_badcolumns"] == True:
BiasCombImg = effects.BadColumns(BiasCombImg-float(config["ins_effects"]["bias_level"])+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(config["ins_effects"]["bias_level"])-5
BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.bias_level)-5
BiasCombImg = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
nsecy = 2, nsecx=8,
......@@ -585,7 +622,7 @@ class Chip(FocalPlane):
BiasCombImg.quantize()
BiasCombImg = galsim.ImageUS(BiasCombImg)
# BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
datetime_obs = datetime.fromtimestamp(timestamp_obs)
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60
......@@ -608,7 +645,7 @@ class Chip(FocalPlane):
self.logger.info(" Output N frame Flat-Field files")
else:
print(" Output N frame Flat-Field files", flush=True)
NFlat = int(config["ins_effects"]["NFlat"])
NFlat = int(config["output_setting"]["NFlat"])
if config["ins_effects"]["add_bias"] == True:
biaslevel = self.bias_level
overscan = biaslevel-2
......@@ -689,7 +726,7 @@ class Chip(FocalPlane):
FlatCombImg.quantize()
FlatCombImg = galsim.ImageUS(FlatCombImg)
# FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
datetime_obs = datetime.fromtimestamp(timestamp_obs)
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60
......@@ -717,7 +754,7 @@ class Chip(FocalPlane):
self.logger.info(" Output N frame Dark Current files")
else:
print(" Output N frame Dark Current files", flush=True)
NDark = int(config["ins_effects"]["NDark"])
NDark = int(config["output_setting"]["NDark"])
if config["ins_effects"]["add_bias"] == True:
biaslevel = self.bias_level
overscan = biaslevel-2
......@@ -746,7 +783,7 @@ class Chip(FocalPlane):
cr_map[cr_map < 0] = 0
crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
del cr_map
datetime_obs = datetime.fromtimestamp(timestamp_obs)
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal(
......@@ -815,7 +852,7 @@ class Chip(FocalPlane):
DarkCombImg.quantize()
DarkCombImg = galsim.ImageUS(DarkCombImg)
# DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
datetime_obs = datetime.fromtimestamp(timestamp_obs)
datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60
......
ObservationSim/Instrument/Filter.py
View file @ fd6c3108
......@@ -17,8 +17,9 @@ class Filter(object):
self.filter_id = filter_id
self.filter_type = filter_type
self.ccd_bandpass = ccd_bandpass
# Load basic filter parameters.
# Assume t_exp = 150s for limiting/saturation magnitudes.
self._getParam(filter_param, filter_type)
# if self.filter_dir is not None:
self.bandpass_full, self.bandpass_sub_list = self._get_bandpasses()
self.survey_type = self._getSurveyType()
......@@ -37,7 +38,6 @@ class Filter(object):
self.sky_background = filter_param.param[filter_type][5]
self.mag_saturation = filter_param.param[filter_type][6]
self.mag_limiting = filter_param.param[filter_type][7]
# self.filter_dir = filter_param.filter_dir
def is_too_bright(self, mag, margin=-2.5):
return mag <= self.mag_saturation + margin
......@@ -48,9 +48,6 @@ class Filter(object):
def _get_bandpasses(self, filter_dir=None, unit='A'):
if self.filter_id < 7: # Photometric
# Get full-bandpass
# filter_file = os.path.join(filter_dir, self.filter_type+".dat")
# bandpass_full = galsim.Bandpass(filter_file, wave_type=unit)
try:
with pkg_resources.files('ObservationSim.Instrument.data.filters').joinpath(self.filter_type.lower() + '.txt') as filter_file:
self.filter_bandpass = galsim.Bandpass(str(filter_file), wave_type=unit)
......@@ -127,3 +124,5 @@ class Filter(object):
throughput_file = self.filter_type.lower() + '_throughput.txt'
self.mag_limiting, self.mag_saturation = calculateLimitMag(psf_fwhm=psf_fwhm, pixelSize=pix_scale, throughputFn=throughput_file, readout=5.0, skyFn=skyFn, darknoise=dark_noise, exTime=exptime, fw=full_well)
print("for filter %s: mag_limiting: %.3f, mag_saturation: %.3f"%(self.filter_type, self.mag_limiting, self.mag_saturation))
ObservationSim/Instrument/FilterParam.py
View file @ fd6c3108
......@@ -33,6 +33,7 @@ class FilterParam(object):
"i": [7653.2, 1588.1, 6760.0, 8550.0, 0.4960, 0.212, 16.7, 25.9],
"z": [9600.6, 2490.5, 8240.0, 11000.0, 0.2000, 0.123, 15.7, 25.2],
"y": [10051.0, 1590.6, 9130.0, 11000.0, 0.0960, 0.037, 14.4, 24.4],
"FGS": [5000.0, 8000.0, 3000.0, 11000.0, 0.6500, 0.164, 0., 30.], # [TODO]
"GU": [0.0, 0.0, 2550.0, 4200.0, 1.0, 0.037, 14.0, 26.0],
"GV": [0.0, 0.0, 4000.0, 6500.0, 1.0, 0.037, 14.0, 26.0],
......
ObservationSim/Instrument/FocalPlane.py
View file @ fd6c3108
......@@ -4,8 +4,8 @@ import numpy as np
class FocalPlane(object):
def __init__(self, config=None, survey_type='Photometric', bad_chips=None):
"""Get the focal plane layout
Note: (If applicable) make sure the config file have the most basic info, e.g.: "nchip_x"...
"""
self.nchips = 42
if bad_chips == None:
self.bad_chips = []
else:
......@@ -13,22 +13,24 @@ class FocalPlane(object):
self.ignore_chips = []
if survey_type == 'Photometric':
# for i in range(6):
# self.ignore_chips.append(i+1)
# self.ignore_chips.append(i+25)
for i in range(5):
self.ignore_chips.append(i+1)
self.ignore_chips.append(i+26)
self.ignore_chips.append(10)
self.ignore_chips.append(21)
for i in range(31, 43):
self.ignore_chips.append(i)
elif survey_type == 'Spectroscopic':
# for i in range(18):
# self.ignore_chips.append(i+7)
for i in range(6, 26):
if i == 10 or i == 21:
continue
else:
self.ignore_chips.append(i)
for i in range(31, 43):
self.ignore_chips.append(i)
elif survey_type == 'FGS':
for i in range(1, 31):
self.ignore_chips.append(i)
if config is not None:
self.nchip_x = config["nchip_x"]
......@@ -42,12 +44,6 @@ class FocalPlane(object):
if "bad_chips" in config:
self.bad_chips = config["bad_chips"]
else:
# self.nchip_x = 5
# self.nchip_y = 6
# self.npix_tot_x = 49752
# self.npix_tot_y = 59516
# self.npix_gap_x = 898
# self.npix_gap_y = (534, 1309)
self.nchip_x = 6
self.nchip_y = 5
self.npix_tot_x = 59516
......@@ -58,8 +54,6 @@ class FocalPlane(object):
self._getCenter()
def _getCenter(self):
# self.cen_pix_x = 0.5 * (self.npix_tot_x + 1.0)
# self.cen_pix_y = 0.5 * (self.npix_tot_y + 1.0)
self.cen_pix_x = 0
self.cen_pix_y = 0
......
ObservationSim/Instrument/Telescope.py
View file @ fd6c3108
......@@ -41,7 +41,4 @@ class Telescope(object):
columns = line.split()
opticsEff[str(columns[0])] = float(columns[2])
f.close()
return opticsEff
def _get_effCurve(self, effCurve_path):
pass
\ No newline at end of file
return opticsEff
\ No newline at end of file
ObservationSim/Instrument/_util.py
View file @ fd6c3108
......@@ -14,6 +14,17 @@ VC_A = 2.99792458e+18 # speed of light: A/s
VC_M = 2.99792458e+8 # speed of light: m/s
H_PLANK = 6.626196e-27 # Plank constant: erg s
def rotate_conterclockwise(x0, y0, x, y, angle):
"""
Rotate a point counterclockwise by a given angle around a given origin.
The angle should be given in radians.
"""
angle = np.deg2rad(angle)
qx = x0 + np.cos(angle)*(x - x0) - np.sin(angle) * (y - y0)
qy = y0 + np.sin(angle)*(x - x0) + np.cos(angle) * (y - y0)
return qx, qy
def photonEnergy(lambd):
""" The energy of photon at a given wavelength
......@@ -26,23 +37,22 @@ def photonEnergy(lambd):
eph = H_PLANK * nu
return eph
'''
description:
param {*} aperture: unit m, default 2 m
param {*} psf_fwhm: psf fwhm, default 0.1969"
param {*} pixelSize: pixel size, default 0.074"
param {*} pmRation: the ratio of souce flux in the limit mag calculation
param {*} throughputFn: throuput file name
param {*} readout: unit, e-/pixel
param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
param {*} darknoise: unit, e-/pixel/s
param {*} exTime: exposure time one time, default 150s
param {*} exNum: exposure number, defautl 1
param {*} fw, full well value( or saturation value),default 90000e-/pixel
return {*} limit mag and saturation mag
'''
def calculateLimitMag(aperture = 2.0, psf_fwhm = 0.1969,pixelSize = 0.074, pmRation = 0.8, throughputFn = 'i_throughput.txt', readout = 5.0, skyFn= 'sky_emiss_hubble_50_50_A.dat', darknoise = 0.02,exTime = 150, exNum = 1, fw = 90000):
'''
description:
param {*} aperture: unit m, default 2 m
param {*} psf_fwhm: psf fwhm, default 0.1969"
param {*} pixelSize: pixel size, default 0.074"
param {*} pmRation: the ratio of souce flux in the limit mag calculation
param {*} throughputFn: throuput file name
param {*} readout: unit, e-/pixel
param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
param {*} darknoise: unit, e-/pixel/s
param {*} exTime: exposure time one time, default 150s
param {*} exNum: exposure number, defautl 1
param {*} fw, full well value( or saturation value),default 90000e-/pixel
return {*} limit mag and saturation mag
'''
try:
with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(throughputFn) as data_file:
throughput_f = np.loadtxt(data_file)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment