rearrange the pipeline structure

README.md

-Pipeline for source injection
\ No newline at end of file
+## Pipeline for source injection
\ No newline at end of file

config_injection.yaml → config/config_injection.yaml

@@ -5,15 +5,21 @@
 #      Last modified: 2022/06/19
 #
 ###############################################
-n_objects: 20
-rotate_objs: False
+# n_objects: 500
+rotate_objs: NO
+use_mpi: YES
+run_name: "test_20230509"
+
 pos_sampling:
   # type: "HexGrid"
-  type: "RectGrid"
-  grid_spacing: 74 # arcsec (~1000 pixels)
+  # type: "RectGrid"
+  # grid_spacing: 18.5 # arcsec (~250 pixels)
+  type: "uniform"
+  object_density: 50 # arcmin^-2
+
+output_img_dir: "./workspace"
+input_img_list: "/share/home/fangyuedong/injection_pipeline/input_L1_img_MSC_0000000.list"
 
-# output_img_name: "test_HexGrid_20220628.fits"
-output_img_name: "test_RectGrid_20220628.fits"
 ###############################################
 # PSF setting
 ###############################################
@@ -31,14 +37,14 @@ psf_setting:
 
   # path to PSF data
   # NOTE: only valid for "Interp" PSF
-  psf_dir: "/data/simudata/CSSOSDataProductsSims/data/csstPSFdata/psfCube"
+  psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"
 
 ###############################################
 # Input path setting
 # (NOTE) Used NGP Catalog for testing
 ###############################################
 # Default path settings for NGP footprint simulation
-data_dir: "/data/simudata/CSSOSDataProductsSims/data/"
+data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
 
 input_path:
   cat_dir: "OnOrbitCalibration/CTargets20211231"
@@ -61,13 +67,13 @@ ins_effects:
   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-]
+  # 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-]
 
 ###############################################
 # Random seeds

config/default.config

+# Default configuration file for SExtractor 2.25.0
+# EB 2019-10-03
+#
+ 
+#-------------------------------- Catalog ------------------------------------
+ 
+CATALOG_NAME     test.cat       # name of the output catalog
+CATALOG_TYPE     FITS_1.0     # NONE,ASCII,ASCII_HEAD, ASCII_SKYCAT,
+                                # ASCII_VOTABLE, FITS_1.0 or FITS_LDAC
+PARAMETERS_NAME  default.param  # name of the file containing catalog contents
+ 
+#------------------------------- Extraction ----------------------------------
+ 
+DETECT_TYPE      CCD            # CCD (linear) or PHOTO (with gamma correction)
+DETECT_MINAREA   3              # min. # of pixels above threshold
+DETECT_THRESH    1.5            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+ANALYSIS_THRESH  1.5            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+THRESH_TYPE      RELATIVE
+ 
+FILTER           Y              # apply filter for detection (Y or N)?
+FILTER_NAME      ../measurement_pipeline/data/default.conv   # name of the file containing the filter
+ 
+DEBLEND_NTHRESH  64             # Number of deblending sub-thresholds
+DEBLEND_MINCONT  0.0005          # Minimum contrast parameter for deblending
+ 
+CLEAN            Y              # Clean spurious detections? (Y or N)?
+CLEAN_PARAM      1.0            # Cleaning efficiency
+ 
+MASK_TYPE        CORRECT        # type of detection MASKing: can be one of
+                                # NONE, BLANK or CORRECT
+
+#------------------------------ Flagging -----------------------------------
+#FLAG_TYPE        OR
+
+#------------------------------ Photometry -----------------------------------
+ 
+PHOT_APERTURES   5              # MAG_APER aperture diameter(s) in pixels
+PHOT_AUTOPARAMS  2.5, 3.5       # MAG_AUTO parameters: <Kron_fact>,<min_radius>
+PHOT_PETROPARAMS 2.0, 3.5       # MAG_PETRO parameters: <Petrosian_fact>,
+                                # <min_radius>
+ 
+SATUR_LEVEL      50000.0        # level (in ADUs) at which arises saturation
+SATUR_KEY        SATURATE       # keyword for saturation level (in ADUs)
+ 
+MAG_ZEROPOINT    0.0            # magnitude zero-point
+#MAG_GAMMA        4.0            # gamma of emulsion (for photographic scans)
+GAIN             1.0            # detector gain in e-/ADU
+#GAIN_KEY         GAIN           # keyword for detector gain in e-/ADU
+PIXEL_SCALE      0            # size of pixel in arcsec (0=use FITS WCS info)
+ 
+#------------------------- Star/Galaxy Separation ----------------------------
+ 
+SEEING_FWHM      0.15           # stellar FWHM in arcsec
+STARNNW_NAME     ../measurement_pipeline/data/default.nnw    # Neural-Network_Weight table filename
+ 
+#------------------------------ Background -----------------------------------
+
+BACK_TYPE        AUTO
+BACK_SIZE        400            # Background mesh: <size> or <width>,<height>
+BACK_FILTERSIZE  3              # Background filter: <size> or <width>,<height>
+ 
+BACKPHOTO_TYPE   GLOBAL         # can be GLOBAL or LOCAL
+ 
+#------------------------------ Check Image ----------------------------------
+ 
+CHECKIMAGE_TYPE  NONE           # can be NONE, BACKGROUND, BACKGROUND_RMS,
+                                # MINIBACKGROUND, MINIBACK_RMS, -BACKGROUND,
+                                # FILTERED, OBJECTS, -OBJECTS, SEGMENTATION,
+                                # or APERTURES
+CHECKIMAGE_NAME  ../check.fits  # Filename for the check-image
+ 
+#--------------------- Memory (change with caution!) -------------------------
+ 
+MEMORY_OBJSTACK  3000           # number of objects in stack
+MEMORY_PIXSTACK  300000         # number of pixels in stack
+MEMORY_BUFSIZE   1024           # number of lines in buffer
+ 
+#----------------------------- Miscellaneous ---------------------------------
+ 
+VERBOSE_TYPE     NORMAL         # can be QUIET, NORMAL or FULL
+HEADER_SUFFIX    .head          # Filename extension for additional headers
+WRITE_XML        N              # Write XML file (Y/N)?
+XML_NAME         sex.xml        # Filename for XML output
+
+#------------------------------ Weighting -----------------------------------
+#WEIGHT_TYPE     MAP_WEIGHT
+#RESCALE_WEIGHTS Y
+#WEIGHT_IMAGE    weight.fits
+#WEIGHT_GAIN     N
+#WEIGHT_THRESH   1e-7

config/default.param

+# ID and flags
+NUMBER 
+FLAGS 
+# Basic coordinates
+X_IMAGE
+Y_IMAGE
+XMIN_IMAGE 
+XMAX_IMAGE 
+YMIN_IMAGE 
+YMAX_IMAGE 
+# Accurate coordinates
+XWIN_IMAGE 
+YWIN_IMAGE 
+ERRAWIN_IMAGE
+ERRBWIN_IMAGE
+ERRTHETAWIN_IMAGE
+ALPHAWIN_J2000 
+DELTAWIN_J2000 
+ERRAWIN_WORLD
+ERRBWIN_WORLD
+ERRTHETAWIN_J2000
+# Peak value
+XPEAK_IMAGE 
+YPEAK_IMAGE 
+ALPHAPEAK_J2000 
+DELTAPEAK_J2000 
+# Basic shape info
+A_IMAGE
+B_IMAGE
+THETA_IMAGE
+A_WORLD
+B_WORLD
+THETA_J2000
+ELONGATION
+FLUX_RADIUS
+# Isophotal photometry
+#FLUX_ISO
+#FLUXERR_ISO
+MAG_ISO
+MAGERR_ISO
+ISOAREA_IMAGE
+ISOAREAF_IMAGE
+ISOAREA_WORLD
+ISOAREAF_WORLD
+# Aperture photometry
+FLUX_APER(1)
+FLUXERR_APER(1)
+MAG_APER(1)
+MAGERR_APER(1)
+# Kron photometry
+FLUX_AUTO
+FLUXERR_AUTO
+MAG_AUTO
+MAGERR_AUTO
+KRON_RADIUS
+# Petrosian photometry
+#FLUX_PETRO
+#FLUXERR_PETRO
+MAG_PETRO
+MAGERR_PETRO
+PETRO_RADIUS
+# Basic surface brightnesses
+BACKGROUND
+THRESHOLD
+MU_THRESHOLD
+FLUX_MAX
+MU_MAX
\ No newline at end of file

config/sex.config

+# Default configuration file for SExtractor 2.25.0
+# EB 2019-10-03
+#
+ 
+#-------------------------------- Catalog ------------------------------------
+ 
+CATALOG_NAME     test.cat       # name of the output catalog
+CATALOG_TYPE     ASCII_HEAD     # NONE,ASCII,ASCII_HEAD, ASCII_SKYCAT,
+                                # ASCII_VOTABLE, FITS_1.0 or FITS_LDAC
+PARAMETERS_NAME  ./sex.param_sersic  # name of the file containing catalog contents
+ 
+#------------------------------- Extraction ----------------------------------
+ 
+DETECT_TYPE      CCD            # CCD (linear) or PHOTO (with gamma correction)
+DETECT_MINAREA   5              # min. # of pixels above threshold
+DETECT_THRESH    1.5            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+ANALYSIS_THRESH  1.5            # <sigmas> or <threshold>,<ZP> in mag.arcsec-2
+ 
+FILTER           Y              # apply filter for detection (Y or N)?
+FILTER_NAME      ./default.conv   # name of the file containing the filter
+ 
+DEBLEND_NTHRESH  32             # Number of deblending sub-thresholds
+DEBLEND_MINCONT  0.005          # Minimum contrast parameter for deblending
+ 
+CLEAN            Y              # Clean spurious detections? (Y or N)?
+CLEAN_PARAM      1.0            # Cleaning efficiency
+ 
+MASK_TYPE        CORRECT        # type of detection MASKing: can be one of
+                                # NONE, BLANK or CORRECT
+ 
+#------------------------------ Photometry -----------------------------------
+PHOT_APERTURES  1.85   3.70   5.55   7.41  11.11   14.81   18.52   22.22   25.93   31.00   44.44  66.67  # MAG_APER aperture diameter(s) in pixels
+#PHOT_APERTURES   5              # MAG_APER aperture diameter(s) in pixels
+PHOT_AUTOPARAMS  2.5, 3.5       # MAG_AUTO parameters: <Kron_fact>,<min_radius>
+PHOT_PETROPARAMS 2.0, 3.5       # MAG_PETRO parameters: <Petrosian_fact>,
+                                # <min_radius>
+ 
+SATUR_LEVEL      50000.0        # level (in ADUs) at which arises saturation
+SATUR_KEY        SATURATE       # keyword for saturation level (in ADUs)
+ 
+MAG_ZEROPOINT    0.0            # magnitude zero-point
+MAG_GAMMA        4.0            # gamma of emulsion (for photographic scans)
+GAIN             1.0            # detector gain in e-/ADU
+GAIN_KEY         GAIN           # keyword for detector gain in e-/ADU
+PIXEL_SCALE      1.0            # size of pixel in arcsec (0=use FITS WCS info)
+
+#-----------------Model-Fitting Photometry-------------------------------------
+# PSF_NMAX     1
+ 
+#------------------------- Star/Galaxy Separation ----------------------------
+ 
+SEEING_FWHM      0.15            # stellar FWHM in arcsec
+STARNNW_NAME     ./default.nnw    # Neural-Network_Weight table filename
+ 
+#------------------------------ Background -----------------------------------
+ 
+BACK_SIZE        256             # Background mesh: <size> or <width>,<height>
+BACK_FILTERSIZE  3              # Background filter: <size> or <width>,<height>
+ 
+BACKPHOTO_TYPE   GLOBAL         # can be GLOBAL or LOCAL
+ 
+#------------------------------ Check Image ----------------------------------
+ 
+CHECKIMAGE_TYPE  BACKGROUND           # can be NONE, BACKGROUND, BACKGROUND_RMS,
+                                # MINIBACKGROUND, MINIBACK_RMS, -BACKGROUND,
+                                # FILTERED, OBJECTS, -OBJECTS, SEGMENTATION,
+                                # or APERTURES
+CHECKIMAGE_NAME  check.fits     # Filename for the check-image
+ 
+#--------------------- Memory (change with caution!) -------------------------
+ 
+MEMORY_OBJSTACK  3000           # number of objects in stack
+MEMORY_PIXSTACK  300000         # number of pixels in stack
+MEMORY_BUFSIZE   1024           # number of lines in buffer
+ 
+#----------------------------- Miscellaneous ---------------------------------
+ 
+VERBOSE_TYPE     NORMAL         # can be QUIET, NORMAL or FULL
+HEADER_SUFFIX    .head          # Filename extension for additional headers
+WRITE_XML        N              # Write XML file (Y/N)?

