Merge branch 'sls_rot_wcs' into 'develop' (9bee4eba) · Commits · Liu Dezi / csst_msc_sim

ObservationSim/Config/ChipOutput.py

+7 −2

Original line number	Diff line number	Diff line
		@@ -68,8 +68,13 @@ class ChipOutput(object):
		self.cat.write(self.hdr)

		def cat_add_obj(self, obj, pos_img, pos_shear):
		ximg = pos_img.x - self.chip.bound.xmin + 1.0
		yimg = pos_img.y - self.chip.bound.ymin + 1.0
		# ximg = pos_img.x - self.chip.bound.xmin + 1.0
		# yimg = pos_img.y - self.chip.bound.ymin + 1.0
		# print('-------------debug-----------------')
		# print('from global',ximg, yimg)
		ximg = obj.real_pos.x + 1.0
		yimg = obj.real_pos.y + 1.0
		# print('from loacl',ximg, yimg)

		line = self.fmt%(
		obj.id, int(self.chipLabel), self.filt.filter_type, ximg, yimg, obj.ra, obj.dec, obj.ra_orig, obj.dec_orig, obj.z, obj.getMagFilter(self.filt), obj.type,

ObservationSim/Config/Header/ImageHeader.py

+71 −3

Original line number	Diff line number	Diff line
		@@ -7,7 +7,7 @@ from astropy.io import fits
		import astropy.wcs as pywcs
		from collections import OrderedDict

		from scipy import math
		# from scipy import math
		import random

		import os
		@@ -87,6 +87,17 @@ def rotate_CD_matrix(cd, pa_aper):
		cd_rot = np.dot(mat, cd)
		return cd_rot

		def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232, w = None):
		rad = np.deg2rad(rot_angle)
		mat = np.zeros((2,2))
		mat[0,:] = np.array([np.cos(rad),-np.sin(rad)])
		mat[1,:] = np.array([np.sin(rad),np.cos(rad)])
		center = np.array([xlen/2, ylen/2])
		rot_pix = np.dot(mat, pix_xy-center) + center
		skyCoor = w.wcs_pix2world(np.array([rot_pix]), 1)

		return skyCoor


		# def Header_extention(xlen = 9216, ylen = 9232, gain = 1.0, readout = 5.0, dark = 0.02,saturation=90000, row_num = 1, col_num = 1):
		#
		@@ -153,7 +164,7 @@ def rotate_CD_matrix(cd, pa_aper):


		##9232 9216 898 534 1309 60 -40 -23.4333
		def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1):
		def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'):

		""" Creat a wcs frame for CCST with multiple extensions

		@@ -228,6 +239,63 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r
		r_dat['CD1_2'] = cd_rot[0,1]
		r_dat['CD2_1'] = cd_rot[1,0]
		r_dat['CD2_2'] = cd_rot[1,1]

		if filter in ['GU', 'GV', 'GI']:
		from astropy import wcs

		w = wcs.WCS(naxis=2)
		w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']]
		w.wcs.cd = cd_rot
		w.wcs.crval = [ra_ref, dec_ref]
		w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']]

		# test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)

		sls_rot = 1
		if i > 2:
		sls_rot = -sls_rot

		sn_x = 30
		sn_y = 30
		x_pixs = np.zeros(sn_y * sn_x)
		y_pixs = np.zeros(sn_y * sn_x)
		xpixs_line = np.linspace(1, xlen, sn_x)
		ypixs_line = np.linspace(1, ylen, sn_y)

		sky_coors = []

		for n1, y in enumerate(ypixs_line):
		for n2, x in enumerate(xpixs_line):
		i_pix = n1 * sn_x + n2
		x_pixs[i_pix] = x
		y_pixs[i_pix] = y

		pix_coor = np.array([x, y])
		sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, w=w)
		# print(sc1[0,0],sc1[0,1])
		sky_coors.append((sc1[0, 0], sc1[0, 1]))

		from astropy.coordinates import SkyCoord
		from astropy.wcs.utils import fit_wcs_from_points

		wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]),
		world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN')

		# print(wcs_new)
		# test_center = wcs_new.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1)
		#
		# print(test_center - test_center_o)

		r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0]
		r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1]
		r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0]
		r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1]
		r_dat['CRPIX1'] = wcs_new.wcs.crpix[0]
		r_dat['CRPIX2'] = wcs_new.wcs.crpix[1]

		r_dat['CRVAL1'] = wcs_new.wcs.crval[0]
		r_dat['CRVAL2'] = wcs_new.wcs.crval[1]

		return r_dat


		@@ -347,7 +415,7 @@ def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -2
		h_ext['CCDCHIP'] = 'ccd' + CCDID[k-1].rjust(2,'0')
		h_ext['POS_ANG'] = pa
		header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize,
		row_num=row_num, col_num=col_num)
		row_num=row_num, col_num=col_num, filter = h_ext['FILTER'])

		h_ext['CRPIX1'] = header_wcs['CRPIX1']
		h_ext['CRPIX2'] = header_wcs['CRPIX2']

ObservationSim/Instrument/Chip/Chip.py

+19 −0

Original line number	Diff line number	Diff line
		@@ -82,6 +82,8 @@ class Chip(FocalPlane):
		else:
		self.sensor = galsim.Sensor()

		self.flat_cube = None # for spectroscopic flat field cube simulation

		# def _getChipRowCol(self):
		# self.rowID = (self.chipID - 1) // 5 + 1
		# self.colID = (self.chipID - 1) % 5 + 1
		@@ -846,3 +848,20 @@ class Chip(FocalPlane):
		# sub_img.write("%s/%s_%s.fits" % (chip_output.subdir, img_name_root, rowcoltag))
		# del sub_img
		return img

		def loadSLSFLATCUBE(self, flat_fn='flat_cube.fits'):
		from astropy.io import fits
		try:
		with pkg_resources.files('ObservationSim.Instrument.data').joinpath(flat_fn) as data_path:
		flat_fits = fits.open(data_path, ignore_missing_simple=True)
		except AttributeError:
		with pkg_resources.path('ObservationSim.Instrument.data', flat_fn) as data_path:
		flat_fits = fits.open(data_path, ignore_missing_simple=True)

		fl = len(flat_fits)
		fl_sh = flat_fits[0].data.shape
		assert fl == 4, 'FLAT Field Cube is Not 4 layess!!!!!!!'
		self.flat_cube = np.zeros([fl, fl_sh[0], fl_sh[1]])
		for i in np.arange(0, fl, 1):
		self.flat_cube[i] = flat_fits[i].data

ObservationSim/Instrument/data/flatCube/flat_cube.fits

0 → 100644

+2.54 GiB

File added.

No diff preview for this file type.

View file

ObservationSim/MockObject/Galaxy.py

+97 −30

Original line number	Diff line number	Diff line
		@@ -107,6 +107,18 @@ class Galaxy(MockObject):
		folding_threshold = 5.e-3
		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)

		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)

		real_wcs_local = self.real_wcs.local(self.real_pos)

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

		@@ -150,20 +162,29 @@ class Galaxy(MockObject):
		gal = galsim.Convolve(psf, gal)

		# Use (explicit) stamps to draw
		stamp = gal.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True)
		# stamp = gal.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True)
		# xmax = max(xmax, stamp.xmax)
		# ymax = max(ymax, stamp.ymax)
		# photons = stamp.photons
		# photons.x += self.x_nominal
		# photons.y += self.y_nominal
		# photons_list.append(photons)

		stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=offset, save_photons=True)
		xmax = max(xmax, stamp.xmax)
		ymax = max(ymax, stamp.ymax)
		photons = stamp.photons
		photons.x += self.x_nominal
		photons.y += self.y_nominal
		photons.x += x_nominal
		photons.y += y_nominal
		photons_list.append(photons)

		# print('xmax = %d, ymax = %d '%(xmax, ymax))

		stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1))
		stamp.wcs = self.localWCS
		stamp.setCenter(self.x_nominal, self.y_nominal)
		bounds = stamp.bounds & chip.img.bounds
		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)
		chip.img.setOrigin(0, 0)
		stamp[bounds] = chip.img[bounds]

		if not big_galaxy:
		@@ -180,10 +201,36 @@ class Galaxy(MockObject):
		else:
		sensor.accumulate(photons_list[i], stamp, resume=True)

		# print(stamp.array.sum())
		# chip.img[bounds] += stamp[bounds]
		chip.img[bounds] = stamp[bounds]
		# print("nphotons_sum = ", nphotons_sum)

		chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)



		# stamp = galsim.ImageF(int(xmax1.1), int(ymax1.1))
		# stamp.wcs = self.localWCS
		# stamp.setCenter(self.x_nominal, self.y_nominal)
		# bounds = stamp.bounds & chip.img.bounds
		# stamp[bounds] = chip.img[bounds]
		#
		# if not big_galaxy:
		# for i in range(len(photons_list)):
		# if i == 0:
		# chip.sensor.accumulate(photons_list[i], stamp)
		# else:
		# chip.sensor.accumulate(photons_list[i], stamp, resume=True)
		# else:
		# sensor = galsim.Sensor()
		# for i in range(len(photons_list)):
		# if i == 0:
		# sensor.accumulate(photons_list[i], stamp)
		# else:
		# sensor.accumulate(photons_list[i], stamp, resume=True)
		#
		# # print(stamp.array.sum())
		# # chip.img[bounds] += stamp[bounds]
		# chip.img[bounds] = stamp[bounds]
		# # print("nphotons_sum = ", nphotons_sum)
		del photons_list
		del stamp
		return True, pos_shear
		@@ -201,6 +248,19 @@ class Galaxy(MockObject):
		normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T,
		names=('WAVELENGTH', 'FLUX'))

		self.real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y,
		img_real_wcs=self.real_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)

		real_wcs_local = self.real_wcs.local(self.real_pos)


		big_galaxy = False
		if self.hlr_disk > 3.0: # Very big galaxy
		big_galaxy = True
		@@ -211,8 +271,11 @@ class Galaxy(MockObject):
		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 = chip.bound.xmin + grating_split_pos
		grating_split_pos_chip = 0 + grating_split_pos
		for i in range(len(bandpass_list)):
		bandpass = bandpass_list[i]

		@@ -236,10 +299,10 @@ class Galaxy(MockObject):
		gal = gal.shear(gal_shear)
		gal = galsim.Convolve(psf, gal)

		starImg = gal.drawImage(wcs=self.localWCS)
		starImg = gal.drawImage(wcs=real_wcs_local)

		origin_star = [self.y_nominal - (starImg.center.y - starImg.ymin),
		self.x_nominal - (starImg.center.x - starImg.xmin)]
		origin_star = [y_nominal - (starImg.center.y - starImg.ymin),
		x_nominal - (starImg.center.x - starImg.xmin)]

		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]
		@@ -251,52 +314,56 @@ class Galaxy(MockObject):
		subImg_p1 = starImg.array[:, 0:subSlitPos]
		star_p1 = galsim.Image(subImg_p1)
		origin_p1 = origin_star
		xcenter_p1 = min(self.x_nominal,grating_split_pos_chip-1) - chip.bound.xmin
		ycenter_p1 = self.y_nominal-chip.bound.ymin
		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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
		conf=chip.sls_conf[0],
		isAlongY=0)
		isAlongY=0,
		flat_cube=flat_cube)

		self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus)
		self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)

		subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]]
		star_p2 = galsim.Image(subImg_p2)
		origin_p2 = [origin_star[0], grating_split_pos_chip]
		xcenter_p2 = max(self.x_nominal, grating_split_pos_chip - 1) - chip.bound.xmin
		ycenter_p2 = self.y_nominal - chip.bound.ymin
		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=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
		conf=chip.sls_conf[1],
		isAlongY=0)
		isAlongY=0,
		flat_cube=flat_cube)

		self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus)
		self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)

		del sdp_p1
		del sdp_p2
		elif grating_split_pos_chip<=gal_origin[1]:
		sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin,
		ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star,
		sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
		ycenter=y_nominal - 0, origin=origin_star,
		tar_spec=normalSED,
		band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
		conf=chip.sls_conf[1],
		isAlongY=0)
		self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus)
		isAlongY=0,
		flat_cube=flat_cube)
		self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
		del sdp
		elif grating_split_pos_chip>=gal_end[1]:
		sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin,
		ycenter=self.y_nominal - chip.bound.ymin, origin=origin_star,
		sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0,
		ycenter=y_nominal - 0, origin=origin_star,
		tar_spec=normalSED,
		band_start=bandpass.blue_limit * 10, band_end=bandpass.red_limit * 10,
		conf=chip.sls_conf[0],
		isAlongY=0)
		self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus)
		isAlongY=0,
		flat_cube=flat_cube)
		self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local)
		del sdp

		# print(self.y_nominal, starImg.center.y, starImg.ymin)