add notebook for photometry evaluation

.gitignore

@@ -6,4 +6,5 @@
 *.pyc
 *.so
 *.out
-*pnodes
\ No newline at end of file
+*pnodes
+*-checkpoint.ipynb
\ No newline at end of file

injection_pipeline/Catalog/C5_SimCat.py

-import os
-import numpy as np
-import h5py as h5
-import random
-import galsim
-import astropy.constants as cons
-from astropy.table import Table
-from scipy import interpolate
-
-from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
-from ObservationSim.MockObject._util import seds, sed_assign, extAv, tag_sed, getObservedSED
-
-try:
-    import importlib.resources as pkg_resources
-except ImportError:
-    # Try backported to PY<37 'importlib_resources'
-    import importlib_resources as pkg_resources
-
-class SimCat(CatalogBase):
-    def __init__(self, config, chip, nobjects=None):
-        super().__init__()
-        self.cat_dir = os.path.join(config["data_dir"], config["input_path"]["cat_dir"])
-        self.config = config
-        self.chip = chip
-        self.seed_Av = config["random_seeds"]["seed_Av"]
-
-        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))
-        if "star_cat" in config["input_path"] and config["input_path"]["star_cat"]:
-            star_file = config["input_path"]["star_cat"]
-            star_SED_file = config["SED_templates_path"]["star_SED"]
-            self.star_path = os.path.join(self.cat_dir, star_file)
-            self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
-            self._load_SED_lib_star()
-        if "galaxy_cat" in config["input_path"] and config["input_path"]["galaxy_cat"]:
-            galaxy_file = config["input_path"]["galaxy_cat"]
-            self.galaxy_path = os.path.join(self.cat_dir, galaxy_file)
-            self.galaxy_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["galaxy_SED"])
-            self._load_SED_lib_gals()
-
-        self._load(nobjects=nobjects)
-    
-    def _load_SED_lib_star(self):
-        self.tempSED_star = h5.File(self.star_SED_path,'r')
-
-    def _load_SED_lib_gals(self):
-        self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir=self.galaxy_SED_path)
-
-    def load_norm_filt(self, obj):
-        if obj.type == "star":
-            return self.normF_star
-        elif obj.type == "galaxy" or obj.type == "quasar":
-            return self.normF_galaxy
-        else:
-            return None
-
-    def _load_gals(self, gals, pix_id=None, nobjects=None):
-        # Load how mnay objects?
-        max_ngals = len(gals)
-        if nobjects is None:
-            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 i in range(ngals):
-            param = self.initialize_param()
-            igals = i % max_ngals
-            param['ra'] = gals['ra_true'][igals]
-            param['dec'] = gals['dec_true'][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'] = 21
-            
-            param['z'] = gals['redshift_true'][igals]
-            param['model_tag'] = 'None'
-            param['gamma1'] = 0
-            param['gamma2'] = 0
-            param['kappa'] = 0
-            param['delta_ra'] = 0
-            param['delta_dec'] = 0
-
-            hlrMajB = gals['size_bulge_true'][igals]
-            hlrMinB = gals['size_minor_bulge_true'][igals]
-
-            hlrMajD = gals['size_disk_true'][igals]
-            hlrMinD = gals['size_minor_disk_true'][igals]
-            aGal = gals['size_true'][igals]
-            bGal = gals['size_minor_true'][igals]
-            param['bfrac'] = gals['bulge_to_total_ratio_i'][igals]
-            param['theta'] = gals['position_angle_true'][igals]
-            param['hlr_bulge'] = np.sqrt(hlrMajB * hlrMinB)
-            param['hlr_disk'] = np.sqrt(hlrMajD * hlrMinD)
-            param['ell_bulge'] = (hlrMajB - hlrMinB)/(hlrMajB + hlrMinB)
-            param['ell_disk'] = (hlrMajD - hlrMinD)/(hlrMajD + hlrMinD)
-            param['ell_tot'] = (aGal - bGal) / (aGal + bGal)
-
-            # Assign each galaxy a template SED
-            param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
-            param['redden'] = self.tempRed_gal[param['sed_type']]
-            param['av'] = self.avGal[int(self.ud()*self.nav)]
-            if param['sed_type'] <= 5:
-                param['av'] = 0.0
-                param['redden'] = 0
-            param['star'] = 0   # Galaxy
-            if param['sed_type'] >= 29:
-                param['av'] = 0.6 * param['av'] / 3.0 # for quasar, av=[0, 0.2], 3.0=av.max-av.im
-                param['star'] = 2 # Quasar
-            
-            param['id'] = gals['galaxyID'][igals]
-            
-            if param['star'] == 0:
-                obj = Galaxy(param)
-            if param['star'] == 2:
-                obj = Quasar(param)
-            
-            self.objs.append(obj)
-    
-    def _load(self, nobjects=None):
-        # (TEST) use objects in healpix:
-        pix = 48656
-        self.nav = 15005
-        self.avGal = extAv(self.nav, seed=self.seed_Av)
-        self.objs = []
-
-        gals_cat = h5.File(self.galaxy_path, 'r')['galaxies']
-        gals = gals_cat[str(pix)]
-        self._load_gals(gals, pix_id=pix, nobjects=nobjects)
-        del gals
-
-    def load_sed(self, obj, **kwargs):
-        if obj.type == 'star':
-            _, wave, flux = tag_sed(
-                h5file=self.tempSED_star,
-                model_tag=obj.param['model_tag'],
-                teff=obj.param['teff'],
-                logg=obj.param['logg'],
-                feh=obj.param['feh']
-            )
-        elif obj.type == 'galaxy' or obj.type == 'quasar':
-            sed_data = getObservedSED(
-                sedCat=self.tempSed_gal[obj.sed_type],
-                redshift=obj.z,
-                av=obj.param["av"],
-                redden=obj.param["redden"]
-            )
-            wave, flux = sed_data[0], sed_data[1]
-        else:
-            raise ValueError("Object type not known")
-        speci = interpolate.interp1d(wave, flux)
-        lamb = np.arange(2000, 18001 + 0.5, 0.5)
-        y = speci(lamb)
-        # erg/s/cm2/A --> photo/s/m2/A
-        all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
-        sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
-        del wave
-        del flux
-        return sed
\ No newline at end of file

measurement_pipeline/run_csst_msc_mbi_photometry.py

@@ -9,6 +9,7 @@ from typing import Optional
 
 def run_csst_msc_mbi_photometry(image_path,
                                 output_dir,
+                                flag_weight_dir: str = None,
                                 weight_path: str = None,
                                 flag_path: str = None,
                                 psf_ccds_path: Optional[str] = None,
@@ -18,11 +19,17 @@ def run_csst_msc_mbi_photometry(image_path,
     input_dir = os.path.dirname(image_path)
     img_filename = os.path.basename(image_path)
 
+    if flag_weight_dir is None:
+        flag_weight_dir = input_dir
+
     if weight_path is None:
+        weight_file_name = img_filename.replace("_injected", "")
         weight_path = os.path.join(
-            input_dir, img_filename.replace("img", "wht"))
+            flag_weight_dir, weight_file_name.replace("img", "wht"))
     if flag_path is None:
-        flag_path = os.path.join(input_dir, img_filename.replace("img", "flg"))
+        flag_file_name = img_filename.replace("_injected", "")
+        flag_path = os.path.join(
+            flag_weight_dir, flag_file_name.replace("img", "flg"))
 
     psf_local_path = os.path.join(
         output_dir, img_filename.replace("img", "psf"))
@@ -65,12 +72,12 @@ def genearte_path_list_for_one_pointing(input_dir,
     pointing_dir = os.path.join(input_dir, pointing_label)
     if chip_label_list is None:
         image_path_list = glob(
-            pointing_dir + '/CSST_MSC_MS_SCIE_*_' + '*_img_*')
+            pointing_dir + '/CSST_MSC_MS_SCIE_*_' + '*_img_*.fits')
     else:
         image_path_list = []
         for chip_label in chip_label_list:
             image_path = glob(pointing_dir + '/CSST_MSC_MS_SCIE_*_' +
-                              chip_label + '_img_*')[0]
+                              chip_label + '_img_*.fits')[0]
             image_path_list.append(image_path)
 
     return image_path_list
@@ -79,9 +86,11 @@ def genearte_path_list_for_one_pointing(input_dir,
 def run_pointing_list(input_dir,
                       pointing_label_list,
                       output_dir,
-                      chip_label_list=None):
+                      chip_label_list=None,
+                      flag_weight_dir=None):
     image_path_list = []
     output_path_list = []
+    fw_dir_list = []
     try:
         if not os.path.exists(output_dir):
             os.makedirs(output_dir)
@@ -90,6 +99,7 @@ def run_pointing_list(input_dir,
 
     for pointing_label in pointing_label_list:
         output_pointing_dir = os.path.join(output_dir, pointing_label)
+        fw_dir = os.path.join(flag_weight_dir, pointing_label)
         try:
             if not os.path.exists(output_pointing_dir):
                 os.makedirs(output_pointing_dir)
@@ -101,6 +111,7 @@ def run_pointing_list(input_dir,
         image_path_list = image_path_list + temp_img_path_list
         output_path_list = output_path_list + \
             [output_pointing_dir] * len(temp_img_path_list)
+        fw_dir_list = fw_dir_list + [fw_dir] * len(temp_img_path_list)
 
     comm = MPI.COMM_WORLD
     ind_thread = comm.Get_rank()
@@ -111,21 +122,27 @@ def run_pointing_list(input_dir,
             continue
         image_path = image_path_list[i]
         output_path = output_path_list[i]
+        fw_dir = fw_dir_list[i]
         run_csst_msc_mbi_photometry(image_path=image_path,
-                                    output_dir=output_path)
+                                    output_dir=output_path,
+                                    flag_weight_dir=fw_dir)
 
 
 if __name__ == "__main__":
-    input_dir = "/public/home/fangyuedong/project/50sqDeg_L1_outputs"
+    # input_dir = "/public/home/fangyuedong/project/50sqDeg_L1_outputs"
+    input_dir = "/public/home/fangyuedong/project/injected_50sqDeg_L1_outputs"
     pointing_label_list = ["MSC_0000000", "MSC_0000001",
                            "MSC_0000002", "MSC_0000003", "MSC_0000004", "MSC_0000005"]
     chip_label_list = None
-    output_dir = "/public/home/fangyuedong/project/test_photometry"
+    # output_dir = "/public/home/fangyuedong/project/test_photometry"
+    output_dir = "/public/home/fangyuedong/project/processed_injected_50sqDeg_L1_outputs"
+    flag_weight_dir = "/public/home/fangyuedong/project/50sqDeg_L1_outputs"
 
     # pointing_label_list = ["MSC_0000000"]
     # chip_label_list = ["08"]
     # output_dir = "/public/home/fangyuedong/project/test_photometry"
     run_pointing_list(input_dir=input_dir,
+                      flag_weight_dir=flag_weight_dir,
                       pointing_label_list=pointing_label_list,
                       output_dir=output_dir,
                       chip_label_list=chip_label_list)