From fedb2da1e093dd4f64a84f7e861cfa1b463488ad Mon Sep 17 00:00:00 2001 From: xin Date: Tue, 27 Sep 2022 15:54:09 +0800 Subject: [PATCH] local wcs test --- ObservationSim/MockObject/MockObject.py | 74 +++++++++++++++++++++---- ObservationSim/ObservationSim.py | 16 +++++- setup.py | 24 ++++---- 3 files changed, 89 insertions(+), 25 deletions(-) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 4bdb2a9..99a72ed 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,15 @@ class MockObject(object): folding_threshold = 5.e-3 gsp = galsim.GSParams(folding_threshold=folding_threshold) + real_pos = self.getRealPos(chip.img, global_x=self.posImg.x, global_y=self.posImg.y, img_real_wcs=self.real_wcs) + + x, y = real_pos.x + 0.5, 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) + for i in range(len(bandpass_list)): bandpass = bandpass_list[i] try: @@ -170,21 +195,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=self.real_wcs, 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 = self.real_wcs.local(real_pos) + stamp.setCenter(x_nominal, y_nominal) + bounds = stamp.bounds & galsim.BoundsD(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 +226,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 diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 4ca0669..a661b22 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") + h_ext = generateExtensionHeader( + xlen=chip.npix_x, + ylen=chip.npix_y, + ra=pointing.ra, + dec=pointing.dec, + 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) 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 7db1250..d3bcf71 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 +) -- GitLab