From b546bca9092421c0ff976a85025dadd17e6da906 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Sun, 30 Oct 2022 08:58:36 +0800 Subject: [PATCH 1/5] test WCS using sheared PSF --- Catalog/NGPCatalog.py | 3 ++- ObservationSim/PSF/PSFInterp.py | 13 ++++++++++--- config/config_NGP_dev.yaml | 8 ++++---- run_NGP.pbs | 4 ++-- 4 files changed, 18 insertions(+), 10 deletions(-) diff --git a/Catalog/NGPCatalog.py b/Catalog/NGPCatalog.py index 15e195c..4e6c4ab 100644 --- a/Catalog/NGPCatalog.py +++ b/Catalog/NGPCatalog.py @@ -244,7 +244,8 @@ class NGPCatalog(CatalogBase): param['parallax'] = parallax_arr[istars] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue - param['mag_use_normal'] = stars['app_sdss_g'][istars] + # param['mag_use_normal'] = stars['app_sdss_g'][istars] + param['mag_use_normal'] = 18. if param['mag_use_normal'] >= 26.5: continue self.ids += 1 diff --git a/ObservationSim/PSF/PSFInterp.py b/ObservationSim/PSF/PSFInterp.py index 058acca..fe90489 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -350,12 +350,19 @@ 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) + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).rotate(pointing_pa*galsim.degrees) + self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp, wcs=chip.img.wcs.local(pos_img)).shear(g1=0.8, g2=0.) + 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..7a11526 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: "sheared_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 f1d83af8a781860275404d38a5fd427dad797699 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Wed, 2 Nov 2022 04:42:19 +0800 Subject: [PATCH 2/5] projected interpolated PSF back to sky coordinates using focal plane wcs before convolving with objects. For testing --- ObservationSim/PSF/PSFInterp.py | 4 ++-- config/config_NGP_dev.yaml | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ObservationSim/PSF/PSFInterp.py b/ObservationSim/PSF/PSFInterp.py index fe90489..b35fb40 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -355,8 +355,8 @@ class PSFInterp(PSFModel): 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) - self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) - # self.psf = galsim.InterpolatedImage(img, gsparams=gsp, wcs=chip.img.wcs.local(pos_img)).shear(g1=0.8, g2=0.) + # self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.) + 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)) diff --git a/config/config_NGP_dev.yaml b/config/config_NGP_dev.yaml index 7a11526..9f5ffdb 100644 --- a/config/config_NGP_dev.yaml +++ b/config/config_NGP_dev.yaml @@ -12,7 +12,7 @@ # work_dir: "/public/home/fangyuedong/sim_code_release/CSST/test/" work_dir: "/share/home/fangyuedong/test_psf_rot/csst-simulation/workplace/" data_dir: "/share/simudata/CSSOSDataProductsSims/data/" -run_name: "sheared_PSF_TEST" +run_name: "PSF_TEST" # (Optional) a file of point list # if you just want to run default pointing: -- GitLab From 7e11e7785c0cf5a058a36abb0e18fd4e330c52d1 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 15 Nov 2022 10:23:52 +0800 Subject: [PATCH 3/5] fix the PSF orientation issue by projecting the interpolated PSF to sky coordinates via focal plane wcs, before convolving with objects --- 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 | 4 +++- 6 files changed, 39 insertions(+), 11 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..e66f925 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 b35fb40..45bf211 100644 --- a/ObservationSim/PSF/PSFInterp.py +++ b/ObservationSim/PSF/PSFInterp.py @@ -354,8 +354,10 @@ class PSFInterp(PSFModel): 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) + ############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) -- GitLab From 343a003a8f9988d272c06167748d63f75ea521fc Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 15 Nov 2022 10:31:13 +0800 Subject: [PATCH 4/5] change the magnitudes back --- Catalog/NGPCatalog.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Catalog/NGPCatalog.py b/Catalog/NGPCatalog.py index 4e6c4ab..15e195c 100644 --- a/Catalog/NGPCatalog.py +++ b/Catalog/NGPCatalog.py @@ -244,8 +244,7 @@ class NGPCatalog(CatalogBase): param['parallax'] = parallax_arr[istars] if not self.chip.isContainObj(ra_obj=param['ra'], dec_obj=param['dec'], margin=200): continue - # param['mag_use_normal'] = stars['app_sdss_g'][istars] - param['mag_use_normal'] = 18. + param['mag_use_normal'] = stars['app_sdss_g'][istars] if param['mag_use_normal'] >= 26.5: continue self.ids += 1 -- GitLab From c4cca012d44acf221c6fdbd09cfffa7f95e3e747 Mon Sep 17 00:00:00 2001 From: fangyuedong Date: Tue, 15 Nov 2022 10:33:14 +0800 Subject: [PATCH 5/5] remove debugging prints --- ObservationSim/MockObject/MockObject.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ObservationSim/MockObject/MockObject.py b/ObservationSim/MockObject/MockObject.py index e66f925..a8316ab 100755 --- a/ObservationSim/MockObject/MockObject.py +++ b/ObservationSim/MockObject/MockObject.py @@ -272,7 +272,7 @@ class MockObject(object): ######################################################### # DEBUG ######################################################### - print("DEBUG: orig_off is", orig_off) + # 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 -- GitLab