From d4b41c3682cf9d09d48ac2077a75829655e153b5 Mon Sep 17 00:00:00 2001 From: weichengliang Date: Mon, 8 Jan 2024 16:07:24 +0800 Subject: [PATCH] update fits module --- Catalog/C6_fits_Catalog.py | 26 ++- ObservationSim/MockObject/Stamp.py | 266 +++++++++++++++++++++++------ config/config_C6_fits.yaml | 88 ++++++---- 3 files changed, 282 insertions(+), 98 deletions(-) diff --git a/Catalog/C6_fits_Catalog.py b/Catalog/C6_fits_Catalog.py index 80af1cb..0f03dc4 100644 --- a/Catalog/C6_fits_Catalog.py +++ b/Catalog/C6_fits_Catalog.py @@ -60,6 +60,9 @@ class Catalog(CatalogBase): with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path: self.normF_star = Table.read(str(filter_path)) + with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path: + self.normF_galaxy = Table.read(str(filter_path)) + self.config = config self.chip = chip @@ -138,8 +141,8 @@ class Catalog(CatalogBase): return None ###mock_stamp_START elif obj.type == "stamp": - #return self.normF_galaxy ###normalize_filter for stamp - return None + return self.normF_galaxy ###normalize_filter for stamp + #return None ###mock_stamp_END else: return None @@ -197,8 +200,8 @@ class Catalog(CatalogBase): for igals in range(ngals): # # (TEST) - # if igals > 100: - # break + if igals > 2000: + break param = self.initialize_param() param['ra'] = ra_arr[igals] @@ -321,8 +324,8 @@ class Catalog(CatalogBase): ) for istars in range(nstars): # # (TEST) - # if istars > 100: - # break + if istars > 100: + break param = self.initialize_param() param['ra'] = ra_arr[istars] @@ -389,6 +392,9 @@ class Catalog(CatalogBase): input_time_str=time_str ) for iAGNs in range(nAGNs): + if iAGNs > 100: + break + param = self.initialize_param() param['ra'] = ra_arr[iAGNs] param['dec'] = dec_arr[iAGNs] @@ -425,13 +431,15 @@ class Catalog(CatalogBase): self.ud = galsim.UniformDeviate(pix_id) for istamp in range(nstamps): + print('DEBUG:::istamp=', istamp) + fitsfile = os.path.join(self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8')) hdu=fitsio.open(fitsfile) param = self.initialize_param() param['id'] = hdu[0].header['index'] #istamp param['star'] = 3 # Stamp type in .cat file - param['lensGalaxyID'] = hdu[0].header['lensGID'] + ###param['lensGalaxyID'] = hdu[0].header['lensGID'] param['ra'] = hdu[0].header['ra'] param['dec']= hdu[0].header['dec'] param['pixScale']= hdu[0].header['pixScale'] @@ -440,9 +448,9 @@ class Catalog(CatalogBase): #param['PA']= hdu[0].header['PA'] #param['bfrac']= hdu[0].header['bfrac'] #param['z']= hdu[0].header['z'] - param['mag_use_normal'] = 22 #hdu[0].header['m_normal'] #gals['mag_true_g_lsst'] + param['mag_use_normal'] = 20 #hdu[0].header['m_normal'] #gals['mag_true_g_lsst'] - assert(stamps['lensGID'][istamp] == param['lensGalaxyID']) + ###assert(stamps['lensGID'][istamp] == param['lensGalaxyID']) # Apply astrometric modeling # in C3 case only aberration diff --git a/ObservationSim/MockObject/Stamp.py b/ObservationSim/MockObject/Stamp.py index 3e985d6..1da5cd4 100644 --- a/ObservationSim/MockObject/Stamp.py +++ b/ObservationSim/MockObject/Stamp.py @@ -1,22 +1,17 @@ import os, sys -import random import numpy as np +import galsim import astropy.constants as cons from astropy.table import Table from scipy import interpolate -import astropy.io.fits as fitsio - -import galsim -import gc from ObservationSim.MockObject.MockObject import MockObject - -from ObservationSim.MockObject._util import magToFlux,VC_A +from ObservationSim.MockObject.SpecDisperser import SpecDisperser from ObservationSim.MockObject._util import eObs, integrate_sed_bandpass, getNormFactorForSpecWithABMAG, getObservedSED, getABMAG,convolveGaussXorders class Stamp(MockObject): - def __init__(self, param): - super().__init__(param) + def __init__(self, param, logger=None): + super().__init__(param, logger=logger) def unload_SED(self): """(Test) free up SED memory @@ -35,9 +30,9 @@ class Stamp(MockObject): self.logger.error(e) return False - nphotons_sum = 0 - photons_list = [] - xmax, ymax = 0, 0 + #nphotons_sum = 0 + #photons_list = [] + #xmax, ymax = 0, 0 if self.getMagFilter(filt) <= 15: folding_threshold = 5.e-4 @@ -46,7 +41,7 @@ class Stamp(MockObject): 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) + img_real_wcs=self.chip_wcs) x, y = self.real_pos.x + 0.5, self.real_pos.y + 0.5 x_nominal = int(np.floor(x + 0.5)) @@ -55,72 +50,237 @@ class Stamp(MockObject): dy = y - y_nominal offset = galsim.PositionD(dx, dy) - real_wcs_local = self.real_wcs.local(self.real_pos) + chip_wcs_local = self.chip_wcs.local(self.real_pos) + is_updated = 0 + + if fd_shear: + g1 += fd_shear.g1 + g2 += fd_shear.g2 + gal_shear = galsim.Shear(g1=g1, g2=g2) for i in range(len(bandpass_list)): bandpass = bandpass_list[i] - try: sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass) except Exception as e: print(e) - self.logger.error(e) - # return False + if self.logger: + self.logger.error(e) continue - ratio = sub/full if not (ratio == -1 or (ratio != ratio)): nphotons = ratio * nphotons_tot else: - # return False continue - nphotons_sum += nphotons + #nphotons_sum += nphotons psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) _gal = self.param['image'] - galIm = galsim.ImageF(_gal, scale=self.param['pixScale']) - gal = galsim.InterpolatedImage(galIm) - gal = gal.withFlux(nphotons) - #gal_shear = galsim.Shear(g1=g1, g2=g2) - #gal = gal.shear(gal_shear) + galImg= galsim.ImageF(_gal, scale=self.param['pixScale']) + gal_temp= galsim.InterpolatedImage(galImg) + gal_temp= gal_temp.shear(gal_shear) + gal_temp= gal_temp.withFlux(nphotons) - gal = galsim.Convolve(psf, gal) - if fd_shear is not None: - gal = gal.shear(fd_shear) + gal_temp= galsim.Convolve(psf, gal_temp) - stamp = gal.drawImage(wcs=real_wcs_local, method='phot', offset=self.offset, save_photons=True) - - xmax = max(xmax, stamp.xmax - stamp.xmin) - ymax = max(ymax, stamp.ymax - stamp.ymin) - - photons = stamp.photons - photons.x += x_nominal - photons.y += y_nominal - photons_list.append(photons) - del gal + if i == 0: + gal = gal_temp + else: + gal = gal + gal_temp - # print('xmax = %d, ymax = %d '%(xmax, ymax)) + stamp = gal.drawImage(wcs=chip_wcs_local, offset=offset) + if np.sum(np.isnan(stamp.array)) > 0: + # ERROR happens + return 2, pos_shear - 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) if bounds.area() > 0: 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] + is_updated = 1 + chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + del stamp - chip.img[bounds] = stamp[bounds] + if is_updated == 0: + print("fits obj %s missed"%(self.id)) + if self.logger: + self.logger.info("fits obj %s missed"%(self.id)) + return 0, pos_shear - chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin) + return 1, pos_shear + + + def drawObj_slitless(self, tel, pos_img, psf_model, bandpass_list, filt, chip, nphotons_tot=None, g1=0, g2=0, + exptime=150., normFilter=None, grating_split_pos=3685, fd_shear=None): + if normFilter is not None: + norm_thr_rang_ids = normFilter['SENSITIVITY'] > 0.001 + sedNormFactor = getNormFactorForSpecWithABMAG(ABMag=self.param['mag_use_normal'], spectrum=self.sed, + norm_thr=normFilter, + sWave=np.floor(normFilter[norm_thr_rang_ids][0][0]), + eWave=np.ceil(normFilter[norm_thr_rang_ids][-1][0])) + if sedNormFactor == 0: + return 2, None + else: + sedNormFactor = 1. + 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.chip_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) + + chip_wcs_local = self.chip_wcs.local(self.real_pos) + + + if self.getMagFilter(filt) <= 15: + folding_threshold = 5.e-4 + else: + folding_threshold = 5.e-3 + gsp = galsim.GSParams(folding_threshold=folding_threshold) + # nphotons_sum = 0 + + flat_cube = chip.flat_cube + + xOrderSigPlus = {'A':1.3909419820029296,'B':1.4760376591236062,'C':4.035447379743442,'D':5.5684364343742825,'E':16.260021029735388} + grating_split_pos_chip = 0 + grating_split_pos + + branges = np.zeros([len(bandpass_list), 2]) + + # print(hasattr(psf_model, 'bandranges')) + + if hasattr(psf_model, 'bandranges'): + if psf_model.bandranges is None: + return 2, None + if len(psf_model.bandranges) != len(bandpass_list): + return 2, None + branges = psf_model.bandranges + else: + for i in range(len(bandpass_list)): + branges[i, 0] = bandpass_list[i].blue_limit * 10 + branges[i, 1] = bandpass_list[i].red_limit * 10 + + for i in range(len(bandpass_list)): + # bandpass = bandpass_list[i] + brange = branges[i] + + # psf, pos_shear = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, folding_threshold=folding_threshold) + + _gal = self.param['image'] + galImg= galsim.ImageF(_gal, scale=self.param['pixScale']) + gal = galsim.InterpolatedImage(galImg) + + # (TEST) Random knots + # knots = galsim.RandomKnots(npoints=100, profile=disk) + # kfrac = np.random.random()*(1.0 - self.bfrac) + # gal = self.bfrac * bulge + (1.0 - self.bfrac - kfrac) * disk + kfrac * knots + + gal = gal.withFlux(tel.pupil_area * exptime) + if fd_shear: + g1 += fd_shear.g1 + g2 += fd_shear.g2 + gal_shear = galsim.Shear(g1=g1, g2=g2) + gal = gal.shear(gal_shear) + # gal = galsim.Convolve(psf, gal) + + # if not big_galaxy: # Not apply PSF for very big galaxy + # gal = galsim.Convolve(psf, gal) + # # if fd_shear is not None: + # # gal = gal.shear(fd_shear) + + starImg = gal.drawImage(wcs=chip_wcs_local, offset=offset,method = 'real_space') + + origin_star = [y_nominal - (starImg.center.y - starImg.ymin), + x_nominal - (starImg.center.x - starImg.xmin)] + starImg.setOrigin(0, 0) + 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] + + if gal_origin[1] < grating_split_pos_chip < gal_end[1]: + subSlitPos = int(grating_split_pos_chip - gal_origin[1] + 1) + ## part img disperse + + subImg_p1 = starImg.array[:, 0:subSlitPos] + star_p1 = galsim.Image(subImg_p1) + star_p1.setOrigin(0, 0) + origin_p1 = origin_star + xcenter_p1 = min(x_nominal,grating_split_pos_chip-1) - 0 + ycenter_p1 = y_nominal-0 + + sdp_p1 = SpecDisperser(orig_img=star_p1, xcenter=xcenter_p1, + ycenter=ycenter_p1, origin=origin_p1, + tar_spec=normalSED, + band_start=brange[0], band_end=brange[1], + conf=chip.sls_conf[0], + isAlongY=0, + flat_cube=flat_cube) + + # self.addSLStoChipImage(sdp=sdp_p1, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p1, chip=chip, pos_img_local=[xcenter_p1, ycenter_p1], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) + + subImg_p2 = starImg.array[:, subSlitPos+1:starImg.array.shape[1]] + star_p2 = galsim.Image(subImg_p2) + star_p2.setOrigin(0, 0) + 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 + + sdp_p2 = SpecDisperser(orig_img=star_p2, xcenter=xcenter_p2, + ycenter=ycenter_p2, origin=origin_p2, + tar_spec=normalSED, + band_start=brange[0], band_end=brange[1], + conf=chip.sls_conf[1], + isAlongY=0, + flat_cube=flat_cube) + + # self.addSLStoChipImage(sdp=sdp_p2, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp_p2, chip=chip, pos_img_local=[xcenter_p2, ycenter_p2], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) + + 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, + tar_spec=normalSED, + band_start=brange[0], band_end=brange[1], + conf=chip.sls_conf[1], + isAlongY=0, + flat_cube=flat_cube) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) + 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, + tar_spec=normalSED, + band_start=brange[0], band_end=brange[1], + conf=chip.sls_conf[0], + isAlongY=0, + flat_cube=flat_cube) + # self.addSLStoChipImage(sdp=sdp, chip=chip, xOrderSigPlus = xOrderSigPlus, local_wcs=chip_wcs_local) + pos_shear = self.addSLStoChipImageWithPSF(sdp=sdp, chip=chip, pos_img_local=[x_nominal, y_nominal], + psf_model=psf_model, bandNo=i + 1, + grating_split_pos=grating_split_pos, + local_wcs=chip_wcs_local, pos_img = pos_img) + del sdp - del photons_list - del stamp - gc.collect() - return True, pos_shear + # print(self.y_nominal, starImg.center.y, starImg.ymin) + # del psf + return 1, pos_shear diff --git a/config/config_C6_fits.yaml b/config/config_C6_fits.yaml index 074a627..1def928 100644 --- a/config/config_C6_fits.yaml +++ b/config/config_C6_fits.yaml @@ -9,9 +9,11 @@ # Base diretories and naming setup # 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: "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C6/" +work_dir: "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C8/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "C6_fits_testRun" +run_name: "testRun_fits" +project_cycle: 6 +run_counter: 1 # Whether to use MPI run_option: @@ -46,7 +48,7 @@ catalog_options: #stamp_SED: # Only simulate stars? - star_only: NO + star_only: YES # Only simulate galaxies? galaxy_only: NO @@ -63,13 +65,23 @@ catalog_options: # Observation setting ############################################### obs_setting: + # Options for survey types: # "Photometric": simulate photometric chips only # "Spectroscopic": simulate slitless spectroscopic chips only # "FGS": simulate FGS chips only (31-42) # "All": simulate full focal plane - survey_type: "Photometric" - + # "CALIBRATION": falt, bias, dark with or without postflash + survey_type: "All" + # "LED": ['LED1','LED2','LED3','LED4','LED5','LED6','LED7','LED8','LED9','LED10','LED11','LED12','LED13','LED14'] or null + # 'LED1': '275', 'LED2': '310', 'LED3': '430', 'LED4': '505', 'LED5': '545', 'LED6': '590', 'LED7': '670', + # 'LED8': '760', 'LED9': '940', 'LED10': '940', 'LED11': '1050', 'LED12': '1050','LED13': '340', 'LED14': '365' + # LED_TYPE: ['LED5'] + # LED_TIME: [1.] + # # unit e- ,flat level + # FLAT_LEVEL: 20000 + # FLAT_LEVEL_FIL: 'g' + # Exposure time [seconds] exp_time: 150. @@ -79,19 +91,17 @@ obs_setting: # Default Pointing [degrees] # Note: NOT valid when a pointing list file is specified - ra_center: 244.972773 - dec_center: 39.895901 + ra_center: 192.8595 + dec_center: 27.1283 # Image rotation [degree] - image_rot: 109.59 + image_rot: -113.4333 - # (Optional) a file of point list + # (Optional) a file of point list # if you just want to run default pointing: # - pointing_dir: null # - pointing_file: null pointing_dir: "/share/simudata/CSSOSDataProductsSims/data/" pointing_file: "pointing_radec_246.5_40.dat" - # pointing_dir: null - # pointing_file: null # Number of calibration pointings np_cal: 0 @@ -101,21 +111,25 @@ obs_setting: # - run all pointings: null # Note: only valid when a pointing list is specified run_pointings: [0] - + # Run specific chip(s): # - give a list of indexes of chips: [ip_1, ip_2...] # - run all chips: null # Note: for all pointings - run_chips: [7] + run_chips: [1, 2, 3, 7, 8, 9] # Whether to enable astrometric modeling - enable_astrometric_model: False + enable_astrometric_model: True + + # Whether to enable straylight model + enable_straylight_model: True # Cut by saturation magnitude in which band? cut_in_band: "z" # saturation magnitude margin - mag_sat_margin: -2.5 + # mag_sat_margin: -2.5 + mag_sat_margin: -15. # limiting magnitude margin mag_lim_margin: +1.0 @@ -138,7 +152,7 @@ psf_setting: # path to PSF data # NOTE: only valid for "Interp" PSF psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1" - + psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/" ############################################### # Shear setting ############################################### @@ -160,25 +174,27 @@ ins_effects: # switches # Note: bias_16channel, gain_16channel, and shutter_effect # is currently not applicable to "FGS" observations - field_dist: OFF # Whether to add field distortions - add_back: OFF # Whether to add sky background - add_dark: OFF # Whether to add dark noise - add_readout: OFF # Whether to add read-out (Gaussian) noise - add_bias: OFF # Whether to add bias-level to images - bias_16channel: OFF # Whether to add different biases for 16 channels - gain_16channel: OFF # Whether to make different gains for 16 channels - shutter_effect: OFF # Whether to add shutter effect - flat_fielding: OFF # Whether to add flat-fielding effect - prnu_effect: OFF # Whether to add PRNU effect - non_linear: OFF # Whether to add non-linearity - cosmic_ray: OFF # Whether to add cosmic-ray - cray_differ: OFF # Whether to generate different cosmic ray maps CAL and MS output - cte_trail: OFF # Whether to simulate CTE trails - saturbloom: OFF # Whether to simulate Saturation & Blooming - add_badcolumns: OFF # Whether to add bad columns - add_hotpixels: OFF # Whether to add hot pixels - add_deadpixels: OFF # Whether to add dead(dark) pixels - bright_fatter: OFF # Whether to simulate Brighter-Fatter (also diffusion) effect + field_dist: YES # Whether to add field distortions + add_back: YES # Whether to add sky background + add_dark: YES # Whether to add dark noise + add_readout: YES # Whether to add read-out (Gaussian) noise + add_bias: YES # Whether to add bias-level to images + bias_16channel: YES # Whether to add different biases for 16 channels + gain_16channel: YES # Whether to make different gains for 16 channels + shutter_effect: YES # Whether to add shutter effect + flat_fielding: YES # Whether to add flat-fielding effect + prnu_effect: YES # Whether to add PRNU effect + non_linear: YES # Whether to add non-linearity + cosmic_ray: YES # Whether to add cosmic-ray + cray_differ: YES # Whether to generate different cosmic ray maps CAL and MS output + cte_trail: NO # Whether to simulate CTE trails, CTI_lgl_v0.3.tar.gz + saturbloom: YES # Whether to simulate Saturation & Blooming + add_badcolumns: YES # Whether to add bad columns + add_hotpixels: YES # Whether to add hot pixels + add_deadpixels: YES # Whether to add dead(dark) pixels + bright_fatter: YES # Whether to simulate Brighter-Fatter (also diffusion) effect + add_prescan: NO # Whether to add pre/over-scan + format_output: NO ##1*16 output # Values: # default values have been defined individually for each chip in: @@ -220,4 +236,4 @@ random_seeds: seed_badcolumns: 20240309 # Seed for bad columns seed_defective: 20210304 # Seed for defective (bad) pixels seed_readout: 20210601 # Seed for read-out gaussian noise -... \ No newline at end of file +... -- GitLab