add new psf convolve for single source; fix header bug (18700088) · Commits · csst-sims / csst_msc_sim

observation_sim/config/header/ImageHeader.py

+1 −1

Original line number	Diff line number	Diff line
		@@ -434,7 +434,7 @@ def generatePrimaryHeader(xlen=9216, ylen=9232, pointing_id='00000001', pointing
		co = coord.SkyCoord(ra, dec, unit='deg')

		ra_hms = format(co.ra.hms.h, '02.0f') + format(co.ra.hms.m,
		'02.0f') + format(co.ra.hms.s, '02.1f')
		'02.0f') + format(co.ra.hms.s, '04.1f')
		dec_hms = format(co.dec.dms.d, '02.0f') + format(abs(co.dec.dms.m),
		'02.0f') + format(abs(co.dec.dms.s), '02.0f')
		if dec >= 0:

observation_sim/mock_objects/MockObject.py

+69 −38

Original line number	Diff line number	Diff line
		@@ -120,7 +120,6 @@ class MockObject(object):
		if self.logger:
		self.logger.error(e)
		return 2, None

		# Set Galsim Parameters
		if self.getMagFilter(filt) <= 15:
		folding_threshold = 5.e-4
		@@ -241,11 +240,10 @@ class MockObject(object):

		def addSLStoChipImageWithPSF(self, sdp=None, chip=None, pos_img_local=[1, 1], psf_model=None, bandNo=1, grating_split_pos=3685, local_wcs=None, pos_img=None):
		spec_orders = sdp.compute_spec_orders()

		pos_shear = galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
		if chip.slsPSFOptim:
		for k, v in spec_orders.items():
		img_s = v[0]
		pos_shear = galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)

		nan_ids = np.isnan(img_s)
		if img_s[nan_ids].shape[0] > 0:
		@@ -280,47 +278,80 @@ class MockObject(object):

		else:
		for k, v in spec_orders.items():
		img_s = v[0]
		# print(bandNo,k)
		try:
		psf, pos_shear = psf_model.get_PSF(
		chip, pos_img_local=pos_img_local, bandNo=bandNo, galsimGSObject=True, g_order=k, grating_split_pos=grating_split_pos)
		except:
		psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)
		# img_s = v[0]
		# # print(bandNo,k)
		# try:
		# psf, pos_shear = psf_model.get_PSF(
		# chip, pos_img_local=pos_img_local, bandNo=bandNo, galsimGSObject=True, g_order=k, grating_split_pos=grating_split_pos)
		# except:
		# psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)

		# psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs)

		# psf_img_m = psf_img.array

		# #########################################################
		# # DEBUG
		# #########################################################
		# # ids_p = psf_img_m < 0
		# # psf_img_m[ids_p] = 0

		# # from astropy.io import fits
		# # fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m)

		# # print("DEBUG: orig_off is", orig_off)
		# nan_ids = np.isnan(img_s)
		# if img_s[nan_ids].shape[0] > 0:
		# img_s[nan_ids] = 0
		# print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
		# #########################################################
		# img_s, orig_off = convolveImg(img_s, psf_img_m)
		# origin_order_x = v[1] - orig_off[0]
		# origin_order_y = v[2] - orig_off[1]

		psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs)
		# specImg = galsim.ImageF(img_s)
		# # photons = galsim.PhotonArray.makeFromImage(specImg)
		# # photons.x += origin_order_x
		# # photons.y += origin_order_y

		psf_img_m = psf_img.array
		# # xlen_imf = int(specImg.xmax - specImg.xmin + 1)
		# # ylen_imf = int(specImg.ymax - specImg.ymin + 1)
		# # stamp = galsim.ImageF(xlen_imf, ylen_imf)
		# # stamp.wcs = local_wcs
		# # stamp.setOrigin(origin_order_x, origin_order_y)

		#########################################################
		# DEBUG
		#########################################################
		# ids_p = psf_img_m < 0
		# psf_img_m[ids_p] = 0
		# specImg.wcs = local_wcs
		# specImg.setOrigin(origin_order_x, origin_order_y)

		# from astropy.io import fits
		# fits.writeto(str(bandNo) + '_' + str(k) + '_psf.fits', psf_img_m)
		print('DEBUG: BEGIN -----------',bandNo,k)

		img_s = v[0]

		# print("DEBUG: orig_off is", orig_off)
		nan_ids = np.isnan(img_s)
		if img_s[nan_ids].shape[0] > 0:
		img_s[nan_ids] = 0
		print("DEBUG: specImg nan num is", img_s[nan_ids].shape[0])
		#########################################################
		origin_order_x = v[1]
		origin_order_y = v[2]
		specImg = galsim.ImageF(img_s)

		specImg.wcs = local_wcs
		specImg.setOrigin(origin_order_x, origin_order_y)
		try:
		specImg = psf_model.get_PSF_AND_convolve_withsubImg(chip, cutImg=specImg, bandNo=bandNo, g_order=k, grating_split_pos=grating_split_pos)
		except:
		psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img)

		psf_img = psf.drawImage(nx=100, ny=100, wcs=local_wcs)

		psf_img_m = psf_img.array

		img_s, orig_off = convolveImg(img_s, psf_img_m)
		origin_order_x = v[1] - orig_off[0]
		origin_order_y = v[2] - orig_off[1]

		specImg = galsim.ImageF(img_s)
		# photons = galsim.PhotonArray.makeFromImage(specImg)
		# photons.x += origin_order_x
		# photons.y += origin_order_y

		# xlen_imf = int(specImg.xmax - specImg.xmin + 1)
		# ylen_imf = int(specImg.ymax - specImg.ymin + 1)
		# stamp = galsim.ImageF(xlen_imf, ylen_imf)
		# stamp.wcs = local_wcs
		# stamp.setOrigin(origin_order_x, origin_order_y)

		specImg.wcs = local_wcs
		specImg.setOrigin(origin_order_x, origin_order_y)

observation_sim/psf/PSFInterpSLS.py

+93 −0

Original line number	Diff line number	Diff line
		@@ -481,6 +481,99 @@ class PSFInterpSLS(PSFModel):

		return PSF_int_trans, PSF_int

		def get_PSF_AND_convolve_withsubImg(self, chip, cutImg=None, bandNo=1, g_order='A', grating_split_pos=3685):
		"""
		Get the PSF at a given image position

		Parameters:
		chip: A 'Chip' object representing the chip we want to extract PSF from.
		pos_img: A 'galsim.Position' object representing the image position.
		bandpass: A 'galsim.Bandpass' object representing the wavelength range.
		pixSize: The pixels size of psf matrix
		findNeighMode: 'treeFind' or 'hoclistFind'
		Returns:
		PSF: A 'galsim.GSObject'.
		"""
		order_IDs = {'A': '1', 'B': '0', 'C': '0', 'D': '0', 'E': '0'}
		contam_order_sigma = {'C': 0.28032344707964174,
		'D': 0.39900182912061344, 'E': 1.1988309797685412} # arcsec
		x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
		y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
		# print(pos_img.x - x_start)
		# pos_img_x = pos_img_local[0] + x_start
		# pos_img_y = pos_img_local[1] + y_start
		# pos_img = galsim.PositionD(pos_img_x, pos_img_y)
		centerPos_local = cutImg.ncol/2.
		if centerPos_local < grating_split_pos:
		psf_data = self.grating1_data
		else:
		psf_data = self.grating2_data

		grating_order = order_IDs[g_order]
		# if grating_order in ['-2','-1','2']:
		# grating_order = '1'

		# if grating_order in ['0', '1']:
		psf_order = psf_data['order'+grating_order]
		psf_order_b = psf_order['band'+str(bandNo)]
		psf_b_dat = psf_order_b['band_data']
		# pos_p = psf_b_dat[1].data
		pos_p = psf_b_dat[1].data/chip.pix_size - np.array([y_start, x_start])

		pc_coeff = psf_b_dat[2].data
		pcs = psf_b_dat[0].data

		npc = 10
		m_size = int(pcs.shape[0]**0.5)
		tmp_img = cutImg*0
		for j in np.arange(npc):
		X_ = jnp.hstack((pos_p[:,1].flatten()[:, None], pos_p[:,0].flatten()[:, None]),dtype=np.float32)
		Z_ = (pc_coeff[j].astype(np.float32)).flatten()
		# print(pc_coeff[j].shape[0], pos_p[:,1].shape[0], pos_p[:,0].shape[0])
		cx_len = int(chip.npix_x)
		cy_len = int(chip.npix_y)
		n_x = jnp.arange(0, cx_len, 1, dtype = int)
		n_y = jnp.arange(0, cy_len, 1, dtype = int)

		M, N = jnp.meshgrid(n_x, n_y)
		t1=datetime.datetime.now()
		# U = interpolate.griddata(X_, Z_, (M[0:cy_len, 0:cx_len],N[0:cy_len, 0:cx_len]),
		# method='nearest',fill_value=1.0)
		ys = cutImg.ymin
		if ys < 0:
		ys = 0
		ye = cutImg.ymin+cutImg.nrow
		if ye >= cy_len-1:
		ye = cy_len-1
		if ye - ys <=0:
		continue
		xs = cutImg.xmin
		if xs < 0:
		xs = 0
		xe = cutImg.xmin+cutImg.ncol
		if xe >= cx_len-1:
		xe = cx_len-1
		if xe - xs <=0:
		continue
		U = interpolate.griddata(X_, Z_, (M[ys:ye, xs:xe],N[ys:ye, xs:xe]),
		method='nearest',fill_value=1.0)
		t2=datetime.datetime.now()

		print("time interpolate:", t2-t1)

		if U.shape != cutImg.array.shape:
		print('DEBUG:SHAPE',cutImg.ncol,cutImg.nrow,cutImg.xmin, cutImg.ymin)
		continue
		img_tmp = cutImg.array*U
		psf = pcs[:, j].reshape(m_size, m_size)
		tmp_img = tmp_img + signal.fftconvolve(img_tmp, psf, mode='same', axes=None)

		t3=datetime.datetime.now()
		print("time convole:", t3-t2)
		del U
		return tmp_img


		def convolveFullImgWithPCAPSF(self, chip, folding_threshold=5.e-3):
		keys_L1= chip_utils.getChipSLSGratingID(chip.chipID)
		# keys_L2 = ['order-2','order-1','order0','order1','order2']