Commit 4c4aafe6 authored by Fang Yuedong's avatar Fang Yuedong
Browse files

1. add project_cycle, run_counter to config file. Modified obsid to reflect...

1. add project_cycle, run_counter to config file. Modified obsid to reflect this change. Applied the same obsid to ".cat" and ".log" files.
2. limiting mags are now based on the field types (WIDE or DEEP) which can be specified in pointing list.
3. Stars are now drawn directly as PSF.
4. Reorgnaized lists of filter types in Instrument._util.
parent 898bdd0e
Loading
Loading
Loading
Loading
+2 −2
Original line number Original line Diff line number Diff line
@@ -344,8 +344,8 @@ class Catalog(CatalogBase):
            )
            )
        for istars in range(nstars):
        for istars in range(nstars):
            # (TEST)
            # (TEST)
            if istars > 100:
            # if istars > 100:
                break
            #     break


            param = self.initialize_param()
            param = self.initialize_param()
            param['ra'] = ra_arr[istars]
            param['ra'] = ra_arr[istars]
+8 −4
Original line number Original line Diff line number Diff line
import os
import os
import logging
import logging
import ObservationSim.Config._util as _util


class ChipOutput(object):
class ChipOutput(object):
    def __init__(self, config, focal_plane, chip, filt, imgKey0="", imgKey1="", imgKey2="", exptime=150., mjdTime="", ra_cen=None, dec_cen=None, pointing_type='MS', pointing_ID='0', subdir="./", prefix=""):
    def __init__(self, config, focal_plane, chip, filt, imgKey0="", imgKey1="", imgKey2="", exptime=150., mjdTime="", ra_cen=None, dec_cen=None, pointing_type='SCI', pointing_ID='0', subdir="./", prefix=""):
        self.focal_plane = focal_plane
        self.focal_plane = focal_plane
        self.chip = chip
        self.chip = chip
        self.filt = filt
        self.filt = filt
@@ -20,12 +21,15 @@ class ChipOutput(object):
            self.dec_cen = config["obs_setting"]["dec_center"]
            self.dec_cen = config["obs_setting"]["dec_center"]


        self.chipLabel = focal_plane.getChipLabel(chip.chipID)
        self.chipLabel = focal_plane.getChipLabel(chip.chipID)
        self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat"
        obs_id = _util.get_obs_id(img_type=pointing_type, project_cycle=config["project_cycle"], run_counter=config["run_counter"], pointing_num=pointing_ID)
        # self.cat_name = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".cat"
        self.cat_name = "MSC_%s_chip_%s_filt_%s"%(obs_id, self.chipLabel, filt.filter_type) + ".cat"


        self.subdir = subdir
        self.subdir = subdir


        # Setup logger for each chip
        # Setup logger for each chip
        logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".log"
        # logger_filename = "MSC_1%s_chip_%s_filt_%s"%(str(pointing_ID).rjust(8, '0'), self.chipLabel, filt.filter_type) + ".log"
        logger_filename = "MSC_%s_chip_%s_filt_%s"%(obs_id, self.chipLabel, filt.filter_type) + ".log"
        self.logger = logging.getLogger()
        self.logger = logging.getLogger()
        fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
        fh = logging.FileHandler(os.path.join(self.subdir, logger_filename), mode='w+', encoding='utf-8')
        fh.setLevel(logging.DEBUG)
        fh.setLevel(logging.DEBUG)
@@ -56,7 +60,7 @@ class ChipOutput(object):
        self.hdr += additional_column_names
        self.hdr += additional_column_names
    
    
    def create_output_file(self):
    def create_output_file(self):
        if self.pointing_type == 'MS':
        if self.pointing_type == 'SCI':
            self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
            self.cat = open(os.path.join(self.subdir, self.cat_name), "w")
            self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
            self.logger.info("Creating catalog file %s ...\n"%(os.path.join(self.subdir, self.cat_name)))
            if not self.hdr.endswith("\n"):
            if not self.hdr.endswith("\n"):
+11 −7
Original line number Original line Diff line number Diff line
@@ -17,6 +17,7 @@ from astropy.coordinates import SkyCoord
from astropy.wcs.utils import fit_wcs_from_points
from astropy.wcs.utils import fit_wcs_from_points
from astropy.time import Time
from astropy.time import Time
from astropy import wcs
from astropy import wcs
from ObservationSim.Config._util import get_obs_id, get_file_type


