diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index 36da1803de00642ced3cbaa94fba02cf79737715..d9253fd03936e848b026d2fe065a86611c909ec6 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -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, diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 75229b718b932cc73623ddb0b50cc02ac3ec2ff9..0281dcddae182eafd299637289f9cf2c3b0ceead 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -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'] diff --git a/ObservationSim/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 5b2cb4488a2689a60e6717df862c1a8deed6e22e..d026560f5873a9d86649db5efa85da97ca52d6a3 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -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 = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.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) + 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(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[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 @@ -212,7 +272,7 @@ class Galaxy(MockObject): gsp = galsim.GSParams(folding_threshold=folding_threshold) # nphotons_sum = 0 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 +296,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,8 +311,8 @@ 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, @@ -261,13 +321,13 @@ class Galaxy(MockObject): conf=chip.sls_conf[0], isAlongY=0) - 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, @@ -276,27 +336,27 @@ class Galaxy(MockObject): conf=chip.sls_conf[1], isAlongY=0) - 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) + 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) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) del sdp # print(self.y_nominal, starImg.center.y, starImg.ymin) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 4bdb2a9a9859f32590c7d9140e5c044109dc316e..5b6dffb9d3582753fe5e36576aa671de533dd428 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -75,7 +75,7 @@ class MockObject(object): dec = self.param["dec"] return galsim.CelestialCoord(ra=ra*galsim.degrees,dec=dec*galsim.degrees) - def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True): + def getPosImg_Offset_WCS(self, img, fdmodel=None, chip=None, verbose=True, img_header=None): self.posImg = img.wcs.toImage(self.getPosWorld()) self.localWCS = img.wcs.local(self.posImg) if (fdmodel is not None) and (chip is not None): @@ -93,7 +93,23 @@ class MockObject(object): dx = x - self.x_nominal dy = y - self.y_nominal self.offset = galsim.PositionD(dx, dy) - return self.posImg, self.offset, self.localWCS + + from astropy import wcs + + if img_header is not None: + self.real_wcs = galsim.FitsWCS(header=img_header) + else: + self.real_wcs = None + + return self.posImg, self.offset, self.localWCS, self.real_wcs + + def getRealPos(self,img, global_x = 0., global_y = 0., img_real_wcs = None): + img_global_pos = galsim.PositionD(global_x, global_y) + cel_pos = img.wcs.toWorld(img_global_pos) + realPos = img_real_wcs.toImage(cel_pos) + + return realPos + def drawObject(self, img, final, flux=None, filt=None, tel=None, exptime=150.): """ Draw (point like) object on img. @@ -146,6 +162,17 @@ class MockObject(object): 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] try: @@ -170,21 +197,29 @@ class MockObject(object): star = star.withFlux(nphotons) star = galsim.Convolve(psf, star) - stamp = star.drawImage(wcs=self.localWCS, method='phot', offset=self.offset, save_photons=True) + # stamp = star.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 = star.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) - # Test stamp size - # print(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 = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.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) + chip.img.setOrigin(0, 0) stamp[bounds] = chip.img[bounds] for i in range(len(photons_list)): if i == 0: @@ -193,6 +228,23 @@ class MockObject(object): chip.sensor.accumulate(photons_list[i], stamp, resume=True) chip.img[bounds] = stamp[bounds] + + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + # Test stamp size + # print(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[bounds] = chip.img[bounds] + # 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) + # + # chip.img[bounds] = stamp[bounds] # print(chip.img.array.sum()) # print("nphotons_sum = ", nphotons_sum) del photons_list @@ -200,7 +252,7 @@ class MockObject(object): return True, pos_shear - def addSLStoChipImage(self, sdp = None, chip = None, xOrderSigPlus = None): + def addSLStoChipImage(self, sdp = None, chip = None, xOrderSigPlus = None, local_wcs=None): spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] @@ -215,15 +267,17 @@ class MockObject(object): 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 = self.localWCS + stamp.wcs = local_wcs stamp.setOrigin(origin_order_x, origin_order_y) - bounds = stamp.bounds & chip.img.bounds + bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x-1, 0, chip.npix_y-1) if bounds.area() == 0: continue + chip.img.setOrigin(0, 0) stamp[bounds] = chip.img[bounds] chip.sensor.accumulate(photons, stamp) chip.img[bounds] = stamp[bounds] + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) del stamp del spec_orders @@ -249,19 +303,30 @@ class MockObject(object): 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) + 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] psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) star = galsim.DeltaFunction(gsparams=gsp) star = star.withFlux(tel.pupil_area * exptime) star = galsim.Convolve(psf, star) - starImg = star.drawImage(nx=100, ny=100, wcs=self.localWCS) + starImg = star.drawImage(nx=100, ny=100, 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] @@ -273,8 +338,8 @@ class MockObject(object): 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, @@ -283,13 +348,13 @@ class MockObject(object): conf=chip.sls_conf[0], isAlongY=0) - 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, @@ -298,27 +363,27 @@ class MockObject(object): conf=chip.sls_conf[1], isAlongY=0) - 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) + 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) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) del sdp del psf return True, pos_shear diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 4ca0669e0065792472db544b399870cc5095f18d..5e79ec38bb51437f32599d6fce48a53d1129bc7a 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -258,8 +258,22 @@ class Observation(object): raise ValueError("Unknown shear input") # chip_output.logger.info("debug point #4") + header_wcs = generateExtensionHeader( + xlen=chip.npix_x, + ylen=chip.npix_y, + ra=ra_cen, + dec=dec_cen, + pa=pointing.img_pa.deg, + gain=chip.gain, + readout=chip.read_noise, + dark=chip.dark_noise, + saturation=90000, + psize=chip.pix_scale, + row_num=chip.rowID, + col_num=chip.colID, + extName='raw') - pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False) + pos_img, offset, local_wcs, real_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=header_wcs) if pos_img.x == -1 or pos_img.y == -1: # Exclude object which is outside the chip area (after field distortion) # print("obj missed!!") diff --git a/setup.py b/setup.py index 7db1250b91c1c868f3c79bf139fdf0cd55e02df3..d3bcf71401a26dd881e684c3470aa421eb1d64c2 100644 --- a/setup.py +++ b/setup.py @@ -29,17 +29,17 @@ setup(name='CSSTSim', version='1.0.4', packages=find_packages(), install_requires=[ - 'numpy>=1.18.5', - 'galsim>=2.2.4', - 'pyyaml>=5.3.1', - 'astropy>=4.0.1', - 'scipy>=1.5.0', - 'mpi4py>=3.0.3', - 'sep>=1.0.3', - 'healpy>=1.14.0', - 'h5py>=2.10.0', - 'Cython>=0.29.21', - 'numba>=0.50.1' + # 'numpy>=1.18.5', + # 'galsim>=2.2.4', + # 'pyyaml>=5.3.1', + # 'astropy>=4.0.1', + # 'scipy>=1.5.0', + # 'mpi4py>=3.0.3', + # 'sep>=1.0.3', + # 'healpy>=1.14.0', + # 'h5py>=2.10.0', + # 'Cython>=0.29.21', + # 'numba>=0.50.1' ], package_data = { 'ObservationSim.Astrometry.lib': ['libshao.so'], @@ -52,4 +52,4 @@ setup(name='CSSTSim', 'Catalog.data': ['*.fits'], }, ext_modules = cythonize(extensions), -) \ No newline at end of file +)