Commit b4eb2513 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

add photometry pipeline, warpper for L1_instrument and injection pipelines

parent 5617dc34
......@@ -5,3 +5,5 @@
*.png
*.pyc
*.so
*.out
*pnodes
\ No newline at end of file
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_06_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_07_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_09_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_11_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_12_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_13_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_14_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_15_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_16_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_17_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_18_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_19_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_20_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_22_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_23_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_24_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_25_flg_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_06_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_07_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_09_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_11_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_12_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_13_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_14_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_15_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_16_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_17_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_18_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_19_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_20_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_22_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_23_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_24_wht_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_25_wht_L1.fits
......@@ -7,18 +7,19 @@
###############################################
# n_objects: 500
rotate_objs: NO
use_mpi: YES
run_name: "test_20230517"
run_name: "test_20240421"
project_cycle: 9
run_counter: 1
pos_sampling:
# type: "HexGrid"
type: "HexGrid"
# type: "RectGrid"
# grid_spacing: 18.5 # arcsec (~250 pixels)
type: "uniform"
object_density: 37 # arcmin^-2
grid_spacing: 15 # arcsec (~500 pixels)
# type: "uniform"
# object_density: 37 # arcmin^-2
output_img_dir: "/share/home/fangyuedong/injection_pipeline/workspace"
input_img_list: "/share/home/fangyuedong/injection_pipeline/input_L1_img_MSC_0000000.list"
# output_img_dir: "/public/home/fangyuedong/project/injection_pipeline/workspace"
###############################################
# PSF setting
......@@ -37,23 +38,23 @@ psf_setting:
# path to PSF data
# NOTE: only valid for "Interp" PSF
psf_dir: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"
psf_pho_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/dataC6/psfCube1"
psf_sls_dir: "/share/simudata/CSSOSDataProductsSims/data/SLS_PSF_PCA_fp/"
###############################################
# Input path setting
# (NOTE) Used NGP Catalog for testing
###############################################
# Default path settings for NGP footprint simulation
data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
catalog_options:
input_path:
cat_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/"
galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/"
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:
galaxy_SED: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/sedlibs/"
SED_templates_path:
star_SED: "Catalog_20210126/SpecLib.hdf5"
galaxy_SED: "Templates/Galaxy/"
# rotate galaxy ellipticity
rotateEll: 0. # [degree]
###############################################
......@@ -64,7 +65,7 @@ SED_templates_path:
###############################################
ins_effects:
# switches
bright_fatter: ON # Whether to simulate Brighter-Fatter (also diffusion) effect
# bright_fatter: ON # Whether to simulate Brighter-Fatter (also diffusion) effect
# values
# dark_exptime: 300 # Exposure time for dark current frames [seconds]
......@@ -78,8 +79,8 @@ ins_effects:
###############################################
# Random seeds
###############################################
random_seeds:
seed_Av: 121212 # Seed for generating random intrinsic extinction
# random_seeds:
# seed_Av: 121212 # Seed for generating random intrinsic extinction
###############################################
# Measurement setting
......
---
###############################################
#
# Configuration file for CSST object injection
# Last modified: 2022/06/19
#
###############################################
# n_objects: 500
rotate_objs: NO
use_mpi: YES
run_name: "test_20231203"
project_cycle: 6
run_counter: 1
pos_sampling:
type: "HexGrid"
# type: "RectGrid"
# grid_spacing: 18.5 # arcsec (~250 pixels)
grid_spacing: 15 # arcsec (~500 pixels)
# type: "uniform"
# object_density: 37 # arcmin^-2
output_img_dir: "/share/home/fangyuedong/injection_pipeline/workspace"
input_img_list: "/share/home/fangyuedong/injection_pipeline/input_L1_IMG_20231203.list"
###############################################
# 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: "/share/simudata/CSSOSDataProductsSims/data/psfCube1"
###############################################
# Input path setting
# (NOTE) Used NGP Catalog for testing
###############################################
# Default path settings for NGP footprint simulation
data_dir: "/share/simudata/CSSOSDataProductsSims/data/"
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/"
###############################################
# Instrumental effects setting
# (NOTE) Here only used to construct
# ObservationSim.Instrument.Chip object
# (TODO) Should readout from header
###############################################
ins_effects:
# switches
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-]
###############################################
# Random seeds
###############################################
random_seeds:
seed_Av: 121212 # Seed for generating random intrinsic extinction
###############################################
# Measurement setting
###############################################
measurement_setting:
input_img_list: "/share/home/fangyuedong/injection_pipeline/L1_INJECTED_20231203.list"
# input_img_list: "/share/home/fangyuedong/injection_pipeline/input_L1_img_MSC_0000000.list"
input_wht_list: "/share/home/fangyuedong/injection_pipeline/L1_WHT_20231203.list"
input_flg_list: "/share/home/fangyuedong/injection_pipeline/L1_FLG_20231203.list"
# input_psf_list: "/share/home/fangyuedong/injection_pipeline/psf_img_MSC_0000000.list"
sex_config: "/share/home/fangyuedong/injection_pipeline/config/default.config"
sex_param: "/share/home/fangyuedong/injection_pipeline/config/default.param"
n_jobs: 18
output_dir: "/share/home/fangyuedong/injection_pipeline/workspace"
\ No newline at end of file
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_06_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_07_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_09_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_11_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_12_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_13_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_14_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_15_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_16_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_17_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_18_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_19_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_20_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_22_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_23_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_24_img_L1_injected.fits
/share/home/fangyuedong/injection_pipeline/workspace/test_20230517/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_25_img_L1_injected.fits
......@@ -6,6 +6,7 @@ import h5py as h5
import healpy as hp
import astropy.constants as cons
import traceback
from itertools import cycle
from astropy.coordinates import spherical_to_cartesian
from astropy.table import Table
from scipy import interpolate
......@@ -29,40 +30,44 @@ except ImportError:
# CONSTANTS
NSIDE = 128
bundle_file_list = ['galaxies_C6_bundle000199.h5', 'galaxies_C6_bundle000200.h5', 'galaxies_C6_bundle000241.h5', 'galaxies_C6_bundle000242.h5', 'galaxies_C6_bundle000287.h5', 'galaxies_C6_bundle000288.h5', 'galaxies_C6_bundle000714.h5', 'galaxies_C6_bundle000715.h5', 'galaxies_C6_bundle000778.h5', 'galaxies_C6_bundle000779.h5', 'galaxies_C6_bundle000842.h5', 'galaxies_C6_bundle000843.h5', 'galaxies_C6_bundle002046.h5', 'galaxies_C6_bundle002110.h5', 'galaxies_C6_bundle002111.h5',
'galaxies_C6_bundle002173.h5', 'galaxies_C6_bundle002174.h5', 'galaxies_C6_bundle002238.h5', 'galaxies_C6_bundle002596.h5', 'galaxies_C6_bundle002597.h5', 'galaxies_C6_bundle002656.h5', 'galaxies_C6_bundle002657.h5', 'galaxies_C6_bundle002711.h5', 'galaxies_C6_bundle002712.h5', 'galaxies_C6_bundle002844.h5', 'galaxies_C6_bundle002845.h5', 'galaxies_C6_bundle002884.h5', 'galaxies_C6_bundle002885.h5', 'galaxies_C6_bundle002921.h5', 'galaxies_C6_bundle002922.h5']
def get_bundleIndex(healpixID_ring, bundleOrder=4, healpixOrder=7):
assert NSIDE == 2**healpixOrder
shift = healpixOrder - bundleOrder
shift = 2*shift
nside_bundle = 2**bundleOrder
nside_healpix= 2**healpixOrder
nside_healpix = 2**healpixOrder
healpixID_nest= hp.ring2nest(nside_healpix, healpixID_ring)
healpixID_nest = hp.ring2nest(nside_healpix, healpixID_ring)
bundleID_nest = (healpixID_nest >> shift)
bundleID_ring = hp.nest2ring(nside_bundle, bundleID_nest)
return bundleID_ring
class SimCat(CatalogBase):
def __init__(self, config, chip, nobjects=None):
def __init__(self, config, chip, filt, nobjects=None, logger=None):
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.seed_Av = config["catalog_options"]["seed_Av"]
self.cat_dir = config["catalog_options"]["input_path"]["cat_dir"]
self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
self.normF_star = Table.read(str(filter_path))
self.config = config
self.chip = chip
self.filt = filt
self.logger = logger
galaxy_dir = config["catalog_options"]["input_path"]["galaxy_cat"]
self.galaxy_path = os.path.join(self.cat_dir, galaxy_dir)
self.galaxy_SED_path = os.path.join(config["data_dir"], config["catalog_options"]["SED_templates_path"]["galaxy_SED"])
self.galaxy_SED_path = config["catalog_options"]["SED_templates_path"]["galaxy_SED"]
self._load_SED_lib_gals()
if "rotateEll" in config["catalog_options"]:
self.rotation = float(int(config["catalog_options"]["rotateEll"]/45.))
self.rotation = float(
int(config["catalog_options"]["rotateEll"]/45.))
else:
self.rotation = 0.
......@@ -70,7 +75,8 @@ class SimCat(CatalogBase):
self._load(nobjects=nobjects)
def _get_healpix_list(self):
self.sky_coverage = self.chip.getSkyCoverageEnlarged(self.chip.img.wcs, margin=0.2)
self.sky_coverage = self.chip.getSkyCoverageEnlarged(
self.chip.img.wcs, margin=0.2)
ra_min, ra_max, dec_min, dec_max = self.sky_coverage.xmin, self.sky_coverage.xmax, self.sky_coverage.ymin, self.sky_coverage.ymax
ra = np.deg2rad(np.array([ra_min, ra_max, ra_max, ra_min]))
dec = np.deg2rad(np.array([dec_max, dec_max, dec_min, dec_min]))
......@@ -86,11 +92,12 @@ class SimCat(CatalogBase):
print("HEALPix List: ", self.pix_list)
def load_norm_filt(self, obj):
if obj.type == "star":
return self.normF_star
elif obj.type == "galaxy" or obj.type == "quasar":
return None
else:
# if obj.type == "star":
# return self.normF_star
# elif obj.type == "galaxy" or obj.type == "quasar":
# return None
# else:
# return None
return None
def _load_SED_lib_gals(self):
......@@ -107,12 +114,16 @@ class SimCat(CatalogBase):
if remain == 0:
break
param = self.initialize_param()
param['ra'] = ra_arr[igals]
param['dec'] = dec_arr[igals]
param['ra'] = gals['ra'][igals]
param['dec'] = gals['dec'][igals]
param['ra_orig'] = gals['ra'][igals]
param['dec_orig'] = gals['dec'][igals]
# [TODO]
param['mag_use_normal'] = gals['mag_csst_%s'%(self.filt.filter_type)][igals]
if self.filt.filter_type == 'NUV':
param['mag_use_normal'] = gals['mag_csst_nuv'][igals]
else:
param['mag_use_normal'] = gals['mag_csst_%s' %
(self.filt.filter_type)][igals]
# if self.filt.is_too_dim(mag=param['mag_use_normal'], margin=self.config["obs_setting"]["mag_lim_margin"]):
# continue
......@@ -126,7 +137,12 @@ class SimCat(CatalogBase):
# For shape calculation
param['ell_total'] = np.sqrt(param['e1']**2 + param['e2']**2)
# For shape calculation
param['e1'], param['e2'], param['ell_total'] = self.rotate_ellipticity(
e1=gals['ellipticity_true'][igals][0],
e2=gals['ellipticity_true'][igals][1],
rotation=self.rotation,
unit='radians')
if param['ell_total'] > 0.9:
continue
......@@ -136,7 +152,6 @@ class SimCat(CatalogBase):
param['e1_bulge'] = param['e1']
param['e2_bulge'] = param['e2']
param['delta_ra'] = 0
param['delta_dec'] = 0
......@@ -145,15 +160,14 @@ class SimCat(CatalogBase):
param['diskmass'] = gals['diskmass'][igals]
param['size'] = gals['size'][igals]
if param['size'] > self.max_size:
self.max_size = param['size']
# Sersic index
param['disk_sersic_idx'] = 1.
param['bulge_sersic_idx'] = 4.
# Sizes
param['bfrac'] = param['bulgemass']/(param['bulgemass'] + param['diskmass'])
param['bfrac'] = param['bulgemass'] / \
(param['bulgemass'] + param['diskmass'])
if param['bfrac'] >= 0.6:
param['hlr_bulge'] = param['size']
param['hlr_disk'] = param['size'] * (1. - param['bfrac'])
......@@ -178,16 +192,16 @@ class SimCat(CatalogBase):
# TEMP
self.ids += 1
# param['id'] = self.ids
param['id'] = '%06d'%(int(pix_id)) + '%06d'%(cat_id) + '%08d'%(igals)
param['id'] = '%06d' % (int(pix_id)) + \
'%06d' % (cat_id) + '%08d' % (igals)
if param['star'] == 0:
obj = Galaxy(param, self.rotation, logger=self.logger)
obj = Galaxy(param, logger=self.logger)
self.objs.append(obj)
return remain
def _load(self, nobjects=1):
from itertools import cycle
self.objs = []
self.ids = 0
to_be_read_in = nobjects
......@@ -197,10 +211,12 @@ class SimCat(CatalogBase):
if to_be_read_in == 0:
break
bundleID = get_bundleIndex(pix)
file_path = os.path.join(self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
file_path = os.path.join(
self.galaxy_path, "galaxies_C6_bundle{:06}.h5".format(bundleID))
gals_cat = h5.File(file_path, 'r')['galaxies']
gals = gals_cat[str(pix)]
to_be_read_in = self._load_gals(gals, pix_id=pix, cat_id=bundleID, n_objects=to_be_read_in)
to_be_read_in = self._load_gals(
gals, pix_id=pix, cat_id=bundleID, nobjects=to_be_read_in)
del gals
except Exception as e:
traceback.print_exc()
......
......@@ -2,10 +2,12 @@ import math
import numpy as np
import matplotlib.pyplot as plt
class BaseGrid(object):
_valid_grid_types = ['RectGrid', 'HexGrid']
_valid_mixed_types = ['MixedGrid']
class Grid(BaseGrid):
def __init__(self, grid_spacing, wcs, Npix_x=10000, Npix_y=10000, pixelscale=0.074, rot_angle=None, pos_offset=None, angle_unit='rad'):
self.grid_spacing = grid_spacing
......@@ -28,13 +30,17 @@ class Grid(BaseGrid):
theta = np.deg2rad(rot_angle)
else:
theta = rot_angle
self.startx = (0.-dx) * np.cos(theta) - (Npix_y-dy) * np.sin(theta) + dx
self.endx = (Npix_x-dx) * np.cos(theta) - (0.-dy) * np.sin(theta) + dx
self.starty = (0.-dx) * np.cos(theta) + (0.-dy) * np.sin(theta) + dx
self.endy = (Npix_x-dx) * np.cos(theta) + (Npix_y-dy) * np.sin(theta) + dx
self.startx = (0.-dx) * np.cos(theta) - \
(Npix_y-dy) * np.sin(theta) + dx
self.endx = (Npix_x-dx) * np.cos(theta) - \
(0.-dy) * np.sin(theta) + dx
self.starty = (0.-dx) * np.cos(theta) + \
(0.-dy) * np.sin(theta) + dx
self.endy = (Npix_x-dx) * np.cos(theta) + \
(Npix_y-dy) * np.sin(theta) + dx
else:
self.startx, self.endx= 0., Npix_x
self.starty, self.endy= 0., Npix_y
self.startx, self.endx = 0., Npix_x
self.starty, self.endy = 0., Npix_y
def rotate_grid(self, theta, offset=None, angle_unit='rad'):
......@@ -44,20 +50,22 @@ class Grid(BaseGrid):
raise ValueError('`angle_unit` can only be `deg` or `rad`! ' +
'Passed unit of {}'.format(angle_unit))
if not offset: offset = [0., 0.]
if not offset:
offset = [0., 0.]
c, s = np.cos(theta), np.sin(theta)
R = np.array(((c,-s), (s, c)))
R = np.array(((c, -s), (s, c)))
offset_grid = np.array([self.im_ra - offset[0], self.im_dec - offset[1]])
offset_grid = np.array(
[self.im_ra - offset[0], self.im_dec - offset[1]])
translate = np.empty_like(offset_grid)
translate[0,:] = offset[0]
translate[1,:] = offset[1]
translate[0, :] = offset[0]
translate[1, :] = offset[1]
rotated_grid = np.dot(R, offset_grid) + translate
self.im_pos = rotated_grid.T
self.im_ra, self.im_dec = self.im_pos[0,:], self.im_pos[1,:]
self.im_ra, self.im_dec = self.im_pos[0, :], self.im_pos[1, :]
def cut2buffer(self):
'''
......@@ -66,16 +74,17 @@ class Grid(BaseGrid):
possible rotations.
'''
b = self.im_gs
in_region = np.where( (self.im_pos[:,0]>b) & (self.im_pos[:,0]<self.Npix_x-b) &
(self.im_pos[:,1]>b) & (self.im_pos[:,1]<self.Npix_y-b) )
in_region = np.where((self.im_pos[:, 0] > b) & (self.im_pos[:, 0] < self.Npix_x-b) &
(self.im_pos[:, 1] > b) & (self.im_pos[:, 1] < self.Npix_y-b))
self.im_pos = self.im_pos[in_region]
self.im_ra = self.im_pos[:,0]
self.im_dec = self.im_pos[:,1]
self.im_ra = self.im_pos[:, 0]
self.im_dec = self.im_pos[:, 1]
# Get all image coordinate pairs
self.pos = self.wcs.wcs_pix2world(self.im_pos, 1)
self.ra = self.pos[:,0]
self.dec = self.pos[:,1]
self.ra = self.pos[:, 0]
self.dec = self.pos[:, 1]
class RectGrid(Grid):
def __init__(self, grid_spacing, wcs, Npix_x=10000, Npix_y=10000, pixelscale=0.074,
......@@ -96,8 +105,8 @@ class RectGrid(Grid):
# Get all image coordinate pairs
self.im_pos = np.array(np.meshgrid(self.im_ra, self.im_dec)).T.reshape(
-1, 2)
self.im_ra = self.im_pos[:,0]
self.im_dec = self.im_pos[:,1]
self.im_ra = self.im_pos[:, 0]
self.im_dec = self.im_pos[:, 1]
if self.rot_angle:
self.rotate_grid(self.rot_angle, angle_unit=self.angle_unit,
......@@ -105,6 +114,7 @@ class RectGrid(Grid):
self.cut2buffer()
class HexGrid(Grid):
def __init__(self, grid_spacing, wcs, Npix_x=10000, Npix_y=10000, pixelscale=0.074,
rot_angle=None, pos_offset=None, angle_unit='rad'):
......@@ -118,7 +128,8 @@ class HexGrid(Grid):
po = self.pos_offset
im_po = [p / self.pixelscale for p in po]
self.im_pos = HexGrid.calc_hex_coords(self.startx, self.starty, self.endx, self.endy, im_gs)
self.im_pos = HexGrid.calc_hex_coords(
self.startx, self.starty, self.endx, self.endy, im_gs)
self.im_ra = self.im_pos[:, 0]
self.im_dec = self.im_pos[:, 1]
......@@ -142,18 +153,20 @@ class HexGrid(Grid):
ys = []
while startx < endx:
x = [startx, startx, startx + r, startx + dx, startx + dx, startx + r, startx + r]
x = [startx, startx, startx + r, startx +
dx, startx + dx, startx + r, startx + r]
xs.append(x)
startx += dx
while starty < endy:
y = [starty + p, starty + 3*p, starty + h, starty + 3*p, starty + p, starty, starty + dy]
y = [starty + p, starty + 3*p, starty + h,
starty + 3*p, starty + p, starty, starty + dy]
ys.append(y)
starty += 2*p
row += 1
print(xs)
print(ys)
# print(xs)
# print(ys)
polygons = [zip(x, y) for x in xs for y in ys]
polygons = [np.column_stack((x, y)) for x in xs for y in ys]
......@@ -161,25 +174,31 @@ class HexGrid(Grid):
hexgrid = cls.polygons2coords(polygons)
# Some hexagonal elements go beyond boundary; cut these out
indx = np.where( (hexgrid[:,0]<endx) & (hexgrid[:,1]<endy) )
indx = np.where((hexgrid[:, 0] < endx) & (hexgrid[:, 1] < endy))
return hexgrid[indx]
@classmethod
def polygons2coords(HexGrid, p):
print(p)
# print(p)
s = np.shape(p)
print(s)
# print(s)
L = s[0]*s[1]
pp = np.array(p).reshape(L,2)
c = np.vstack({tuple(row) for row in pp})
pp = np.array(p).reshape(L, 2)
# print(pp)
# print(pp.shape)
# c = np.vstack({tuple(row) for row in pp})
c = np.vstack([tuple(row) for row in pp])
# Some of the redundant coordinates are offset by ~1e-10 pixels
return np.unique(c.round(decimals=6), axis=0)
def _build_grid(grid_type, **kwargs):
if grid_type in GRID_TYPES:
return GRID_TYPES[grid_type](**kwargs)
else:
raise ValueError('There is not yet an implemnted default Grid of type {}'.format(grid_type))
raise ValueError(
'There is not yet an implemnted default Grid of type {}'.format(grid_type))
GRID_TYPES = {
'RectGrid': RectGrid,
......
......@@ -5,30 +5,38 @@ import galsim
import traceback
from astropy import wcs
from astropy.io import fits
from astropy.time import Time
from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
# from ObservationSim.Config.Header import generatePrimaryHeader, generateExtensionHeader
from ObservationSim.Instrument import Chip, Filter, FilterParam, FocalPlane, Telescope
from ObservationSim.PSF import PSFGauss, PSFInterp, FieldDistortion
from ObservationSim.Config import ChipOutput
from ObservationSim.Config import Pointing
class SingleEpochImage(object):
def __init__(self, config, input_img_path, output_img_path):
def __init__(self, config, input_img_path, output_dir):
# Read the origianl image file
self.config = config
self.read_initial_image(input_img_path)
self._get_wcs(self.header_img)
self._get_wcs(hdr_0=self.header0, hdr_img=self.header_img)
self._determine_unique_area(config)
self.output_img_fname = output_img_path
self.output_dir = os.path.dirname(output_img_path)
output_img_name = input_img_path.split(
'/')[-1].split('.')[0] + '_injected.fits'
output_img_path = os.path.join(output_dir, output_img_name)
self.output_img_path = output_img_path
self.output_dir = output_dir
self.pointing.output_dir = self.output_dir
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 = 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 * self.config["pos_sampling"]['object_density'])
print("number of injected obj = ",self.objs_per_real)
self.objs_per_real = round(
self.u_area * self.config["pos_sampling"]['object_density'])
print("number of injected obj = ", self.objs_per_real)
else:
# Grid types: calculate nobjects later
self.objs_per_real = None
......@@ -41,12 +49,6 @@ class SingleEpochImage(object):
else:
self.fd_model = None
# Load PSF model
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()
......@@ -59,20 +61,29 @@ class SingleEpochImage(object):
self.header_ext = self.header_img
# Load PSF model
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":
if self.filt.survey_type == "photometric":
self.psf_model = PSFInterp(chip=self.chip, npsf=self.chip.n_psf_samples,
PSF_data_file=self.config["psf_setting"]["psf_pho_dir"])
else:
self.psf_model = PSFInterp(chip=self.chip, npsf=self.chip.n_psf_samples,
PSF_data_file=self.config["psf_setting"]["psf_sls_dir"])
# (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,
config=config,
chip=self.chip,
filt=self.filt,
pointing=self.pointing,
logger_filename=self.output_img_path.split(
'/')[-1].split('.')[0] + '.log'
)
# 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.cat_name = self.output_img_path.split(
'/')[-1].split('.')[0] + '.cat'
self.chip_output.create_output_file()
def setup_image_for_injection(self):
......@@ -99,7 +110,7 @@ class SingleEpochImage(object):
self.chip_ID = int(self.header_img["CHIPID"])
# Construnct Chip object
self.chip = Chip(chipID=self.chip_ID, config=self.config)
self.exp_time = float(self.header0 ['EXPTIME'])
self.exp_time = float(self.header0['EXPTIME'])
# self.chip.gain = float(self.header_img["GAIN1"])
self.chip.gain = float(self.header_img["GAIN01"])
......@@ -119,7 +130,8 @@ class SingleEpochImage(object):
# Process L1 FLAG image
# flag_img_path = os.path.join(data_dir, img_name.replace("img_L1", "flg_L1"))
flag_img_path = os.path.join(data_dir, img_name.replace("IMG", "FLG"))
# flag_img_path = os.path.join(data_dir, img_name.replace("IMG", "FLG"))
flag_img_path = os.path.join(data_dir, img_name.replace("img", "flg"))
flag_img = fits.getdata(flag_img_path)
self.image[flag_img > 0] = 0.
......@@ -133,16 +145,27 @@ class SingleEpochImage(object):
# fname = "nullwt_image_for_injection.fits"
# hdu1.writeto(fname, output_verify='ignore', overwrite=True)
def _get_wcs(self, header):
def _get_wcs(self, hdr_0, hdr_img):
# self.pos_ang = float(header['POS_ANG'])
self.wcs = wcs.WCS(header)
self.wcs = wcs.WCS(hdr_img)
self.pixel_scale = 0.074
self.Npix_x = int(header['NAXIS1'])
self.Npix_y = int(header['NAXIS2'])
self.Npix_x = int(hdr_img['NAXIS1'])
self.Npix_y = int(hdr_img['NAXIS2'])
self.pointing = Pointing(
ra=hdr_0['RA_PNT0'],
dec=hdr_0['DEC_PNT0'],
img_pa=hdr_0["POS_ANG0"],
timestamp=Time(hdr_0['DATE-OBS']).unix,
exp_time=hdr_0['EXPTIME'],
pointing_id=hdr_0['OBSID'],
pointing_type='SCI',
pointing_type_code='101')
def _determine_unique_area(self, config):
coners = np.array([(1, 1), (1, self.Npix_y), (self.Npix_x, 1), (self.Npix_x, self.Npix_y)])
coners = np.array(
[(1, 1), (1, self.Npix_y), (self.Npix_x, 1), (self.Npix_x, self.Npix_y)])
coners = self.wcs.wcs_pix2world(coners, 1)
ra_coners = coners[:, 0]
dec_coners = coners[:, 1]
......@@ -177,7 +200,7 @@ class SingleEpochImage(object):
try:
sed_data = cat.load_sed(obj)
norm_filt = cat.load_norm_filt(obj)
obj.sed, obj.param["mag_%s"%self.filt.filter_type.lower()], obj.param["flux_%s"%self.filt.filter_type.lower()] = cat.convert_sed(
obj.sed, obj.param["mag_%s" % self.filt.filter_type.lower()], obj.param["flux_%s" % self.filt.filter_type.lower()] = cat.convert_sed(
mag=obj.param["mag_use_normal"],
sed=sed_data,
target_filt=self.filt,
......@@ -188,9 +211,15 @@ class SingleEpochImage(object):
# Update object position to a point on grid
obj.param['ra'], obj.param['dec'] = pos[i][0], pos[i][1]
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, chip_wcs=chip_wcs)
self.chip_output.Log_info("pos_img_x = %.3f, pos_img_y = %.3f"%(pos_img.x, pos_img.y))
obj.ra, obj.dec = pos[i][0], pos[i][1]
obj.ra_orig, obj.dec_orig = pos[i][0], pos[i][1]
# 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, chip_wcs=chip_wcs)
# 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,
......@@ -201,7 +230,8 @@ class SingleEpochImage(object):
chip=self.chip,
g1=obj.g1,
g2=obj.g2,
exptime=self.exp_time)
exptime=self.exp_time,
fd_shear=fd_shear)
if isUpdated:
# TODO: add up stats
# self.chip_output.Log_info("updating output catalog...")
......@@ -231,7 +261,5 @@ class SingleEpochImage(object):
hdu1 = fits.PrimaryHDU(header=self.header0)
hdu2 = fits.ImageHDU(self.chip.img.array, header=self.header_img)
hdu1 = fits.HDUList([hdu1, hdu2])
fname = self.output_img_fname
fname = self.output_img_path
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
from astropy.io import fits
from astropy import wcs
import mpi4py.MPI as MPI
import galsim
import os
import yaml
from glob import glob
from SingleEpochImage import SingleEpochImage
from InjectionCatalog import InjectionCatalog
from Catalog.C5_SimCat import SimCat
from Grid import RectGrid
from Catalog.C6_SimCat import SimCat
from config import parse_args
class InjectionPipeline(object):
def __init__(self):
def __init__(self, config_file=None):
# Load configuration
if config_file is None:
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:
config_file = args.config_file
with open(config_file, "r") as stream:
try:
self.config = yaml.safe_load(stream)
for key, value in self.config.items():
print (key + " : " + str(value))
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"]):
def genearte_path_list_for_one_pointing(self,
input_dir,
pointing_label,
chip_label_list=None):
"""_summary_
Args:
input_dir (_type_): _description_
pointing_label (_type_): _description_
chip_label_list (_type_, optional): _description_. Defaults to None.
Returns:
_type_: _description_
"""
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_*')
else:
image_path_list = []
for chip_label in chip_label_list:
image_path = glob(pointing_dir + '/CSST_MSC_MS_SCIE_*_' +
chip_label + '_img_*')[0]
image_path_list.append(image_path)
return image_path_list
def run_pointing_list(self,
input_dir,
pointing_label_list,
output_dir,
Catalog,
chip_label_list=None):
image_path_list = []
output_path_list = []
try:
os.makedirs(self.config["output_img_dir"])
if not os.path.exists(output_dir):
os.makedirs(output_dir)
except OSError:
pass
self.output_dir = os.path.join(self.config["output_img_dir"], self.config["run_name"])
if not os.path.exists(self.output_dir):
for pointing_label in pointing_label_list:
output_pointing_dir = os.path.join(output_dir, pointing_label)
try:
os.makedirs(self.output_dir)
if not os.path.exists(output_pointing_dir):
os.makedirs(output_pointing_dir)
except OSError:
pass
temp_img_path_list = self.genearte_path_list_for_one_pointing(
input_dir=input_dir,
pointing_label=pointing_label,
chip_label_list=chip_label_list)
image_path_list = image_path_list + temp_img_path_list
output_path_list = output_path_list + \
[output_pointing_dir] * len(temp_img_path_list)
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"]:
for i in range(len(image_path_list)):
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)
image_path = image_path_list[i]
output_path = output_path_list[i]
image = SingleEpochImage(
config=self.config, input_img_path=image_path, output_dir=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)
print("number of galaxies to be injected: %d" %
(inject_cat.nobjects))
input_cat = Catalog(config=self.config,
chip=image.chip,
filt=image.filt,
nobjects=inject_cat.nobjects,
logger=image.chip_output.logger)
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
input_dir = "/public/home/fangyuedong/project/50sqDeg_L1_outputs"
pointing_label_list = ["MSC_0000000", "MSC_0000001",
"MSC_0000002", "MSC_0000003", "MSC_0000004", "MSC_0000005"]
chip_label_list = None
# pointing_label_list = ["MSC_0000000"]
# chip_label_list = ["08"]
output_dir = "/public/home/fangyuedong/project/injected_50sqDeg_L1_outputs"
config_file = "/public/home/fangyuedong/project/injection_pipeline/config/config_injection.yaml"
pipeline = InjectionPipeline(config_file=config_file)
pipeline.run_pointing_list(input_dir=input_dir,
pointing_label_list=pointing_label_list,
output_dir=output_dir,
chip_label_list=chip_label_list,
Catalog=SimCat)
#! /bin/bash
#SBATCH -J INJECT
#SBATCH -N 1
#SBATCH --ntasks-per-node=36
#SBATCH -p batch
#SBATCH --mem=240G
module load mpi/openmpi/4.0.2/gcc-7.3.1
date
srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes
mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 36 --map-by node python injection_pipeline.py
\ No newline at end of file
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_06_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_07_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_08_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_09_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_11_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_12_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_13_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_14_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_15_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_16_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_17_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_18_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_19_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_20_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_22_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_23_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_24_img_L1.fits
/share/L1Result/C6_2sq_sheared_photometry/MSC_0000000/CSST_MSC_MS_SCI_20240817060512_20240817060742_100000000_25_img_L1.fits
"""
Identifier: csst_msc_mbi_photometry/__init__.py
Name: __init__.py
Description: photometry of CSST MS-MBI
Author: Hu Zou (zouhu@nao.cas.cn)
Created: 2022-11-11
Modified-History:
2022-11-11, created
2023-11-23, new API and format
"""
import os
from .csst_photometry import core_msc_l1_mbi_phot
__version__ = "1.4.1"
PACKAGE_PATH = os.path.dirname(__file__)
__all__ = ["core_msc_l1_mbi_phot"]
\ No newline at end of file
"""
Identifier: KSC-SJ1-MSC-MBI/stats.py
Name: stats.py
Description: some stats programs
Author: Hu Zou (zouhu@nao.cas.cn)
Created: 2022-11-11
Modified-History:
2022-11-11, created
2023-12-04, add docstring
v1.1.2 first version of csst photometry pipeline based on sextractor
v1.2.0 change weight and flag names due to the CSST naming
add FLAGS_ISO to remove flagged objects in PSFEx
v1.2.1 change gain to exposure time
add options of output images
change catalog columns to the defined L2 data format
change config directory
v1.2.2 add pick_psfstars, select PSF stars
add photometry types: PSF or MODEL
v1.2.3 add computing time
v1.2.4 bug in pick_psfstars with missing max_elp
bug in psf name if outdir is not provided
add saturate for phot
v1.2.5 bug in pick_psfstars, the match_dist should be set
add PHOT_FLUXFRAC to 0.2,0.5,0.8 to give the radii of containing different fraction of the total flux
v1.2.6 rename the catalog and clean columns; remove sersic bulge; only use exp disk + dev buldge
bug in adding multi-PSF fitting, cleaning catalog
v1.3.0 add ID for single-epoch images defined as obsid+ccdno+objid;
combine two image header into one primary header in the catalog
add more annotations
v1.3.1
fixing version for the CSST C6
v1.3.2 coding style fixed
add return status and file recorder in the do_phot main program (rubish API design)
v1.3.3 fixed some output defaults and change default bad values
v1.3.4 add some keywords:VER_, STM_, STA; change output file postfix to upper case; remove PAErr, change ccdno
v1.3.5 change APCOR? keywords
change unittest
change unit test for C6.2
shrink code
v1.3.6 CHAGE CODE style
v1.3.7 CCDCHIP keywords change to CHIPID and value from ccd+number to number
chang CsstMbiDataManager to CsstMsDataManager
change all APIs from DM to function
v1.4.0 new cycle C7:
get_psf: new limits on FWHM, max_area, PSF size, degree
FWHM from [1.5,20] to [1.2,5]
max_area from 0 (no limit) to 150
psf cutout: from 101 to 31
psf size from 101 to 25
degree from 3 to 2
new strategy for picking psf stars: remove pollution -
select good stars - keep max stars if overdensity
v1.4.1 new API
v1.4.2 docstring issues
v1.4.3 add getpsf_flag in API and units in catalog
add filter in catalog
v1.4.4 new sextractor background
"""
from __future__ import (absolute_import, division,
print_function, unicode_literals)
import numpy as np
import matplotlib.pyplot as plt
import argparse
from astropy import wcs as awcs
from astropy.io import fits
from astropy import table
from astropy.table import Table
from astropy.time import Time
import re
import os
import sys
import time
from scipy.interpolate import UnivariateSpline
from astropy.coordinates import SkyCoord
from shutil import which
import astropy.units as u
from .stats import sigmaclip_limitsig, weighted_mean, closest_match
from .mag_flux_convert import fluxerr2magerr
# from csst_common import CsstResult, CsstStatus
# from typing import Optional
prog_dir = os.path.dirname(__file__)
config_path = os.path.join(prog_dir, 'data/')
__version__ = "1.4.3"
def pick_psfstars(psffile, remove_polution=False, min_nstar=10, max_nstar=1500, class_star=0.7, match_dist=10,
min_sn=20, max_elp=0.3, fwhm_range=[1.5, 20]):
""" pick PSF stars (remove PSF stars with neighbour objects, limit maximum number of stars
Parameters
----------
psffile str
catalog used for selecting PSF stars
remove_polution: bool
wether to remove objects with pollution from nearby objects
min_nstar: int
mininum stars for constructing PSF model
max_nstar: int
maximum stars for constructing PSF model
max_elp: float
maximum ellpticity
class_star: float
class threshhold for PSF stars
match_dist: float
matching distance for checking polution in coordinates pixel
min_sn: float
minimum S/N for PSF stars
Returns
-------
bool
whether succeed to pick PSF stars
"""
# read catalog and keep coordinates
hdulist = fits.open(psffile, mode='update')
cat = hdulist[2].data
if len(cat) < min_nstar:
print("too few psf stars...")
hdulist.close()
return False
if remove_polution:
coord = np.array((cat['XWIN_IMAGE'], cat['YWIN_IMAGE'])).transpose()
coord1 = coord.copy()
index1, index2 = closest_match(coord, coord1, min_dist=match_dist)
index1 = np.array(index1)
mask = np.array([len(ind) == 1 for ind in index2])
index = index1[mask]
if len(index) >= min_nstar:
cat = cat[index]
print("selecting isolated stars:", len(cat))
mask = (cat['flags'] == 0) & (cat['nimaflags_iso'] < 3) & (cat['CLASS_STAR'] > class_star) & (cat['ELLIPTICITY'] < max_elp) & (
cat['SNR_WIN'] > min_sn) & (cat['FWHM_IMAGE'] > fwhm_range[0]) & (cat['FWHM_IMAGE'] < fwhm_range[1])
cat = cat[mask]
if len(cat) < min_nstar:
print("too few psf stars...")
hdulist.close()
return False
if len(cat) > max_nstar:
print('stars in the catalog, needed:', len(cat), max_nstar)
if len(cat) > max_nstar:
indsort = np.argsort(cat['flux_aper'])
cat = cat[indsort[-max_nstar:]]
hdulist[2].data = cat
hdulist.flush()
hdulist.close()
print('final picked stars:', len(cat))
return True
# if save_starpos:
# # f=open(fitsname+'-psfxy.txt','w')
# # for i in range(len(cat)):
# # f.write('%7.2f %7.2f\n' % (cat['XWIN_IMAGE'][i],cat['YWIN_IMAGE'][i]))
# # f.close()
def get_psf(imgfile, whtfile, flgfile, psffile, psf_size=101, degree=3, variability=0.3, fwhm_range=[2.0, 20.0],
max_elp=0.3, sampling=0, min_sn=20.0, detect_thresh=5.0, detect_minarea=5, seeing=0.15, pixel_scale=0.075,
filter_name='gauss_4.0_7x7.conv', phot_aper=10.0, back_size='400,400', nthread=1, **kwd):
# save_starpos=False
""" get PSF profile for a specified image
Parameters
----------
imgfile: str
input image
outdir: str
output directory
psf_size: int
size of the PSF image
degree: int
order of spatially varied PSF
variablility: float
allowed FWHM variability
fwhm_range: array with two elements
FWHM range of selected stars
max_elp: float
max (A-B)/(A+B) for selected stars
sampling: float
sampling step in pixels (0 = Auto)
min_sn: minimum S/N for souces
detect_thresh: float
threshold for detections
detect_minarea: float
min pixel area for detections
seeing: float
PSF fwhm in arcsec
pixel_scale: float
pixel scale in arcsec
filter_name: str
convolve kernel file
phot_aper: float
aperture used for photometry
back_size: array
background grid size
check_plots: bool
output the check plots
nthread: int
number of threads to get PSF
**kwd: keywords for PICK_PSFSTARS
Returns
-------
bool, whether succeed to get PSF profile
"""
# if not system.cmd_exists('sex') or not system.cmd_exists('psfex'):
if which('sex') is None or which('psfex') is None:
raise OSError('No sex or psfex SHELL command found! Please install')
# rootname,_=os.path.splitext(imgfile)
# indpos=str.rfind(rootname,'_img')
# flagfile=rootname[:indpos]+'_flg.fits'
# weightfile=rootname[:indpos]+'_wht.fits'
# imgfile = dm.l1_detector(detector=detector, post="IMG.fits")
# whtfile = dm.l1_detector(detector=detector, post="WHT.fits")
# flgfile = dm.l1_detector(detector=detector, post="FLG.fits")
# psffile = dm.l1_detector(detector=detector, post="PSF.fits")
# fitsname = dm.l1_detector(detector=detector, post="")
head = fits.getheader(imgfile, 0)
gain = head['exptime']
saturate = 50000.0 / gain
# fitsname,_=os.path.splitext(imgfile)
# fitsname=rootname[:indpos]
# if outdir is not None:
# # _,filename=os.path.split(rootname[:indpos])
# # fitsname=os.path.join(outdir,filename)
sexfile = os.path.join(config_path, 'csst_psfex.sex')
covfile = os.path.join(config_path, filter_name)
nnwfile = os.path.join(config_path, 'default.nnw')
parfile = os.path.join(config_path, 'csst_psfex.param')
pexfile = os.path.join(config_path, 'csst.psfex')
command = 'sex ' + imgfile + ' -c ' + sexfile + ' -CATALOG_NAME ' + psffile + ' -PARAMETERS_NAME ' + parfile \
+ ' -DETECT_MINAREA ' + str(detect_minarea) + ' -DETECT_THRESH ' + str(detect_thresh) + ' -FILTER_NAME ' \
+ covfile + ' -WEIGHT_IMAGE ' + whtfile + ' -FLAG_IMAGE ' + flgfile + ' -PHOT_APERTURES ' \
+ str(phot_aper) + ' -GAIN ' + str(gain) + ' -PIXEL_SCALE ' + str(pixel_scale) + ' -SEEING_FWHM ' \
+ str(seeing) + ' -SATUR_LEVEL ' + str(saturate) + ' -STARNNW_NAME ' + nnwfile + ' -BACK_SIZE ' \
+ back_size + ' -NTHREADS ' + str(nthread)
print(command)
os.system(command)
pickflag = pick_psfstars(psffile, min_sn=min_sn,
fwhm_range=fwhm_range, **kwd)
if pickflag is False:
if os.path.isfile(psffile):
os.remove(psffile)
return False
check_plots = False
if not check_plots:
command = 'psfex ' + psffile + ' -c ' + pexfile + ' -PSF_SIZE ' + str(psf_size) + ',' + str(psf_size) \
+ ' -SAMPLE_MINSN ' + str(min_sn) + \
' -CHECKPLOT_DEV NULL -CHECKPLOT_TYPE NONE -CHECKIMAGE_TYPE NONE -PSFVAR_DEGREES ' + str(degree) \
+ ' -SAMPLE_VARIABILITY ' + str(variability) + ' -SAMPLE_FWHMRANGE ' + str(fwhm_range[0]) + ',' \
+ str(fwhm_range[1]) + ' -PSF_SAMPLING ' + str(sampling) + ' -SAMPLE_MAXELLIP ' + str(max_elp) \
+ ' -NTHREADS ' + str(nthread)
# else:
# checknames=fitsname+'psffwhm.png,'+fitsname+'psfelp.png'
# command='psfex '+psffile+' -c '+pexfile+' -PSF_SIZE '+str(psf_size)+','+str(psf_size)+' -SAMPLE_MINSN '+str(min_sn)+' -PSFVAR_DEGREES '+str(degree)+' -SAMPLE_VARIABILITY '+str(variability)+' -SAMPLE_FWHMRANGE '+str(fwhm_range[0])+','+str(fwhm_range[1])+' -PSF_SAMPLING '+str(sampling) + ' -SAMPLE_MAXELLIP '+str(max_elp)+' -NTHREADS '+str(nthread)+' -CHECKPLOT_TYPE FWHM,ELLIPTICITY -CHECKPLOT_DEV PNG -CHECKPLOT_NAME '+checknames
print(command)
os.system(command)
return True
"""
Identifier: KSC-SJ1-MSC-MBI/magfluxconvert.py
Name: stats.py
Description: conversion between magnitude and flux
Author: Hu Zou (zouhu@nao.cas.cn)
Created: 2022-11-11
Modified-History:
2022-11-11, created
2023-12-04, add docstring
"""
import numpy as np
def fluxerr2magerr(flux, fluxerr, asinh=False, filter='u', zp=22.5):
"""
convert flux and flux error to mag and mag error (in pogson or asinh form)
Parameters:
flux, fluxerr in nanamaggie
return mag and magerr
"""
flux = np.array(flux)
fluxerr = np.array(fluxerr)
# f0 = 1.0e9
f0 = 10**(zp/2.5)
nn = flux.size
mag = np.array(np.ones_like(flux)*99.0)
magerr = np.array(np.ones_like(flux)*99.0)
if not asinh:
mask = flux > 0
if mask.any():
mag[mask] = -2.5*np.log10(flux[mask]/f0)
magerr[mask] = 2.5/np.log(10.0)*fluxerr[mask]/flux[mask]
else:
bs = {'u': 1.4e-10, 'g': 0.9e-10,
'r': 1.2e-10, 'i': 1.8e-10, 'z': 7.4e-10}
b = bs[filter]
mag = -(2.5/np.log(10.0))*(np.arcsinh((flux/f0)/(2.0*b))+np.log(b))
magerr = 2.5/np.log(10.0)*(fluxerr/f0)/(2.0*b) / \
np.sqrt(1.0+((flux/f0)/(2.0*b))**2)
return mag, magerr
"""
Identifier: KSC-SJ1-MSC-MBI/stats.py
Name: stats.py
Description: some stats programs
Author: Hu Zou (zouhu@nao.cas.cn)
Created: 2022-11-11
Modified-History:
2022-11-11, created
2023-12-04, add docstring
"""
import astropy.stats as ast
import numpy as np
def valid_coordinates(coordinates, size=None):
""" convert tuple, list or np.ndarray of cooridinates to a 2d ndarray and check the coordinates within an image size
Parameters
----------
coordinates: 2d array for coordinates
size: size of an image
Returns:
----------
coord: reshape coordinates
indcoord: index of coordinates in a size range
"""
if isinstance(coordinates, (list, tuple, np.ndarray)):
coord = np.atleast_2d(coordinates)
if coord.shape[0] != 2 and coord.shape[1] != 2:
raise ValueError(
"coordinates should have at least one axis with 2 elements")
if coord.shape[1] != 2 and coord.shape[0] == 2:
coord = coord.transpose()
else:
raise TypeError(
"coordinates should be list or array of (x,y) pixel positions")
if size is None:
return coord
else:
if len(size) != 2:
raise ValueError("size should have 2 elements")
nx, ny = size
x = coord[:, 0]
y = coord[:, 1]
indcoord = np.arange(coord.shape[0])
good = (x >= 0.5) & (x < nx+0.5) & (y >= 0.5) & (y < ny+0.5)
if np.any(good):
indcoord = indcoord[good]
else:
raise ValueError('coordinates are not in the image range')
return coord, indcoord
def closest_match(coord1, coord2, min_dist=1.0):
""" find closest pairs between two sets of coordinates
Parameters
----------
coord1: 2d array
coordinates to be matched
coord2: 2d array
coordinates matched to
min_dist: float
separation tolerance
Returns
-------
idx1: 1d array
matched index for coord1
idx2: 1d array
matched index for coord2
"""
coord1 = valid_coordinates(coord1)
coord2 = valid_coordinates(coord2)
n1 = len(coord1)
n2 = len(coord2)
index1 = []
index2 = []
x2 = coord2[:, 0]
y2 = coord2[:, 1]
for i in range(n1):
ix1, iy1 = coord1[i]
index = np.where((np.abs(x2-ix1) < min_dist) &
(np.abs(y2-iy1) < min_dist))[0]
nmatch = len(index)
if nmatch < 1:
continue
if nmatch > 1:
x2tmp = x2[index]
y2tmp = y2[index]
dist = ((x2tmp-ix1)**2+(y2tmp-iy1)**2)
indsort = np.argsort(dist)
index1.append(i)
index2.append(index[indsort])
else:
index1.append(i)
index2.append(index)
return index1, index2
def sigmaclip_limitsig(data, error=None, sig_limit=None, **kwd):
""" sigma-clipping the data with upper limits
Parameters
data: 1d array
error: error corresponding to the data
sig_limit: the upper limit of the std of the clipped data
**kwd: keywords for astropy.stats.sigma_clip
OUPUTS:
mdata: masked data
"""
data = np.array(data)
mdata = ast.sigma_clip(data, **kwd)
if sig_limit is not None:
while True:
med = np.ma.median(mdata)
sig = np.ma.std(mdata)
if sig < sig_limit:
break
index = np.ma.argmax(np.ma.abs(mdata-med))
mdata.mask[index] = True
return mdata
def weighted_mean(x, sigmax, weight_square=True):
"""
Calculate the mean and estimated errors for a set of data points
DESCRIPTION:
This routine is adapted from Program 5-1, XFIT, from "Data Reduction
and Error Analysis for the Physical Sciences", p. 76, by Philip R.
Bevington, McGraw Hill. This routine computes the weighted mean using
Instrumental weights (w=1/sigma^2).
Parameters:
x - Array of data points
sigmax - array of standard deviations for data points
weight_square - if True, weight is invariance, else the reciprocal of the error
Returns:
xmean - weighted mean
sigmam - standard deviation of mean
stdm - standard deviation of data
"""
x = np.atleast_1d(x).copy()
sigmax = np.atleast_1d(sigmax).copy()
if len(x) == 1:
xmean = x
sigmam = sigmax
stdm = sigmax
else:
weight = 1.0/sigmax**2
weight1 = weight
if not weight_square:
weight1 = 1.0/sigmax
wsum = weight.sum()
xmean = (weight1*x).sum()/weight1.sum()
sigmam = np.sqrt(1.0/wsum)
stdm = np.sqrt(np.sum((x-xmean)**2)/len(x))
return xmean, sigmam, stdm
"""_summary_
"""
import os
from glob import glob
import mpi4py.MPI as MPI
# from .L1_pipeline.csst_msc_instrument.csst_msc_mbi_instrument import core_msc_l1_mbi_instrument
from L1_pipeline.csst_msc_instrument.csst_msc_mbi_instrument import core_msc_l1_mbi_instrument
# core_msc_l1_mbi_instrument(
# image_path=r"D:\Zhou\Desktop\data\data_09.fits",
# bias_path=r"D:\Zhou\Desktop\data\bias_09.fits",
# dark_path=r"D:\Zhou\Desktop\data\dark_09.fits",
# flat_path=r"D:\Zhou\Desktop\data\flat_09.fits",
# shutter_path=r"D:\Zhou\Desktop\data\csst_msc_ms_shutter_09_000001.fits",
# image_output_path=r"D:\Zhou\Desktop\data\result\image_test_output.fits",
# weight_output_path=r"D:\Zhou\Desktop\data\result\weight_test_output.fits",
# flag_output_path=r"D:\Zhou\Desktop\data\result\flag_test_output.fits",
# )
\ No newline at end of file
def run_csst_msc_instrument(image_path,
output_dir,
calib_data_path="/public/home/fangyuedong/project/calib_data/"):
"""_summary_
Args:
image_path (_type_): _description_
output_dir (_type_): _description_
calib_data_path (str, optional): _description_. Defaults to "/public/home/fangyuedong/project/calib_data/".
"""
img_filename = os.path.splitext(os.path.basename(image_path))[0]
chip_label = img_filename.split("_")[-3]
# Get paths to corresponding calibration files
bias_path = os.path.join(calib_data_path, "bias_%s.fits" % (chip_label))
dark_path = os.path.join(calib_data_path, "dark_%s.fits" % (chip_label))
flat_path = os.path.join(calib_data_path, "flat_%s.fits" % (chip_label))
output_img_types = ['img', 'wht', 'flg']
output_paths = []
for i in range(len(output_img_types)):
output_paths.append(os.path.join(output_dir,
img_filename[:-6] + output_img_types[i] + "_L1.fits"))
core_msc_l1_mbi_instrument(
image_path=image_path,
bias_path=bias_path,
dark_path=dark_path,
flat_path=flat_path,
shutter_path="/dummy/csst_msc_ms_shutter_09_000001.fits", # Didn't use at this moment
image_output_path=output_paths[0],
weight_output_path=output_paths[1],
flag_output_path=output_paths[2],
)
def genearte_path_list_for_one_pointing(input_dir,
pointing_label,
chip_label_list=None):
"""_summary_
Args:
input_dir (_type_): _description_
pointing_label (_type_): _description_
chip_label_list (_type_, optional): _description_. Defaults to None.
Returns:
_type_: _description_
"""
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_*_' + '*_*')
else:
image_path_list = []
for chip_label in chip_label_list:
image_path = glob(pointing_dir + '/CSST_MSC_MS_SCIE_*_' +
chip_label + '_*')[0]
image_path_list.append(image_path)
return image_path_list
def run_pointing_list(input_dir,
pointing_label_list,
output_dir,
calib_data_path="/public/home/fangyuedong/project/calib_data/",
chip_label_list=None):
"""_summary_
Args:
input_dir (_type_): _description_
pointing_label_list (_type_): _description_
output_dir (_type_): _description_
calib_data_path (str, optional): _description_. Defaults to "/public/home/fangyuedong/project/calib_data/".
chip_label_list (_type_, optional): _description_. Defaults to None.
"""
image_path_list = []
output_path_list = []
try:
if not os.path.exists(output_dir):
os.makedirs(output_dir)
except OSError:
pass
for pointing_label in pointing_label_list:
output_pointing_dir = os.path.join(output_dir, pointing_label)
try:
if not os.path.exists(output_pointing_dir):
os.makedirs(output_pointing_dir)
except OSError:
pass
temp_img_path_list = genearte_path_list_for_one_pointing(input_dir=input_dir,
pointing_label=pointing_label,
chip_label_list=chip_label_list)
image_path_list = image_path_list + temp_img_path_list
output_path_list = output_path_list + \
[output_pointing_dir] * len(temp_img_path_list)
comm = MPI.COMM_WORLD
ind_thread = comm.Get_rank()
num_thread = comm.Get_size()
for i in range(len(image_path_list)):
if i % num_thread != ind_thread:
continue
image_path = image_path_list[i]
output_path = output_path_list[i]
run_csst_msc_instrument(image_path=image_path,
output_dir=output_path,
calib_data_path=calib_data_path)
if __name__ == "__main__":
input_dir = "/public/share/yangxuliu/CSSOSDataProductsSims/outputs_50sqDeg/50sqDeg_Photo_W1/"
pointing_label_list = ["MSC_0000000", "MSC_0000001",
"MSC_0000002", "MSC_0000003", "MSC_0000004", "MSC_0000005"]
# chip_label_list = ["08"]
chip_label_list = None
output_dir = "/public/home/fangyuedong/project/50sqDeg_L1_outputs"
run_pointing_list(input_dir=input_dir,
pointing_label_list=pointing_label_list,
output_dir=output_dir,
chip_label_list=chip_label_list)
"""_summary_
"""
import os
from glob import glob
from astropy.io import fits
from L1_pipeline.ref_combine import combine_images
ref_path = "/public/share/yangxuliu/CSSOSDataProductsSims/outputs_cali/"
output_path = "/public/home/fangyuedong/project/calib_data"
def combine_ref_func(ref_path, output_path, num="01"):
bias_path_list = glob(ref_path + '*/CSST_MSC_MS_BIAS_*_' + num + '_*')
dark_path_list = glob(ref_path + '*/CSST_MSC_MS_DARK_*_' + num + '_*')
flat_path_list = glob(ref_path + '*/CSST_MSC_MS_FLAT_*_' + num + '_*')
def combine_ref_func(
ref_path: str,
output_path: str,
chip_id: str = "01"
) -> None:
"""_summary_
Args:
ref_path (str): _description_
output_path (str): _description_
chip_id (str, optional): _description_. Defaults to "01".
"""
bias_path_list = glob(ref_path + '*/CSST_MSC_MS_BIAS_*_' + chip_id + '_*')
dark_path_list = glob(ref_path + '*/CSST_MSC_MS_DARK_*_' + chip_id + '_*')
flat_path_list = glob(ref_path + '*/CSST_MSC_MS_FLAT_*_' + chip_id + '_*')
bias, dark, flat = combine_images(b_p_lst=bias_path_list,
d_p_lst=dark_path_list,
f_p_lst=flat_path_list, )
bias_out_path = os.path.join(output_path, "bias_" + str(num) + ".fits")
dark_out_path = os.path.join(output_path, "dark_" + str(num) + ".fits")
flat_out_path = os.path.join(output_path, "flat_" + str(num) + ".fits")
bias_out_path = os.path.join(output_path, "bias_" + str(chip_id) + ".fits")
dark_out_path = os.path.join(output_path, "dark_" + str(chip_id) + ".fits")
flat_out_path = os.path.join(output_path, "flat_" + str(chip_id) + ".fits")
hdu = fits.PrimaryHDU(bias)
hdu.writeto(bias_out_path)
hdu.writeto(bias_out_path, overwrite=True)
hdu = fits.PrimaryHDU(dark)
hdu.writeto(dark_out_path)
hdu.writeto(dark_out_path, overwrite=True)
hdu = fits.PrimaryHDU(flat)
hdu.writeto(flat_out_path)
hdu.writeto(flat_out_path, overwrite=True)
if __name__ == "__main__":
ref_path = "/public/share/yangxuliu/CSSOSDataProductsSims/outputs_cali/"
output_path = "/public/home/fangyuedong/project/calib_data"
num = '08'
REFERENCE_PATH = "/public/share/yangxuliu/CSSOSDataProductsSims/outputs_cali/"
OUTPUT_DIR = "/public/home/fangyuedong/project/calib_data"
# CHIP_IDS = '08'
CHIP_IDS = ['06', '07', '08', '09', '11', '12', '13', '14', '15',
'16', '17', '18', '19', '20', '22', '23', '24', '25']
for CHIP_ID in CHIP_IDS:
combine_ref_func(
ref_path=ref_path,
output_path=output_path,
num=num
ref_path=REFERENCE_PATH,
output_path=OUTPUT_DIR,
chip_id=CHIP_ID
)
#! /bin/bash
#SBATCH -J L1_INST
#SBATCH -N 1
#SBATCH --ntasks-per-node=24
#SBATCH -p batch
#SBATCH --mem=240G
module load mpi/openmpi/4.0.2/gcc-7.3.1
date
srun hostname -s | sort -n | awk -F"-" '{print $2}' | uniq > pnodes
mpirun -mca pml ucx -x UCX_NET_DEVICES=mlx5_0:1 -machinefile pnodes -np 24 --map-by node python run_csst_msc_instrument.py
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment