Merge remote-tracking branch 'origin/new_sim' into sim_scheduler (33a3e5e3) · Commits · Liu Dezi / 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/MockObject.py

+87 −0

Original line number	Diff line number	Diff line
		import os
		import galsim
		import numpy as np
		import astropy.constants as cons
		from astropy import wcs
		from astropy.table import Table
		import astropy.io.fits as fitsio

		from ObservationSim.MockObject._util import magToFlux, VC_A, convolveGaussXorders, convolveImg
		from ObservationSim.MockObject._util import integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, \
		@@ -444,3 +446,88 @@ class MockObject(object):
		sig_obj = np.std(stamp.array)
		snr_obj = img_flux / sig_obj
		return snr_obj



		def drawObj_PSF(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0,
		exptime=150., fd_shear=None, chip_output=None):
		if nphotons_tot == None:
		nphotons_tot = self.getElectronFluxFilt(filt, tel, exptime)
		# print("nphotons_tot = ", nphotons_tot)

		try:
		full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
		except Exception as e:
		print(e)
		if self.logger:
		self.logger.error(e)
		return 2, None

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

		# Get real image position of object (deal with chip rotation w.r.t its center)
		self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
		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))
		y_nominal = int(np.floor(y + 0.5))
		dx = x - x_nominal
		dy = y - y_nominal
		offset = galsim.PositionD(dx, dy)
		# Get real local wcs of object (deal with chip rotation w.r.t its center)
		chip_wcs_local = self.chip_wcs.local(self.real_pos)
		is_updated = 0

		# Loop over all sub-bandpasses
		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)
		continue
		ratio = sub / full
		if not (ratio == -1 or (ratio != ratio)):
		nphotons = ratio * nphotons_tot
		else:
		continue

		# Get PSF model
		psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass,
		folding_threshold=folding_threshold)
		star_temp = psf.withFlux(nphotons)

		if i==0:
		star = star_temp
		else:
		star = star+star_temp

		pixelScale = 0.074
		stamp = star.drawImage(wcs=chip_wcs_local, offset=offset)
		#stamp = star.drawImage(nx=256, ny=256, scale=pixelScale)
		if np.sum(np.isnan(stamp.array)) > 0:
		return None


		fn = chip_output.subdir + "/psfIDW"
		os.makedirs(fn, exist_ok=True)
		fn = fn + "/ccd_{:}".format(chip.chipID)+"_psf_"+str(self.param['id'])+".fits"
		if fn != None:
		if os.path.exists(fn):
		os.remove(fn)
		hdu = fitsio.PrimaryHDU()
		hdu.data = stamp.array
		hdu.header.set('name', self.type)
		hdu.header.set('pixScale', pixelScale)
		hdu.header.set('objID', self.param['id'])
		hdu.writeto(fn)

		del stamp
		return None

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

ObservationSim/ObservationSim.py

+14 −0

Original line number	Diff line number	Diff line
		@@ -304,6 +304,20 @@ class Observation(object):
		normFilter=norm_filt,
		fd_shear=fd_shear)

		if isUpdated == 1 and self.config["run_option"]["out_psf"]:
		obj.drawObj_PSF(
		tel=self.tel,
		pos_img=pos_img,
		psf_model=psf_model,
		bandpass_list=filt.bandpass_sub_list,
		filt=filt,
		chip=chip,
		g1=obj.g1,
		g2=obj.g2,
		exptime=pointing.exp_time,
		fd_shear=fd_shear,
		chip_output=chip_output)

		if isUpdated == 1:
		# TODO: add up stats
		chip_output.cat_add_obj(obj, pos_img, pos_shear)

config/config_C6_fits.yaml

+56 −37

File changed.

Preview size limit exceeded, changes collapsed.