config/sex.param_sersic

+# ID and flags
+NUMBER 
+FLAGS 
+# Basic coordinates
+X_IMAGE
+Y_IMAGE
+XMIN_IMAGE 
+XMAX_IMAGE 
+YMIN_IMAGE 
+YMAX_IMAGE 
+# Accurate coordinates
+XWIN_IMAGE 
+YWIN_IMAGE 
+ERRAWIN_IMAGE
+ERRBWIN_IMAGE
+ERRTHETAWIN_IMAGE
+ALPHAWIN_J2000 
+DELTAWIN_J2000 
+ERRAWIN_WORLD
+ERRBWIN_WORLD
+ERRTHETAWIN_J2000
+# Peak value
+XPEAK_IMAGE 
+YPEAK_IMAGE 
+ALPHAPEAK_J2000 
+DELTAPEAK_J2000 
+# Basic shape info
+A_IMAGE
+B_IMAGE
+THETA_IMAGE
+A_WORLD
+B_WORLD
+THETA_J2000
+ELONGATION
+FLUX_RADIUS
+# Isophotal photometry
+#FLUX_ISO
+#FLUXERR_ISO
+MAG_ISO
+MAGERR_ISO
+ISOAREA_IMAGE
+ISOAREAF_IMAGE
+ISOAREA_WORLD
+ISOAREAF_WORLD
+# Aperture photometry
+FLUX_APER(12)
+FLUXERR_APER(12)
+MAG_APER(12)
+MAGERR_APER(12)
+# Kron photometry
+FLUX_AUTO
+FLUXERR_AUTO
+MAG_AUTO
+MAGERR_AUTO
+KRON_RADIUS
+# Petrosian photometry
+#FLUX_PETRO
+#FLUXERR_PETRO
+MAG_PETRO
+MAGERR_PETRO
+PETRO_RADIUS
+#hybrid magnitudes
+MAG_HYBRID
+MAGERR_HYBRID
+# Basic surface brightnesses
+BACKGROUND
+THRESHOLD
+MU_THRESHOLD
+FLUX_MAX
+MU_MAX
+# PSF fitting
+#XPSF_IMAGE 
+#YPSF_IMAGE 
+#ERRAPSF_IMAGE
+#ERRBPSF_IMAGE
+#ERRTHETAPSF_IMAGE
+#ALPHAPSF_J2000 
+#DELTAPSF_J2000 
+#ERRAPSF_WORLD
+#ERRBPSF_WORLD
+#ERRTHETAPSF_J2000
+#FLUX_PSF
+#FLUXERR_PSF
+#MAG_PSF
+#MAGERR_PSF
+#NITER_PSF 
+#CHI2_PSF
+# Model fitting coordinates
+XMODEL_IMAGE 
+YMODEL_IMAGE 
+ERRAMODEL_IMAGE
+ERRBMODEL_IMAGE
+ERRTHETAMODEL_IMAGE
+ALPHAMODEL_J2000 
+DELTAMODEL_J2000 
+ERRAMODEL_WORLD
+ERRBMODEL_WORLD
+ERRTHETAMODEL_J2000
+# Global shape measurements from model fitting
+AMODEL_IMAGE
+BMODEL_IMAGE
+THETAMODEL_IMAGE
+AMODEL_WORLD
+BMODEL_WORLD
+THETAMODEL_J2000
+# Shear measurements from model fitting
+ELLIP1MODEL_IMAGE
+ELLIP2MODEL_IMAGE
+ELLIP1MODEL_WORLD
+ELLIP2MODEL_WORLD
+# Global model-fitting photometry
+FLUX_MODEL
+FLUXERR_MODEL
+MAG_MODEL
+MAGERR_MODEL
+#FLUX_MAX_MODEL
+#FLUX_EFF_MODEL
+#FLUX_MEAN_MODEL
+MU_MAX_MODEL
+MU_EFF_MODEL
+MU_MEAN_MODEL
+# Spheroid model component
+#FLUX_SPHEROID
+#FLUXERR_SPHEROID
+MAG_SPHEROID
+MAGERR_SPHEROID
+SPHEROID_REFF_IMAGE
+SPHEROID_REFFERR_IMAGE
+SPHEROID_REFF_WORLD
+SPHEROID_REFFERR_WORLD
+SPHEROID_ASPECT_IMAGE
+SPHEROID_ASPECTERR_IMAGE
+SPHEROID_ASPECT_WORLD
+SPHEROID_ASPECTERR_WORLD
+SPHEROID_THETA_IMAGE
+SPHEROID_THETAERR_IMAGE
+SPHEROID_THETA_J2000
+SPHEROID_THETAERR_WORLD
+SPREAD_MODEL
+SPREADERR_MODEL
+CLASS_STAR
+# Model fitting info
+CHI2_MODEL
+NITER_MODEL 
+#new ones
+FLAGS_WEIGHT
+NLOWWEIGHT_ISO 
+NLOWDWEIGHT_ISO
+#
+SPHEROID_SERSICN
+SPHEROID_SERSICNERR
+#FWHMPSF_IMAGE
+#FWHMPSF_WORLD
+
+
+# added 08/30/2011 per new coadd_objects table def
+ELLIP1ERRMODEL_WORLD
+ELLIP2ERRMODEL_WORLD
+ELLIPCORRMODEL_WORLD
+ELLIPTICITY
+ERRA_IMAGE
+ERRB_IMAGE
+ERRTHETA_IMAGE
+ERRX2_WORLD
+ERRXY_WORLD
+ERRY2_WORLD
+FLAGS_MODEL
+FWHM_WORLD
+X2_WORLD
+XY_WORLD
+Y2_WORLD
+# added on feb 27 2015 after SE cataloging is also done like in dual image mode
+MAG_DETMODEL
+MAGERR_DETMODEL
\ No newline at end of file

evaluation/cross_match_catalogs.py

+import numpy as np
+import astropy.units as u
+import matplotlib.pyplot as plt
+from astropy.coordinates import SkyCoord
+from astropy.io import fits
+from astropy.io import ascii
+from sklearn.neighbors import BallTree
+
+# TU_catalog = "test_RectGrid_20220628.cat"
+# source_catalog = "extracted_test_RectGrid_20220628.fits"
+TU_catalog = "injected_bkgsub_img.cat"
+source_catalog = "extracted_injected_bkgsub_img.fits"
+
+def convert_catalog(catname):
+    text_file = ascii.read(catname)
+    text_file.write('test_ascii_to_fits.fits', overwrite=True)
+
+def read_catalog(catname, ext_num=1, ra_name='ra', dec_name='dec', col_list=[]):
+    hdu = fits.open(catname)
+    hdr = hdu[ext_num].header
+    # print(hdr)
+    data = hdu[ext_num].data
+    ra = data[ra_name]
+    dec = data[dec_name]
+    col_other = []
+    if len(col_list) > 0:
+        for col in col_list:
+            col_other.append(data[col])
+    # print(ra, dec)
+    return ra, dec, col_other
+
+def match_catalogs_sky(ra1, dec1, ra2, dec2, max_dist=0.6, others1=[], others2=[], thresh=[]):
+    cat1 = SkyCoord(ra=ra1*u.degree, dec=dec1*u.degree)
+    cat2 = SkyCoord(ra=ra2*u.degree, dec=dec2*u.degree)
+    idx1, idx2, d2d, d3d = cat2.search_around_sky(cat1, 1*u.deg)
+    # print(idx1)
+    # print(idx2)
+    # print(np.shape(idx1))
+    # print(np.shape(idx2))
+
+def match_catalogs_img(x1, y1, x2, y2, max_dist=0.5, others1=[], others2=[], thresh=[]):
+    cat1 = np.array([(x, y) for x,y in zip(x1, y1)])
+    cat2 = np.array([(x, y) for x,y in zip(x2, y2)])
+    print(np.shape(cat1))
+    print(np.shape(cat2))
+    tree = BallTree(cat2)
+    idx1 = tree.query_radius(cat1, r = max_dist)
+    tree = BallTree(cat1)
+    idx2 = tree.query_radius(cat2, r = max_dist)
+    # print(np.shape(idx1))
+    tot = 0
+    print(idx1)
+    for idx in idx1:
+        if len(idx) == 0:
+            continue
+        tot += 1
+    print(tot)
+    return idx1, idx2
+
+def validation_hist(val1, idx1, name="val1", nbins=10):
+    counts, bins = np.histogram(val1)
+    plt.stairs(bins, counts)
+    plt.xlabel(name, size='x-large')
+    plt.savefig("detection_completeness.png")
+    # plt.show()
+
+if __name__=="__main__":
+    convert_catalog(TU_catalog)
+    ra_TU, dec_TU, _ = read_catalog('test_ascii_to_fits.fits', ext_num=1, ra_name="ra", dec_name="dec")
+    x_TU, y_TU, col_list = read_catalog('test_ascii_to_fits.fits', ext_num=1, ra_name="xImage", dec_name="yImage", col_list=["mag"])
+    mag_TU = col_list[0]
+    # ra_source, dec_source, _ = read_catalog(source_catalog, ext_num=2, ra_name="ALPHAPEAK_J2000", dec_name="DELTAPEAK_J2000")
+    x_source, y_source, _ = read_catalog(source_catalog, ext_num=1, ra_name="X_IMAGE", dec_name="Y_IMAGE")
+    # match_catalogs_sky(ra1=ra_TU, dec1=dec_TU, ra2=ra_source, dec2=dec_source)
+    idx1, idx2, = match_catalogs_img(x1=x_TU, y1=y_TU, x2=x_source, y2=y_source)
+    # print(ra_TU, dec_TU)
+    validation_hist(mag_TU, idx1, name="mag_injected")
\ No newline at end of file

injection_pipeline.py

-from astropy.io import fits
-from astropy import wcs
-import galsim
-import os
-import yaml
-
-from SingleEpochImage import SingleEpochImage
-from InjectionCatalog import InjectionCatalog
-from InputCatalogs import SimCat
-from Grid import RectGrid
-from config import parse_args
-
-# img_file = 'CSST_MSC_MS_SCI_20240617065639_20240617065909_100000000_24_L0_1.fits'
-img_file = 'CSST_MSC_MS_SCI_20240617065639_20240617065909_100000000_23_L0_1.fits' # null-weight test image
-
-args = parse_args()
-if args.config_dir is None:
-    args.config_dir = ''
-args.config_dir = os.path.abspath(args.config_dir)
-args.config_file = os.path.join(args.config_dir, args.config_file)
-print(args.config_file)
-with open(args.config_file, "r") as stream:
-    try:
-        config = yaml.safe_load(stream)
-        for key, value in config.items():
-            print (key + " : " + str(value))
-    except yaml.YAMLError as exc:
-        print(exc)
-
-image = SingleEpochImage(config=config, filepath=img_file)
-# print(image.ramin, image.ramax, image.u_area, image.ra_boundary_cross)
-inject_cat = InjectionCatalog(image=image)
-inject_cat.generate_positions(config=config)
-input_cat = SimCat(config=config, chip=image.chip, nobjects=inject_cat.nobjects)
-print(inject_cat.pos, inject_cat.nobjects)
-image.inject_objects(pos=inject_cat.pos, cat=input_cat)
-image.save_injected_img()
\ No newline at end of file

InputCatalogs.py → injection_pipeline/Catalog/C5_SimCat.py

@@ -58,23 +58,25 @@ class SimCat(CatalogBase):
 
     def _load_gals(self, gals, pix_id=None, nobjects=None):
         # Load how mnay objects?
+        max_ngals = len(gals)
         if nobjects is None:
-            ngals = 5000
+            ngals = max_ngals
         else:
             ngals = nobjects
         self.rng_sedGal = random.Random()
         self.rng_sedGal.seed(pix_id) # Use healpix index as the random seed
         self.ud = galsim.UniformDeviate(pix_id)
 
-        for igals in range(ngals):
+        for i in range(ngals):
             param = self.initialize_param()
+            igals = ngals % max_ngals
             param['ra'] = gals['ra_true'][igals]
             param['dec'] = gals['dec_true'][igals]
 
-            # param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
+            param['mag_use_normal'] = gals['mag_true_g_lsst'][igals]
             # (TEST) use same magnitude 
             # (there will be slight difference due to randomness in SED)
-            param['mag_use_normal'] = 18
+            # param['mag_use_normal'] = 21
             
             param['z'] = gals['redshift_true'][igals]
             param['model_tag'] = 'None'

InjectionCatalog.py → injection_pipeline/InjectionCatalog.py

@@ -7,12 +7,13 @@ import mathutil as util
 class InjectionCatalog(object):
 
     def __init__(self, image):
-        self.rmin, self.rmax = image.ramin, image.ramax
+        self.ramin, self.ramax = image.ramin, image.ramax
         self.decmin, self.decmax = image.decmin, image.decmax
         self.ra_boundary_cross = image.ra_boundary_cross
         self.Npix_x, self.Npix_y = image.Npix_x, image.Npix_y
         self.pixel_scale = image.pixel_scale
         self.wcs = image.wcs
+        print(self.wcs)
 
         # This is also set for a given image, as it may depend on its area
         self.objs_per_real = image.objs_per_real
@@ -119,6 +120,3 @@ class InjectionCatalog(object):
                         angle_unit=self.angle_unit,
                         pos_offset=self.grid_offset)
         return grid_kwargs
-
-    def get_truth_outfile(self):
-        pass

SingleEpochImage.py → injection_pipeline/SingleEpochImage.py

 import numpy as np
 import copy
+import os
+import galsim
+import traceback
 from astropy import wcs
 from astropy.io import fits
-import galsim
 
+from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
 from ObservationSim.Instrument import Chip, Filter, FilterParam, FocalPlane, Telescope
-from ObservationSim.PSF import PSFGauss, PSFInterp
+from ObservationSim.PSF import PSFGauss, PSFInterp, FieldDistortion
+from ObservationSim.Config import ChipOutput
 
 class SingleEpochImage(object):
     
-    def __init__(self, config, filepath):
-        self.header0, self.header_img, self.img = self.read_initial_image(filepath)
+    def __init__(self, config, input_img_path, output_img_path):
+        # Read the origianl image file
+        self.config = config
+        self.read_initial_image(input_img_path)
         self._get_wcs(self.header_img)
         self._determine_unique_area(config)
-        self.output_img_fname = config['output_img_name']
+        self.output_img_fname = output_img_path
+        self.output_dir = os.path.dirname(output_img_path)
 
-        if config['n_objects'] is not None:
+        if 'n_objects' in self.config["pos_sampling"] and self.config["pos_sampling"]['n_objects'] is not None:
             # Fixed number of objects per image
-            self.objs_per_real = config['n_objects']
-        elif config['object_density'] is not None:
+            self.objs_per_real = self.config["pos_sampling"]['n_objects']
+        elif 'object_density' in self.config["pos_sampling"] and self.config["pos_sampling"]['object_density'] is not None:
             # Fixed number density of objects
-            self.objs_per_real = round(self.u_area * config['object_density'])
+            self.objs_per_real = round(self.u_area * self.config["pos_sampling"]['object_density'])
         else:
             # Grid types: calculate nobjects later
             self.objs_per_real = None
 
         self.tel = Telescope()
-        # Determine which CCD
-        self.chip_ID = int(self.header0['DETECTOR'][-2:])
-        # Determine epxosure time
-        self.exp_time = float(self.header0['EXPTIME'])
-        config["obs_setting"]={}
-        config["obs_setting"]["exp_time"] = self.exp_time
-        # Construnct Chip object
-        self.chip = Chip(chipID=self.chip_ID, config=config)
+
+        # Setup field distortion model:
+        if ("field_dist" in self.config) and (self.config["ins_effects"]["field_dist"] == True):
+            self.fd_model = FieldDistortion(chip=chip)
+        else:
+            self.fd_model = None
+
         # Load PSF model
-        if config["psf_setting"]["psf_model"] == "Gauss":
-            self.psf_model = PSFGauss(chip=self.chip)
-        elif config["psf_setting"]["psf_model"] == "Interp":
-            self.psf_model = PSFInterp(chip=self.chip, PSF_data_file=config["psf_setting"]["psf_dir"])
+        if self.config["psf_setting"]["psf_model"] == "Gauss":
+            self.psf_model = PSFGauss(chip=self.chip, psfRa=self.config["psf_setting"]["psf_rcont"])
+        elif self.config["psf_setting"]["psf_model"] == "Interp":
+            self.psf_model = PSFInterp(chip=self.chip, npsf=self.chip.n_psf_samples, PSF_data_file=self.config["psf_setting"]["psf_dir"])
 
+        # Setup Filter
         filter_id, filter_type = self.chip.getChipFilter()
         filter_param = FilterParam()
         self.filt = Filter(filter_id=filter_id,
@@ -49,63 +56,89 @@ class SingleEpochImage(object):
         
         self.setup_image_for_injection()
 
+        self.header_ext = generateExtensionHeader(
+            chip=self.chip,
+            xlen=self.chip.npix_x,
+            ylen=self.chip.npix_y,
+            ra=self.ra_cen,
+            dec=self.dec_cen,
+            pa=self.pos_ang,
+            gain=self.chip.gain, # {TODO}
+            readout=self.chip.read_noise,
+            dark=self.chip.dark_noise,
+            saturation=90000,
+            pixel_scale=self.chip.pix_scale, 
+            pixel_size=self.chip.pix_size,
+            xcen=self.chip.x_cen,
+            ycen=self.chip.y_cen,
+            extName = 'raw')
+
+        # (TODO) specify sub-directory
+        self.chip_output = ChipOutput(
+            config = config,
+            focal_plane = self.focal_plane,
+            chip = self.chip,
+            filt = self.filt,
+            ra_cen = self.ra_cen,
+            dec_cen = self.dec_cen,
+            exptime = self.exp_time,
+            subdir = self.output_dir,
+        )
+        # temp_str = config["output_img_name"].split('.')[0] + ".cat"
+        temp_str = self.output_img_fname.split('/')[-1].split('.')[0] + '.cat'
+        self.chip_output.cat_name = temp_str
+        self.chip_output.create_output_file()
+
     def setup_image_for_injection(self):
-        ra_cen = self.wcs.wcs.crval[0]
-        dec_cen = self.wcs.wcs.crval[1]
-        self.wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, self.pos_ang*galsim.degrees, self.pixel_scale)
-        # self.inj_img = galsim.ImageF(self.chip.npix_x, self.chip.npix_y)
-        self.chip.img = galsim.Image(self.img, copy=True)
+        self.ra_cen = self.wcs.wcs.crval[0]
+        self.dec_cen = self.wcs.wcs.crval[1]
+        self.wcs_fp = self.focal_plane.getTanWCS(self.ra_cen, self.dec_cen, self.pos_ang*galsim.degrees, self.pixel_scale)
+        self.chip.img = galsim.Image(self.image, copy=True)
         self.chip.img.setOrigin(self.chip.bound.xmin, self.chip.bound.ymin)
         self.chip.img.wcs = self.wcs_fp
-        print(self.chip.img.array)
 
     def read_initial_image(self, filepath):
+        data_dir = os.path.dirname(filepath)
+        img_name = os.path.basename(filepath)
         data = fits.open(filepath)
-        header0 = data[0].header
-        header1 = data[1].header
-        image = fits.getdata(filepath)
-
-        # (TEMP)
-        image = np.float64(image)
-        image *= 1.1
-        image -= 500.
-
-        temp_img = galsim.Image(image, copy=True)
-
-        temp_img.array[temp_img.array > 65535] = 65535
-        temp_img.replaceNegative(replace_value=0)
-        temp_img.quantize()
-        temp_img = galsim.Image(temp_img.array, dtype=np.uint16)
-        # self.chip.img = galsim.Image(self.chip.img.array, dtype=np.int32)
-        hdu1 = fits.PrimaryHDU(header=header0)
-        hdu2 = fits.ImageHDU(temp_img.array, header=header1)
-        hdu1 = fits.HDUList([hdu1, hdu2])
-        fname = "nullwt_image_for_injection.fits"
-        hdu1.writeto(fname, output_verify='ignore', overwrite=True)
-        return header0, header1, image
-
-    def _get_wcs(self, header):
-        crpix1 = float(header['CRPIX1'])
-        crpix2 = float(header['CRPIX2'])
+        self.header0 = data[0].header
+        self.header_img = data[1].header
+        self.image = fits.getdata(filepath)
         
-        crval1 = float(header['CRVAL1'])
-        crval2 = float(header['CRVAL2'])
-
-        ctype1 = str(header['CTYPE1'])
-        ctype2 = str(header['CTYPE2'])
+        # Determine which CCD
+        self.chip_ID = int(self.header0['DETECTOR'][-2:])
+        # Construnct Chip object
+        self.chip = Chip(chipID=self.chip_ID, config=self.config)
+        self.exp_time = float(self.header0 ['EXPTIME'])
+
+        # Process L1 image
+        self.image = self.image * self.exp_time / self.chip.gain
+
+        # Process L1 SKY image
+        # [TODO]
+        # sky_img_path = "CSST_MSC_MS_SCI_20240617065639_20240617065909_100000020_08_sky.fits"
+        # sky_img = fits.getdata(sky_img_path)
+        # sky_img = sky_img * self.exp_time / self.chip.gain
+        # image -= sky_img
+
+        # Process L1 FLAG image
+        flag_img_path = os.path.join(data_dir, img_name.replace("img_L1", "flg_L1"))
+        flag_img = fits.getdata(flag_img_path)
+        self.image[flag_img > 0] = 0.
+
+        # Process L1 WEIGHT image
+        # [TODO]
+        
+        # Save null weight image
+        # hdu1 = fits.PrimaryHDU(header=header0)
+        # hdu2 = fits.ImageHDU(image, header=header1)
+        # hdu1 = fits.HDUList([hdu1, hdu2])
+        # fname = "nullwt_image_for_injection.fits"
+        # hdu1.writeto(fname, output_verify='ignore', overwrite=True)
 
-        cd1_1 = float(header['CD1_1'])
-        cd1_2 = float(header['CD1_2'])
-        cd2_1 = float(header['CD2_1'])
-        cd2_2 = float(header['CD2_2'])
+    def _get_wcs(self, header):
         self.pos_ang = float(header['POS_ANG'])
-
-        # Create WCS object
-        self.wcs = wcs.WCS()
-        self.wcs.wcs.crpix = [crpix1, crpix2]
-        self.wcs.wcs.crval = [crval1, crval2]
-        self.wcs.wcs.ctype = [ctype1, ctype2]
-        self.wcs.wcs.cd = [[cd1_1, cd1_2], [cd2_1, cd2_2]]
+        self.wcs = wcs.WCS(header)
 
         self.pixel_scale = 0.074
         self.Npix_x = int(header['NAXIS1'])
@@ -151,14 +184,14 @@ class SingleEpochImage(object):
                         target_filt=self.filt, 
                         norm_filt=norm_filt)
             except Exception as e:
-                print(e)
+                traceback.print_exc()
                 continue
 
             # Update object position to a point on grid
             obj.param['ra'], obj.param['dec'] = pos[i][0], pos[i][1]
-
-            pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=self.chip.img)
-            print(pos_img.x, pos_img.y)
+            self.chip_output.Log_info("ra = %.3f, dec = %.3f"%(obj.param['ra'], obj.param['dec']))
+            pos_img, offset, local_wcs, real_wcs, fd_shear = obj.getPosImg_Offset_WCS(img=self.chip.img, fdmodel=self.fd_model, chip=self.chip, verbose=False, img_header=self.header_ext)
+            self.chip_output.Log_info("pos_img_x = %.3f, pos_img_y = %.3f"%(pos_img.x, pos_img.y))
             try:
                 isUpdated, pos_shear = obj.drawObj_multiband(
                     tel=self.tel,
@@ -172,30 +205,32 @@ class SingleEpochImage(object):
                     exptime=self.exp_time)
                 if isUpdated:
                     # TODO: add up stats
-                    # print("updating output catalog...")
-                    print('Updated')
+                    # self.chip_output.Log_info("updating output catalog...")
+                    self.chip_output.cat_add_obj(obj, pos_img, pos_shear)
                     pass
                 else:
-                    # print("object omitted", flush=True)
+                    self.chip_output.Log_info("object omitted")
                     continue
             except Exception as e:
-                print(e)
+                traceback.print_exc()
                 pass
             # Unload SED:
             obj.unload_SED()
             del obj
+    
+    def add_to_truth_cat(self):
+        pass
 
     def save_injected_img(self):
-        self.chip.img.array[self.chip.img.array > 65535] = 65535
-        self.chip.img.replaceNegative(replace_value=0)
-        self.chip.img.quantize()
-        self.chip.img = galsim.Image(self.chip.img.array, dtype=np.uint16)
-        # self.chip.img = galsim.Image(self.chip.img.array, dtype=np.int32)
+        self.chip.img = galsim.Image(self.chip.img.array, dtype=np.float32)
+        
+        # [TODO] TEST
+        self.chip.img *= self.gain
+        self.chip.img /= self.exp_time
+
         hdu1 = fits.PrimaryHDU(header=self.header0)
         hdu2 = fits.ImageHDU(self.chip.img.array, header=self.header_img)
         hdu1 = fits.HDUList([hdu1, hdu2])
-        # fname = 'test_inject.fits'
-        # fname = '20220621_test_injection.fits'
         fname = self.output_img_fname
         hdu1.writeto(fname, output_verify='ignore', overwrite=True)
 

injection_pipeline/injection_pipeline.py

+from astropy.io import fits
+from astropy import wcs
+import mpi4py.MPI as MPI
+import galsim
+import os
+import yaml
+
+from SingleEpochImage import SingleEpochImage
+from InjectionCatalog import InjectionCatalog
+from Catalog.C5_SimCat import SimCat
+from Grid import RectGrid
+from config import parse_args
+
+class InjectionPipeline(object):
+
+    def __init__(self):
+        # Load configuration
+        args = parse_args()
+        if args.config_dir is None:
+            args.config_dir = ''
+        args.config_dir = os.path.abspath(args.config_dir)
+        args.config_file = os.path.join(args.config_dir, args.config_file)
+        with open(args.config_file, "r") as stream:
+            try:
+                self.config = yaml.safe_load(stream)
+                for key, value in self.config.items():
+                    print (key + " : " + str(value))
+            except yaml.YAMLError as exc:
+                print(exc)
+        # Read the input image paths
+        with open(self.config["input_img_list"]) as input_list:
+            self.img_list = [line.rstrip() for line in input_list]
+
+        # Prepare the output directory
+        if not os.path.exists(self.config["output_img_dir"]):
+            os.makedirs(self.config["output_img_dir"])
+        
+        self.output_dir = os.path.join(self.config["output_img_dir"], self.config["run_name"])
+        if not os.path.exists(self.output_dir):
+            os.makedirs(self.output_dir)
+
+    def run_injection_img_list(self, Catalog):
+        if self.config["use_mpi"]:
+            comm = MPI.COMM_WORLD
+            ind_thread = comm.Get_rank()
+            num_thread = comm.Get_size()
+        
+        for i in range(len(self.img_list)):
+            if self.config["use_mpi"]:
+                if i % num_thread != ind_thread:
+                    continue
+            img_path = self.img_list[i]
+            # Prepare output image path
+            img_name = img_path.split('/')[-1].split('.')[0]
+            output_name = img_name + '_injected.fits'
+            output_path = os.path.join(self.output_dir, output_name)
+
+            image = SingleEpochImage(config=self.config, input_img_path=img_path, output_img_path=output_path)
+            inject_cat = InjectionCatalog(image=image)
+            inject_cat.generate_positions(config=self.config)
+            print("number of galaxies to be injected: %d"%(inject_cat.nobjects))
+            input_cat = Catalog(config=self.config, chip=image.chip, nobjects=inject_cat.nobjects)
+            image.inject_objects(pos=inject_cat.pos, cat=input_cat)
+            image.save_injected_img()
+    
+if __name__ == "__main__":
+    pipeline = InjectionPipeline()
+    pipeline.run_injection_img_list(Catalog=SimCat)
\ No newline at end of file

mathutil.py → injection_pipeline/mathutil.py

@@ -48,9 +48,9 @@ def sample_uniform_dec(d1, d2, N=None, unit='deg'):
     delta = np.arcsin(P * (np.sin(d2) - np.sin(d1)) +np.sin(d1))
 
     # Return correct unit
-    if unit is 'deg' or unit is 'degree':
+    if unit == 'deg' or unit == 'degree':
         return np.rad2deg(delta)
-    elif unit is 'rad' or unit is 'radian':
+    elif unit == 'rad' or unit == 'radian':
         return delta
     else:
         raise TypeError('Only deg/degree or rad/radian are allowed as inputs!')

measurement_pipeline/data/default.conv

+CONV NORM
+# 3x3 ``all-ground'' convolution mask with FWHM = 2 pixels.
+1 2 1
+2 4 2
+1 2 1

measurement_pipeline/data/default.nnw

+NNW
+# Neural Network Weights for the SExtractor star/galaxy classifier (V1.3)
+# inputs:	9 for profile parameters + 1 for seeing.
+# outputs:	``Stellarity index'' (0.0 to 1.0)
+# Seeing FWHM range: from 0.025 to 5.5'' (images must have 1.5 < FWHM < 5 pixels)
+# Optimized for Moffat profiles with 2<= beta <= 4.
+
+ 3 10 10  1
+
+-1.56604e+00 -2.48265e+00 -1.44564e+00 -1.24675e+00 -9.44913e-01 -5.22453e-01  4.61342e-02  8.31957e-01  2.15505e+00  2.64769e-01
+ 3.03477e+00  2.69561e+00  3.16188e+00  3.34497e+00  3.51885e+00  3.65570e+00  3.74856e+00  3.84541e+00  4.22811e+00  3.27734e+00
+
+-3.22480e-01 -2.12804e+00  6.50750e-01 -1.11242e+00 -1.40683e+00 -1.55944e+00 -1.84558e+00 -1.18946e-01  5.52395e-01 -4.36564e-01 -5.30052e+00
+ 4.62594e-01 -3.29127e+00  1.10950e+00 -6.01857e-01  1.29492e-01  1.42290e+00  2.90741e+00  2.44058e+00 -9.19118e-01  8.42851e-01 -4.69824e+00
+-2.57424e+00  8.96469e-01  8.34775e-01  2.18845e+00  2.46526e+00  8.60878e-02 -6.88080e-01 -1.33623e-02  9.30403e-02  1.64942e+00 -1.01231e+00
+ 4.81041e+00  1.53747e+00 -1.12216e+00 -3.16008e+00 -1.67404e+00 -1.75767e+00 -1.29310e+00  5.59549e-01  8.08468e-01 -1.01592e-02 -7.54052e+00
+ 1.01933e+01 -2.09484e+01 -1.07426e+00  9.87912e-01  6.05210e-01 -6.04535e-02 -5.87826e-01 -7.94117e-01 -4.89190e-01 -8.12710e-02 -2.07067e+01
+-5.31793e+00  7.94240e+00 -4.64165e+00 -4.37436e+00 -1.55417e+00  7.54368e-01  1.09608e+00  1.45967e+00  1.62946e+00 -1.01301e+00  1.13514e-01
+ 2.20336e-01  1.70056e+00 -5.20105e-01 -4.28330e-01  1.57258e-03 -3.36502e-01 -8.18568e-02 -7.16163e+00  8.23195e+00 -1.71561e-02 -1.13749e+01
+ 3.75075e+00  7.25399e+00 -1.75325e+00 -2.68814e+00 -3.71128e+00 -4.62933e+00 -2.13747e+00 -1.89186e-01  1.29122e+00 -7.49380e-01  6.71712e-01
+-8.41923e-01  4.64997e+00  5.65808e-01 -3.08277e-01 -1.01687e+00  1.73127e-01 -8.92130e-01  1.89044e+00 -2.75543e-01 -7.72828e-01  5.36745e-01
+-3.65598e+00  7.56997e+00 -3.76373e+00 -1.74542e+00 -1.37540e-01 -5.55400e-01 -1.59195e-01  1.27910e-01  1.91906e+00  1.42119e+00 -4.35502e+00
+
+-1.70059e+00 -3.65695e+00  1.22367e+00 -5.74367e-01 -3.29571e+00  2.46316e+00  5.22353e+00  2.42038e+00  1.22919e+00 -9.22250e-01 -2.32028e+00
+
+
+ 0.00000e+00 
+ 1.00000e+00