update fits module (d4b41c36) · Commits · csst-sims / csst_msc_sim

Catalog/C6_fits_Catalog.py

+17 −9

Original line number	Diff line number	Diff line
		@@ -60,6 +60,9 @@ class Catalog(CatalogBase):

		with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
		self.normF_star = Table.read(str(filter_path))
		with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path:
		self.normF_galaxy = Table.read(str(filter_path))


		self.config = config
		self.chip = chip
		@@ -138,8 +141,8 @@ class Catalog(CatalogBase):
		return None
		###mock_stamp_START
		elif obj.type == "stamp":
		#return self.normF_galaxy ###normalize_filter for stamp
		return None
		return self.normF_galaxy ###normalize_filter for stamp
		#return None
		###mock_stamp_END
		else:
		return None
		@@ -197,8 +200,8 @@ class Catalog(CatalogBase):

		for igals in range(ngals):
		# # (TEST)
		# if igals > 100:
		# break
		if igals > 2000:
		break

		param = self.initialize_param()
		param['ra'] = ra_arr[igals]
		@@ -321,8 +324,8 @@ class Catalog(CatalogBase):
		)
		for istars in range(nstars):
		# # (TEST)
		# if istars > 100:
		# break
		if istars > 100:
		break

		param = self.initialize_param()
		param['ra'] = ra_arr[istars]
		@@ -389,6 +392,9 @@ class Catalog(CatalogBase):
		input_time_str=time_str
		)
		for iAGNs in range(nAGNs):
		if iAGNs > 100:
		break

		param = self.initialize_param()
		param['ra'] = ra_arr[iAGNs]
		param['dec'] = dec_arr[iAGNs]
		@@ -425,13 +431,15 @@ class Catalog(CatalogBase):
		self.ud = galsim.UniformDeviate(pix_id)

		for istamp in range(nstamps):
		print('DEBUG:::istamp=', istamp)

		fitsfile = os.path.join(self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8'))
		hdu=fitsio.open(fitsfile)

		param = self.initialize_param()
		param['id'] = hdu[0].header['index'] #istamp
		param['star'] = 3 # Stamp type in .cat file
		param['lensGalaxyID'] = hdu[0].header['lensGID']
		###param['lensGalaxyID'] = hdu[0].header['lensGID']
		param['ra'] = hdu[0].header['ra']
		param['dec']= hdu[0].header['dec']
		param['pixScale']= hdu[0].header['pixScale']
		@@ -440,9 +448,9 @@ class Catalog(CatalogBase):
		#param['PA']= hdu[0].header['PA']
		#param['bfrac']= hdu[0].header['bfrac']
		#param['z']= hdu[0].header['z']
		param['mag_use_normal'] = 22 #hdu[0].header['m_normal'] #gals['mag_true_g_lsst']
		param['mag_use_normal'] = 20 #hdu[0].header['m_normal'] #gals['mag_true_g_lsst']

		assert(stamps['lensGID'][istamp] == param['lensGalaxyID'])
		###assert(stamps['lensGID'][istamp] == param['lensGalaxyID'])

		# Apply astrometric modeling
		# in C3 case only aberration

ObservationSim/MockObject/Stamp.py

+213 −53

Original line number	Diff line number	Diff line
		import os, sys
		import random
		import numpy as np
		import galsim
		import astropy.constants as cons
		from astropy.table import Table
		from scipy import interpolate

		import astropy.io.fits as fitsio

		import galsim
		import gc
		from ObservationSim.MockObject.MockObject import MockObject

		from ObservationSim.MockObject._util import magToFlux,VC_A
		from ObservationSim.MockObject.SpecDisperser import SpecDisperser
		from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders

		class Stamp(MockObject):
		def __init__(self, param):
		super().__init__(param)
		def __init__(self, param, logger=None):
		super().__init__(param, logger=logger)

		def unload_SED(self):
		"""(Test) free up SED memory
		@@ -35,9 +30,9 @@ class Stamp(MockObject):
		self.logger.error(e)
		return False

		nphotons_sum = 0
		photons_list = []
		xmax, ymax = 0, 0
		#nphotons_sum = 0
		#photons_list = []
		#xmax, ymax = 0, 0

		if self.getMagFilter(filt) <= 15:
		folding_threshold = 5.e-4
		@@ -46,7 +41,7 @@ class Stamp(MockObject):
		gsp = galsim.GSParams(folding_threshold=folding_threshold)

		self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
		img_real_wcs=self.real_wcs)
		img_real_wcs=self.chip_wcs)

		x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
		x_nominal = int(np.floor(x + 0.5))
		@@ -55,72 +50,237 @@ class Stamp(MockObject):
		dy = y - y_nominal
		offset = galsim.PositionD(dx, dy)

		real_wcs_local = self.real_wcs.local(self.real_pos)
		chip_wcs_local = self.chip_wcs.local(self.real_pos)
		is_updated = 0

		if fd_shear:
		g1 += fd_shear.g1
		g2 += fd_shear.g2
		gal_shear = galsim.Shear(g1=g1, g2=g2)

		for i in range(len(bandpass_list)):
		bandpass = bandpass_list[i]

		try:
		sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
		except Exception as e:
		print(e)
		if self.logger:
		self.logger.error(e)
		# return False
		continue

		ratio = sub/full
		if not (ratio == -1 or (ratio != ratio)):
		nphotons = ratio * nphotons_tot
		else:
		# return False
		continue
		nphotons_sum += nphotons
		#nphotons_sum += nphotons

		psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)

		_gal = self.param['image']
		galIm = galsim.ImageF(_gal, scale=self.param['pixScale'])
		gal = galsim.InterpolatedImage(galIm)
		gal = gal.withFlux(nphotons)
		#gal_shear = galsim.Shear(g1=g1, g2=g2)
		#gal = gal.shear(gal_shear)

		gal = galsim.Convolve(psf, gal)
		if fd_shear is not None:
		gal = gal.shear(fd_shear)
		galImg= galsim.ImageF(_gal, scale=self.param['pixScale'])
		gal_temp= galsim.InterpolatedImage(galImg)
		gal_temp= gal_temp.shear(gal_shear)
		gal_temp= gal_temp.withFlux(nphotons)

		stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=self.offset, save_photons=True)
		gal_temp= galsim.Convolve(psf, gal_temp)

		xmax = max(xmax, stamp.xmax - stamp.xmin)
		ymax = max(ymax, stamp.ymax - stamp.ymin)

		photons = stamp.photons
		photons.x += x_nominal
		photons.y += y_nominal
		photons_list.append(photons)
		del gal
		if i == 0:
		gal = gal_temp
		else:
		gal = gal + gal_temp

		# print('xmax = %d, ymax = %d '%(xmax, ymax))
		stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset)
		if np.sum(np.isnan(stamp.array)) > 0:
		# ERROR happens
		return 2, pos_shear

		stamp = galsim.ImageF(int(xmax1.1), int(ymax1.1))
		stamp.wcs = real_wcs_local
		stamp.setCenter(x_nominal, y_nominal)
		bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x - 1, 0, chip.npix_y - 1)

		if bounds.area() > 0:
		chip.img.setOrigin(0, 0)
		stamp[bounds] = chip.img[bounds]
		for i in range(len(photons_list)):
		if i == 0:
		chip.sensor.accumulate(photons_list[i], stamp)
		chip.img[bounds] += stamp[bounds]
		is_updated = 1
		chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
		del stamp

		if is_updated == 0:
		print("fits obj %s missed"%(self.id))
		if self.logger:
		self.logger.info("fits obj %s missed"%(self.id))
		return 0, pos_shear

		return 1, pos_shear


		def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
		exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None):
		if normFilter is not None:
		norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001
		sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed,
		norm_thr=normFilter,
		sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]),
		eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0]))
		if sedNormFactor == 0:
		return 2, None
		else:
		chip.sensor.accumulate(photons_list[i], stamp, resume=True)
		sedNormFactor = 1.
		normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
		names=('WAVELENGTH', 'FLUX'))

		chip.img[bounds] = stamp[bounds]
		self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
		img_real_wcs=self.chip_wcs)

		chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
		x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5
		x_nominal = int(np.floor(x + 0.5))
		y_nominal = int(np.floor(y + 0.5))
		dx = x - x_nominal
		dy = y - y_nominal
		offset = galsim.PositionD(dx, dy)

		del photons_list
		del stamp
		gc.collect()
		return True, pos_shear
		chip_wcs_local = self.chip_wcs.local(self.real_pos)


		if self.getMagFilter(filt) <= 15:
		folding_threshold = 5.e-4
		else:
		folding_threshold = 5.e-3
		gsp = galsim.GSParams(folding_threshold=folding_threshold)
		# nphotons_sum = 0

		flat_cube = chip.flat_cube

		xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388}
		grating_split_pos_chip = 0 + grating_split_pos

		branges = np.zeros([len(bandpass_list), 2])

		# print(hasattr(psf_model, 'bandranges'))

		if hasattr(psf_model, 'bandranges'):
		if psf_model.bandranges is None:
		return 2, None
		if len(psf_model.bandranges) != len(bandpass_list):
		return 2, None
		branges = psf_model.bandranges
		else:
		for i in range(len(bandpass_list)):
		branges[i, 0] = bandpass_list[i].blue_limit * 10
		branges[i, 1] = bandpass_list[i].red_limit * 10

		for i in range(len(bandpass_list)):
		# bandpass = bandpass_list[i]
		brange = branges[i]

		# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold)

		_gal = self.param['image']
		galImg= galsim.ImageF(_gal, scale=self.param['pixScale'])
		gal = galsim.InterpolatedImage(galImg)

		# (TEST) Random knots
		# knots = galsim.RandomKnots(npoints=100, profile=disk)
		# kfrac = np.random.random()*(1.0 - self.bfrac)
		# gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots

		gal = gal.withFlux(tel.pupil_area * exptime)
		if fd_shear:
		g1 += fd_shear.g1
		g2 += fd_shear.g2
		gal_shear = galsim.Shear(g1=g1, g2=g2)
		gal = gal.shear(gal_shear)
		# gal = galsim.Convolve(psf, gal)

		# if not big_galaxy: # Not apply PSF for very big galaxy
		# gal = galsim.Convolve(psf, gal)
		# # if fd_shear is not None:
		# # gal = gal.shear(fd_shear)

		starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space')

		origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
		x_nominal - (starImg.center.x - starImg.xmin)]
		starImg.setOrigin(0, 0)
		gal_origin = [origin_star[0], origin_star[1]]
		gal_end = [origin_star[0] + starImg.array.shape[0] - 1, origin_star[1] + starImg.array.shape[1] - 1]

		if gal_origin[1] < grating_split_pos_chip < gal_end[1]:
		subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1)
		## part img disperse

		subImg_p1 = starImg.array[:, 0:subSlitPos]
		star_p1 = galsim.Image(subImg_p1)
		star_p1.setOrigin(0, 0)
		origin_p1 = origin_star
		xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0
		ycenter_p1 = y_nominal-0

		sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1,
		ycenter=ycenter_p1, origin=origin_p1,
		tar_spec=normalSED,
		band_start=brange[0], band_end=brange[1],
		conf=chip.sls_conf[0],
		isAlongY=0,
		flat_cube=flat_cube)

		# self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
		pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1],
		psf_model=psf_model, bandNo=i + 1,
		grating_split_pos=grating_split_pos,
		local_wcs=chip_wcs_local, pos_img = pos_img)

		subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
		star_p2 = galsim.Image(subImg_p2)
		star_p2.setOrigin(0, 0)
		origin_p2 = [origin_star[0], grating_split_pos_chip]
		xcenter_p2 = max(x_nominal, grating_split_pos_chip - 1) - 0
		ycenter_p2 = y_nominal - 0

		sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2,
		ycenter=ycenter_p2, origin=origin_p2,
		tar_spec=normalSED,
		band_start=brange[0], band_end=brange[1],
		conf=chip.sls_conf[1],
		isAlongY=0,
		flat_cube=flat_cube)

		# self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
		pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2],
		psf_model=psf_model, bandNo=i + 1,
		grating_split_pos=grating_split_pos,
		local_wcs=chip_wcs_local, pos_img = pos_img)

		del sdp_p1
		del sdp_p2
		elif grating_split_pos_chip<=gal_origin[1]:
		sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
		ycenter=y_nominal - 0, origin=origin_star,
		tar_spec=normalSED,
		band_start=brange[0], band_end=brange[1],
		conf=chip.sls_conf[1],
		isAlongY=0,
		flat_cube=flat_cube)
		# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
		pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
		psf_model=psf_model, bandNo=i + 1,
		grating_split_pos=grating_split_pos,
		local_wcs=chip_wcs_local, pos_img = pos_img)
		del sdp
		elif grating_split_pos_chip>=gal_end[1]:
		sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
		ycenter=y_nominal - 0, origin=origin_star,
		tar_spec=normalSED,
		band_start=brange[0], band_end=brange[1],
		conf=chip.sls_conf[0],
		isAlongY=0,
		flat_cube=flat_cube)
		# self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local)
		pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal],
		psf_model=psf_model, bandNo=i + 1,
		grating_split_pos=grating_split_pos,
		local_wcs=chip_wcs_local, pos_img = pos_img)
		del sdp

		# print(self.y_nominal, starImg.center.y, starImg.ymin)
		# del psf
		return 1, pos_shear

config/config_C6_fits.yaml

+52 −36

Original line number	Diff line number	Diff line
		@@ -9,9 +9,11 @@
		# Base diretories and naming setup
		# Can add some of the command-line arguments here as well;
		# OK to pass either way or both, as long as they are consistent
		work_dir: "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C6/"
		work_dir: "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C8/"
		data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
		run_name: "C6_fits_testRun"
		run_name: "testRun_fits"
		project_cycle: 6
		run_counter: 1

		# Whether to use MPI
		run_option:
		@@ -46,7 +48,7 @@ catalog_options:
		#stamp_SED:

		# Only simulate stars?
		star_only: NO
		star_only: YES

		# Only simulate galaxies?
		galaxy_only: NO
		@@ -63,12 +65,22 @@ catalog_options:
		# Observation setting
		###############################################
		obs_setting:

		# Options for survey types:
		# "Photometric": simulate photometric chips only
		# "Spectroscopic": simulate slitless spectroscopic chips only
		# "FGS": simulate FGS chips only (31-42)
		# "All": simulate full focal plane
		survey_type: "Photometric"
		# "CALIBRATION": falt, bias, dark with or without postflash
		survey_type: "All"
		# "LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null
		# 'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670',
		# 'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050','LED13': '340', 'LED14': '365'
		# LED_TYPE: ['LED5']
		# LED_TIME: [1.]
		# # unit e- ,flat level
		# FLAT_LEVEL: 20000
		# FLAT_LEVEL_FIL: 'g'

		# Exposure time [seconds]
		exp_time: 150.
		@@ -79,10 +91,10 @@ obs_setting:

		# Default Pointing [degrees]
		# Note: NOT valid when a pointing list file is specified
		ra_center: 244.972773
		dec_center: 39.895901
		ra_center: 192.8595
		dec_center: 27.1283
		# Image rotation [degree]
		image_rot: 109.59
		image_rot: -113.4333

		# (Optional) a file of point list
		# if you just want to run default pointing:
		@@ -90,8 +102,6 @@ obs_setting:
		# - pointing_file: null
		pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/"
		pointing_file: "pointing_radec_246.5_40.dat"
		# pointing_dir: null
		# pointing_file: null

		# Number of calibration pointings
		np_cal: 0
		@@ -106,16 +116,20 @@ obs_setting:
		# - give a list of indexes of chips: [ip_1, ip_2...]
		# - run all chips: null
		# Note: for all pointings
		run_chips: [7]
		run_chips: [1, 2, 3, 7, 8, 9]

		# Whether to enable astrometric modeling
		enable_astrometric_model: False
		enable_astrometric_model: True

		# Whether to enable straylight model
		enable_straylight_model: True

		# Cut by saturation magnitude in which band?
		cut_in_band: "z"

		# saturation magnitude margin
		mag_sat_margin: -2.5
		# mag_sat_margin: -2.5
		mag_sat_margin: -15.

		# limiting magnitude margin
		mag_lim_margin: +1.0
		@@ -138,7 +152,7 @@ psf_setting:
		# path to PSF data
		# NOTE: only valid for "Interp" PSF
		psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"

		psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/"
		###############################################
		# Shear setting
		###############################################
		@@ -160,25 +174,27 @@ ins_effects:
		# switches
		# Note: bias_16channel, gain_16channel, and shutter_effect
		# is currently not applicable to "FGS" observations
		field_dist: OFF # Whether to add field distortions
		add_back: OFF # Whether to add sky background
		add_dark: OFF # Whether to add dark noise
		add_readout: OFF # Whether to add read-out (Gaussian) noise
		add_bias: OFF # Whether to add bias-level to images
		bias_16channel: OFF # Whether to add different biases for 16 channels
		gain_16channel: OFF # Whether to make different gains for 16 channels
		shutter_effect: OFF # Whether to add shutter effect
		flat_fielding: OFF # Whether to add flat-fielding effect
		prnu_effect: OFF # Whether to add PRNU effect
		non_linear: OFF # Whether to add non-linearity
		cosmic_ray: OFF # Whether to add cosmic-ray
		cray_differ: OFF # Whether to generate different cosmic ray maps CAL and MS output
		cte_trail: OFF # Whether to simulate CTE trails
		saturbloom: OFF # Whether to simulate Saturation & Blooming
		add_badcolumns: OFF # Whether to add bad columns
		add_hotpixels: OFF # Whether to add hot pixels
		add_deadpixels: OFF # Whether to add dead(dark) pixels
		bright_fatter: OFF # Whether to simulate Brighter-Fatter (also diffusion) effect
		field_dist: YES # Whether to add field distortions
		add_back: YES # Whether to add sky background
		add_dark: YES # Whether to add dark noise
		add_readout: YES # Whether to add read-out (Gaussian) noise
		add_bias: YES # Whether to add bias-level to images
		bias_16channel: YES # Whether to add different biases for 16 channels
		gain_16channel: YES # Whether to make different gains for 16 channels
		shutter_effect: YES # Whether to add shutter effect
		flat_fielding: YES # Whether to add flat-fielding effect
		prnu_effect: YES # Whether to add PRNU effect
		non_linear: YES # Whether to add non-linearity
		cosmic_ray: YES # Whether to add cosmic-ray
		cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output
		cte_trail: NO # Whether to simulate CTE trails, CTI_lgl_v0.3.tar.gz
		saturbloom: YES # Whether to simulate Saturation & Blooming
		add_badcolumns: YES # Whether to add bad columns
		add_hotpixels: YES # Whether to add hot pixels
		add_deadpixels: YES # Whether to add dead(dark) pixels
		bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect
		add_prescan: NO # Whether to add pre/over-scan
		format_output: NO ##1*16 output

		# Values:
		# default values have been defined individually for each chip in: