local wcs test

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

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!!")

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
+)