init tools/getPSF.py

import os

import numpy as np

import ObservationSim.PSF.PSFInterp as PSFInterp

from ObservationSim.Instrument import Chip, Filter, FilterParam

import yaml

import galsim

import astropy.io.fits as fitsio

# Setup PATH

SIMPATH = "/share/simudata/CSSOSDataProductsSims/data/CSSTSimImage_C8/testRun_FGS"

config_filename= SIMPATH+"/config_C6_fits.yaml"

cat_filename   = SIMPATH+"/MSC_00000000/MSC_10106100000000_chip_40_filt_FGS.cat"

# Read cat file

catFn = open(cat_filename,"r")

line = catFn.readline()

print(cat_filename,'\n',line)

imgPos = []

chipID = -1

for line in catFn:

    line = line.strip()

    columns = line.split()

    if chipID == -1:

        chipID = int(columns[1])

    else:

        assert chipID == int(columns[1])

    ximg = float(columns[3])

    yimg = float(columns[4])

    imgPos.append([ximg, yimg])

imgPos = np.array(imgPos)

nobj = imgPos.shape[0]

print('chipID, nobj::', chipID, nobj)

# Read config file

with open(config_filename, "r") as stream:

    try:

        config = yaml.safe_load(stream)

        for key, value in config.items():

            print (key + " : " + str(value))

    except yaml.YAMLError as exc:

        print(exc)

# Setup Chip

chip = Chip(chipID=chipID, config=config)

print('chip.bound::', chip.bound.xmin, chip.bound.xmax, chip.bound.ymin, chip.bound.ymax)

for iobj in range(nobj):

    print("\nget psf for iobj-", iobj, '\t', 'bandpass:', end=" ", flush=True)

    # Setup Position on focalplane

    x, y = imgPos[iobj, :] # try get the PSF at some location (1234, 1234) on the chip

    x = x+chip.bound.xmin

    y = y+chip.bound.ymin

    pos_img = galsim.PositionD(x, y)

    # Setup sub-bandpass

    # (There are 4 sub-bandpasses for each PSF sample)

    filter_param = FilterParam()

    filter_id, filter_type = chip.getChipFilter()

    filt = Filter(

        filter_id=filter_id, 

        filter_type=filter_type, 

        filter_param=filter_param, 

        ccd_bandpass=chip.effCurve)

    bandpass_list = filt.bandpass_sub_list

    for i in range(len(bandpass_list)):

        print(i, end=" ", flush=True)

        bandpass = bandpass_list[i] # say you want to access the PSF for the sub-bandpass at the blue end for that chip

        # Get corresponding PSF model

        psf_model = PSFInterp(chip=chip, npsf=100, PSF_data_file=config["psf_setting"]["psf_dir"])

        psf  = psf_model.get_PSF(chip=chip, pos_img=pos_img, bandpass=bandpass, galsimGSObject=False)

        if True:

            fn = "psf_{:}.{:}.{:}.fits".format(chipID, iobj, i)

            if fn != None:

                if os.path.exists(fn):

                    os.remove(fn)

            hdu = fitsio.PrimaryHDU()

            hdu.data = psf

            hdu.header.set('pixScale', 5)

            hdu.writeto(fn)