From 264eb57f7786d7023c5731990de3347e48b37007 Mon Sep 17 00:00:00 2001 From: xin Date: Wed, 7 Sep 2022 16:48:44 +0800 Subject: [PATCH 01/16] chip of slitless spectoscory rotate 1 deg around chip center, add these feature to code --- ObservationSim/Config/Header/ImageHeader.py | 72 ++++++++++++++++++++- 1 file changed, 70 insertions(+), 2 deletions(-) diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 75229b7..0281dcd 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'] -- GitLab From fedb2da1e093dd4f64a84f7e861cfa1b463488ad Mon Sep 17 00:00:00 2001 From: xin Date: Tue, 27 Sep 2022 15:54:09 +0800 Subject: [PATCH 02/16] 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 From 2468ff7c6dacf991e9373782edfdb8f307a66982 Mon Sep 17 00:00:00 2001 From: xin Date: Tue, 27 Sep 2022 23:43:25 +0800 Subject: [PATCH 03/16] fix bug no parameter --- ObservationSim/MockObject/MockObject.py | 6 ++++-- ObservationSim/ObservationSim.py | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 99a72ed..bbb1b6e 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -171,6 +171,8 @@ class MockObject(object): dy = y - y_nominal offset = galsim.PositionD(dx, dy) + real_wcs_local = self.real_wcs.local(real_pos) + for i in range(len(bandpass_list)): bandpass = bandpass_list[i] try: @@ -205,7 +207,7 @@ class MockObject(object): - stamp = star.drawImage(wcs=self.real_wcs, method='phot', offset=offset, save_photons=True) + 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 @@ -214,7 +216,7 @@ class MockObject(object): photons_list.append(photons) stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) - stamp.wcs = self.real_wcs.local(real_pos) + stamp.wcs = self.real_wcs_local 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) diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index a661b22..d80b0dc 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -273,7 +273,7 @@ class Observation(object): col_num=chip.colID, extName='raw') - pos_img, offset, local_wcs, real_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=h_ext) if pos_img.x == -1 or pos_img.y == -1: # Exclude object which is outside the chip area (after field distortion) # print("obj missed!!") -- GitLab From 9129df8f51b94535cea7ef704c95739af597acae Mon Sep 17 00:00:00 2001 From: xin Date: Wed, 28 Sep 2022 09:43:51 +0800 Subject: [PATCH 04/16] upload code to server --- ObservationSim/MockObject/MockObject.py | 64 ++++++++++++------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index bbb1b6e..ca9d41f 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -197,46 +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) - # 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) + 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 += x_nominal - photons.y += y_nominal + photons.x += self.x_nominal + photons.y += self.y_nominal photons_list.append(photons) - stamp = galsim.ImageF(int(xmax * 1.1), int(ymax * 1.1)) - stamp.wcs = self.real_wcs_local - 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: - chip.sensor.accumulate(photons_list[i], stamp) - else: - 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 = 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 += x_nominal + # photons.y += y_nominal + # photons_list.append(photons) - # 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: @@ -245,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 -- GitLab From 364e4c76df9e8b409ca017c83349c24726e27bab Mon Sep 17 00:00:00 2001 From: xin Date: Wed, 28 Sep 2022 11:38:52 +0800 Subject: [PATCH 05/16] local wcs sls --- ObservationSim/MockObject/MockObject.py | 113 +++++++++++++----------- 1 file changed, 63 insertions(+), 50 deletions(-) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index ca9d41f..153f097 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -197,29 +197,46 @@ 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) + 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: + chip.sensor.accumulate(photons_list[i], stamp) + else: + chip.sensor.accumulate(photons_list[i], stamp, resume=True) + chip.img[bounds] = stamp[bounds] - # 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 += x_nominal - # photons.y += y_nominal - # photons_list.append(photons) + 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 = 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 = 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: @@ -228,23 +245,6 @@ 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 @@ -252,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] @@ -267,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 @@ -301,9 +303,20 @@ class MockObject(object): normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T, names=('WAVELENGTH', 'FLUX')) + 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) + + real_wcs_local = self.real_wcs.local(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) @@ -312,8 +325,8 @@ class MockObject(object): star = galsim.Convolve(psf, star) starImg = star.drawImage(nx=100, ny=100, wcs=self.localWCS) - 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] @@ -325,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, @@ -335,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, @@ -350,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 -- GitLab From cd1351bc3191a615b4d8164413b91691d0cdc7b0 Mon Sep 17 00:00:00 2001 From: xin Date: Wed, 28 Sep 2022 12:50:47 +0800 Subject: [PATCH 06/16] upload server --- ObservationSim/MockObject/MockObject.py | 49 +++++++++---------------- 1 file changed, 18 insertions(+), 31 deletions(-) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 153f097..271bce8 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -252,7 +252,7 @@ class MockObject(object): return True, pos_shear - def addSLStoChipImage(self, sdp = None, chip = None, xOrderSigPlus = None, local_wcs=None): + def addSLStoChipImage(self, sdp = None, chip = None, xOrderSigPlus = None): spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] @@ -267,17 +267,15 @@ 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 = local_wcs + stamp.wcs = self.localWCS stamp.setOrigin(origin_order_x, origin_order_y) - bounds = stamp.bounds & galsim.BoundsI(0, chip.npix_x-1, 0, chip.npix_y-1) + bounds = stamp.bounds & chip.img.bounds 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 @@ -303,20 +301,9 @@ class MockObject(object): normalSED = Table(np.array([self.sed['WAVELENGTH'], self.sed['FLUX'] * sedNormFactor]).T, names=('WAVELENGTH', 'FLUX')) - 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) - - real_wcs_local = self.real_wcs.local(real_pos) - xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} - grating_split_pos_chip = 0 + grating_split_pos + grating_split_pos_chip = chip.bound.xmin + 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) @@ -325,8 +312,8 @@ class MockObject(object): star = galsim.Convolve(psf, star) starImg = star.drawImage(nx=100, ny=100, wcs=self.localWCS) - origin_star = [y_nominal - (starImg.center.y - starImg.ymin), - x_nominal - (starImg.center.x - starImg.xmin)] + origin_star = [self.y_nominal - (starImg.center.y - starImg.ymin), + self.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] @@ -338,8 +325,8 @@ class MockObject(object): subImg_p1 = starImg.array[:, 0:subSlitPos] star_p1 = galsim.Image(subImg_p1) origin_p1 = origin_star - xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0 - ycenter_p1 = y_nominal-0 + xcenter_p1 = min(self.x_nominal,grating_split_pos_chip-1) - chip.bound.xmin + ycenter_p1 = self.y_nominal-chip.bound.ymin sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, ycenter=ycenter_p1, origin=origin_p1, @@ -348,13 +335,13 @@ class MockObject(object): conf=chip.sls_conf[0], isAlongY=0) - self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus) 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(x_nominal, grating_split_pos_chip - 1) - 0 - ycenter_p2 = y_nominal - 0 + xcenter_p2 = max(self.x_nominal, grating_split_pos_chip - 1) - chip.bound.xmin + ycenter_p2 = self.y_nominal - chip.bound.ymin sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, ycenter=ycenter_p2, origin=origin_p2, @@ -363,27 +350,27 @@ class MockObject(object): conf=chip.sls_conf[1], isAlongY=0) - self.addSLStoChipImage(sdp=sdp_p2, chip=chip,xOrderSigPlus = xOrderSigPlus, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp_p2, chip=chip,xOrderSigPlus = xOrderSigPlus) del sdp_p1 del sdp_p2 elif grating_split_pos_chip<=gal_origin[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, - ycenter=y_nominal - 0, origin=origin_star, + sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin, + ycenter=self.y_nominal - chip.bound.ymin, 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, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus) del sdp elif grating_split_pos_chip>=gal_end[1]: - sdp = SpecDisperser(orig_img=starImg, xcenter=x_nominal - 0, - ycenter=y_nominal - 0, origin=origin_star, + sdp = SpecDisperser(orig_img=starImg, xcenter=self.x_nominal - chip.bound.xmin, + ycenter=self.y_nominal - chip.bound.ymin, 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, local_wcs=real_wcs_local) + self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus) del sdp del psf return True, pos_shear -- GitLab From 1b917613a8a277c3b171ef5e59af51c653a6651b Mon Sep 17 00:00:00 2001 From: xin Date: Wed, 28 Sep 2022 15:25:25 +0800 Subject: [PATCH 07/16] updata galaxy, mockobject wcs --- ObservationSim/Config/ChipOutput.py | 5 ++ ObservationSim/MockObject/Galaxy.py | 110 ++++++++++++++++++------ ObservationSim/MockObject/MockObject.py | 57 +++++++----- 3 files changed, 125 insertions(+), 47 deletions(-) diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index 36da180..fae606d 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -70,6 +70,11 @@ class ChipOutput(object): 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 + 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/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 5b2cb44..65da1f4 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 @@ -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 271bce8..5b6dffb 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -162,16 +162,16 @@ 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) + 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 = real_pos.x + 0.5, real_pos.y + 0.5 + 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(real_pos) + real_wcs_local = self.real_wcs.local(self.real_pos) for i in range(len(bandpass_list)): bandpass = bandpass_list[i] @@ -252,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] @@ -267,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 @@ -301,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] @@ -325,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, @@ -335,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, @@ -350,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 -- GitLab From a93829b08d9eb2dbef87a74c6e98846c78aca000 Mon Sep 17 00:00:00 2001 From: xin Date: Fri, 30 Sep 2022 10:00:08 +0800 Subject: [PATCH 08/16] galay spec bug, split pos bug --- ObservationSim/Config/ChipOutput.py | 10 +++++----- ObservationSim/MockObject/Galaxy.py | 2 +- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/ObservationSim/Config/ChipOutput.py b/ObservationSim/Config/ChipOutput.py index fae606d..d9253fd 100755 --- a/ObservationSim/Config/ChipOutput.py +++ b/ObservationSim/Config/ChipOutput.py @@ -68,13 +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 - print('-------------debug-----------------') - print('from global',ximg, yimg) + # 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) + # 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/MockObject/Galaxy.py b/ObservationSim/MockObject/Galaxy.py index 65da1f4..d026560 100755 --- a/ObservationSim/MockObject/Galaxy.py +++ b/ObservationSim/MockObject/Galaxy.py @@ -272,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] -- GitLab From 6904f2fcc511338b9be839cf902ceb20c8ed4b29 Mon Sep 17 00:00:00 2001 From: xin Date: Fri, 30 Sep 2022 10:35:36 +0800 Subject: [PATCH 09/16] fix pointing bug in wcs for astrometry --- ObservationSim/ObservationSim.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index d80b0dc..5e79ec3 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -258,11 +258,11 @@ class Observation(object): raise ValueError("Unknown shear input") # chip_output.logger.info("debug point #4") - h_ext = generateExtensionHeader( + header_wcs = generateExtensionHeader( xlen=chip.npix_x, ylen=chip.npix_y, - ra=pointing.ra, - dec=pointing.dec, + ra=ra_cen, + dec=dec_cen, pa=pointing.img_pa.deg, gain=chip.gain, readout=chip.read_noise, @@ -273,7 +273,7 @@ class Observation(object): col_num=chip.colID, extName='raw') - pos_img, offset, local_wcs, real_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False, img_header=h_ext) + 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!!") -- GitLab From adaef7ac07e9fc2df609586e346e82dbc5725f4d Mon Sep 17 00:00:00 2001 From: Fang Yuedong Date: Tue, 15 Nov 2022 02:45:18 +0000 Subject: [PATCH 10/16] fix PSF orientation --- ObservationSim/Config/Header/ImageHeader.py | 2 +- ObservationSim/MockObject/MockObject.py | 19 +++++++++++++++++++ ObservationSim/MockObject/_util.py | 5 ++--- ObservationSim/ObservationSim.py | 10 +++++++--- ObservationSim/PSF/PSFGauss.py | 10 +++++++--- ObservationSim/PSF/PSFInterp.py | 15 ++++++++++++--- config/config_NGP_dev.yaml | 8 ++++---- run_NGP.pbs | 4 ++-- 8 files changed, 54 insertions(+), 19 deletions(-) diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 0281dcd..86a46cb 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -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 diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index 5b6dffb..a8316ab 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -256,9 +256,28 @@ class MockObject(object): spec_orders = sdp.compute_spec_orders() for k, v in spec_orders.items(): img_s = v[0] + ######################################################### + # DEBUG + ######################################################### + # print("before convolveGaussXorders, img_s:", img_s) + nan_ids = np.isnan(img_s) + if img_s[nan_ids].shape[0] > 0: + # img_s[nan_ids] = 0 + print("DEBUG: before convolveGaussXorders specImg nan num is", img_s[nan_ids].shape[0]) + ######################################################### img_s, orig_off = convolveGaussXorders(img_s, xOrderSigPlus[k]) origin_order_x = v[1] - orig_off origin_order_y = v[2] - orig_off + + ######################################################### + # DEBUG + ######################################################### + # 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]) + ######################################################### specImg = galsim.ImageF(img_s) photons = galsim.PhotonArray.makeFromImage(specImg) photons.x += origin_order_x diff --git a/ObservationSim/MockObject/_util.py b/ObservationSim/MockObject/_util.py index ffb978b..c6cb573 100755 --- a/ObservationSim/MockObject/_util.py +++ b/ObservationSim/MockObject/_util.py @@ -545,10 +545,9 @@ def tag_sed(h5file, model_tag, teff=5000, logg=2, feh=0): flux = np.array(h5file["sed"][path][()]).ravel() return path, wave, flux -from astropy.modeling.models import Gaussian2D -from scipy import signal def convolveGaussXorders(img=None, sigma = 1): - + from astropy.modeling.models import Gaussian2D + from scipy import signal offset = int(np.ceil(sigma*10)) g_size = 2*offset+1 diff --git a/ObservationSim/ObservationSim.py b/ObservationSim/ObservationSim.py index 5e79ec3..4426177 100755 --- a/ObservationSim/ObservationSim.py +++ b/ObservationSim/ObservationSim.py @@ -7,6 +7,8 @@ import psutil from astropy.io import fits from datetime import datetime +import traceback + from ObservationSim.Config import config_dir, ChipOutput from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip @@ -77,7 +79,7 @@ class Observation(object): chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::') if self.config["psf_setting"]["psf_model"] == "Gauss": - psf_model = PSFGauss(chip=chip) + psf_model = PSFGauss(chip=chip, psfRa=self.config["psf_setting"]["psf_rcont"]) elif self.config["psf_setting"]["psf_model"] == "Interp": psf_model = PSFInterp(chip=chip, PSF_data_file=self.path_dict["psf_dir"]) else: @@ -207,7 +209,8 @@ class Observation(object): ) except Exception as e: - print(e) + # print(e) + traceback.print_exc() chip_output.logger.error(e) continue @@ -325,7 +328,8 @@ class Observation(object): # print("object omitted", flush=True) continue except Exception as e: - print(e) + # print(e) + traceback.print_exc() chip_output.logger.error(e) pass # Unload SED: diff --git a/ObservationSim/PSF/PSFGauss.py b/ObservationSim/PSF/PSFGauss.py index 5083880..7d09c3c 100644 --- a/ObservationSim/PSF/PSFGauss.py +++ b/ObservationSim/PSF/PSFGauss.py @@ -6,12 +6,16 @@ from scipy.interpolate import interp1d from ObservationSim.PSF.PSFModel import PSFModel class PSFGauss(PSFModel): - def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=0.15): + def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=None): self.pix_size = chip.pix_scale self.chip = chip - self.fwhm = fwhm + if psfRa is None: + self.fwhm = fwhm + self.sigGauss = 0.15 + else: + self.fwhm = self.fwhmGauss(r=psfRa) + self.sigGauss = psfRa # 80% light radius self.sigSpin = sigSpin - self.sigGauss = psfRa # 80% light radius self.psf = galsim.Gaussian(flux=1.0,fwhm=fwhm) def perfGauss(self, r, sig): diff --git a/ObservationSim/PSF/PSFInterp.py b/ObservationSim/PSF/PSFInterp.py index 058acca..45bf211 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -350,12 +350,21 @@ class PSFInterp(PSFModel): if galsimGSObject: imPSFt = np.zeros([257,257]) imPSFt[0:256, 0:256] = imPSF + # imPSFt[120:130, 0:256] = 1. img = galsim.ImageF(imPSFt, scale=pixSize) gsp = galsim.GSParams(folding_threshold=folding_threshold) - self.psf = galsim.InterpolatedImage(img, gsparams=gsp).rotate(pointing_pa*galsim.degrees) - - return self.PSFspin(x=px/0.01, y=py/0.01) + ############TEST: START + # Use sheared PSF to test the PSF orientation + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + ############TEST: END + self.psf = galsim.InterpolatedImage(img, gsparams=gsp) + wcs = chip.img.wcs.local(pos_img) + scale = galsim.PixelScale(0.074) + self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img)) + + # return self.PSFspin(x=px/0.01, y=py/0.01) + return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians) return imPSF def PSFspin(self, x, y): diff --git a/config/config_NGP_dev.yaml b/config/config_NGP_dev.yaml index 27bcee0..9f5ffdb 100644 --- a/config/config_NGP_dev.yaml +++ b/config/config_NGP_dev.yaml @@ -10,9 +10,9 @@ # Can add some of the command-line arguments here as well; # OK to pass either way or both, as long as they are consistent # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/" -work_dir: "/share/home/fangyuedong/csst-simulation/workplace/" +work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "header_TEST" +run_name: "PSF_TEST" # (Optional) a file of point list # if you just want to run default pointing: @@ -34,7 +34,7 @@ run_option: out_cat_only: NO # Only simulate stars? - star_only: NO + star_only: YES # Only simulate galaxies? galaxy_only: NO @@ -48,7 +48,7 @@ obs_setting: # "Photometric": simulate photometric chips only # "Spectroscopic": simulate slitless spectroscopic chips only # "All": simulate full focal plane - survey_type: "Photometric" + survey_type: "All" # Exposure time [seconds] exp_time: 150. diff --git a/run_NGP.pbs b/run_NGP.pbs index 2763d98..f4d181e 100755 --- a/run_NGP.pbs +++ b/run_NGP.pbs @@ -1,7 +1,7 @@ #!/bin/bash #PBS -N SIMS -#PBS -l nodes=wcl-1:ppn=32+wcl-6:ppn=32 +#PBS -l nodes=wcl-1:ppn=8+wcl-2:ppn=8+wcl-3:ppn=8+wcl-4:ppn=8+wcl-5:ppn=8+wcl-6:ppn=8 #PBS -u fangyuedong ###PBS -j oe @@ -12,5 +12,5 @@ date echo $NP # mpirun -np $NP python /public/home/fangyuedong/test/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/test/CSST/config -mpirun -np $NP python3 /share/home/fangyuedong/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/csst-simulation/config +mpirun -np $NP python3 /share/home/fangyuedong/test_psf_rot/csst-simulation/run_sim.py config_NGP_dev.yaml -c /share/home/fangyuedong/test_psf_rot/csst-simulation/config -- GitLab From 98d946960c9d15d1176deecfab7d6c85b08c4c8e Mon Sep 17 00:00:00 2001 From: xin Date: Wed, 16 Nov 2022 14:31:22 +0800 Subject: [PATCH 11/16] add the simulaiton of flat-field cube for slitless spectroscopic sim --- ObservationSim/Instrument/Chip/Chip.py | 19 +++ .../Instrument/data/flatCube/flat_cube.fits | Bin 0 -> 2722648320 bytes ObservationSim/MockObject/Galaxy.py | 15 ++- ObservationSim/MockObject/MockObject.py | 97 +++++++------- ObservationSim/MockObject/SkybackgroundMap.py | 14 +- .../MockObject/SpecDisperser/SpecDisperser.py | 120 +++++++++++------- .../SpecDisperser/disperse_c/disperse.pyx | 43 +++++-- ObservationSim/ObservationSim.py | 20 +-- setup.py | 1 + 9 files changed, 209 insertions(+), 120 deletions(-) create mode 100644 ObservationSim/Instrument/data/flatCube/flat_cube.fits diff --git a/ObservationSim/Instrument/Chip/Chip.py b/ObservationSim/Instrument/Chip/Chip.py index 69d0c10..96aa78f 100755 --- a/ObservationSim/Instrument/Chip/Chip.py +++ b/ObservationSim/Instrument/Chip/Chip.py @@ -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 + diff --git a/ObservationSim/Instrument/data/flatCube/flat_cube.fits b/ObservationSim/Instrument/data/flatCube/flat_cube.fits new file mode 100644 index 0000000000000000000000000000000000000000..cca7b145abfca4b8f012a728f588d318d3b734f4 GIT binary patch