Merge branch 'customizable_catalog'

ObservationSim/MockObject/MockObject.py

@@ -44,6 +44,7 @@ class MockObject(object):
         self.hlr_bulge = self.param["hlr_bulge"]
         self.e1_disk, self.e2_disk = 0., 0.
         self.e1_bulge, self.e2_bulge = 0., 0.
+        self.additional_output_str = ""
 
         self.logger = logger
 
@@ -126,6 +127,7 @@ class MockObject(object):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            self.logger.error(e)
             return False
 
         nphotons_sum = 0
@@ -145,6 +147,7 @@ class MockObject(object):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                self.logger.error(e)
                 # return False
                 continue
         

ObservationSim/MockObject/Star.py

@@ -37,6 +37,7 @@ class Star(MockObject):
             full = integrate_sed_bandpass(sed=self.sed, bandpass=filt.bandpass_full)
         except Exception as e:
             print(e)
+            self.logger.error(e)
             return -1
 
         for i in range(len(bandpass_list)):
@@ -46,6 +47,7 @@ class Star(MockObject):
                 sub = integrate_sed_bandpass(sed=self.sed, bandpass=bandpass)
             except Exception as e:
                 print(e)
+                self.logger.error(e)
                 return -1
         
             ratio = sub/full

ObservationSim/ObservationSim.py

@@ -25,6 +25,7 @@ class Observation(object):
         self.filter_param = FilterParam() 
         self.chip_list = []
         self.filter_list = []
+        self.all_filter = []
         self.Catalog = Catalog
 
         # if we want to apply field distortion?
@@ -37,8 +38,6 @@ class Observation(object):
         nchips = self.focal_plane.nchip_x*self.focal_plane.nchip_y
         for i in range(nchips):
             chipID = i + 1
-            if self.focal_plane.isIgnored(chipID=chipID):
-                continue
 
             # Make Chip & Filter lists
             chip = Chip(
@@ -48,8 +47,10 @@ class Observation(object):
             filt = Filter(filter_id=filter_id, 
                 filter_type=filter_type, 
                 filter_param=self.filter_param)
-            self.chip_list.append(chip)
-            self.filter_list.append(filt)
+            if not self.focal_plane.isIgnored(chipID=chipID):
+                self.chip_list.append(chip)
+                self.filter_list.append(filt)
+            self.all_filter.append(filt)
 
         # Read catalog and shear(s)
         self.g1_field, self.g2_field, self.nshear = get_shear_field(config=self.config)
@@ -75,9 +76,6 @@ class Observation(object):
         chip_output.logger.info('Chip : %d' % chip.chipID)
         chip_output.logger.info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
 
-        # Update the limiting magnitude using exposure time in pointing
-        filt.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)
-
         if self.config["psf_setting"]["psf_model"] == "Gauss":
             psf_model = PSFGauss(chip=chip)
         elif self.config["psf_setting"]["psf_model"] == "Interp":
@@ -160,9 +158,19 @@ class Observation(object):
 
         if pointing.pointing_type == 'MS':
             # Load catalogues and templates
-            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, logger=chip_output.logger)
+            self.cat = self.Catalog(config=self.config, chip=chip, pointing=pointing, cat_dir=cat_dir, sed_dir=sed_dir, chip_output=chip_output)
+            chip_output.create_output_file()
             self.nobj = len(self.cat.objs)
 
+            for ifilt in range(len(self.all_filter)):
+                temp_filter = self.all_filter[ifilt]
+                # Update the limiting magnitude using exposure time in pointing
+                temp_filter.update_limit_saturation_mags(exptime=pointing.exp_time, chip=chip)
+
+                # Select cutting band filter for saturation/limiting magnitude
+                if temp_filter.filter_type.lower() == self.config["obs_setting"]["cut_in_band"].lower():
+                    cut_filter = temp_filter
+
             # Loop over objects
             missed_obj = 0
             bright_obj = 0
@@ -181,7 +189,7 @@ class Observation(object):
                 elif obj.type == 'quasar' and self.config["run_option"]["star_only"]:
                     continue
 
-                # load SED
+                # load and convert SED; also caculate object's magnitude in all CSST bands
                 try:
                     sed_data = self.cat.load_sed(obj)
                     norm_filt = self.cat.load_norm_filt(obj)
@@ -191,26 +199,39 @@ class Observation(object):
                         target_filt=filt, 
                         norm_filt=norm_filt,
                     )
+                    _, obj.param["mag_%s"%cut_filter.filter_type] = self.cat.convert_sed(
+                        mag=obj.param["mag_use_normal"],
+                        sed=sed_data,
+                        target_filt=cut_filter, 
+                        norm_filt=norm_filt,
+                    )
+
                 except Exception as e:
-                    # print(e)
+                    print(e)
                     chip_output.logger.error(e)
                     continue
+                
+                # chip_output.logger.info("debug point #1")
 
                 # Exclude very bright/dim objects (for now)
                 # if filt.is_too_bright(mag=obj.getMagFilter(filt)):
-                if filt.is_too_bright(mag=obj.mag_use_normal):
+                # if filt.is_too_bright(mag=obj.mag_use_normal):
+                if cut_filter.is_too_bright(mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()]):
                     # print("obj too birght!!", flush=True)
                     if obj.type != 'galaxy':
                         bright_obj += 1
                         obj.unload_SED()
                         continue
                 if filt.is_too_dim(mag=obj.getMagFilter(filt)):
+                # if cut_filter.is_too_dim(mag=obj.param["mag_%s"%self.config["obs_setting"]["cut_in_band"].lower()]):
                     # print("obj too dim!!", flush=True)
                     dim_obj += 1
                     obj.unload_SED()
                     # print(obj.getMagFilter(filt))
                     continue
 
+                # chip_output.logger.info("debug point #2")
+
                 if self.config["shear_setting"]["shear_type"] == "constant":
                     if obj.type == 'star':
                         obj.g1, obj.g2 = 0., 0.
@@ -224,11 +245,15 @@ class Observation(object):
                         # print("failed to load external shear.")
                         chip_output.logger.error("failed to load external shear.")
                         pass
+
+                    # chip_output.logger.info("debug point #3")
                 elif self.config["shear_setting"]["shear_type"] == "catalog":
                     pass
                 else:
                     chip_output.logger.error("Unknown shear input")
                     raise ValueError("Unknown shear input")
+                
+                # chip_output.logger.info("debug point #4")
 
                 pos_img, offset, local_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:
@@ -238,8 +263,12 @@ class Observation(object):
                     obj.unload_SED()
                     continue
 
+                # chip_output.logger.info("debug point #5")
+
                 # Draw object & update output catalog
                 try:
+                    # chip_output.logger.info("debug point #6")
+                    # chip_output.logger.info("current filter type: %s"%filt.filter_type)
                     if self.config["run_option"]["out_cat_only"]:
                         isUpdated = True
                         pos_shear = 0.
@@ -268,6 +297,7 @@ class Observation(object):
                             exptime=pointing.exp_time,
                             normFilter=norm_filt,
                             )
+                    # chip_output.logger.info("debug point #7")
                     if isUpdated:
                         # TODO: add up stats
                         # print("updating output catalog...")
@@ -277,7 +307,7 @@ class Observation(object):
                         # print("object omitted", flush=True)
                         continue
                 except Exception as e:
-                    # print(e)
+                    print(e)
                     chip_output.logger.error(e)
                     pass
                 # Unload SED:

ObservationSim/_util.py

@@ -13,7 +13,7 @@ def parse_args():
     '''
     parser = argparse.ArgumentParser()
     parser.add_argument('config_file', help='.yaml config file for simulation settings.')
-    parser.add_argument('-c', '--config_dir', help='Directory that houses the ,yaml config file.')
+    parser.add_argument('-c', '--config_dir', help='Directory that houses the .yaml config file.')
     parser.add_argument('-d', '--data_dir', help='Directory that houses the input data.')
     parser.add_argument('-w', '--work_dir', help='The path for output.')
     return parser.parse_args()

README.md

@@ -2,6 +2,7 @@
 
 ## 重要更新或问题修复：
 
+* 2022.04.13: 修复了天测模块中的调用错误
 * 2022.04.08: 修复了仅输出星表选项("out_cat_only")存在的问题
 
 ## 使用方法及相关说明

config/config_C3.yaml

@@ -84,6 +84,9 @@ obs_setting:
   # astrometric_lib: "libshao.so"
   enable_astrometric_model: True
 
+  # Cut by saturation/limiting magnitude in which band?
+  cut_in_band: "g"
+
 ###############################################
 # Input path setting
 ###############################################

config/config_NGP.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: "/public/home/fangyuedong/temp/CSST/workplace/"
+work_dir: "/public/home/fangyuedong/test/CSST/workplace/"
 data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
-run_name: "NGP_Astrometry-on"
+run_name: "TEST_Speed"
 
 # (Optional) a file of point list 
 # if you just want to run default pointing:
@@ -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.
@@ -73,18 +73,21 @@ obs_setting:
   # - give a list of indexes of pointings: [ip_1, ip_2...]
   # - run all pointings: null
   # Note: only valid when a pointing list is specified
-  run_pointings: [ 5,  7, 11, 14, 19, 60, 70, 82, 88]
+  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: null
+  run_chips: [24]
 
   # Whether to enable astrometric modeling
   # astrometric_lib: "libshao.so"
   enable_astrometric_model: True
 
+  # Cut by saturation/limiting magnitude in which band?
+  cut_in_band: "g"
+
 ###############################################
 # Input path setting
 ###############################################
@@ -165,7 +168,7 @@ ins_effects:
   non_linear:     ON # Whether to add non-linearity
   cosmic_ray:     ON # Whether to add cosmic-ray
   cray_differ:    ON # Whether to generate different cosmic ray maps CAL and MS output
-  cte_trail:      OFF # Whether to simulate CTE trails
+  cte_trail:      ON # Whether to simulate CTE trails
   saturbloom:     ON # Whether to simulate Saturation & Blooming
   add_badcolumns: ON # Whether to add bad columns
   add_hotpixels:  ON # Whether to add hot pixels

config/config_NGP_dev.yaml

+---
+###############################################
+#
+#  Configuration file for CSST simulation
+#  CSST-Sim Group, 2021/10/07
+#
+###############################################
+
+# 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: "/public/home/fangyuedong/sim_code_release/CSST/test/"
+work_dir: "/public/home/fangyuedong/test/CSST/workplace/"
+data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
+run_name: "TEST_16channel_off"
+
+# (Optional) a file of point list 
+# if you just want to run default pointing:
+# - pointing_dir: null
+# - pointing_file: null
+pointing_dir: "/data/simudata/CSSOSDataProductsSims/data/"
+pointing_file: "pointing_test_NGP_2.17.dat"
+
+# Whether to use MPI
+run_option:
+  use_mpi: YES
+  # NOTE: "n_threads" paramters is currently not used in the backend
+  # simulation codes. It should be implemented later in the web frontend
+  # in order to config the number of threads to request from NAOC cluster
+  n_threads: 80
+
+  # Output catalog only?
+  # If yes, no imaging simulation will run
+  out_cat_only: NO
+
+  # Only simulate stars?
+  star_only: NO
+
+  # Only simulate galaxies?
+  galaxy_only: NO
+
+###############################################
+# Observation setting
+###############################################
+obs_setting:
+
+  # Options for survey types:
+  # "Photometric": simulate photometric chips only
+  # "Spectroscopic": simulate slitless spectroscopic chips only
+  # "All": simulate full focal plane
+  survey_type: "All"
+  
+  # Exposure time [seconds]
+  exp_time: 150.
+
+  # Observation starting date & time
+  # (Subject to change)
+  date_obs: "210525" # [yymmdd]
+  time_obs: "120000" # [hhmmss]
+
+  # Default Pointing [degrees]
+  # Note: NOT valid when a pointing list file is specified
+  ra_center: 192.8595
+  dec_center:  27.1283
+
+  # Image rotation [degree]
+  image_rot: -113.4333
+
+  # Number of calibration pointings
+  np_cal: 0
+
+  # Run specific pointing(s):
+  # - give a list of indexes of pointings: [ip_1, ip_2...]
+  # - 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: [24]
+
+  # Whether to enable astrometric modeling
+  # astrometric_lib: "libshao.so"
+  enable_astrometric_model: True
+
+  # Cut by saturation/limiting magnitude in which band?
+  cut_in_band: "g"
+
+###############################################
+# Input path setting
+###############################################
+
+# Default path settings for WIDE survey simulation
+input_path:
+  cat_dir: "OnOrbitCalibration/CTargets20211231"
+  star_cat: "CT-NGP_r1.8_G28.hdf5"
+  galaxy_cat: "galaxyCats_r_3.0_healpix_shift_192.859500_27.128300.hdf5"
+
+SED_templates_path:
+  star_SED: "Catalog_20210126/SpecLib.hdf5"
+  galaxy_SED: "Templates/Galaxy/"
+
+###############################################
+# PSF setting
+###############################################
+psf_setting:
+  
+  # Which PSF model to use:
+  # "Gauss": simple gaussian profile
+  # "Interp": Interpolated PSF from sampled ray-tracing data
+  psf_model: "Interp"
+
+  # PSF size [arcseconds]
+  # radius of 80% energy encircled
+  # NOTE: only valid for "Gauss" PSF
+  psf_rcont: 0.15
+
+  # path to PSF data
+  # NOTE: only valid for "Interp" PSF
+  psf_dir: "/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube"
+
+  # path to field-distortion model
+  # Note: only valid when ins_effects: field_dist is "ON"
+  fd_path: "FieldDistModelGlobal_v1.0.pickle"
+  
+  # sigma_spin: 0.0 # psf spin?
+
+###############################################
+# Shear setting
+###############################################
+
+shear_setting:
+  # Options to generate mock shear field:
+  # "constant": all galaxies are assigned a constant reduced shear
+  # "catalog": from catalog (not available yet)
+  # "extra": from seprate file
+  shear_type: "constant"
+
+  # For constant shear filed
+  reduced_g1: 0.026
+  reduced_g2: 0.015
+
+  # Representation of the shear vector?
+  reShear: "E"
+
+  # rotate galaxy ellipticity
+  rotateEll: 0. # [degree]
+
+  # Extra shear catalog
+  # (currently not used)
+  # shear_cat: "mockShear.cat"
+
+###############################################
+# Instrumental effects setting
+###############################################
+ins_effects:
+  # switches
+  field_dist:     ON  # Whether to add field distortions
+  add_back:       ON  # Whether to add sky background
+  add_dark:       ON  # Whether to add dark noise
+  add_readout:    ON # Whether to add read-out (Gaussian) noise
+  add_bias:       ON # Whether to add bias-level to images
+  bias_16channel: ON  # Whether to add different biases for 16 channels
+  gain_16channel: ON  # Whether to make different gains for 16 channels
+  shutter_effect: ON # Whether to add shutter effect
+  flat_fielding:  ON # Whether to add flat-fielding effect
+  prnu_effect:    ON # Whether to add PRNU effect
+  non_linear:     ON # Whether to add non-linearity
+  cosmic_ray:     ON # Whether to add cosmic-ray
+  cray_differ:    ON # Whether to generate different cosmic ray maps CAL and MS output
+  cte_trail:      ON # Whether to simulate CTE trails
+  saturbloom:     ON # Whether to simulate Saturation & Blooming
+  add_badcolumns: ON # Whether to add bad columns
+  add_hotpixels:  ON # Whether to add hot pixels
+  add_deadpixels: ON # Whether to add dead(dark) pixels
+  bright_fatter:  ON # Whether to simulate Brighter-Fatter (also diffusion) effect
+
+  # values
+  dark_exptime:   300   # Exposure time for dark current frames [seconds]
+  flat_exptime:   150   # Exposure time for flat-fielding frames [seconds]
+  readout_time:   40    # The read-out time for each channel [seconds]
+  df_strength:    2.3   # Sillicon sensor diffusion strength
+  bias_level:     500   # bias level [e-/pixel]
+  gain:           1.1   # Gain
+  full_well:      90000 # Full well depth [e-]
+  NBias:          1     # Number of bias frames to be exported for each exposure
+  NDark:          1     # Number of dark frames to be exported for each exposure
+  NFlat:          1     # Number of flat frames to be exported for each exposure
+
+###############################################
+# Output options
+###############################################
+output_setting:
+  readout16:      OFF # Whether to export as 16 channels (subimages) with pre- and over-scan
+  shutter_output: OFF # Whether to export shutter effect 16-bit image
+  bias_output:    ON  # Whether to export bias frames
+  dark_output:    ON  # Whether to export the combined dark current files
+  flat_output:    ON  # Whether to export the combined flat-fielding files
+  prnu_output:    OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files
+
+###############################################
+# Random seeds
+###############################################
+random_seeds:
+  seed_Av:              121212    # Seed for generating random intrinsic extinction
+  seed_poisson:         20210601  # Seed for Poisson noise
+  seed_CR:              20210317  # Seed for generating random cosmic ray maps
+  seed_flat:            20210101  # Seed for generating random flat fields
+  seed_prnu:            20210102  # Seed for photo-response non-uniformity
+  seed_gainNonUniform:  20210202  # Seed for gain nonuniformity
+  seed_biasNonUniform:  20210203  # Seed for bias nonuniformity
+  seed_rnNonUniform:    20210204  # Seed for readout-noise nonuniformity
+  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

config/config_example.yaml

@@ -80,6 +80,9 @@ obs_setting:
   # Note: for all pointings
   run_chips: [18]
 
+  # Cut by saturation/limiting magnitude in which band?
+  cut_in_band: "g"
+
 ###############################################
 # Input path setting
 ###############################################

run_NGP.pbs

@@ -12,8 +12,9 @@
 #PBS -q batch
 #PBS -u fangyuedong
 
-NP=30
+NP=80
 date
 echo $NP
 
-mpirun -np $NP python /public/home/fangyuedong/temp/CSST/run_sim.py config_NGP.yaml -c /public/home/fangyuedong/temp/CSST/config
+# 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 python /public/home/fangyuedong/test/CSST/run_sim.py config_NGP_dev.yaml -c /public/home/fangyuedong/test/CSST/config
\ No newline at end of file

run_sim.py

 from ObservationSim.ObservationSim import Observation
 from ObservationSim._util import parse_args, make_run_dirs, generate_pointing_list
+from pkg_resources import get_distribution
 import os
 import yaml
 import shutil
+import datetime
 
 import gc
 gc.enable()
@@ -20,6 +22,12 @@ def run_sim(Catalog):
     ----------
         None
     """
+    # Get version of CSSTSim Package
+    __version__ = get_distribution("CSSTSim").version
+
+    # Get run datetime
+    now = datetime.datetime.now()
+
     args = parse_args()
     if args.config_dir is None:
         args.config_dir = ''
@@ -40,6 +48,12 @@ def run_sim(Catalog):
     if args.work_dir is not None:
         config['work_dir'] = args.work_dir
 
+    # Some default values
+    if "bias_16channel" not in config:
+        config["bias_16channel"] = False
+    if "gain_16channel" not in config:
+        config["gain_16channel"] = False
+
     # Generate lists pointings based on the input pointing list (or default 
     # pointing RA, DEC) and "config["obs_setting"]["run_pointings"]".
     # "config['obs_setting']['np_cal']"" is the number of CAL pointings which will be 
@@ -49,7 +63,19 @@ def run_sim(Catalog):
 
     # Make the main output directories
     run_dir = make_run_dirs(work_dir=config['work_dir'], run_name=config['run_name'], pointing_list=pointing_list)
+    
+    # Copy the config file to output directory & Write Run metadata
     shutil.copy(args.config_file, run_dir)
+    run_meta = os.path.join(run_dir, "run_metadata.yaml")
+    with open(run_meta, "w") as config_out:
+        config_out.write("\n")
+        config_out.write("###############################################\n")
+        config_out.write("CSSTSim_version: \"%s\"\n"%__version__)
+        date_str = datetime.datetime.strftime(now, '%m/%d/%Y')
+        time_str = datetime.datetime.strftime(now, '%H:%M:%S')
+        config_out.write("Run_date: \"%s\"\n"%date_str)
+        config_out.write("Run_time: \"%s\"\n"%time_str)
+        config_out.write("###############################################\n")
 
     # Initialize the simulation
     obs = Observation(config=config, Catalog=Catalog, work_dir=config['work_dir'], data_dir=config['data_dir'])

setup.py

@@ -26,7 +26,7 @@ extensions = [
 
 
 setup(name='CSSTSim', 
-    version='0.3', 
+    version='0.5', 
     packages=find_packages(),
     # install_requires=[
     #     'numpy>=1.18.5',