Skip to content
GitLab
Commit 3226117c authored by Zhang Xin's avatar Zhang Xin
Browse files

Merge branch 'develop' into 'release_v2.0' 

version 2.1.0

See merge request csst_sim/csst-simulation!16
parents 81589f9d f540664f
Show whitespace changes
Inline Side-by-side
Catalog/C6_50sqdeg.py 0 → 100644
View file @ 3226117c
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
from astropy.coordinates import SkyCoord
from astropy.io import fits
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
NSIDE = 128
bundle_file_list = ['galaxies_C6_bundle000199.h5','galaxies_C6_bundle000200.h5','galaxies_C6_bundle000241.h5','galaxies_C6_bundle000242.h5','galaxies_C6_bundle000287.h5','galaxies_C6_bundle000288.h5','galaxies_C6_bundle000714.h5','galaxies_C6_bundle000715.h5','galaxies_C6_bundle000778.h5','galaxies_C6_bundle000779.h5','galaxies_C6_bundle000842.h5','galaxies_C6_bundle000843.h5','galaxies_C6_bundle002046.h5','galaxies_C6_bundle002110.h5','galaxies_C6_bundle002111.h5','galaxies_C6_bundle002173.h5','galaxies_C6_bundle002174.h5','galaxies_C6_bundle002238.h5','galaxies_C6_bundle002596.h5','galaxies_C6_bundle002597.h5','galaxies_C6_bundle002656.h5','galaxies_C6_bundle002657.h5','galaxies_C6_bundle002711.h5','galaxies_C6_bundle002712.h5','galaxies_C6_bundle002844.h5','galaxies_C6_bundle002845.h5','galaxies_C6_bundle002884.h5','galaxies_C6_bundle002885.h5','galaxies_C6_bundle002921.h5','galaxies_C6_bundle002922.h5']
qsosed_file_list = ['quickspeclib_interp1d_run1.fits','quickspeclib_interp1d_run2.fits','quickspeclib_interp1d_run3.fits','quickspeclib_interp1d_run4.fits','quickspeclib_interp1d_run5.fits','quickspeclib_interp1d_run6.fits','quickspeclib_interp1d_run7.fits','quickspeclib_interp1d_run8.fits','quickspeclib_interp1d_run9.fits','quickspeclib_interp1d_run10.fits','quickspeclib_interp1d_run11.fits','quickspeclib_interp1d_run12.fits','quickspeclib_interp1d_run13.fits','quickspeclib_interp1d_run14.fits','quickspeclib_interp1d_run15.fits','quickspeclib_interp1d_run16.fits','quickspeclib_interp1d_run17.fits','quickspeclib_interp1d_run18.fits','quickspeclib_interp1d_run19.fits','quickspeclib_interp1d_run20.fits','quickspeclib_interp1d_run21.fits','quickspeclib_interp1d_run22.fits','quickspeclib_interp1d_run23.fits','quickspeclib_interp1d_run24.fits','quickspeclib_interp1d_run25.fits','quickspeclib_interp1d_run26.fits','quickspeclib_interp1d_run27.fits','quickspeclib_interp1d_run28.fits','quickspeclib_interp1d_run29.fits','quickspeclib_interp1d_run30.fits']
star_file_list = ['C7_Gaia_Galaxia_RA170DECm23_healpix.hdf5', 'C7_Gaia_Galaxia_RA180DECp60_healpix.hdf5', 'C7_Gaia_Galaxia_RA240DECp30_healpix.hdf5', 'C7_Gaia_Galaxia_RA300DECm60_healpix.hdf5', 'C7_Gaia_Galaxia_RA30DECm48_healpix.hdf5']
star_center_list = [(170., -23.), (180., 60.), (240., 30.), (300., -60.), (30., -48.)]
def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
assert NSIDE == 2**healpixOrder
shift = healpixOrder - bundleOrder
shift = 2*shift
nside_bundle = 2**bundleOrder
nside_healpix= 2**healpixOrder
healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
bundleID_nest = (healpixID_nest >> shift)
bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
return bundleID_ring
def get_agnsed_file(bundle_file_name):
return qsosed_file_list[bundle_file_list.index(bundle_file_name)]
def get_star_cat(ra_pointing, dec_pointing):
pointing_c = SkyCoord(ra=ra_pointing*U.deg, dec=dec_pointing*U.deg)
max_dist = 10
return_star_path = None
for star_file, center in zip(star_file_list, star_center_list):
center_c = SkyCoord(ra=center[0]*U.deg, dec=center[1]*U.deg)
dist = pointing_c.separation(center_c).to(U.deg).value
if dist < max_dist:
return_star_path = star_file
max_dist = dist
return return_star_path
class Catalog(CatalogBase):
def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.seed_Av = config["catalog_options"]["seed_Av"]
self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
self.chip_output = chip_output
self.filt = filt
self.logger = chip_output.logger
with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
self.normF_star = Table.read(str(filter_path))
self.config = config
self.chip = chip
self.pointing = pointing
self.max_size = 0.
if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
# Get the cloest star catalog file
star_file_name = get_star_cat(ra_pointing=self.pointing.ra, dec_pointing=self.pointing.dec)
star_path = os.path.join(config["catalog_options"]["input_path"]["star_cat"], star_file_name)
star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
self.star_path = os.path.join(self.cat_dir, star_path)
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()
self.agn_seds = {}
if "AGN_SED" in config["catalog_options"]["SED_templates_path"] and not config["catalog_options"]["star_only"]:
self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
if "rotateEll" in config["catalog_options"]:
self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else:
self.rotation = 0.
# Update output .cat header with catalog specific output columns
self._add_output_columns_header()
self._get_healpix_list()
self._load()
def _add_output_columns_header(self):
self.add_hdr = " model_tag teff logg feh"
self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
self.add_fmt = " %10s %8.4f %8.4f %8.4f"
self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
self.chip_output.update_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]))
self.pix_list = hp.query_polygon(
NSIDE,
hp.ang2vec(np.radians(90.) - dec, ra),
inclusive=True
)
if self.logger is not None:
msg = str(("HEALPix List: ", self.pix_list))
self.logger.info(msg)
else:
print("HEALPix List: ", self.pix_list)
def load_norm_filt(self, obj):
if obj.type == "star":
return self.normF_star
elif obj.type == "galaxy" or obj.type == "quasar":
# return self.normF_galaxy
return None
else:
return None
def _load_SED_lib_star(self):
self.tempSED_star = h5.File(self.star_SED_path,'r')
def _load_SED_lib_gals(self):
pcs = h5.File(os.path.join(self.galaxy_SED_path, "pcs.h5"), "r")
lamb = h5.File(os.path.join(self.galaxy_SED_path, "lamb.h5"), "r")
self.lamb_gal = lamb['lamb'][()]
self.pcs = pcs['pcs'][()]
def _load_gals(self, gals, pix_id=None, cat_id=0, agnsed_file=""):
ngals = len(gals['ra'])
# Apply astrometric modeling
ra_arr = gals['ra'][:]
dec_arr = gals['dec'][:]
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = ra_arr.tolist()
dec_list = dec_arr.tolist()
pmra_list = np.zeros(ngals).tolist()
pmdec_list = np.zeros(ngals).tolist()
rv_list = np.zeros(ngals).tolist()
parallax_list = [1e-9] * ngals
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
input_ra_list=ra_list,
input_dec_list=dec_list,
input_pmra_list=pmra_list,
input_pmdec_list=pmdec_list,
input_rv_list=rv_list,
input_parallax_list=parallax_list,
input_nstars=ngals,
input_x=self.pointing.sat_x,
input_y=self.pointing.sat_y,
input_z=self.pointing.sat_z,
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2000",
input_date_str=date_str,
input_time_str=time_str
)
for igals in range(ngals):
# # (TEST)
# if igals > 100:
# break
param = self.initialize_param()
param['ra'] = ra_arr[igals]
param['dec'] = dec_arr[igals]
param['ra_orig'] = gals['ra'][igals]
param['dec_orig'] = gals['dec'][igals]
# param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
if self.filt.filter_type == 'NUV':
param['mag_use_normal'] = gals['mag_csst_nuv'][igals]
else:
param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
continue
param['z'] = gals['redshift'][igals]
param['model_tag'] = 'None'
param['g1'] = gals['shear'][igals][0]
param['g2'] = gals['shear'][igals][1]
param['kappa'] = gals['kappa'][igals]
param['e1'] = gals['ellipticity_true'][igals][0]
param['e2'] = gals['ellipticity_true'][igals][1]
# For shape calculation
param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity(
e1=gals['ellipticity_true'][igals][0],
e2=gals['ellipticity_true'][igals][1],
rotation=self.rotation,
unit='radians')
# param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
if param['ell_total'] > 0.9:
continue
# phi_e = cmath.phase(complex(param['e1'], param['e2']))
# param['e1'] = param['ell_total'] * np.cos(phi_e + 2*self.rotation)
# param['e2'] = param['ell_total'] * np.sin(phi_e + 2*self.rotation)
param['e1_disk'] = param['e1']
param['e2_disk'] = param['e2']
param['e1_bulge'] = param['e1']
param['e2_bulge'] = param['e2']
param['delta_ra'] = 0
param['delta_dec'] = 0
# Masses
param['bulgemass'] = gals['bulgemass'][igals]
param['diskmass'] = gals['diskmass'][igals]
param['size'] = gals['size'][igals]
if param['size'] > self.max_size:
self.max_size = param['size']
# Sersic index
param['disk_sersic_idx'] = 1.
param['bulge_sersic_idx'] = 4.
# Sizes
param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
if param['bfrac'] >= 0.6:
param['hlr_bulge'] = param['size']
param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
else:
param['hlr_disk'] = param['size']
param['hlr_bulge'] = param['size'] * param['bfrac']
# SED coefficients
param['coeff'] = gals['coeff'][igals]
param['detA'] = gals['detA'][igals]
# Others
param['galType'] = gals['type'][igals]
param['veldisp'] = gals['veldisp'][igals]
# TEST no redening and no extinction
param['av'] = 0.0
param['redden'] = 0
# Is this an Quasar?
param['qsoindex'] = gals['qsoindex'][igals]
if param['qsoindex'] == -1:
param['star'] = 0 # Galaxy
param['agnsed_file'] = ""
else:
param['star'] = 2 # Quasar
param['agnsed_file'] = agnsed_file
# 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'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0:
obj = Galaxy(param, logger=self.logger)
elif param['star'] == 2:
obj = Quasar(param, logger=self.logger)
# Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
param['bulgemass'], param['diskmass'], param['detA'],
param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
param['galType'], param['veldisp'])
self.objs.append(obj)
def _load_stars(self, stars, pix_id=None):
nstars = len(stars['sourceID'])
# Apply astrometric modeling
ra_arr = stars["RA"][:]
dec_arr = stars["Dec"][:]
pmra_arr = stars['pmra'][:]
pmdec_arr = stars['pmdec'][:]
rv_arr = stars['RV'][:]
parallax_arr = stars['parallax'][:]
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = ra_arr.tolist()
dec_list = dec_arr.tolist()
pmra_list = pmra_arr.tolist()
pmdec_list = pmdec_arr.tolist()
rv_list = rv_arr.tolist()
parallax_list = parallax_arr.tolist()
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
input_ra_list=ra_list,
input_dec_list=dec_list,
input_pmra_list=pmra_list,
input_pmdec_list=pmdec_list,
input_rv_list=rv_list,
input_parallax_list=parallax_list,
input_nstars=nstars,
input_x=self.pointing.sat_x,
input_y=self.pointing.sat_y,
input_z=self.pointing.sat_z,
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2000",
input_date_str=date_str,
input_time_str=time_str
)
for istars in range(nstars):
# # (TEST)
# if istars > 100:
# break
param = self.initialize_param()
param['ra'] = ra_arr[istars]
param['dec'] = dec_arr[istars]
param['ra_orig'] = stars["RA"][istars]
param['dec_orig'] = stars["Dec"][istars]
param['pmra'] = pmra_arr[istars]
param['pmdec'] = pmdec_arr[istars]
param['rv'] = rv_arr[istars]
param['parallax'] = parallax_arr[istars]
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
param['mag_use_normal'] = stars['app_sdss_g'][istars]
self.ids += 1
param['id'] = stars['sourceID'][istars]
param['sed_type'] = stars['sourceID'][istars]
param['model_tag'] = stars['model_tag'][istars]
param['teff'] = stars['teff'][istars]
param['logg'] = stars['grav'][istars]
param['feh'] = stars['feh'][istars]
param['z'] = 0.0
param['star'] = 1 # Star
obj = Star(param, logger=self.logger)
# Append additional output columns to the .cat file
obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
self.objs.append(obj)
def _load(self, **kwargs):
self.objs = []
self.ids = 0
if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
star_cat = h5.File(self.star_path, 'r')['catalog']
for pix in self.pix_list:
try:
stars = star_cat[str(pix)]
self._load_stars(stars, pix_id=pix)
del stars
except Exception as e:
self.logger.error(str(e))
print(e)
if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
for pix in self.pix_list:
try:
bundleID = get_bundleIndex(pix)
bundle_file = "galaxies_C6_bundle{:06}.h5".format(bundleID)
file_path = os.path.join(self.galaxy_path, bundle_file)
gals_cat = h5.File(file_path, 'r')['galaxies']
gals = gals_cat[str(pix)]
# Get corresponding AGN SED file
agnsed_file = get_agnsed_file(bundle_file)
agnsed_path = os.path.join(self.AGN_SED_path, agnsed_file)
self.agn_seds[agnsed_file] = fits.open(agnsed_path)[0].data
self._load_gals(gals, pix_id=pix, cat_id=bundleID, agnsed_file=agnsed_file)
del gals
except Exception as e:
traceback.print_exc()
self.logger.error(str(e))
print(e)
if self.logger is not None:
self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
else:
print("number of objects in catalog: ", len(self.objs))
def load_sed(self, obj, **kwargs):
if obj.type == 'star':
_, wave, flux = tag_sed(
h5file=self.tempSED_star,
model_tag=obj.param['model_tag'],
teff=obj.param['teff'],
logg=obj.param['logg'],
feh=obj.param['feh']
)
elif obj.type == 'galaxy' or obj.type == 'quasar':
factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
if obj.type == 'galaxy':
flux = np.matmul(self.pcs, obj.coeff) * factor
# if np.any(flux < 0):
# raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
flux[flux < 0] = 0.
sedcat = np.vstack((self.lamb_gal, flux)).T
sed_data = getObservedSED(
sedCat=sedcat,
redshift=obj.z,
av=obj.param["av"],
redden=obj.param["redden"]
)
wave, flux = sed_data[0], sed_data[1]
elif obj.type == 'quasar':
flux = self.agn_seds[obj.agnsed_file][int(obj.qsoindex)] * 1e-17
flux[flux < 0] = 0.
wave = self.lamb_gal
else:
raise ValueError("Object type not known")
speci = interpolate.interp1d(wave, flux)
lamb = np.arange(2000, 11001+0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
if obj.type == 'quasar':
# integrate to get the magnitudes
sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
# mag = getABMAG(interFlux, self.filt.bandpass_full)
# print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
del wave
del flux
return sed
Catalog/C6_Catalog.py
View file @ 3226117c
...@@ -83,7 +83,7 @@ class Catalog(CatalogBase): ...@@ -83,7 +83,7 @@ class Catalog(CatalogBase):
self._load_SED_lib_AGN() self._load_SED_lib_AGN()
if "rotateEll" in config["catalog_options"]: if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else: else:
self.rotation = 0. self.rotation = 0.
...@@ -259,7 +259,7 @@ class Catalog(CatalogBase): ...@@ -259,7 +259,7 @@ class Catalog(CatalogBase):
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0: if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger) obj = Galaxy(param, logger=self.logger)
# Need to deal with additional output columns # Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
......
Catalog/C6_fits_Catalog.py
View file @ 3226117c
...@@ -94,7 +94,7 @@ class Catalog(CatalogBase): ...@@ -94,7 +94,7 @@ class Catalog(CatalogBase):
###mock_stamp_END ###mock_stamp_END
if "rotateEll" in config["catalog_options"]: if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else: else:
self.rotation = 0. self.rotation = 0.
...@@ -272,7 +272,7 @@ class Catalog(CatalogBase): ...@@ -272,7 +272,7 @@ class Catalog(CatalogBase):
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0: if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger) obj = Galaxy(param, logger=self.logger)
# Need to deal with additional output columns # Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
......
Catalog/FGS_Catalog.py
View file @ 3226117c
...@@ -84,7 +84,7 @@ class Catalog(CatalogBase): ...@@ -84,7 +84,7 @@ class Catalog(CatalogBase):
self._load_SED_lib_AGN() self._load_SED_lib_AGN()
if "rotateEll" in config["catalog_options"]: if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else: else:
self.rotation = 0. self.rotation = 0.
...@@ -260,7 +260,7 @@ class Catalog(CatalogBase): ...@@ -260,7 +260,7 @@ class Catalog(CatalogBase):
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0: if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger) obj = Galaxy(param, logger=self.logger)
# Need to deal with additional output columns # Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
......
Catalog/NGPCatalog.py
View file @ 3226117c
...@@ -52,7 +52,7 @@ class Catalog(CatalogBase): ...@@ -52,7 +52,7 @@ class Catalog(CatalogBase):
self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["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() self._load_SED_lib_gals()
if "rotateEll" in config["catalog_options"]: if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else: else:
self.rotation = 0. self.rotation = 0.
...@@ -191,7 +191,7 @@ class Catalog(CatalogBase): ...@@ -191,7 +191,7 @@ class Catalog(CatalogBase):
param['id'] = gals['galaxyID'][igals] param['id'] = gals['galaxyID'][igals]
if param['star'] == 0: if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger) obj = Galaxy(param, logger=self.logger)
if param['star'] == 2: if param['star'] == 2:
obj = Quasar(param, logger=self.logger) obj = Quasar(param, logger=self.logger)
......
Catalog/fd_test_C6.py 0 → 100644
View file @ 3226117c
import os
import galsim
import random
import numpy as np
import h5py as h5
import healpy as hp
import astropy.constants as cons
import traceback
from astropy.coordinates import spherical_to_cartesian
from astropy.table import Table
from scipy import interpolate
from datetime import datetime
from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
# (TEST)
from astropy.cosmology import FlatLambdaCDM
from astropy import constants
from astropy import units as U
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
NSIDE = 128
def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
assert NSIDE == 2**healpixOrder
shift = healpixOrder - bundleOrder
shift = 2*shift
nside_bundle = 2**bundleOrder
nside_healpix= 2**healpixOrder
healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
bundleID_nest = (healpixID_nest >> shift)
bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
return bundleID_ring
class Catalog(CatalogBase):
def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.seed_Av = config["catalog_options"]["seed_Av"]
self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
self.chip_output = chip_output
self.filt = filt
self.logger = chip_output.logger
with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
self.normF_star = Table.read(str(filter_path))
self.config = config
self.chip = chip
self.pointing = pointing
self.max_size = 0.
if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
star_file = config["catalog_options"]["input_path"]["star_cat"]
star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
self.star_path = os.path.join(self.cat_dir, star_file)
self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
self._load_SED_lib_star()
if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
self._load_SED_lib_gals()
if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
self._load_SED_lib_AGN()
if "rotateEll" in config["catalog_options"]:
self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else:
self.rotation = 0.
# Update output .cat header with catalog specific output columns
self._add_output_columns_header()
self._get_healpix_list()
self._load()
def _add_output_columns_header(self):
self.add_hdr = " model_tag teff logg feh"
self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
self.add_fmt = " %10s %8.4f %8.4f %8.4f"
self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
self.chip_output.update_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):
for igals in range(0, ngals, 5):
# # (TEST)
# if igals > 1000:
# break
param = self.initialize_param()
param['ra'] = ra_arr[igals]
param['dec'] = dec_arr[igals]
param['ra_orig'] = gals['ra'][igals]
param['dec_orig'] = gals['dec'][igals]
# param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
param['mag_use_normal'] = 20.
if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
continue
param['z'] = gals['redshift'][igals]
param['model_tag'] = 'None'
# param['g1'] = gals['shear'][igals][0]
# param['g2'] = gals['shear'][igals][1]
param['g1'] = 0.
param['g2'] = 0.
param['kappa'] = gals['kappa'][igals]
param['e1'] = gals['ellipticity_true'][igals][0]
param['e2'] = gals['ellipticity_true'][igals][1]
# For shape calculation
param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
if param['ell_total'] > 0.9:
continue
# param['e1_disk'] = param['e1']
# param['e2_disk'] = param['e2']
# param['e1_bulge'] = param['e1']
# param['e2_bulge'] = param['e2']
param['e1_disk'] = 0.
param['e2_disk'] = 0.
param['e1_bulge'] = 0.
param['e2_bulge'] = 0.
param['delta_ra'] = 0
param['delta_dec'] = 0
# Masses
param['bulgemass'] = gals['bulgemass'][igals]
param['diskmass'] = gals['diskmass'][igals]
# param['size'] = gals['size'][igals]
param['size'] = 1.
if param['size'] > self.max_size:
self.max_size = param['size']
# Sersic index
param['disk_sersic_idx'] = 1.
param['bulge_sersic_idx'] = 4.
# Sizes
# param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
param['bfrac'] = 0.
if param['bfrac'] >= 0.6:
param['hlr_bulge'] = param['size']
param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
else:
param['hlr_disk'] = param['size']
param['hlr_bulge'] = param['size'] * param['bfrac']
# SED coefficients
param['coeff'] = gals['coeff'][igals]
param['detA'] = gals['detA'][igals]
# Others
param['galType'] = gals['type'][igals]
param['veldisp'] = gals['veldisp'][igals]
# TEST no redening and no extinction
param['av'] = 0.0
param['redden'] = 0
param['star'] = 0 # Galaxy
# NOTE: this cut cannot be put before the SED type has been assigned
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
# TEMP
self.ids += 1
# param['id'] = self.ids
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0:
obj = Galaxy(param, logger=self.logger)
# Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
param['bulgemass'], param['diskmass'], param['detA'],
param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
param['galType'], param['veldisp'])
self.objs.append(obj)
def _load_stars(self, stars, pix_id=None):
nstars = len(stars['sourceID'])
# Apply astrometric modeling
ra_arr = stars["RA"][:]
dec_arr = stars["Dec"][:]
pmra_arr = stars['pmra'][:]
pmdec_arr = stars['pmdec'][:]
rv_arr = stars['RV'][:]
parallax_arr = stars['parallax'][:]
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = ra_arr.tolist()
dec_list = dec_arr.tolist()
pmra_list = pmra_arr.tolist()
pmdec_list = pmdec_arr.tolist()
rv_list = rv_arr.tolist()
parallax_list = parallax_arr.tolist()
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
input_ra_list=ra_list,
input_dec_list=dec_list,
input_pmra_list=pmra_list,
input_pmdec_list=pmdec_list,
input_rv_list=rv_list,
input_parallax_list=parallax_list,
input_nstars=nstars,
input_x=self.pointing.sat_x,
input_y=self.pointing.sat_y,
input_z=self.pointing.sat_z,
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2000",
input_date_str=date_str,
input_time_str=time_str
)
for istars in range(nstars):
# # (TEST)
# if istars > 100:
# break
param = self.initialize_param()
param['ra'] = ra_arr[istars]
param['dec'] = dec_arr[istars]
param['ra_orig'] = stars["RA"][istars]
param['dec_orig'] = stars["Dec"][istars]
param['pmra'] = pmra_arr[istars]
param['pmdec'] = pmdec_arr[istars]
param['rv'] = rv_arr[istars]
param['parallax'] = parallax_arr[istars]
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
param['mag_use_normal'] = stars['app_sdss_g'][istars]
# if param['mag_use_normal'] >= 26.5:
# continue
self.ids += 1
param['id'] = stars['sourceID'][istars]
param['sed_type'] = stars['sourceID'][istars]
param['model_tag'] = stars['model_tag'][istars]
param['teff'] = stars['teff'][istars]
param['logg'] = stars['grav'][istars]
param['feh'] = stars['feh'][istars]
param['z'] = 0.0
param['star'] = 1 # Star
obj = Star(param, logger=self.logger)
# Append additional output columns to the .cat file
obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
self.objs.append(obj)
def _load_AGNs(self):
data = Table.read(self.AGN_path)
ra_arr = data['ra']
dec_arr = data['dec']
nAGNs = len(data)
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = ra_arr.tolist()
dec_list = dec_arr.tolist()
pmra_list = np.zeros(nAGNs).tolist()
pmdec_list = np.zeros(nAGNs).tolist()
rv_list = np.zeros(nAGNs).tolist()
parallax_list = [1e-9] * nAGNs
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
input_ra_list=ra_list,
input_dec_list=dec_list,
input_pmra_list=pmra_list,
input_pmdec_list=pmdec_list,
input_rv_list=rv_list,
input_parallax_list=parallax_list,
input_nstars=nAGNs,
input_x=self.pointing.sat_x,
input_y=self.pointing.sat_y,
input_z=self.pointing.sat_z,
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2000",
input_date_str=date_str,
input_time_str=time_str
)
for iAGNs in range(nAGNs):
param = self.initialize_param()
param['ra'] = ra_arr[iAGNs]
param['dec'] = dec_arr[iAGNs]
param['ra_orig'] = data['ra'][iAGNs]
param['dec_orig'] = data['dec'][iAGNs]
param['z'] = data['z'][iAGNs]
param['appMag'] = data['appMag'][iAGNs]
param['absMag'] = data['absMag'][iAGNs]
# NOTE: this cut cannot be put before the SED type has been assigned
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
# TEST no redening and no extinction
param['av'] = 0.0
param['redden'] = 0
param['star'] = 2 # Quasar
param['id'] = data['igmlos'][iAGNs]
if param['star'] == 2:
obj = Quasar(param, logger=self.logger)
# Append additional output columns to the .cat file
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
self.objs.append(obj)
def _load(self, **kwargs):
self.objs = []
self.ids = 0
# if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
# star_cat = h5.File(self.star_path, 'r')['catalog']
# for pix in self.pix_list:
# try:
# stars = star_cat[str(pix)]
# self._load_stars(stars, pix_id=pix)
# del stars
# except Exception as e:
# self.logger.error(str(e))
# print(e)
if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
for pix in self.pix_list:
try:
bundleID = get_bundleIndex(pix)
file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
gals_cat = h5.File(file_path, 'r')['galaxies']
gals = gals_cat[str(pix)]
self._load_gals(gals, pix_id=pix, cat_id=bundleID)
del gals
except Exception as e:
traceback.print_exc()
self.logger.error(str(e))
print(e)
# if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
# try:
# self._load_AGNs()
# except Exception as e:
# traceback.print_exc()
# self.logger.error(str(e))
# print(e)
if self.logger is not None:
self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
else:
print("number of objects in catalog: ", len(self.objs))
def load_sed(self, obj, **kwargs):
if obj.type == 'star':
_, wave, flux = tag_sed(
h5file=self.tempSED_star,
model_tag=obj.param['model_tag'],
teff=obj.param['teff'],
logg=obj.param['logg'],
feh=obj.param['feh']
)
elif obj.type == 'galaxy' or obj.type == 'quasar':
factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
if obj.type == 'galaxy':
flux = np.matmul(self.pcs, obj.coeff) * factor
# if np.any(flux < 0):
# raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
flux[flux < 0] = 0.
sedcat = np.vstack((self.lamb_gal, flux)).T
sed_data = getObservedSED(
sedCat=sedcat,
redshift=obj.z,
av=obj.param["av"],
redden=obj.param["redden"]
)
wave, flux = sed_data[0], sed_data[1]
elif obj.type == 'quasar':
flux = self.SED_AGN[int(obj.id)] * 1e-17
# if np.any(flux < 0):
# raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
flux[flux < 0] = 0.
# sedcat = np.vstack((self.lamb_AGN, flux)).T
wave = self.lamb_AGN
# print("sed (erg/s/cm2/A) = ", sed_data)
# np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
else:
raise ValueError("Object type not known")
speci = interpolate.interp1d(wave, flux)
lamb = np.arange(2000, 11001+0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
if obj.type == 'quasar':
# integrate to get the magnitudes
sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
# if obj.param['mag_use_normal'] >= 30:
# print("obj ID = %d"%obj.id)
# print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
# print("integrated flux = %.7f"%(interFlux))
# print("app mag = %.3f"%obj.param['appMag'])
# np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
# print("obj ID = %d"%obj.id)
# print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
# print("integrated flux = %.7f"%(interFlux))
# print("app mag = %.3f"%obj.param['appMag'])
# print("abs mag = %.3f"%obj.param['absMag'])
# mag = getABMAG(interFlux, self.filt.bandpass_full)
# print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
del wave
del flux
return sed
Catalog/testCat_galaxy.py
View file @ 3226117c
...@@ -81,7 +81,7 @@ class Catalog(CatalogBase): ...@@ -81,7 +81,7 @@ class Catalog(CatalogBase):
self._load_SED_lib_AGN() self._load_SED_lib_AGN()
if "rotateEll" in config["catalog_options"]: if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.)) self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else: else:
self.rotation = 0. self.rotation = 0.
...@@ -253,7 +253,7 @@ class Catalog(CatalogBase): ...@@ -253,7 +253,7 @@ class Catalog(CatalogBase):
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals) param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0: if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger) obj = Galaxy(param, logger=self.logger)
# Need to deal with additional output columns # Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0., obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
......
Catalog/wcs_test_C6.py 0 → 100644
View file @ 3226117c
import os
import galsim
import random
import numpy as np
import h5py as h5
import healpy as hp
import astropy.constants as cons
import traceback
from astropy.coordinates import spherical_to_cartesian
from astropy.table import Table
from scipy import interpolate
from datetime import datetime
from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
# (TEST)
from astropy.cosmology import FlatLambdaCDM
from astropy import constants
from astropy import units as U
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
NSIDE = 128
def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
assert NSIDE == 2**healpixOrder
shift = healpixOrder - bundleOrder
shift = 2*shift
nside_bundle = 2**bundleOrder
nside_healpix= 2**healpixOrder
healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
bundleID_nest = (healpixID_nest >> shift)
bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
return bundleID_ring
class Catalog(CatalogBase):
def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.seed_Av = config["catalog_options"]["seed_Av"]
self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
self.chip_output = chip_output
self.filt = filt
self.logger = chip_output.logger
with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
self.normF_star = Table.read(str(filter_path))
self.config = config
self.chip = chip
self.pointing = pointing
self.max_size = 0.
if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"] and not config["catalog_options"]["galaxy_only"]:
star_file = config["catalog_options"]["input_path"]["star_cat"]
star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
self.star_path = os.path.join(self.cat_dir, star_file)
self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
self._load_SED_lib_star()
if "galaxy_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["galaxy_cat"] and not config["catalog_options"]["star_only"]:
galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
self._load_SED_lib_gals()
if "AGN_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["AGN_cat"] and not config["catalog_options"]["star_only"]:
AGN_dir = config["catalog_options"]["input_path"]["AGN_cat"]
self.AGN_path = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["AGN_cat"])
self.AGN_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["AGN_SED"])
self.AGN_SED_wave_path = os.path.join(config['data_dir'], config["catalog_options"]["SED_templates_path"]["AGN_SED_WAVE"])
self._load_SED_lib_AGN()
if "rotateEll" in config["catalog_options"]:
self.rotation = np.radians(float(config["catalog_options"]["rotateEll"]))
else:
self.rotation = 0.
# Update output .cat header with catalog specific output columns
self._add_output_columns_header()
self._get_healpix_list()
self._load()
def _add_output_columns_header(self):
self.add_hdr = " model_tag teff logg feh"
self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
self.add_fmt = " %10s %8.4f %8.4f %8.4f"
self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
self.chip_output.update_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]
# param['mag_use_normal'] = 20.
if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
continue
param['z'] = gals['redshift'][igals]
param['model_tag'] = 'None'
# param['g1'] = gals['shear'][igals][0]
# param['g2'] = gals['shear'][igals][1]
param['g1'] = 0.
param['g2'] = 0.
param['kappa'] = gals['kappa'][igals]
# param['e1'] = gals['ellipticity_true'][igals][0]
# param['e2'] = gals['ellipticity_true'][igals][1]
param['e1'] = 0.
param['e2'] = 0.
# For shape calculation
param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
if param['ell_total'] > 0.9:
continue
param['e1_disk'] = param['e1']
param['e2_disk'] = param['e2']
param['e1_bulge'] = param['e1']
param['e2_bulge'] = param['e2']
param['delta_ra'] = 0
param['delta_dec'] = 0
# Masses
param['bulgemass'] = gals['bulgemass'][igals]
param['diskmass'] = gals['diskmass'][igals]
param['size'] = gals['size'][igals]
if param['size'] > self.max_size:
self.max_size = param['size']
# Sersic index
param['disk_sersic_idx'] = 1.
param['bulge_sersic_idx'] = 4.
# Sizes
param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
if param['bfrac'] >= 0.6:
param['hlr_bulge'] = param['size']
param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
else:
param['hlr_disk'] = param['size']
param['hlr_bulge'] = param['size'] * param['bfrac']
# SED coefficients
param['coeff'] = gals['coeff'][igals]
param['detA'] = gals['detA'][igals]
# Others
param['galType'] = gals['type'][igals]
param['veldisp'] = gals['veldisp'][igals]
# TEST no redening and no extinction
param['av'] = 0.0
param['redden'] = 0
param['star'] = 0 # Galaxy
# NOTE: this cut cannot be put before the SED type has been assigned
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
# TEMP
self.ids += 1
# param['id'] = self.ids
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
if param['star'] == 0:
obj = Galaxy(param, logger=self.logger)
# Need to deal with additional output columns
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
param['bulgemass'], param['diskmass'], param['detA'],
param['e1'], param['e2'], param['kappa'], param['g1'], param['g2'], param['size'],
param['galType'], param['veldisp'])
self.objs.append(obj)
def _load_stars(self, stars, pix_id=None):
nstars = len(stars['sourceID'])
# Apply astrometric modeling
ra_arr = stars["RA"][:]
dec_arr = stars["Dec"][:]
pmra_arr = stars['pmra'][:]
pmdec_arr = stars['pmdec'][:]
rv_arr = stars['RV'][:]
parallax_arr = stars['parallax'][:]
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = ra_arr.tolist()
dec_list = dec_arr.tolist()
pmra_list = pmra_arr.tolist()
pmdec_list = pmdec_arr.tolist()
rv_list = rv_arr.tolist()
parallax_list = parallax_arr.tolist()
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
input_ra_list=ra_list,
input_dec_list=dec_list,
input_pmra_list=pmra_list,
input_pmdec_list=pmdec_list,
input_rv_list=rv_list,
input_parallax_list=parallax_list,
input_nstars=nstars,
input_x=self.pointing.sat_x,
input_y=self.pointing.sat_y,
input_z=self.pointing.sat_z,
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2000",
input_date_str=date_str,
input_time_str=time_str
)
for istars in range(nstars):
# # (TEST)
# if istars > 100:
# break
param = self.initialize_param()
param['ra'] = ra_arr[istars]
param['dec'] = dec_arr[istars]
param['ra_orig'] = stars["RA"][istars]
param['dec_orig'] = stars["Dec"][istars]
param['pmra'] = pmra_arr[istars]
param['pmdec'] = pmdec_arr[istars]
param['rv'] = rv_arr[istars]
param['parallax'] = parallax_arr[istars]
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
# param['mag_use_normal'] = stars['app_sdss_g'][istars]
param['mag_use_normal'] = 20.
# if param['mag_use_normal'] >= 26.5:
# continue
self.ids += 1
param['id'] = stars['sourceID'][istars]
param['sed_type'] = stars['sourceID'][istars]
param['model_tag'] = stars['model_tag'][istars]
param['teff'] = stars['teff'][istars]
param['logg'] = stars['grav'][istars]
param['feh'] = stars['feh'][istars]
param['z'] = 0.0
param['star'] = 1 # Star
obj = Star(param, logger=self.logger)
# Append additional output columns to the .cat file
obj.additional_output_str = self.add_fmt%(param["model_tag"], param['teff'], param['logg'], param['feh'],
0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
self.objs.append(obj)
def _load_AGNs(self):
data = Table.read(self.AGN_path)
ra_arr = data['ra']
dec_arr = data['dec']
nAGNs = len(data)
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = ra_arr.tolist()
dec_list = dec_arr.tolist()
pmra_list = np.zeros(nAGNs).tolist()
pmdec_list = np.zeros(nAGNs).tolist()
rv_list = np.zeros(nAGNs).tolist()
parallax_list = [1e-9] * nAGNs
dt = datetime.utcfromtimestamp(self.pointing.timestamp)
date_str = dt.date().isoformat()
time_str = dt.time().isoformat()
ra_arr, dec_arr = on_orbit_obs_position(
input_ra_list=ra_list,
input_dec_list=dec_list,
input_pmra_list=pmra_list,
input_pmdec_list=pmdec_list,
input_rv_list=rv_list,
input_parallax_list=parallax_list,
input_nstars=nAGNs,
input_x=self.pointing.sat_x,
input_y=self.pointing.sat_y,
input_z=self.pointing.sat_z,
input_vx=self.pointing.sat_vx,
input_vy=self.pointing.sat_vy,
input_vz=self.pointing.sat_vz,
input_epoch="J2000",
input_date_str=date_str,
input_time_str=time_str
)
for iAGNs in range(nAGNs):
param = self.initialize_param()
param['ra'] = ra_arr[iAGNs]
param['dec'] = dec_arr[iAGNs]
param['ra_orig'] = data['ra'][iAGNs]
param['dec_orig'] = data['dec'][iAGNs]
param['z'] = data['z'][iAGNs]
param['appMag'] = data['appMag'][iAGNs]
param['absMag'] = data['absMag'][iAGNs]
# NOTE: this cut cannot be put before the SED type has been assigned
if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200):
continue
# TEST no redening and no extinction
param['av'] = 0.0
param['redden'] = 0
param['star'] = 2 # Quasar
param['id'] = data['igmlos'][iAGNs]
if param['star'] == 2:
obj = Quasar(param, logger=self.logger)
# Append additional output columns to the .cat file
obj.additional_output_str = self.add_fmt%("n", 0., 0., 0.,
0., 0., 0., 0., 0., 0., 0., 0., 0., -1, 0.)
self.objs.append(obj)
def _load(self, **kwargs):
self.objs = []
self.ids = 0
if "star_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["star_cat"] and not self.config["catalog_options"]["galaxy_only"]:
star_cat = h5.File(self.star_path, 'r')['catalog']
for pix in self.pix_list:
try:
stars = star_cat[str(pix)]
self._load_stars(stars, pix_id=pix)
del stars
except Exception as e:
self.logger.error(str(e))
print(e)
if "galaxy_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["galaxy_cat"] and not self.config["catalog_options"]["star_only"]:
for pix in self.pix_list:
try:
bundleID = get_bundleIndex(pix)
file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
gals_cat = h5.File(file_path, 'r')['galaxies']
gals = gals_cat[str(pix)]
self._load_gals(gals, pix_id=pix, cat_id=bundleID)
del gals
except Exception as e:
traceback.print_exc()
self.logger.error(str(e))
print(e)
if "AGN_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["AGN_cat"] and not self.config["catalog_options"]["star_only"]:
try:
self._load_AGNs()
except Exception as e:
traceback.print_exc()
self.logger.error(str(e))
print(e)
if self.logger is not None:
self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
else:
print("number of objects in catalog: ", len(self.objs))
def load_sed(self, obj, **kwargs):
if obj.type == 'star':
_, wave, flux = tag_sed(
h5file=self.tempSED_star,
model_tag=obj.param['model_tag'],
teff=obj.param['teff'],
logg=obj.param['logg'],
feh=obj.param['feh']
)
elif obj.type == 'galaxy' or obj.type == 'quasar':
factor = 10**(-.4 * self.cosmo.distmod(obj.z).value)
if obj.type == 'galaxy':
flux = np.matmul(self.pcs, obj.coeff) * factor
# if np.any(flux < 0):
# raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
flux[flux < 0] = 0.
sedcat = np.vstack((self.lamb_gal, flux)).T
sed_data = getObservedSED(
sedCat=sedcat,
redshift=obj.z,
av=obj.param["av"],
redden=obj.param["redden"]
)
wave, flux = sed_data[0], sed_data[1]
elif obj.type == 'quasar':
flux = self.SED_AGN[int(obj.id)] * 1e-17
# if np.any(flux < 0):
# raise ValueError("Glaxy %s: negative SED fluxes"%obj.id)
flux[flux < 0] = 0.
# sedcat = np.vstack((self.lamb_AGN, flux)).T
wave = self.lamb_AGN
# print("sed (erg/s/cm2/A) = ", sed_data)
# np.savetxt(os.path.join(self.config["work_dir"], "%s_sed.txt"%(obj.id)), sedcat)
else:
raise ValueError("Object type not known")
speci = interpolate.interp1d(wave, flux)
lamb = np.arange(2000, 11001+0.5, 0.5)
y = speci(lamb)
# erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
if obj.type == 'quasar':
# integrate to get the magnitudes
sed_photon = np.array([sed['WAVELENGTH'], sed['FLUX']]).T
sed_photon = galsim.LookupTable(x=np.array(sed_photon[:, 0]), f=np.array(sed_photon[:, 1]), interpolant='nearest')
sed_photon = galsim.SED(sed_photon, wave_type='A', flux_type='1', fast=False)
interFlux = integrate_sed_bandpass(sed=sed_photon, bandpass=self.filt.bandpass_full)
obj.param['mag_use_normal'] = getABMAG(interFlux, self.filt.bandpass_full)
# if obj.param['mag_use_normal'] >= 30:
# print("obj ID = %d"%obj.id)
# print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
# print("integrated flux = %.7f"%(interFlux))
# print("app mag = %.3f"%obj.param['appMag'])
# np.savetxt('./AGN_SED_test/sed_objID_%d.txt'%obj.id, np.transpose([self.lamb_AGN, self.SED_AGN[int(obj.id)]]))
# print("obj ID = %d"%obj.id)
# print("mag_use_normal = %.3f"%obj.param['mag_use_normal'])
# print("integrated flux = %.7f"%(interFlux))
# print("app mag = %.3f"%obj.param['appMag'])
# print("abs mag = %.3f"%obj.param['absMag'])
# mag = getABMAG(interFlux, self.filt.bandpass_full)
# print("mag diff = %.3f"%(mag - obj.param['mag_use_normal']))
del wave
del flux
return sed
ObservationSim/Config/Header/ImageHeader.py
View file @ 3226117c
...@@ -18,6 +18,10 @@ from astropy.wcs.utils import fit_wcs_from_points ...@@ -18,6 +18,10 @@ from astropy.wcs.utils import fit_wcs_from_points
from astropy.time import Time from astropy.time import Time
from astropy import wcs from astropy import wcs
from datetime import datetime
# import socket
import platform
def chara2digit(char): def chara2digit(char):
""" Function to judge and convert characters to digitals """ Function to judge and convert characters to digitals
...@@ -287,7 +291,8 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r ...@@ -287,7 +291,8 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
r_dat['CRPIX1'] = -x1 r_dat['CRPIX1'] = -x1
r_dat['CRPIX2'] = -y1 r_dat['CRPIX2'] = -y1
cd = np.array([[ pixel_scale, 0], [0, pixel_scale]])/3600.*flag_x[1] # cd = np.array([[ pixel_scale, 0], [0, pixel_scale]])/3600.*flag_x[1]
cd = np.array([[ pixel_scale, 0], [0, -pixel_scale]])/3600.
cd_rot = rotate_CD_matrix(cd, pa_aper) cd_rot = rotate_CD_matrix(cd, pa_aper)
r_dat['CD1_1'] = cd_rot[0,0] r_dat['CD1_1'] = cd_rot[0,0]
r_dat['CD1_2'] = cd_rot[0,1] r_dat['CD1_2'] = cd_rot[0,1]
...@@ -381,7 +386,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec ...@@ -381,7 +386,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
h_prim['EXPTIME'] = exptime h_prim['EXPTIME'] = exptime
# Define file types # Define file types
file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'cosmic_ray', 'CRD':'cosmic_ray','CALS':'CALS','CALF':'CALF'} file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
h_prim['FILETYPE'] = file_type[im_type] h_prim['FILETYPE'] = file_type[im_type]
co = coord.SkyCoord(ra, dec, unit='deg') co = coord.SkyCoord(ra, dec, unit='deg')
...@@ -417,7 +422,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec ...@@ -417,7 +422,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
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] 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 + '_' + OBS_id + '_' + CCDID[ # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
# k - 1].rjust(2, '0') + '_L0_V01' # 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['FILENAME'] = 'CSST_MSC_MS_' + file_type[im_type] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
h_prim['POSI0_X'] = sat_pos[0] h_prim['POSI0_X'] = sat_pos[0]
...@@ -433,12 +438,15 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec ...@@ -433,12 +438,15 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
# Get version of CSSTSim Package # Get version of CSSTSim Package
from pkg_resources import get_distribution 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 currentDateAndTime = datetime.now()
compute_name = platform.node()
h_prim['FITSCREA'] = get_distribution("CSSTSim").version +'_' + currentDateAndTime.strftime("%Y%m%d") + '_' +compute_name
h_prim['EPOCH'] = round((Time(h_prim['EXPSTART'], format='mjd', scale='tcb')).jyear, 1)
return h_prim return h_prim
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, 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): extName='SCI', row_num = None, col_num = None, xcen=None, ycen=None, timestamp = 1621915200,exptime = 150., readoutTime = 40.):
e_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/extension_header.header' 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') f = open(os.path.split(os.path.realpath(__file__))[0] + '/filter.lst')
...@@ -460,47 +468,60 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p ...@@ -460,47 +468,60 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
# h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0') # h_ext['CCDCHIP'] = CCDID[k - 1].rjust(2, '0')
# h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1] # h_ext['CCDLABEL'] = filters[k-1] + '-' + filterID[k-1]
# h_ext['FILTER'] = filters[k-1] # h_ext['FILTER'] = filters[k-1]
h_ext['CCDCHIP'] = chip.chip_name h_ext['CHIPID'] = str(chip.chipID).rjust(2, '0')
h_ext['CCDLABEL'] = chip.filter_type + '-' + str(chip.chipID).rjust(2, '0') h_ext['CHIPLAB'] = chip.chip_name
h_ext['FILTER'] = chip.filter_type h_ext['FILTER'] = chip.filter_type
h_ext['NAXIS1'] = xlen h_ext['NAXIS1'] = xlen
h_ext['NAXIS2'] = ylen h_ext['NAXIS2'] = ylen
h_ext['EXTNAME'] = extName h_ext['EXTNAME'] = extName
h_ext['GAIN1'] = gain h_ext['GAIN01'] = chip.gain_channel[0]
h_ext['GAIN2'] = gain h_ext['GAIN02'] = chip.gain_channel[1]
h_ext['GAIN3'] = gain h_ext['GAIN03'] = chip.gain_channel[2]
h_ext['GAIN4'] = gain h_ext['GAIN04'] = chip.gain_channel[3]
h_ext['GAIN5'] = gain h_ext['GAIN05'] = chip.gain_channel[4]
h_ext['GAIN6'] = gain h_ext['GAIN06'] = chip.gain_channel[5]
h_ext['GAIN7'] = gain h_ext['GAIN07'] = chip.gain_channel[6]
h_ext['GAIN8'] = gain h_ext['GAIN08'] = chip.gain_channel[7]
h_ext['GAIN9'] = gain h_ext['GAIN09'] = chip.gain_channel[8]
h_ext['GAIN10'] = gain h_ext['GAIN10'] = chip.gain_channel[9]
h_ext['GAIN11'] = gain h_ext['GAIN11'] = chip.gain_channel[10]
h_ext['GAIN12'] = gain h_ext['GAIN12'] = chip.gain_channel[11]
h_ext['GAIN13'] = gain h_ext['GAIN13'] = chip.gain_channel[12]
h_ext['GAIN14'] = gain h_ext['GAIN14'] = chip.gain_channel[13]
h_ext['GAIN15'] = gain h_ext['GAIN15'] = chip.gain_channel[14]
h_ext['GAIN16'] = gain h_ext['GAIN16'] = chip.gain_channel[15]
h_ext['RDNOIS1'] = readout h_ext['RON01'] = readout
h_ext['RDNOIS2'] = readout h_ext['RON02'] = readout
h_ext['RDNOIS3'] = readout h_ext['RON03'] = readout
h_ext['RDNOIS4'] = readout h_ext['RON04'] = readout
h_ext['RDNOIS5'] = readout h_ext['RON05'] = readout
h_ext['RDNOIS6'] = readout h_ext['RON06'] = readout
h_ext['RDNOIS7'] = readout h_ext['RON07'] = readout
h_ext['RDNOIS8'] = readout h_ext['RON08'] = readout
h_ext['RDNOIS9'] = readout h_ext['RON09'] = readout
h_ext['RDNOIS10'] = readout h_ext['RON10'] = readout
h_ext['RDNOIS11'] = readout h_ext['RON11'] = readout
h_ext['RDNOIS12'] = readout h_ext['RON12'] = readout
h_ext['RDNOIS13'] = readout h_ext['RON13'] = readout
h_ext['RDNOIS14'] = readout h_ext['RON14'] = readout
h_ext['RDNOIS15'] = readout h_ext['RON15'] = readout
h_ext['RDNOIS16'] = readout h_ext['RON16'] = readout
h_ext['PIXSCAL1'] = pixel_scale h_ext['PIXSCAL1'] = pixel_scale
h_ext['PIXSCAL2'] = pixel_scale h_ext['PIXSCAL2'] = pixel_scale
h_ext['EXPTIME'] = exptime
h_ext['DARKTIME'] = exptime + 2
datetime_obs = datetime.utcfromtimestamp(timestamp)
tstart = Time(datetime_obs)
tstart_read = Time(tstart.mjd + exptime / 86400., format="mjd")
tend_read = Time(tstart.mjd + (exptime + readoutTime) / 86400., format="mjd")
t_s1 = str(tstart_read.datetime).split()
h_ext['READT0'] = t_s1[0]+'T'+t_s1[1]
t_s2 = str(tend_read.datetime).split()
h_ext['READT1'] = t_s2[0] + 'T' + t_s2[1]
# h_ext['POS_ANG'] = pa # 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, 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,
...@@ -514,8 +535,8 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p ...@@ -514,8 +535,8 @@ def generateExtensionHeader(chip, xlen = 9216, ylen = 9232,ra = 60, dec = -40, p
h_ext['CD1_2'] = header_wcs['CD1_2'] h_ext['CD1_2'] = header_wcs['CD1_2']
h_ext['CD2_1'] = header_wcs['CD2_1'] h_ext['CD2_1'] = header_wcs['CD2_1']
h_ext['CD2_2'] = header_wcs['CD2_2'] h_ext['CD2_2'] = header_wcs['CD2_2']
h_ext['EQUINOX'] = header_wcs['EQUINOX'] # h_ext['EQUINOX'] = header_wcs['EQUINOX']
h_ext['WCSDIM'] = header_wcs['WCSDIM'] # h_ext['WCSDIM'] = header_wcs['WCSDIM']
h_ext['CTYPE1'] = header_wcs['CTYPE1'] h_ext['CTYPE1'] = header_wcs['CTYPE1']
h_ext['CTYPE2'] = header_wcs['CTYPE2'] h_ext['CTYPE2'] = header_wcs['CTYPE2']
......
ObservationSim/Config/Header/extension_header.header
View file @ 3226117c
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 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 ================================================================== CAMERA = 'MS' / camera of main survey DETSN = '12345678' / detector serial number DETNAME = 'CCD' / detector type DETTEMP1= 173.0 / detector temperature at EXPSTART(in Kelvin) DETTEMP2= 173.0 / detector temperature at EXPEND(in Kelvin) DETTEMP3= 173.0 / detector temperature at READT1(in Kelvin) DETSIZE = '9560x9264' / detector size DATASECT= '9216x9232' / data section PIXSCAL1= 0.074 / pixel scale for axis 1 PIXSCAL2= 0.074 / pixel scale for axis 2 PIXSIZE1= 10 / pixel size for axis 1 (in um) PIXSIZE2= 10 / pixel size for axis 2 (in um) COMMENT ================================================================== COMMENT CCD chip information COMMENT ================================================================== CHIPID = '08' / chip ID CHIPLAB = 'y-1' / chip label FILTER = 'y' / filter name NCHAN = 16 / number of readout channels PSCAN1 = 27 / horizontal prescan width, per readout channel PSCAN2 = 8 / vertical prescan width, per readout channel OSCAN1 = 16 / horizontal overscan width,per readout channel OSCAN2 = 16 / vertical 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 CD1_2 = 8.17455836176000E-06 / partial of first axis coordinate w.r.t.y CD2_1 = -8.1745583617600E-06 / partial of second axis coordinate w.r.t.x CD2_2 = 1.88602083707394E-05 / partial of second axis coordinate w.r.t.y OTHERS = '' / COMMENT ================================================================== COMMENT Readout information COMMENT ================================================================== GAINLVL = '01' / gain level GAIN01 = 1.1 / gain (channel 01) GAIN02 = 1.1 / gain (channel 02) GAIN03 = 1.1 / gain (channel 03) GAIN04 = 1.1 / gain (channel 04) GAIN05 = 1.1 / gain (channel 05) GAIN06 = 1.1 / gain (channel 06) GAIN07 = 1.1 / gain (channel 07) GAIN08 = 1.1 / gain (channel 08) GAIN09 = 1.1 / gain (channel 09) GAIN10 = 1.1 / gain (channel 10) GAIN11 = 1.1 / gain (channel 11) GAIN12 = 1.1 / gain (channel 12) GAIN13 = 1.1 / gain (channel 13) GAIN14 = 1.1 / gain (channel 14) GAIN15 = 1.1 / gain (channel 15) GAIN16 = 1.1 / gain (channel 16) RON01 = 5.0 / read noise (channel 01) RON02 = 5.0 / read noise (channel 02) RON03 = 5.0 / read noise (channel 03) RON04 = 5.0 / read noise (channel 04) RON05 = 5.0 / read noise (channel 05) RON06 = 5.0 / read noise (channel 06) RON07 = 5.0 / read noise (channel 07) RON08 = 5.0 / read noise (channel 08) RON09 = 5.0 / read noise (channel 09) RON10 = 5.0 / read noise (channel 10) RON11 = 5.0 / read noise (channel 11) RON12 = 5.0 / read noise (channel 12) RON13 = 5.0 / read noise (channel 13) RON14 = 5.0 / read noise (channel 14) RON15 = 5.0 / read noise (channel 15) RON16 = 5.0 / read noise (channel 16) READT0 = '2024-00-00T00:00:00'/ readout start time(UTC) READT1 = '2024-00-00T00:00:00'/ readout end time(UTC) ROSPEED = 10.0 / readout speed (in MHz) EXPTIME = 150.0 / exposure duration DARKTIME= 150.0 / dark current time COMMENT ================================================================== COMMENT Shutter information COMMENT ================================================================== SHTSTAT = 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 ================================================================== LEDFLAG = 0 / main/backup LED LEDSTAT = '00000000000000' / LED status LEDEXPT = 0.0 / LED flash time (s) LEDTEMP = 173.0 / LED temperature (in K) COMMENT ================================================================== COMMENT Other information COMMENT ================================================================== CHECKSUM= '''abcde''' / HDU checksum updated yyyy-mm-ddTHH:MM:SS DATASUM = '''abcde''' / data unit checksum updated yyyy-mm-ddTHH:MM:SS END
ObservationSim/Config/Header/global_header.header
View file @ 3226117c
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 SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T NEXTEND = 1 / number of array dimensions GROUPS = F / ' ' DATE = '2021-03-04T09:30:00'/ the date on which this file was written FILENAME= 'MSC_MS_210304093000_100000000_06_raw' / file name FILETYPE= 'SCIE ' / observation type TELESCOP= 'CSST ' / telescope used to acquire data INSTRUME= 'MSC ' / instrument used to acquire data RADECSYS= 'ICRS ' / reference coordinates system 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 opt axis at CABSTART MOONANG0= 30.0 / angle between moon and opt axis at CABSTART TEL_ALT0= 20.0 / angle between opt axis and Elimb at CABSTART POS_ANG0= 20.0 / angle between y axis and NP at CABSTART 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 opt axis at CABEND MOONANG1= 30.0 / angle between moon and opt axis at CABEND TEL_ALT1= 20.0 / angle between opt axis and Elimb at CABEND POS_ANG1= 20.0 / angle between y axis and NP 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.0 / coordinate epoch COMMENT Other information COMMENT ================================================================== CHECKSUM= 'abcdefg ' / HDU checksum updated yyyy-mm-ddTHH:MM:SS DATASUM = 'abcdefg ' / data unit checksum updated yyyy-mm-ddTHH:MM:SS END
ObservationSim/Config/Pointing.py
View file @ 3226117c
...@@ -3,16 +3,18 @@ import galsim ...@@ -3,16 +3,18 @@ import galsim
from astropy.time import Time from astropy.time import Time
class Pointing(object): class Pointing(object):
def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='MS'): def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='MS'):
self.id = id self.id = id
self.ra = ra self.ra = ra
self.dec = dec self.dec = dec
self.img_pa = img_pa * galsim.degrees self.img_pa = img_pa * galsim.degrees
self.timestamp = timestamp self.timestamp = timestamp
self.sat_x, self.sat_y, self.sat_z = sat_x, sat_y, sat_z self.sat_x, self.sat_y, self.sat_z = sat_x, sat_y, sat_z
self.sun_x, self.sun_y, self.sun_z = sun_x, sun_y, sun_z
self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz
self.exp_time = exp_time self.exp_time = exp_time
self.pointing_type = pointing_type self.pointing_type = pointing_type
self.jdt = 0.
def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'): def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'):
self.id = id self.id = id
...@@ -24,10 +26,14 @@ class Pointing(object): ...@@ -24,10 +26,14 @@ class Pointing(object):
if col_len > 5: if col_len > 5:
jdt = np.double(columns[5]) jdt = np.double(columns[5])
t_temp = Time(jdt, format='jd') t_temp = Time(jdt, format='jd')
self.jdt = jdt
self.timestamp = t_temp.unix self.timestamp = t_temp.unix
self.sat_x = float(columns[6]) self.sat_x = float(columns[6])
self.sat_y = float(columns[7]) self.sat_y = float(columns[7])
self.sat_z = float(columns[8]) self.sat_z = float(columns[8])
self.sun_x = float(columns[9])
self.sun_y = float(columns[10])
self.sun_z = float(columns[1])
self.sat_vx = float(columns[15]) self.sat_vx = float(columns[15])
self.sat_vy = float(columns[16]) self.sat_vy = float(columns[16])
self.sat_vz = float(columns[17]) self.sat_vz = float(columns[17])
......
ObservationSim/Instrument/Chip/Chip.py
View file @ 3226117c
...@@ -27,7 +27,9 @@ class Chip(FocalPlane): ...@@ -27,7 +27,9 @@ class Chip(FocalPlane):
# self.npix_x = 9216 # self.npix_x = 9216
# self.npix_y = 9232 # self.npix_y = 9232
# self.pix_scale = 0.074 # pixel scale # self.pix_scale = 0.074 # pixel scale
self.nsecy = 2
self.nsecx = 8
self.gain_channel = np.ones(self.nsecy* self.nsecx)
self._set_attributes_from_config(config) self._set_attributes_from_config(config)
self.logger = logger self.logger = logger
...@@ -75,7 +77,8 @@ class Chip(FocalPlane): ...@@ -75,7 +77,8 @@ class Chip(FocalPlane):
self.fdModel = None self.fdModel = None
else: else:
try: try:
with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion: # with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModelGlobal_mainFP_v1.0.pickle") as field_distortion:
with pkg_resources.files('ObservationSim.Instrument.data.field_distortion').joinpath("FieldDistModel_v2.0.pickle") as field_distortion:
with open(field_distortion, "rb") as f: with open(field_distortion, "rb") as f:
self.fdModel = pickle.load(f) self.fdModel = pickle.load(f)
except AttributeError: except AttributeError:
...@@ -115,7 +118,7 @@ class Chip(FocalPlane): ...@@ -115,7 +118,7 @@ class Chip(FocalPlane):
self._getCRdata() self._getCRdata()
# Define the sensor model # Define the sensor model
if config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric": if "bright_fatter" in config["ins_effects"] and config["ins_effects"]["bright_fatter"] == True and self.survey_type == "photometric":
self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func) self.sensor = galsim.SiliconSensor(strength=self.df_strength, treering_func=treering_func)
else: else:
self.sensor = galsim.Sensor() self.sensor = galsim.Sensor()
...@@ -144,18 +147,18 @@ class Chip(FocalPlane): ...@@ -144,18 +147,18 @@ class Chip(FocalPlane):
def _getSurveyType(self): def _getSurveyType(self):
if self.filter_type in ["GI", "GV", "GU"]: if self.filter_type in ["GI", "GV", "GU"]:
return "spectroscopic" return "spectroscopic"
elif self.filter_type in ["nuv", "u", "g", 'r', 'i', 'z', 'y', 'FGS']: elif self.filter_type in ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
return "photometric" return "photometric"
# elif self.filter_type in ["FGS"]: # elif self.filter_type in ["FGS"]:
# return "FGS" # return "FGS"
def _getChipEffCurve(self, filter_type): def _getChipEffCurve(self, filter_type):
# CCD efficiency curves # CCD efficiency curves
if filter_type in ['nuv', 'u', 'GU']: filename = 'UV0.txt' if filter_type in ['NUV', 'u', 'GU']: filename = 'UV0.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 ['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' if filter_type in ['i', 'z', 'y', 'GI']: filename = 'Basic_NIR.txt'
# Mirror efficiency: # Mirror efficiency:
# if filter_type == 'nuv': mirror_eff = 0.54 # if filter_type == 'NUV': mirror_eff = 0.54
# if filter_type == 'u': mirror_eff = 0.68 # 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 ['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 # if filter_type in ['GU', 'GV', 'GI']: mirror_eff = 1. # Not sure if this is right
...@@ -182,7 +185,7 @@ class Chip(FocalPlane): ...@@ -182,7 +185,7 @@ class Chip(FocalPlane):
def getChipFilter(self, chipID=None, filter_layout=None): def getChipFilter(self, chipID=None, filter_layout=None):
"""Return the filter index and type for a given chip #(chipID) """Return the filter index and type for a given chip #(chipID)
""" """
filter_type_list = ["nuv","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"] filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
if filter_layout is not None: if filter_layout is not None:
return filter_layout[chipID][0], filter_layout[chipID][1] return filter_layout[chipID][0], filter_layout[chipID][1]
if chipID == None: if chipID == None:
...@@ -195,7 +198,7 @@ class Chip(FocalPlane): ...@@ -195,7 +198,7 @@ class Chip(FocalPlane):
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 [14, 17]: filter_type = "u"
if chipID in [9, 22]: filter_type = "r" if chipID in [9, 22]: filter_type = "r"
if chipID in [12, 13, 18, 19]: filter_type = "nuv" if chipID in [12, 13, 18, 19]: filter_type = "NUV"
if chipID in [8, 23]: filter_type = "g" if chipID in [8, 23]: filter_type = "g"
if chipID in [1, 10, 21, 30]: filter_type = "GI" if chipID in [1, 10, 21, 30]: filter_type = "GI"
if chipID in [2, 5, 26, 29]: filter_type = "GV" if chipID in [2, 5, 26, 29]: filter_type = "GV"
...@@ -304,7 +307,10 @@ class Chip(FocalPlane): ...@@ -304,7 +307,10 @@ class Chip(FocalPlane):
if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf'] if self.chipID == 30: confFile = ['CSST_GI8.conf', 'CSST_GI7.conf']
return confFile return confFile
def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, exptime=150.): def generateHeader(self, ra_cen, dec_cen, img_rot, im_type, pointing_ID, exptime=150., timestamp = 1621915200):
datetime_obs = datetime.utcfromtimestamp(timestamp)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
h_prim = generatePrimaryHeader( h_prim = generatePrimaryHeader(
xlen=self.npix_x, xlen=self.npix_x,
ylen=self.npix_y, ylen=self.npix_y,
...@@ -333,23 +339,30 @@ class Chip(FocalPlane): ...@@ -333,23 +339,30 @@ class Chip(FocalPlane):
pixel_size=self.pix_size, pixel_size=self.pix_size,
xcen=self.x_cen, xcen=self.x_cen,
ycen=self.y_cen, ycen=self.y_cen,
extName='raw') extName='SCI',
timestamp = timestamp,
exptime = exptime,
readoutTime = 40.)
return h_prim, h_ext return h_prim, h_ext
def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, date_obs, time_obs, output_dir, exptime=150.): def outputCal(self, img, ra_cen, dec_cen, img_rot, im_type, pointing_ID, output_dir, exptime=150., timestamp = 1621915200):
h_prim, h_ext = self.generateHeader( h_prim, h_ext = self.generateHeader(
ra_cen=ra_cen, ra_cen=ra_cen,
dec_cen=dec_cen, dec_cen=dec_cen,
img_rot=img_rot, img_rot=img_rot,
im_type=im_type, im_type=im_type,
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs, exptime=exptime,
time_obs=time_obs, timestamp = timestamp)
exptime=exptime)
hdu1 = fits.PrimaryHDU(header=h_prim) hdu1 = fits.PrimaryHDU(header=h_prim)
hdu1.add_checksum() hdu1.add_checksum()
hdu1.header.comments['CHECKSUM'] = 'HDU checksum'
hdu1.header.comments['DATASUM'] = 'data unit checksum'
hdu2 = fits.ImageHDU(img.array, header=h_ext) hdu2 = fits.ImageHDU(img.array, header=h_ext)
hdu2.add_checksum() hdu2.add_checksum()
hdu2.header.comments['XTENSION'] = 'extension type'
hdu2.header.comments['CHECKSUM'] = 'HDU checksum'
hdu2.header.comments['DATASUM'] = 'data unit checksum'
hdu1 = fits.HDUList([hdu1, hdu2]) hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits') fname = os.path.join(output_dir, h_prim['FILENAME']+'.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True) hdu1.writeto(fname, output_verify='ignore', overwrite=True)
...@@ -456,9 +469,9 @@ class Chip(FocalPlane): ...@@ -456,9 +469,9 @@ class Chip(FocalPlane):
del cr_map del cr_map
# crmap_gsimg.write("%s/CosmicRay_%s_1.fits" % (chip_output.subdir, self.chipID)) # 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)) # crmap_gsimg.write("%s/CosmicRay_%s.fits" % (chip_output.subdir, self.chipID))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal( self.outputCal(
img=crmap_gsimg, img=crmap_gsimg,
ra_cen=ra_cen, ra_cen=ra_cen,
...@@ -466,10 +479,9 @@ class Chip(FocalPlane): ...@@ -466,10 +479,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='CRS', im_type='CRS',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=exptime) exptime=exptime,
timestamp=timestamp_obs)
del crmap_gsimg del crmap_gsimg
# Apply PRNU effect and output PRNU flat file: # Apply PRNU effect and output PRNU flat file:
...@@ -503,21 +515,6 @@ class Chip(FocalPlane): ...@@ -503,21 +515,6 @@ class Chip(FocalPlane):
if config["ins_effects"]["add_badcolumns"] == True: if config["ins_effects"]["add_badcolumns"] == True:
img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) img = effects.BadColumns(img, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger)
# Add Bias level
if config["ins_effects"]["add_bias"] == True:
if self.logger is not None:
self.logger.info(" Adding Bias level and 16-channel non-uniformity")
else:
print(" Adding Bias level and 16-channel non-uniformity")
if config["ins_effects"]["bias_16channel"] == True:
img = effects.AddBiasNonUniform16(img,
bias_level=float(self.bias_level),
nsecy = 2, nsecx=8,
seed=SeedBiasNonuni+self.chipID,
logger=self.logger)
elif config["ins_effects"]["bias_16channel"] == False:
img += self.bias_level
# Apply Nonlinearity on the chip image # Apply Nonlinearity on the chip image
if config["ins_effects"]["non_linear"] == True: if config["ins_effects"]["non_linear"] == True:
if self.logger is not None: if self.logger is not None:
...@@ -542,6 +539,21 @@ class Chip(FocalPlane): ...@@ -542,6 +539,21 @@ class Chip(FocalPlane):
print(" Apply CTE Effect") print(" Apply CTE Effect")
img = effects.CTE_Effect(GSImage=img, threshold=27) img = effects.CTE_Effect(GSImage=img, threshold=27)
# Add Bias level
if config["ins_effects"]["add_bias"] == True:
if self.logger is not None:
self.logger.info(" Adding Bias level and 16-channel non-uniformity")
else:
print(" Adding Bias level and 16-channel non-uniformity")
if config["ins_effects"]["bias_16channel"] == True:
img = effects.AddBiasNonUniform16(img,
bias_level=float(self.bias_level),
nsecy = 2, nsecx=8,
seed=SeedBiasNonuni+self.chipID,
logger=self.logger)
elif config["ins_effects"]["bias_16channel"] == False:
img += self.bias_level
# Add Read-out Noise # Add Read-out Noise
if config["ins_effects"]["add_readout"] == True: if config["ins_effects"]["add_readout"] == True:
seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID seed = int(config["random_seeds"]["seed_readout"]) + pointing_ID*30 + self.chipID
...@@ -555,7 +567,7 @@ class Chip(FocalPlane): ...@@ -555,7 +567,7 @@ class Chip(FocalPlane):
else: else:
print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True) print(" Applying Gain (and 16 channel non-uniformity) & Quantization", flush=True)
if config["ins_effects"]["gain_16channel"] == True: if config["ins_effects"]["gain_16channel"] == True:
img = effects.ApplyGainNonUniform16( img, self.gain_channel = effects.ApplyGainNonUniform16(
img, gain=self.gain, img, gain=self.gain,
nsecy = 2, nsecx=8, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
...@@ -610,7 +622,7 @@ class Chip(FocalPlane): ...@@ -610,7 +622,7 @@ class Chip(FocalPlane):
if config["ins_effects"]["add_badcolumns"] == True: if config["ins_effects"]["add_badcolumns"] == True:
BiasCombImg = effects.BadColumns(BiasCombImg-float(self.bias_level)+5, seed=SeedBadColumns, chipid=self.chipID, logger=self.logger) + float(self.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, BiasCombImg, self.gain_channel = effects.ApplyGainNonUniform16(BiasCombImg, gain=self.gain,
nsecy = 2, nsecx=8, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
logger=self.logger) logger=self.logger)
...@@ -622,9 +634,9 @@ class Chip(FocalPlane): ...@@ -622,9 +634,9 @@ class Chip(FocalPlane):
BiasCombImg.quantize() BiasCombImg.quantize()
BiasCombImg = galsim.ImageUS(BiasCombImg) BiasCombImg = galsim.ImageUS(BiasCombImg)
# BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1)) # BiasCombImg.write("%s/BiasImg_%s_%s_%s.fits" % (chip_output.subdir, BiasTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60 timestamp_obs += 10 * 60
self.outputCal( self.outputCal(
img=BiasCombImg, img=BiasCombImg,
...@@ -633,10 +645,9 @@ class Chip(FocalPlane): ...@@ -633,10 +645,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='BIAS', im_type='BIAS',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=0.0) exptime=0.0,
timestamp=timestamp_obs)
del BiasCombImg del BiasCombImg
# Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file # Export combined (ncombine, Vignetting + PRNU) & single vignetting flat-field file
...@@ -717,7 +728,7 @@ class Chip(FocalPlane): ...@@ -717,7 +728,7 @@ class Chip(FocalPlane):
readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
FlatCombImg.addNoise(readout_noise) FlatCombImg.addNoise(readout_noise)
FlatCombImg = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain, FlatCombImg, self.gain_channel = effects.ApplyGainNonUniform16(FlatCombImg, gain=self.gain,
nsecy = 2, nsecx=8, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
logger=self.logger) logger=self.logger)
...@@ -726,9 +737,9 @@ class Chip(FocalPlane): ...@@ -726,9 +737,9 @@ class Chip(FocalPlane):
FlatCombImg.quantize() FlatCombImg.quantize()
FlatCombImg = galsim.ImageUS(FlatCombImg) FlatCombImg = galsim.ImageUS(FlatCombImg)
# FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1)) # FlatCombImg.write("%s/FlatImg_%s_%s_%s.fits" % (chip_output.subdir, FlatTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60 timestamp_obs += 10 * 60
self.outputCal( self.outputCal(
img=FlatCombImg, img=FlatCombImg,
...@@ -737,10 +748,9 @@ class Chip(FocalPlane): ...@@ -737,10 +748,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='FLAT', im_type='FLAT',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=self.flat_exptime) exptime=self.flat_exptime,
timestamp=timestamp_obs)
del FlatCombImg, FlatSingle, prnu_img del FlatCombImg, FlatSingle, prnu_img
# flat_img.replaceNegative(replace_value=0) # flat_img.replaceNegative(replace_value=0)
...@@ -783,9 +793,9 @@ class Chip(FocalPlane): ...@@ -783,9 +793,9 @@ class Chip(FocalPlane):
cr_map[cr_map < 0] = 0 cr_map[cr_map < 0] = 0
crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16) crmap_gsimg = galsim.Image(cr_map, dtype=np.uint16)
del cr_map del cr_map
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
self.outputCal( self.outputCal(
img=crmap_gsimg, img=crmap_gsimg,
ra_cen=ra_cen, ra_cen=ra_cen,
...@@ -793,10 +803,9 @@ class Chip(FocalPlane): ...@@ -793,10 +803,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='CRD', im_type='CRD',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=self.dark_exptime) exptime=self.dark_exptime,
timestamp=timestamp_obs)
del crmap_gsimg del crmap_gsimg
# Non-Linearity for Dark # Non-Linearity for Dark
...@@ -839,7 +848,7 @@ class Chip(FocalPlane): ...@@ -839,7 +848,7 @@ class Chip(FocalPlane):
readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise) readout_noise = galsim.GaussianNoise(rng=rng_readout, sigma=self.read_noise)
DarkCombImg.addNoise(readout_noise) DarkCombImg.addNoise(readout_noise)
DarkCombImg = effects.ApplyGainNonUniform16( DarkCombImg, self.gain_channel = effects.ApplyGainNonUniform16(
DarkCombImg, gain=self.gain, DarkCombImg, gain=self.gain,
nsecy = 2, nsecx=8, nsecy = 2, nsecx=8,
seed=SeedGainNonuni+self.chipID, seed=SeedGainNonuni+self.chipID,
...@@ -852,9 +861,9 @@ class Chip(FocalPlane): ...@@ -852,9 +861,9 @@ class Chip(FocalPlane):
DarkCombImg.quantize() DarkCombImg.quantize()
DarkCombImg = galsim.ImageUS(DarkCombImg) DarkCombImg = galsim.ImageUS(DarkCombImg)
# DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1)) # DarkCombImg.write("%s/DarkImg_%s_%s_%s.fits" % (chip_output.subdir, DarkTag, self.chipID, i+1))
datetime_obs = datetime.utcfromtimestamp(timestamp_obs) # datetime_obs = datetime.utcfromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d") # date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S") # time_obs = datetime_obs.strftime("%H%M%S")
timestamp_obs += 10 * 60 timestamp_obs += 10 * 60
self.outputCal( self.outputCal(
img=DarkCombImg, img=DarkCombImg,
...@@ -863,10 +872,9 @@ class Chip(FocalPlane): ...@@ -863,10 +872,9 @@ class Chip(FocalPlane):
img_rot=img_rot, img_rot=img_rot,
im_type='DARK', im_type='DARK',
pointing_ID=pointing_ID, pointing_ID=pointing_ID,
date_obs=date_obs,
time_obs=time_obs,
output_dir=chip_output.subdir, output_dir=chip_output.subdir,
exptime=self.dark_exptime) exptime=self.dark_exptime,
timestamp = timestamp_obs)
del DarkCombImg del DarkCombImg
# img = galsim.ImageUS(img) # img = galsim.ImageUS(img)
......
ObservationSim/Instrument/Chip/Effects.py
View file @ 3226117c
...@@ -152,6 +152,7 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg ...@@ -152,6 +152,7 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg
rg = Generator(PCG64(int(seed))) rg = Generator(PCG64(int(seed)))
Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1 # sigma~1% Random16 = (rg.random(nsecy*nsecx)-0.5)*0.04+1 # sigma~1%
Gain16 = Random16.reshape((nsecy,nsecx))/gain Gain16 = Random16.reshape((nsecy,nsecx))/gain
gain_array = np.ones(nsecy*nsecx)*gain
if logger is not None: if logger is not None:
msg = str("Gain of 16 channels: " + str(Gain16)) msg = str("Gain of 16 channels: " + str(Gain16))
logger.info(msg) logger.info(msg)
...@@ -163,7 +164,8 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg ...@@ -163,7 +164,8 @@ def ApplyGainNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logg
for rowi in range(nsecy): for rowi in range(nsecy):
for coli in range(nsecx): for coli in range(nsecx):
GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli] GSImage.array[rowi*secsize_y:(rowi+1)*secsize_y,coli*secsize_x:(coli+1)*secsize_x] *= Gain16[rowi,coli]
return GSImage gain_array[rowi*nsecx+coli] = 1/Gain16[rowi,coli]
return GSImage, gain_array
def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None): def GainsNonUniform16(GSImage, gain=1, nsecy = 2, nsecx=8, seed=202102, logger=None):
......
ObservationSim/Instrument/Filter.py
View file @ 3226117c
...@@ -2,9 +2,11 @@ import galsim ...@@ -2,9 +2,11 @@ import galsim
import pylab as pl import pylab as pl
import os import os
import numpy as np import numpy as np
import gc
from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag from ObservationSim.Instrument._util import photonEnergy, calculateLimitMag
from ObservationSim.Instrument.FilterParam import FilterParam from ObservationSim.Instrument.FilterParam import FilterParam
from ObservationSim.Straylight import Straylight
try: try:
import importlib.resources as pkg_resources import importlib.resources as pkg_resources
...@@ -38,6 +40,7 @@ class Filter(object): ...@@ -38,6 +40,7 @@ class Filter(object):
self.sky_background = filter_param.param[filter_type][5] self.sky_background = filter_param.param[filter_type][5]
self.mag_saturation = filter_param.param[filter_type][6] self.mag_saturation = filter_param.param[filter_type][6]
self.mag_limiting = filter_param.param[filter_type][7] self.mag_limiting = filter_param.param[filter_type][7]
self.zodical_spec = None
def is_too_bright(self, mag, margin=-2.5): def is_too_bright(self, mag, margin=-2.5):
return mag <= self.mag_saturation + margin return mag <= self.mag_saturation + margin
...@@ -108,6 +111,24 @@ class Filter(object): ...@@ -108,6 +111,24 @@ class Filter(object):
def getSkyNoise(self, exptime, gain=1.): def getSkyNoise(self, exptime, gain=1.):
return self.sky_background * exptime / gain return self.sky_background * exptime / gain
def setFilterStrayLightPixel(self,jtime = 2460843., sat_pos = np.array([0,0,0]), pointing_radec = np.array([0,0]), sun_pos = np.array([0,0,0])):
sl = Straylight(jtime=jtime, sat_pos=sat_pos, pointing_radec=pointing_radec,sun_pos=sun_pos)
if self.filter_type in ["GU","GV","GI"]:
s_pix, spec = sl.calculateStrayLightGrating(grating = self.filter_type.upper())
if s_pix>0.8:
s_pix = 0.8
self.sky_background = s_pix
self.zodical_spec = spec
elif self.filter_type.lower() in ["nuv","u","g","r","i","z","y"]:
s_pix = sl.calculateStrayLightFilter(filter=self.filter_type.lower())
if s_pix>1:
s_pix = 1
self.sky_background = s_pix
self.zodical_spec = None
del sl
gc.collect()
def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None): def update_limit_saturation_mags(self, exptime=150., psf_fwhm=0.1969, skyFn='sky_emiss_hubble_50_50_A.dat', chip=None):
if self.filter_type in ["GI", "GV", "GU"]: if self.filter_type in ["GI", "GV", "GU"]:
return return
......
ObservationSim/Instrument/FilterParam.py
View file @ 3226117c
...@@ -26,7 +26,7 @@ class FilterParam(object): ...@@ -26,7 +26,7 @@ class FilterParam(object):
# 8) dim end magnitude # 8) dim end magnitude
if filter_param == None: if filter_param == None:
filtP = { filtP = {
"nuv": [2867.7, 705.4, 2470.0, 3270.0, 0.1404, 0.004, 15.7, 25.4], "NUV": [2867.7, 705.4, 2470.0, 3270.0, 0.1404, 0.004, 15.7, 25.4],
"u": [3601.1, 852.1, 3120.0, 4090.0, 0.2176, 0.021, 16.1, 25.4], "u": [3601.1, 852.1, 3120.0, 4090.0, 0.2176, 0.021, 16.1, 25.4],
"g": [4754.5, 1569.8, 3900.0, 5620.0, 0.4640, 0.164, 17.2, 26.3], "g": [4754.5, 1569.8, 3900.0, 5620.0, 0.4640, 0.164, 17.2, 26.3],
"r": [6199.8, 1481.2, 5370.0, 7030.0, 0.5040, 0.207, 17.0, 26.0], "r": [6199.8, 1481.2, 5370.0, 7030.0, 0.5040, 0.207, 17.0, 26.0],
......
ObservationSim/Instrument/FocalPlane.py
View file @ 3226117c
...@@ -86,15 +86,16 @@ class FocalPlane(object): ...@@ -86,15 +86,16 @@ class FocalPlane(object):
if (xcen == None) or (ycen == None): if (xcen == None) or (ycen == None):
xcen = self.cen_pix_x xcen = self.cen_pix_x
ycen = self.cen_pix_y ycen = self.cen_pix_y
# dudx = -np.cos(img_rot.rad) * pix_scale
# dudy = -np.sin(img_rot.rad) * pix_scale
# dvdx = -np.sin(img_rot.rad) * pix_scale
# dvdy = +np.cos(img_rot.rad) * pix_scale
dudx = -np.cos(img_rot.rad) * pix_scale dudx = -np.cos(img_rot.rad) * pix_scale
dudy = -np.sin(img_rot.rad) * pix_scale dudy = +np.sin(img_rot.rad) * pix_scale
dvdx = -np.sin(img_rot.rad) * pix_scale dvdx = -np.sin(img_rot.rad) * pix_scale
dvdy = +np.cos(img_rot.rad) * pix_scale dvdy = -np.cos(img_rot.rad) * pix_scale
# dudx = +np.sin(img_rot.rad) * pix_scale
# dudy = +np.cos(img_rot.rad) * pix_scale
# dvdx = -np.cos(img_rot.rad) * pix_scale
# dvdy = +np.sin(img_rot.rad) * pix_scale
moscen = galsim.PositionD(x=xcen, y=ycen) moscen = galsim.PositionD(x=xcen, y=ycen)
sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees) sky_center = galsim.CelestialCoord(ra=ra*galsim.degrees, dec=dec*galsim.degrees)
affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen) affine = galsim.AffineTransform(dudx, dudy, dvdx, dvdy, origin=moscen)
......
ObservationSim/Instrument/data/ccd/chip_definition.json
View file @ 3226117c
...@@ -468,7 +468,7 @@ ...@@ -468,7 +468,7 @@
"full_well": 90000 "full_well": 90000
}, },
"15": { "15": {
"chip_name": "y-4", "chip_name": "y-2",
"pix_size": 0.01, "pix_size": 0.01,
"pix_scale": 0.074, "pix_scale": 0.074,
"npix_x": 9216, "npix_x": 9216,
...@@ -648,7 +648,7 @@ ...@@ -648,7 +648,7 @@
"full_well": 90000 "full_well": 90000
}, },
"25": { "25": {
"chip_name": "y-2", "chip_name": "y-4",
"pix_size": 0.01, "pix_size": 0.01,
"pix_scale": 0.074, "pix_scale": 0.074,
"npix_x": 9216, "npix_x": 9216,
......
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