Merge branch 'test_sls_rot_wcs' into 'sls_rot_wcs'

fix PSF orientation

See merge request csst_sim/csst-simulation!9

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

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

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
 

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:

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

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

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.

run_NGP.pbs

 #!/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