Skip to content
GitLab
Commit ef8c4609 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

modified fits header etc.

parent 3f932427
Hide whitespace changes
Inline Side-by-side
ObservationSim/ObservationSim.py
View file @ ef8c4609
......@@ -4,8 +4,8 @@ from Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
from MockObject import Catalog, MockObject, Star, Galaxy, Quasar, calculateSkyMap_split_g
from PSF import PSFGauss, PSFInterp, FieldDistortion
from _util import makeSubDir, getShearFiled, makeSubDir_PointingList
from astropy.time import Time as asTime
from astropy.io import fits
from datetime import datetime
import numpy as np
import mpi4py.MPI as MPI
import galsim
......@@ -47,7 +47,7 @@ class Observation(object):
self.g1_field, self.g2_field, self.nshear = getShearFiled(config=self.config)
def runOneChip(self, chip, filt, chip_output, wcs_fp=None, psf_model=None, pointing_ID=0, ra_cen=None, dec_cen=None, img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None, cat_dir=None, sed_dir=None):
def runOneChip(self, chip, filt, chip_output, wcs_fp=None, psf_model=None, pointing_ID=0, ra_cen=None, dec_cen=None, img_rot=None, exptime=150., timestamp_obs=1621915200, pointing_type='MS', input_cat_name=None, shear_cat_file=None, cat_dir=None, sed_dir=None):
if (ra_cen is None) or (dec_cen is None):
ra_cen = self.config["ra_center"]
......@@ -80,169 +80,189 @@ class Observation(object):
skyfile = os.path.join(self.path_dict["data_dir"], 'skybackground/sky_emiss_hubble_50_50_A.dat')
sky_map = calculateSkyMap_split_g(xLen=chip.npix_x, yLen=chip.npix_y, blueLimit=filt.blue_limit, redLimit=filt.red_limit, skyfn=skyfile, conf=chip.sls_conf, pixelSize=chip.pix_scale, isAlongY=0)
# Load catalogues and templates
self.cat = Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir, pRa=ra_cen, pDec=dec_cen, rotation=img_rot, template_dir=self.path_dict["template_dir"])
self.nobj = len(self.cat.objs)
# Loop over objects
missed_obj = 0
bright_obj = 0
dim_obj = 0
for j in range(self.nobj):
# if j >= 20:
# break
obj = self.cat.objs[j]
# Load SED
if obj.type == 'star':
normF = chip.normF_star
try:
if pointing_type == 'MS':
# Load catalogues and templates
self.cat = Catalog(config=self.config, chip=chip, cat_dir=cat_dir, sed_dir=sed_dir, pRa=ra_cen, pDec=dec_cen, rotation=img_rot, template_dir=self.path_dict["template_dir"])
self.nobj = len(self.cat.objs)
# Loop over objects
missed_obj = 0
bright_obj = 0
dim_obj = 0
for j in range(self.nobj):
# if j >= 20:
# break
obj = self.cat.objs[j]
# Load SED
if obj.type == 'star':
normF = chip.normF_star
try:
obj.load_SED(
survey_type=chip.survey_type,
normFilter=normF,
target_filt=filt,
sed_lib=self.cat.tempSED_star)
except Exception as e:
print(e)
continue
elif obj.type == 'galaxy': # or obj.type == quasar
normF = chip.normF_galaxy
obj.load_SED(
sed_path=sed_dir,
survey_type=chip.survey_type,
sed_templates=self.cat.tempSed_gal,
normFilter=normF,
target_filt=filt,
sed_lib=self.cat.tempSED_star)
except Exception as e:
print(e)
continue
elif obj.type == 'galaxy': # or obj.type == quasar
normF = chip.normF_galaxy
obj.load_SED(
sed_path=sed_dir,
survey_type=chip.survey_type,
sed_templates=self.cat.tempSed_gal,
normFilter=normF,
target_filt=filt)
elif obj.type == 'quasar':
normF = chip.normF_galaxy
obj.load_SED(
sed_path=sed_dir,
survey_type=chip.survey_type,
sed_templates=self.cat.tempSed_gal,
normFilter=normF,
target_filt=filt)
# Exclude very bright/dim objects (for now)
if filt.is_too_bright(mag=obj.getMagFilter(filt)):
# print("obj too birght!!", flush=True)
if obj.type != 'galaxy':
bright_obj += 1
target_filt=filt)
elif obj.type == 'quasar':
normF = chip.normF_galaxy
obj.load_SED(
sed_path=sed_dir,
survey_type=chip.survey_type,
sed_templates=self.cat.tempSed_gal,
normFilter=normF,
target_filt=filt)
# Exclude very bright/dim objects (for now)
if filt.is_too_bright(mag=obj.getMagFilter(filt)):
# print("obj too birght!!", flush=True)
if obj.type != 'galaxy':
bright_obj += 1
obj.unload_SED()
continue
if filt.is_too_dim(mag=obj.getMagFilter(filt)):
# print("obj too dim!!", flush=True)
dim_obj += 1
obj.unload_SED()
# print(obj.getMagFilter(filt))
continue
if filt.is_too_dim(mag=obj.getMagFilter(filt)):
# print("obj too dim!!", flush=True)
dim_obj += 1
obj.unload_SED()
# print(obj.getMagFilter(filt))
continue
if self.config["shear_method"] == "constant":
if obj.type == 'star':
g1, g2 = 0, 0
else:
g1, g2 = self.g1_field, self.g2_field
elif self.config["shear_method"] == "extra":
# TODO: every object with individual shear from input catalog(s)
g1, g2 = self.g1_field[j], self.g2_field[j]
pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False)
if pos_img.x == -1 or pos_img.y == -1:
# Exclude object which is outside the chip area (after field distortion)
# print("obj missed!!")
missed_obj += 1
obj.unload_SED()
continue
if self.config["shear_method"] == "constant":
if obj.type == 'star':
g1, g2 = 0, 0
else:
g1, g2 = self.g1_field, self.g2_field
elif self.config["shear_method"] == "extra":
# TODO: every object with individual shear from input catalog(s)
g1, g2 = self.g1_field[j], self.g2_field[j]
pos_img, offset, local_wcs = obj.getPosImg_Offset_WCS(img=chip.img, fdmodel=self.fd_model, chip=chip, verbose=False)
if pos_img.x == -1 or pos_img.y == -1:
# Exclude object which is outside the chip area (after field distortion)
# print("obj missed!!")
missed_obj += 1
obj.unload_SED()
continue
# Draw object & update output catalog
try:
if chip.survey_type == "photometric":
isUpdated, pos_shear = obj.drawObj_multiband(
tel=self.tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=g1,
g2=g2,
exptime=exptime)
elif chip.survey_type == "spectroscopic":
isUpdated, pos_shear = obj.drawObj_slitless(
tel=self.tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=g1,
g2=g2,
exptime=exptime,
normFilter=normF)
if isUpdated:
# TODO: add up stats
chip_output.cat_add_obj(obj, pos_img, pos_shear, g1, g2)
# Draw object & update output catalog
try:
if chip.survey_type == "photometric":
isUpdated, pos_shear = obj.drawObj_multiband(
tel=self.tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=g1,
g2=g2,
exptime=exptime)
elif chip.survey_type == "spectroscopic":
isUpdated, pos_shear = obj.drawObj_slitless(
tel=self.tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=g1,
g2=g2,
exptime=exptime,
normFilter=normF)
if isUpdated:
# TODO: add up stats
chip_output.cat_add_obj(obj, pos_img, pos_shear, g1, g2)
pass
else:
# print("object omitted", flush=True)
continue
except Exception as e:
print(e)
pass
else:
# print("object omitted", flush=True)
continue
except Exception as e:
print(e)
pass
# Unload SED:
obj.unload_SED()
del obj
# Unload SED:
obj.unload_SED()
del obj
del psf_model
del self.cat
del psf_model
del self.cat
print("check running:1: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing_ID, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
# Detector Effects
# ===========================================================
chip.img = chip.addNoise(config=self.config, tel=self.tel, filt=filt, img=chip.img, sky_map=sky_map)
chip.img = chip.addEffects(config=self.config, img=chip.img, chip_output=chip_output, filt=filt, pointing_ID=pointing_ID)
# chip.img = chip.addEffects(config=self.config, img=chip.img, chip_output=chip_output, filt=filt, pointing_ID=pointing_ID)
# whether to output zero, dark, flat calibration images.
chip.img = chip.addEffects(
config=self.config,
img=chip.img,
chip_output=chip_output,
filt=filt,
ra_cen=ra_cen,
dec_cen=dec_cen,
img_rot=img_rot,
pointing_ID=pointing_ID,
timestamp_obs=timestamp_obs,
pointing_type=pointing_type)
h_prim = generatePrimaryHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
pointNum = str(pointing_ID),
ra=ra_cen,
dec=dec_cen,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
date=self.config["date_obs"],
time_obs=self.config["time_obs"],
im_type = 'MS')
h_ext = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra_cen,
dec=dec_cen,
pa=img_rot.deg,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='SCI')
chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
hdu1 = fits.PrimaryHDU(header=h_prim)
hdu2 = fits.ImageHDU(chip.img.array, header=h_ext)
hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
if pointing_type == 'MS':
datetime_obs = datetime.fromtimestamp(timestamp_obs)
date_obs = datetime_obs.strftime("%y%m%d")
time_obs = datetime_obs.strftime("%H%M%S")
h_prim = generatePrimaryHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
pointNum = str(pointing_ID),
ra=ra_cen,
dec=dec_cen,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
# date=self.config["date_obs"],
# time_obs=self.config["time_obs"],
date=date_obs,
time_obs=time_obs,
im_type='MS')
h_ext = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra_cen,
dec=dec_cen,
pa=img_rot.deg,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw')
chip.img = galsim.Image(chip.img.array, dtype=np.uint16)
hdu1 = fits.PrimaryHDU(header=h_prim)
hdu2 = fits.ImageHDU(chip.img.array, header=h_ext)
hdu1 = fits.HDUList([hdu1, hdu2])
fname = os.path.join(chip_output.subdir, h_prim['FILENAME'] + '.fits')
hdu1.writeto(fname, output_verify='ignore', overwrite=True)
print("# objects that are too bright %d out of %d"%(bright_obj, self.nobj))
print("# objects that are too dim %d out of %d"%(dim_obj, self.nobj))
print("# objects that are missed %d out of %d"%(missed_obj, self.nobj))
del chip.img
print("check running:2: pointing-{:} chip-{:} pid-{:} memory-{:6.2}GB".format(pointing_ID, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ), flush=True)
print("# objects that are too bright %d out of %d"%(bright_obj, self.nobj))
print("# objects that are too dim %d out of %d"%(dim_obj, self.nobj))
print("# objects that are missed %d out of %d"%(missed_obj, self.nobj))
def runExposure(self, ra_cen=None, dec_cen=None, pointing_ID=0, img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None, oneChip=None):
def runExposure(self, ra_cen=None, dec_cen=None, timestamp_obs=1621915200, pointing_ID=0, pointing_type='MS', img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None, oneChip=None):
if (ra_cen == None) or (dec_cen == None):
ra_cen = self.config["ra_center"]
......@@ -269,6 +289,7 @@ class Observation(object):
chip=chip,
filt=filt,
exptime=exptime,
pointing_type=pointing_type,
pointing_ID=pointing_ID,
subdir=sub_img_dir,
prefix=prefix)
......@@ -282,11 +303,13 @@ class Observation(object):
dec_cen=dec_cen,
img_rot=img_rot,
exptime=exptime,
timestamp_obs=timestamp_obs,
pointing_type=pointing_type,
cat_dir=self.path_dict["cat_dir"],
sed_dir=self.path_dict["SED_dir"])
print("finished running chip#%d..."%(chip.chipID), flush=True)
def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None):
def runExposure_MPI_PointingList(self, ra_cen=None, dec_cen=None, pRange=None, timestamp_obs=[1621915200], pointing_type=['MS'], img_rot=None, exptime=150., input_cat_name=None, shear_cat_file=None):
comm = MPI.COMM_WORLD
ind_thread = comm.Get_rank()
num_thread = comm.Get_size()
......@@ -302,6 +325,11 @@ class Observation(object):
else:
pStart = 0
# TEMP
if len(timestamp_obs) == 1:
timestamp_obs *= len(pRange)
pointing_type *= len(pRange)
for ipoint in range(len(ra_cen)):
for ichip in range(nchips_per_fp):
i = ipoint*nchips_per_fp + ichip
......@@ -322,6 +350,7 @@ class Observation(object):
chip=chip,
filt=filt,
exptime=exptime,
pointing_type=pointing_type[ipoint],
pointing_ID=pointing_ID,
subdir=sub_img_dir,
prefix=prefix)
......@@ -333,7 +362,9 @@ class Observation(object):
ra_cen=ra_cen[ipoint],
dec_cen=dec_cen[ipoint],
img_rot=img_rot,
exptime=exptime,
exptime=exptime,
timestamp_obs=timestamp_obs[ipoint],
pointing_type=pointing_type[ipoint],
cat_dir=self.path_dict["cat_dir"],
sed_dir=self.path_dict["SED_dir"])
print("finished running chip#%d..."%(chip.chipID), flush=True)
\ No newline at end of file
ObservationSim/Pointing.py
View file @ ef8c4609
import os
import numpy as np
from Config import ReadConfig
from datetime import datetime
work_dir = "/public/home/fangyuedong/CSST_new_framework/test/"
work_dir = "/public/home/fangyuedong/CSST/test/"
data_dir = "/data/simudata/CSSOSDataProductsSims/data/"
config_file = os.path.join(work_dir, "ObservationSim.cfg")
......@@ -25,5 +26,23 @@ f.close()
pRA = np.array(pRA)
pDEC = np.array(pDEC)
# Create calibration pointings
ncal = 1 # Define the number of calibration pointings
pointing_type = ['MS']*len(pRA)
pRA = np.append(pRA[:ncal], pRA)
pDEC = np.append(pDEC[:ncal], pDEC)
pointing_type = ['CAL']*ncal + pointing_type
# Calculate starting time(s)
t0 = datetime(2021, 5, 25, 12, 0, 0)
t = datetime.timestamp(t0)
timestamp_obs = []
for i in range(len(pointing_type)):
timestamp_obs.append(t)
if pointing_type[i] == 'CAL':
t += 3 * 5 * 60 # 3 calibration exposure
elif pointing_type[i] == 'MS':
t += 5 * 60
# Define the range of pointing list
pRange = range(0, 3)
\ No newline at end of file
pRange = range(0, 2)
\ 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