from datetime import datetime
from datetime import datetime
# import socket
# import socket
@@ -341,7 +342,7 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r




#TODO project_cycle is temporary, is not in header defined, delete in future
#TODO project_cycle is temporary, is not in header defined, delete in future
def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, chip_name="01"):
def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec = -40, pixel_scale = 0.074, date='200930', time_obs='120000', im_type = 'MS', exptime=150., sat_pos = [0.,0.,0.], sat_vel = [0., 0., 0.], project_cycle=6, run_counter=0, chip_name="01"):


    # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0
    # array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0


@@ -376,18 +377,21 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
    h_prim['OBJ_RA'] = ra
    h_prim['OBJ_RA'] = ra
    h_prim['OBJ_DEC'] = dec
    h_prim['OBJ_DEC'] = dec


    obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}
    # obs_type = {'SCI': '01', 'BIAS': '03', 'DARK': '07', 'FLAT': '11', 'CRS': '98', 'CRD': '99'}


    OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')
    # # OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + pointNum.rjust(7,'0')
    # OBS_id = '1'+ obs_type[im_type] + str(int(project_cycle)) + str(int(run_counter)).rjust(2, '0') + pointNum.rjust(5,'0')
    OBS_id = get_obs_id(img_type=im_type, project_cycle=project_cycle, run_counter=run_counter, pointing_num=pointNum)


    h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7,'0')
    h_prim['OBJECT'] = str(int(project_cycle)) + pointNum.rjust(7,'0')
    h_prim['OBSID'] = OBS_id
    h_prim['OBSID'] = OBS_id
    # h_prim['TELFOCUS'] = 'f/14'
    # h_prim['TELFOCUS'] = 'f/14'
    h_prim['EXPTIME'] = exptime
    h_prim['EXPTIME'] = exptime
    
    
    # Define file types
    # # Define file types
    file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
    # file_type = {'SCI':'SCIE', 'BIAS':'BIAS', 'DARK':'DARK', 'FLAT':'FLAT', 'CRS':'CRS', 'CRD':'CRD','CALS':'CALS','CALF':'CALF'}
    h_prim['FILETYPE'] = file_type[im_type]
    # h_prim['FILETYPE'] = file_type[im_type]
    h_prim['FILETYPE'] = get_file_type(img_type=im_type)


    co = coord.SkyCoord(ra, dec, unit='deg')
    co = coord.SkyCoord(ra, dec, unit='deg')


@@ -422,7 +426,7 @@ def generatePrimaryHeader(xlen = 9216, ylen = 9232, pointNum = '1', ra = 60, dec
    file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
    file_end_time = end_time_str[0:4] + end_time_str[5:7]+end_time_str[8:10] + end_time_str[11:13] + end_time_str[14:16] + end_time_str[17:19]
    # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
    # h_prim['FILENAME'] = 'CSST_MSC_MS_' + im_type + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + CCDID[
    #     k - 1].rjust(2, '0') + '_L0_V01'
    #     k - 1].rjust(2, '0') + '_L0_V01'
    h_prim['FILENAME'] = 'CSST_MSC_MS_' + file_type[im_type] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'
    h_prim['FILENAME'] = 'CSST_MSC_MS_' + h_prim['FILETYPE'] + '_' + file_start_time + '_' + file_end_time + '_' + OBS_id + '_' + chip_name + '_L0_V01'




    h_prim['POSI0_X'] = sat_pos[0]
    h_prim['POSI0_X'] = sat_pos[0]
+33 −3
Original line number Original line Diff line number Diff line
@@ -2,8 +2,10 @@ import numpy as np
import galsim
import galsim
from astropy.time import Time
from astropy.time import Time


import ObservationSim.Instrument._util as _util

class Pointing(object):
class Pointing(object):
    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='MS'):
    def __init__(self, id=0, ra=0., dec=0., img_pa=0., timestamp=1621915200, sat_x=0., sat_y=0., sat_z=0., sun_x=0., sun_y=0., sun_z=0., sat_vx=0., sat_vy=0., sat_vz=0., exp_time=150., pointing_type='SCI'):
        self.id = id
        self.id = id
        self.ra = ra
        self.ra = ra
        self.dec = dec
        self.dec = dec
@@ -14,9 +16,33 @@ class Pointing(object):
        self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz
        self.sat_vx, self.sat_vy, self.sat_vz = sat_vx, sat_vy, sat_vz
        self.exp_time = exp_time
        self.exp_time = exp_time
        self.pointing_type = pointing_type
        self.pointing_type = pointing_type
        self.survey_field_type = 'WIDE'
        self.jdt = 0. 
        self.jdt = 0. 
    
    
    def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='MS'):
    def get_full_depth_exptime(self, filter_type):
        if self.survey_field_type == 'WIDE':
            if filter_type in _util.SPEC_FILTERS:
                return 150. * 4
            else:
                if filter_type.lower() in ['nuv', 'y']:
                    return 150. * 4
                elif filter_type.lower() in ['u', 'g', 'r', 'i', 'z']:
                    return 150. * 2
                else:
                    return max(150., self.exp_time) # [TODO] for FGS
        elif self.survey_field_type == 'DEEP':
            if filter_type in _util.SPEC_FILTERS:
                return 250. * 4 * 4
            else:
                if filter_type.lower() in ['nuv', 'y']:
                    return 250. * 4 * 4
                elif filter_type.lower() in ['u', 'g', 'r', 'i', 'z']:
                    return 250. * 2 * 4
                else:
                    return max(150., self.exp_time) # [TODO] for FGS


    def read_pointing_columns(self, columns, id=0, t=1621915200, pointing_type='SCI'):
        self.id = id
        self.id = id
        col_len = len(columns)
        col_len = len(columns)
        self.ra = float(columns[0])
        self.ra = float(columns[0])
@@ -38,5 +64,9 @@ class Pointing(object):
            self.sat_vy = float(columns[16])
            self.sat_vy = float(columns[16])
            self.sat_vz = float(columns[17])
            self.sat_vz = float(columns[17])
            self.exp_time = float(columns[18])
            self.exp_time = float(columns[18])
            is_deep = float(columns[19])
            # [TODO] Can also define other survey types
            if is_deep != -1.0:
                self.survey_field_type = "DEEP"
        else:
        else:
            self.timestamp = t
            self.timestamp = t
+17 −6
Original line number Original line Diff line number Diff line
@@ -3,6 +3,7 @@ import os
import numpy as np
import numpy as np
import pickle
import pickle
import json
import json
import ObservationSim.Instrument._util as _util
from astropy.table import Table
from astropy.table import Table
from numpy.random import Generator, PCG64
from numpy.random import Generator, PCG64
from astropy.io import fits
from astropy.io import fits
@@ -142,9 +143,9 @@ class Chip(FocalPlane):
        return rowID, colID
        return rowID, colID


    def _getSurveyType(self):
    def _getSurveyType(self):
        if self.filter_type in ["GI", "GV", "GU"]:
        if self.filter_type in _util.SPEC_FILTERS:
            return "spectroscopic"
            return "spectroscopic"
        elif self.filter_type in ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']:
        elif self.filter_type in _util.PHOT_FILTERS:
            return "photometric"
            return "photometric"
        # elif self.filter_type in ["FGS"]:
        # elif self.filter_type in ["FGS"]:
        #     return "FGS"
        #     return "FGS"
@@ -190,7 +191,7 @@ class Chip(FocalPlane):
    def getChipFilter(self, chipID=None):
    def getChipFilter(self, chipID=None):
        """Return the filter index and type for a given chip #(chipID)
        """Return the filter index and type for a given chip #(chipID)
        """
        """
        filter_type_list = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
        filter_type_list = _util.ALL_FILTERS
        if chipID == None:
        if chipID == None:
            chipID = self.chipID
            chipID = self.chipID


@@ -228,7 +229,7 @@ class Chip(FocalPlane):
        for i in range(4):
        for i in range(4):
            x = xcen + sign_x[i] * self.npix_x / 2.
            x = xcen + sign_x[i] * self.npix_x / 2.
            y = ycen + sign_y[i] * self.npix_y / 2.
            y = ycen + sign_y[i] * self.npix_y / 2.
            x, y = rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
            x, y = _util.rotate_conterclockwise(x0=xcen, y0=ycen, x=x, y=y, angle=self.rotate_angle)
            xmin, xmax = min(xmin, x), max(xmax, x)
            xmin, xmax = min(xmin, x), max(xmax, x)
            ymin, ymax = min(ymin, y), max(ymax, y)
            ymin, ymax = min(ymin, y), max(ymax, y)
        return galsim.BoundsD(xmin, xmax, ymin, ymax)
        return galsim.BoundsD(xmin, xmax, ymin, ymax)
@@ -267,7 +268,7 @@ class Chip(FocalPlane):
        noise = self.dark_noise * exptime + self.read_noise**2
        noise = self.dark_noise * exptime + self.read_noise**2
        return noise
        return noise


    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='MS', sky_map=None, tel=None, logger=None):
    def addEffects(self, config, img, chip_output, filt, ra_cen, dec_cen, img_rot, exptime=150., pointing_ID=0, timestamp_obs=1621915200, pointing_type='SCI', sky_map=None, tel=None, logger=None):
        # Set random seeds
        # Set random seeds
        SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
        SeedGainNonuni=int(config["random_seeds"]["seed_gainNonUniform"])
        SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
        SeedBiasNonuni=int(config["random_seeds"]["seed_biasNonUniform"])
@@ -325,7 +326,7 @@ class Chip(FocalPlane):
        #     img.addNoise(poisson_noise)
        #     img.addNoise(poisson_noise)


        # Add cosmic-rays
        # Add cosmic-rays
        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='MS':
        if config["ins_effects"]["cosmic_ray"] == True and pointing_type=='SCI':
            chip_utils.log_info(msg="  Adding Cosmic-Ray", logger=self.logger)
            chip_utils.log_info(msg="  Adding Cosmic-Ray", logger=self.logger)
            img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime, 
            img, crmap_gsimg, cr_event_num = chip_utils.add_cosmic_rays(img=img, chip=self, exptime=exptime, 
                                                    seed=SeedCosmicRay+pointing_ID*30+self.chipID)
                                                    seed=SeedCosmicRay+pointing_ID*30+self.chipID)
@@ -339,6 +340,8 @@ class Chip(FocalPlane):
                pointing_ID=pointing_ID,
                pointing_ID=pointing_ID,
                output_dir=chip_output.subdir,
                output_dir=chip_output.subdir,
                exptime=exptime,
                exptime=exptime,
                project_cycle=config["project_cycle"],
                run_counter=config["run_counter"],
                timestamp=timestamp_obs)
                timestamp=timestamp_obs)
            del crmap_gsimg
            del crmap_gsimg


@@ -490,6 +493,8 @@ class Chip(FocalPlane):
                    pointing_ID=pointing_ID,
                    pointing_ID=pointing_ID,
                    output_dir=chip_output.subdir,
                    output_dir=chip_output.subdir,
                    exptime=0.0,
                    exptime=0.0,
                    project_cycle=config["project_cycle"],
                    run_counter=config["run_counter"],
                    timestamp=timestamp_obs)
                    timestamp=timestamp_obs)
            del BiasCombImg
            del BiasCombImg


@@ -590,6 +595,8 @@ class Chip(FocalPlane):
                    pointing_ID=pointing_ID,
                    pointing_ID=pointing_ID,
                    output_dir=chip_output.subdir,
                    output_dir=chip_output.subdir,
                    exptime=self.flat_exptime,
                    exptime=self.flat_exptime,
                    project_cycle=config["project_cycle"],
                    run_counter=config["run_counter"],
                    timestamp=timestamp_obs)
                    timestamp=timestamp_obs)


            del FlatCombImg, FlatSingle, prnu_img
            del FlatCombImg, FlatSingle, prnu_img
@@ -643,6 +650,8 @@ class Chip(FocalPlane):
                        pointing_ID=pointing_ID,
                        pointing_ID=pointing_ID,
                        output_dir=chip_output.subdir,
                        output_dir=chip_output.subdir,
                        exptime=self.dark_exptime,
                        exptime=self.dark_exptime,
                        project_cycle=config["project_cycle"],
                        run_counter=config["run_counter"],
                        timestamp=timestamp_obs)
                        timestamp=timestamp_obs)
                    del crmap_gsimg
                    del crmap_gsimg


@@ -709,6 +718,8 @@ class Chip(FocalPlane):
                    pointing_ID=pointing_ID,
                    pointing_ID=pointing_ID,
                    output_dir=chip_output.subdir,
                    output_dir=chip_output.subdir,
                    exptime=self.dark_exptime,
                    exptime=self.dark_exptime,
                    project_cycle=config["project_cycle"],
                    run_counter=config["run_counter"],
                    timestamp = timestamp_obs)
                    timestamp = timestamp_obs)
            del DarkCombImg
            del DarkCombImg
        # img = galsim.ImageUS(img)
        # img = galsim.ImageUS(img)
Loading