Skip to content
GitLab
Commit 36189a3e authored by JX's avatar JX 😵
Browse files

Merge remote-tracking branch 'origin/develop'

parents dd26d370 27646bc4
Pipeline #4509 passed with stage
in 0 seconds
Expand all Hide whitespace changes
Inline Side-by-side
ObservationSim/Config/Header/ImageHeader_1.py deleted 100644 → 0
View file @ dd26d370
"""
generate image header
"""
import numpy as np
from astropy.io import fits
import astropy.wcs as pywcs
from scipy import math
import random
import os
import sys
def chara2digit(char):
""" Function to judge and convert characters to digitals
Parameters
----------
"""
try:
float(char) # for int, long and float
except ValueError:
pass
return char
else:
data = float(char)
return data
def read_header_parameter(filename='global_header.param'):
""" Function to read the header parameters
Parameters
----------
"""
name = []
value = []
description = []
for line in open(filename):
line = line.strip("\n")
arr = line.split('|')
# csvReader = csv.reader(csvDataFile)
# for arr in csvReader:
name.append(arr[0])
value.append(chara2digit(arr[1]))
description.append(arr[2])
# print(value)
return name, value, description
def rotate_CD_matrix(cd, pa_aper):
"""Rotate CD matrix
Parameters
----------
cd: (2,2) array
CD matrix
pa_aper: float
Position angle, in degrees E from N, of y axis of the detector
Returns
-------
cd_rot: (2,2) array
Rotated CD matrix
Comments
--------
`astropy.wcs.WCS.rotateCD` doesn't work for non-square pixels in that it
doesn't preserve the pixel scale! The bug seems to come from the fact
that `rotateCD` assumes a transposed version of its own CD matrix.
"""
rad = np.deg2rad(-pa_aper)
mat = np.zeros((2,2))
mat[0,:] = np.array([np.cos(rad),-np.sin(rad)])
mat[1,:] = np.array([np.sin(rad),np.cos(rad)])
cd_rot = np.dot(mat, cd)
return cd_rot
def Header_extention(xlen = 9232, ylen = 9216, gain = 1.0, readout = 5.0, dark = 0.02,saturation=90000, row_num = 1, col_num = 1):
""" Creat an image frame for CCST with multiple extensions
Parameters
----------
"""
flag_ltm_x = [0,1,-1,1,-1]
flag_ltm_y = [0,1,1,-1,-1]
flag_ltv_x = [0,0,1,0,1]
flag_ltv_y = [0,0,0,1,1]
detector_size_x = int(xlen)
detector_size_y = int(ylen)
data_x = str(int(detector_size_x))
data_y = str(int(detector_size_y))
data_sec = '[1:'+data_x+',1:'+data_y+']'
name = []
value = []
description = []
for k in range(1,2):
# f = open("extension"+str(k)+"_image.param","w")
j = row_num
i = col_num
ccdnum = str((j-1)*5+i)
name = ['EXTNAME',
'BSCALE',
'BZERO',
'OBSID',
'CCDNAME',
'AMPNAME',
'GAIN',
'RDNOISE',
'DARK',
'SATURATE',
'RSPEED',
'CHIPTEMP',
'CCDCHIP',
'DATASEC',
'CCDSUM',
'NSUM',
'LTM1_1',
'LTM2_2',
'LTV1',
'LTV2',
'ATM1_1',
'ATM2_2',
'ATV1',
'ATV2',
'DTV1',
'DTV2',
'DTM1_1',
'DTM2_2']
value = ['IM'+str(k),
1.0,
0.0,
'CSST.20200101T000000',
'ccd' + ccdnum.rjust(2,'0'),
'ccd' + ccdnum.rjust(2,'0') + ':'+str(k),
gain,
readout,
dark,
saturation,
10.0,
-100.0,
'ccd' + ccdnum.rjust(2,'0'),
data_sec,
'1 1',
'1 1',
flag_ltm_x[k],
flag_ltm_y[k],
flag_ltv_x[k]*(detector_size_x-20*2+1),
flag_ltv_y[k]*(detector_size_y+1),
flag_ltm_x[k],
flag_ltm_y[k],
flag_ltv_x[k]*(detector_size_x-20*2+1),
flag_ltv_y[k]*(detector_size_y+1),
0,
0,
1,
1]
description = ['Extension name',
' ',
' ',
'Observation ID',
'CCD name',
'Amplifier name',
'Gain (e-/ADU)',
'Readout noise (e-/pixel)',
'Dark noise (e-/pixel/s)',
'Saturation (e-)',
'Read speed',
'Chip temperature',
'CCD chip ID',
'Data section',
'CCD pixel summing',
'CCD pixel summing',
'CCD to image transformation',
'CCD to image transformation',
'CCD to image transformation',
'CCD to image transformation',
'CCD to amplifier transformation',
'CCD to amplifier transformation',
'CCD to amplifier transformation',
'CCD to amplifier transformation',
'CCD to detector transformatio',
'CCD to detector transformatio',
'CCD to detector transformatio',
'CCD to detector transformatio']
return name, value, description
##9232 9216 898 534 1309 60 -40 -23.4333
def WCS_def(xlen = 9232, ylen = 9216, gapx = 898.0, gapy1 = 534, gapy2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1):
""" Creat a wcs frame for CCST with multiple extensions
Parameters
----------
"""
flag_x = [0, 1, -1, 1, -1]
flag_y = [0, 1, 1, -1, -1]
flag_ext_x = [0,-1,1,-1,1]
flag_ext_y = [0,-1,-1,1,1]
x_num = 5
y_num = 6
detector_num = x_num*y_num
detector_size_x = xlen
detector_size_y = ylen
gap_x = gapx
gap_y = [gapy1,gapy2]
ra_ref = ra
dec_ref = dec
pa_aper = pa
pixel_size = psize
gap_y1_num = 3
gap_y2_num = 2
x_center = (detector_size_x*x_num+gap_x*(x_num-1))/2
y_center = (detector_size_y*y_num+gap_y[0]*gap_y1_num+gap_y[1]*gap_y2_num)/2
gap_y_map = np.array([[0,0,0,0,0],[gap_y[0],gap_y[1],gap_y[1],gap_y[1],gap_y[1]],[gap_y[1],gap_y[0],gap_y[0],gap_y[0],gap_y[0]],[gap_y[0],gap_y[0],gap_y[0],gap_y[0],gap_y[0]],[gap_y[0],gap_y[0],gap_y[0],gap_y[0],gap_y[1]],[gap_y[1],gap_y[1],gap_y[1],gap_y[1],gap_y[0]]])
frame_array = np.empty((5,6),dtype=np.float64)
# print(x_center,y_center)
j = row_num
i = col_num
ccdnum = str((j-1)*5+i)
x_ref, y_ref = (detector_size_x+gap_x)*i-gap_x-detector_size_x/2, detector_size_y*j + sum(gap_y_map[0:j,i-1]) - detector_size_y/2
# print(i,j,x_ref,y_ref,ra_ref,dec_ref)
name = []
value = []
description = []
for k in range(1,2):
cd = np.array([[ pixel_size, 0], [0, pixel_size]])/3600.*flag_x[k]
cd_rot = rotate_CD_matrix(cd, pa_aper)
# f = open("CCD"+ccdnum.rjust(2,'0')+"_extension"+str(k)+"_wcs.param","w")
name = ['EQUINOX',
'WCSDIM',
'CTYPE1',
'CTYPE2',
'CRVAL1',
'CRVAL2',
'CRPIX1',
'CRPIX2',
'CD1_1',
'CD1_2',
'CD2_1',
'CD2_2']
value = [2000.0,
2.0,
'RA---TAN',
'DEC--TAN',
ra_ref,
dec_ref,
flag_ext_x[k]*((x_ref+flag_ext_x[k]*detector_size_x/2)-x_center),
flag_ext_y[k]*((y_ref+flag_ext_y[k]*detector_size_y/2)-y_center),
cd_rot[0,0],
cd_rot[0,1],
cd_rot[1,0],
cd_rot[1,1]]
description = ['Equinox of WCS',
'WCS Dimensionality',
'Coordinate type',
'Coordinate typ',
'Coordinate reference value',
'Coordinate reference value',
'Coordinate reference pixel',
'Coordinate reference pixel',
'Coordinate matrix',
'Coordinate matrix',
'Coordinate matrix',
'Coordinate matrix']
return name, value, description
def generatePrimaryHeader(xlen = 9232, ylen = 9216,pointNum = '1', ra = 60, dec = -40, psize = 0.074, row_num = 1, col_num = 1):
# array_size1, array_size2, flux, sigma = int(argv[1]), int(argv[2]), 1000.0, 5.0
filerParm_fn = os.path.split(os.path.realpath(__file__))[0] + '/filter.lst'
f = open(filerParm_fn)
s = f.readline()
s = s.strip("\n")
filter = s.split(' ')
k = (row_num-1)*5+col_num
ccdnum = str(k)
g_header_fn = os.path.split(os.path.realpath(__file__))[0] + '/global_header.param'
name, value, description = read_header_parameter(g_header_fn)
h_prim = fits.Header()
date = '200930'
time_obs = '120000'
for i in range(len(name)):
if(name[i]=='FILTER'):
value[i] = filter[k-1]
if(name[i]=='FILENAME'):
value[i] = 'CSST_' + date + '_' +time_obs + '_' + pointNum.rjust(6,'0') + '_' +ccdnum.rjust(2,'0')+'_raw'
if(name[i]=='DETSIZE'):
value[i] = '[1:' + str(int(xlen)) + ',1:'+ str(int(ylen)) + ']'
if(name[i]=='PIXSCAL1'):
value[i] = str(psize)
if(name[i]=='PIXSCAL2'):
value[i] = str(psize)
h_prim[name[i]] = (value[i],description[i])
h_prim.add_comment('==================================================================',after='FILETYPE')
h_prim.add_comment('Target information')
h_prim.add_comment('==================================================================')
h_prim.add_comment('==================================================================',after='EQUINOX')
h_prim.add_comment('Exposure information')
h_prim.add_comment('==================================================================')
h_prim.add_comment('==================================================================',after='MJDEND')
h_prim.add_comment('Telescope information')
h_prim.add_comment('==================================================================')
h_prim.add_comment('==================================================================',after='REFFRAME')
h_prim.add_comment('Detector information')
h_prim.add_comment('==================================================================')
h_prim.add_comment('==================================================================',after='FILTER')
h_prim.add_comment('Other information')
h_prim.add_comment('==================================================================')
return h_prim
def generateExtensionHeader(xlen = 9232, ylen = 9216,ra = 60, dec = -40, pa = -23.433, gain = 1.0, readout = 5.0, dark = 0.02, saturation=90000, psize = 0.074, row_num = 1, col_num = 1):
h_ext = fits.Header()
for i in range(1,2):
# NAXIS1:Number of pixels per row; NAXIS2:Number of rows
h_ext['NAXIS1'] = xlen
h_ext['NAXIS2'] = ylen
name, value, description = Header_extention(xlen = xlen, ylen = ylen, gain = gain, readout = readout, dark = dark, saturation=saturation, row_num = row_num, col_num = col_num)
for j in range(len(name)):
h_ext[name[j]] = (value[j],description[j])
name, value, description = WCS_def(xlen = xlen, ylen = ylen, gapx = 898.0, gapy1 = 534, gapy2 = 1309, ra = ra, dec = dec, pa = pa ,psize = psize, row_num = row_num, col_num = col_num)
for j in range(len(name)):
h_ext[name[j]] = (value[j],description[j])
h_ext.add_comment('==================================================================',after='OBSID')
h_ext.add_comment('Readout information')
h_ext.add_comment('==================================================================')
h_ext.add_comment('==================================================================',after='CHIPTEMP')
h_ext.add_comment('Chip information')
h_ext.add_comment('==================================================================')
h_ext.add_comment('==================================================================',after='DTM2_2')
h_ext.add_comment('WCS information')
h_ext.add_comment('==================================================================')
return h_ext
def main(argv):
xlen = int(argv[1])
ylen = int(argv[2])
pointingNum = argv[3]
ra = float(argv[4])
dec = float(argv[5])
pSize = float(argv[6])
ccd_row_num = int(argv[7])
ccd_col_num = int(argv[8])
pa_aper = float(argv[9])
gain = float(argv[10])
readout = float(argv[11])
dark = float(argv[12])
fw = float(argv[13])
h_prim = generatePrimaryHeader(xlen = xlen, ylen = ylen,ra = ra, dec = dec, psize = pSize, row_num = ccd_row_num, col_num = ccd_col_num, pointNum = pointingNum)
h_ext = generateExtensionHeader(xlen = xlen, ylen = ylen,ra = ra, dec = dec, pa = pa_aper, gain = gain, readout = readout, dark = dark, saturation=fw, psize = pSize, row_num = ccd_row_num, col_num = ccd_col_num)
hdu1 = fits.PrimaryHDU(header=h_prim)
hdu2 = fits.ImageHDU(np.zeros([ylen,xlen]),header = h_ext)
hdul = fits.HDUList([hdu1,hdu2])
hdul.writeto(h_prim['FILENAME']+'.fits',output_verify='ignore')
# if __name__ == "__main__":
# main(sys.argv)
ObservationSim/Instrument/Chip/Chip.py deleted 100755 → 0
View file @ dd26d370
This diff is collapsed.
ObservationSim/Instrument/Telescope.py deleted 100755 → 0
View file @ dd26d370
import numpy as np
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
class Telescope(object):
def __init__(self, param=None, optEffCurve_path=None):
self.diameter = 2.0 # in unit of meter
if param is not None:
self.diameter = param["diameter"]
self.pupil_area = np.pi * (0.5 * self.diameter)**2
if optEffCurve_path is not None:
self.efficiency = self._get_efficiency(optEffCurve_path)
else:
try:
with pkg_resources.files('ObservationSim.Instrument.data').joinpath('mirror_ccdnote.txt') as optEffCurve_path:
self.efficiency = self._get_efficiency(optEffCurve_path)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data', 'mirror_ccdnote.txt') as optEffCurve_path:
self.efficiency = self._get_efficiency(optEffCurve_path)
def _get_efficiency(self, effCurve_path):
""" Read in the efficiency of optics
for each band
Parameters:
effCurve_path: the path for efficiency file
Returns:
opticsEff: a dictionary of efficiency (a scalar) for each band
"""
f = open(effCurve_path, 'r')
for _ in range(2):
header = f.readline()
iline = 0
opticsEff = {}
for line in f:
line = line.strip()
columns = line.split()
opticsEff[str(columns[0])] = float(columns[2])
f.close()
return opticsEff
\ No newline at end of file
ObservationSim/Instrument/_util.py deleted 100755 → 0
View file @ dd26d370
import numpy as np
import os
import math
from pylab import *
from scipy import interpolate
try:
import importlib.resources as pkg_resources
except ImportError:
# Try backported to PY<37 'importlib_resources'
import importlib_resources as pkg_resources
VC_A = 2.99792458e+18 # speed of light: A/s
VC_M = 2.99792458e+8 # speed of light: m/s
H_PLANK = 6.626196e-27 # Plank constant: erg s
ALL_FILTERS = ["NUV","u", "g", "r", "i","z","y","GU", "GV", "GI", "FGS"]
PHOT_FILTERS = ["NUV", "u", "g", 'r', 'i', 'z', 'y', 'FGS']
SPEC_FILTERS = ["GI", "GV", "GU"]
def rotate_conterclockwise(x0, y0, x, y, angle):
"""
Rotate a point counterclockwise by a given angle around a given origin.
The angle should be given in radians.
"""
angle = np.deg2rad(angle)
qx = x0 + np.cos(angle)*(x - x0) - np.sin(angle) * (y - y0)
qy = y0 + np.sin(angle)*(x - x0) + np.cos(angle) * (y - y0)
return qx, qy
def photonEnergy(lambd):
""" The energy of photon at a given wavelength
Parameter:
lambd: the wavelength in unit of Angstrom
Return:
eph: energy of photon in unit of erg
"""
nu = VC_A / lambd
eph = H_PLANK * nu
return eph
def calculateLimitMag(aperture = 2.0, psf_fwhm = 0.1969,pixelSize = 0.074, pmRation = 0.8, throughputFn = 'i_throughput.txt', readout = 5.0, skyFn= 'sky_emiss_hubble_50_50_A.dat', darknoise = 0.02,exTime = 150, exNum = 1, fw = 90000):
'''
description:
param {*} aperture: unit m, default 2 m
param {*} psf_fwhm: psf fwhm, default 0.1969"
param {*} pixelSize: pixel size, default 0.074"
param {*} pmRation: the ratio of souce flux in the limit mag calculation
param {*} throughputFn: throuput file name
param {*} readout: unit, e-/pixel
param {*} skyFn: sky sed file name, average of hst, 'sky_emiss_hubble_50_50_A.dat'
param {*} darknoise: unit, e-/pixel/s
param {*} exTime: exposure time one time, default 150s
param {*} exNum: exposure number, defautl 1
param {*} fw, full well value( or saturation value),default 90000e-/pixel
return {*} limit mag and saturation mag
'''
try:
with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(throughputFn) as data_file:
throughput_f = np.loadtxt(data_file)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.throughputs', throughputFn) as data_file:
throughput_f = np.loadtxt(data_file)
thr_i = interpolate.interp1d(throughput_f[:,0]/10, throughput_f[:,1]); # wavelength in anstrom
f_s = 200
f_e = 1100
delt_f = 0.5
data_num = int((f_e-f_s)/delt_f+1)
eff = np.zeros([data_num,2])
eff[:,0] = np.arange(f_s,f_e+delt_f,delt_f)
eff[:,1] = thr_i(eff[:,0])
wave = np.arange(f_s,f_e+delt_f,delt_f)
wavey = np.ones(wave.shape[0])
try:
with pkg_resources.files('ObservationSim.Instrument.data.throughputs').joinpath(skyFn) as data_file:
skydata = np.loadtxt(data_file)
except AttributeError:
with pkg_resources.path('ObservationSim.Instrument.data.throughputs', skyFn) as data_file:
skydata = np.loadtxt(data_file)
skydatai = interpolate.interp1d(skydata[:,0]/10, skydata[:,1]*10)
sky_data = np.zeros([data_num,2])
sky_data[:,0] = np.arange(f_s,f_e+delt_f,delt_f)
sky_data[:,1] = skydatai(sky_data[:,0])
flux_sky = trapz((sky_data[:,1])*eff[:,1],sky_data[:,0])
skyPix = flux_sky*pixelSize*pixelSize*pi*(aperture*aperture/4)
###limit mag
r_pix = psf_fwhm*0.7618080243778568/pixelSize # radius RE80, pixel
cnum = math.pi * r_pix * r_pix
sn = 5
d = skyPix*exTime*exNum*cnum + darknoise*exTime*exNum*cnum+readout*readout*cnum*exNum
a=1
b=-sn*sn
c=-sn*sn*d
flux = (-b+sqrt(b*b-4*a*c))/(2*a)/pmRation
limitMag = -2.5*log10(flux/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)))
### saturation mag
from astropy.modeling.models import Gaussian2D
m_size = int(20 * psf_fwhm/pixelSize)
if m_size%2 == 0:
m_size + 1
m_cen = m_size//2
psf_sigma = psf_fwhm/2.355/pixelSize
gaussShape = Gaussian2D(1, m_cen, m_cen, psf_sigma, psf_sigma)
yp, xp = np.mgrid[0:m_size, 0:m_size]
psfMap = gaussShape(xp, yp)
maxRatio = np.amax(psfMap)/np.sum(psfMap)
# print(maxRatio)
flux_sat = fw/maxRatio*exNum
satMag = -2.5*log10(flux_sat/(54799275581.04437 * trapz(wavey*eff[:,1]/wave,wave, 0.1)*exTime*exNum*pi*(aperture/2)*(aperture/2)));
return limitMag , satMag
\ No newline at end of file
ObservationSim/MockObject/CosmicRay.py deleted 100755 → 0
View file @ dd26d370
from ObservationSim.MockObject.MockObject import MockObject
class CosmicRay(MockObject):
pass
\ No newline at end of file
ObservationSim/PSF/PSFGauss.py deleted 100644 → 0
View file @ dd26d370
import galsim
import sep
import numpy as np
from scipy.interpolate import interp1d
from ObservationSim.PSF.PSFModel import PSFModel
class PSFGauss(PSFModel):
def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=None):
self.pix_size = chip.pix_scale
self.chip = chip
if psfRa is None:
self.fwhm = fwhm
self.sigGauss = 0.15
else:
self.fwhm = self.fwhmGauss(r=psfRa)
self.sigGauss = psfRa # 80% light radius
self.sigSpin = sigSpin
self.psf = galsim.Gaussian(flux=1.0,fwhm=fwhm)
def perfGauss(self, r, sig):
"""
pseudo-error function, i.e. Cumulative distribution function of Gaussian distribution
Parameter:
r: radius
sig: sigma of the Gaussian distribution
Return:
the value of the pseudo CDF
"""
gaussFun = lambda sigma, r: 1.0/(np.sqrt(2.0*np.pi)*sigma) * np.exp(-r**2/(2.0*sigma**2))
nxx = 1000
rArr = np.linspace(0.0,r,nxx)
gauss = gaussFun(sig,rArr)
erf = 2.0*np.trapz(gauss,rArr)
return erf
def fracGauss(self, sig, r=0.15, pscale=None):
"""
For a given Gaussian PSF with sigma=sig,
derive the flux ratio ar the given radius r
Parameters:
sig: sigma of the Gauss PSF Function in arcsec
r: radius in arcsec
pscale: pixel scale
Return: the flux ratio
"""
if pscale == None:
pscale = self.pix_size
gaussx = galsim.Gaussian(flux=1.0,sigma=sig)
gaussImg = gaussx.drawImage(scale=pscale, method='no_pixel')
gaussImg = gaussImg.array
size = np.size(gaussImg,axis=0)
cxy = 0.5*(size-1)
flux, ferr, flag = sep.sum_circle(gaussImg,[cxy],[cxy],[r/pscale],subpix=0)
return flux
def fwhmGauss(self, r=0.15,fr=0.8,pscale=None):
"""
Given a total flux ratio 'fr' within a fixed radius 'r',
estimate the fwhm of the Gaussian function
return the fwhm in arcsec
"""
if pscale == None:
pscale = self.pix_size
err = 1.0e-3
nxx = 100
sig = np.linspace(0.5*pscale,1.0,nxx)
frA = np.zeros(nxx)
for i in range(nxx): frA[i] = self.fracGauss(sig[i],r=r,pscale=pscale)
index = [i for i in range(nxx-1) if (fr-frA[i])*(fr-frA[i+1])<=0.0][0]
while abs(frA[index]-fr)>1.0e-3:
sig = np.linspace(sig[index],sig[index+1],nxx)
for i in range(nxx): frA[i] = self.fracGauss(sig[i],r=r,pscale=pscale)
index = [i for i in range(nxx-1) if (fr-frA[i])*(fr-frA[i+1])<=0.0][0]
fwhm = 2.35482*sig[index]
return fwhm
def get_PSF(self, pos_img, chip=None, bandpass=None, folding_threshold=5.e-3):
dx = pos_img.x - self.chip.cen_pix_x
dy = pos_img.y - self.chip.cen_pix_y
return self.PSFspin(dx, dy)
def PSFspin(self, x, y):
"""
The PSF profile at a given image position relative to the axis center
Parameters:
theta : spin angles in a given exposure in unit of [arcsecond]
dx, dy: relative position to the axis center in unit of [pixels]
Return:
Spinned PSF: g1, g2 and axis ratio 'a/b'
"""
a2Rad = np.pi/(60.0*60.0*180.0)
ff = self.sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
rc = np.sqrt(x*x + y*y)
cpix = rc*(self.sigSpin*a2Rad)
beta = (np.arctan2(y,x) + np.pi/2)
ell = cpix**2/(2.0*ff**2+cpix**2)
#ell *= 10.0
qr = np.sqrt((1.0+ell)/(1.0-ell))
#psfShape = galsim.Shear(e=ell, beta=beta)
#g1, g2 = psfShape.g1, psfShape.g2
#qr = np.sqrt((1.0+ell)/(1.0-ell))
#return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return self.psf.shear(PSFshear), PSFshear
\ No newline at end of file
ObservationSim/PSF/PSFModel.py deleted 100755 → 0
View file @ dd26d370
import galsim
import sep
import numpy as np
from scipy.interpolate import interp1d
import pylab as pl
import os, sys
class PSFModel(object):
def __init__(self, sigSpin=0., psfRa=0.15):
# TODO: what are the nesseary fields in PSFModel class?
pass
def PSFspin(self, psf, sigSpin, sigGauss, dx, dy):
"""
The PSF profile at a given image position relative to the axis center
Parameters:
theta : spin angles in a given exposure in unit of [arcsecond]
dx, dy: relative position to the axis center in unit of [pixels]
Return:
Spinned PSF: g1, g2 and axis ratio 'a/b'
"""
a2Rad = np.pi/(60.0*60.0*180.0)
ff = sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
rc = np.sqrt(dx*dx + dy*dy)
cpix = rc*(sigSpin*a2Rad)
beta = (np.arctan2(dy,dx) + np.pi/2)
ell = cpix**2/(2.0*ff**2+cpix**2)
#ell *= 10.0
qr = np.sqrt((1.0+ell)/(1.0-ell))
#psfShape = galsim.Shear(e=ell, beta=beta)
#g1, g2 = psfShape.g1, psfShape.g2
#qr = np.sqrt((1.0+ell)/(1.0-ell))
#return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return psf.shear(PSFshear), PSFshear
\ No newline at end of file
Catalog/Catalog_example.py → catalog/Catalog_example.py
View file @ 36189a3e
......@@ -4,15 +4,16 @@ import astropy.constants as cons
from astropy.table import Table
from scipy import interpolate
from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar
from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar
class Catalog(CatalogBase):
"""An user customizable class for reading in catalog(s) of objects and SEDs.
NOTE: must inherit the "CatalogBase" abstract class
...
Attributes
----------
cat_dir : str
......@@ -24,7 +25,7 @@ class Catalog(CatalogBase):
objs : list
a list of ObservationSim.MockObject (Star, Galaxy, or Quasar)
NOTE: must have "obj" list when implement your own Catalog
Methods
----------
load_sed(obj, **kwargs):
......@@ -32,9 +33,10 @@ class Catalog(CatalogBase):
load_norm_filt(obj):
load the filter throughput for the input catalog's photometric system.
"""
def __init__(self, config, chip, **kwargs):
"""Constructor method.
Parameters
----------
config : dict
......@@ -44,20 +46,22 @@ class Catalog(CatalogBase):
**kwargs : dict
other needed input parameters (in key-value pairs), please modify corresponding
initialization call in "ObservationSim.py" as you need.
Returns
----------
None
"""
super().__init__()
self.cat_dir = os.path.join(config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.cat_dir = os.path.join(
config["data_dir"], config["catalog_options"]["input_path"]["cat_dir"])
self.chip = chip
if "star_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["star_cat"]:
star_file = config["catalog_options"]["input_path"]["star_cat"]
star_SED_file = config["catalog_options"]["SED_templates_path"]["star_SED"]
self.star_path = os.path.join(self.cat_dir, star_file)
self.star_SED_path = os.path.join(config["data_dir"], star_SED_file)
self.star_SED_path = os.path.join(
config["data_dir"], star_SED_file)
# NOTE: must call _load() method here to read in all objects
self.objs = []
self._load()
......@@ -67,7 +71,7 @@ class Catalog(CatalogBase):
This is a must implemented method which is used to read in all objects, and
then convert them to ObservationSim.MockObject (Star, Galaxy, or Quasar).
Currently,
the model of ObservationSim.MockObject.Star class requires:
param["star"] : int
......@@ -83,7 +87,7 @@ class Catalog(CatalogBase):
NOTE: if that filter is not the corresponding CSST filter, the
load_norm_filt(obj) function must be implemented to load the filter
throughput of that particular photometric system
the model of ObservationSim.MockObject.Galaxy class requires:
param["star"] : int
specify the object type: 0: galaxy, 1: star, 2: quasar
......@@ -173,7 +177,8 @@ class Catalog(CatalogBase):
"""
if obj.type == 'star':
wave = Table.read(self.star_SED_path, path=f"/SED/wave_{obj.model_tag}")
wave = Table.read(self.star_SED_path,
path=f"/SED/wave_{obj.model_tag}")
flux = Table.read(self.star_SED_path, path=f"/SED/{obj.sed_type}")
wave, flux = wave['col0'].data, flux['col0'].data
else:
......@@ -203,4 +208,4 @@ class Catalog(CatalogBase):
norm_filt["WAVELENGTH"] : wavelengthes in Angstroms
norm_filt["SENSITIVITY"] : efficiencies
"""
return None
\ No newline at end of file
return None
Catalog/testCat_fits.py → catalog/testCat_fits.py
View file @ 36189a3e
......@@ -11,12 +11,12 @@ from astropy.table import Table
from scipy import interpolate
from datetime import datetime
from ObservationSim.MockObject import CatalogBase, Star, Galaxy, Quasar, Stamp
from ObservationSim.MockObject._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
from observation_sim.mock_objects import CatalogBase, Star, Galaxy, Quasar, Stamp
from observation_sim.mock_objects._util import tag_sed, getObservedSED, getABMAG, integrate_sed_bandpass, comoving_dist
from observation_sim.astrometry.Astrometry_util import on_orbit_obs_position
import astropy.io.fits as fitsio
from ObservationSim.MockObject._util import seds, sed_assign, extAv
from observation_sim.mock_objects._util import seds, sed_assign, extAv
# (TEST)
from astropy.cosmology import FlatLambdaCDM
......@@ -31,11 +31,12 @@ except ImportError:
NSIDE = 128
class Catalog(CatalogBase):
def __init__(self, config, chip, pointing, chip_output, filt, **kwargs):
super().__init__()
self.cat_dir = config["catalog_options"]["input_path"]["cat_dir"]
self.seed_Av = 121212 #config["catalog_options"]["seed_Av"]
self.seed_Av = 121212 # config["catalog_options"]["seed_Av"]
# (TEST)
self.cosmo = FlatLambdaCDM(H0=67.66, Om0=0.3111)
......@@ -44,11 +45,11 @@ class Catalog(CatalogBase):
self.filt = filt
self.logger = chip_output.logger
with pkg_resources.path('Catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
with pkg_resources.path('catalog.data', 'SLOAN_SDSS.g.fits') as filter_path:
self.normF_star = Table.read(str(filter_path))
with pkg_resources.path('Catalog.data', 'lsst_throuput_g.fits') as filter_path:
with pkg_resources.path('catalog.data', 'lsst_throuput_g.fits') as filter_path:
self.normF_galaxy = Table.read(str(filter_path))
self.config = config
self.chip = chip
self.pointing = pointing
......@@ -58,24 +59,27 @@ class Catalog(CatalogBase):
if "stamp_cat" in config["catalog_options"]["input_path"] and config["catalog_options"]["input_path"]["stamp_cat"] and config["catalog_options"]["stamp_yes"]:
stamp_file = config["catalog_options"]["input_path"]["stamp_cat"]
self.stamp_path = os.path.join(self.cat_dir, stamp_file)
#self.stamp_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["stamp_SED"]) ###shoule be stamp-SED
#self._load_SED_lib_stamps() ###shoule be stamp-SED
self.tempSed_gal, self.tempRed_gal = seds("galaxy.list", seddir="/public/home/chengliang/CSSOSDataProductsSims/testCats/Templates/Galaxy/") #only for test
# self.stamp_SED_path = os.path.join(config["data_dir"], config["SED_templates_path"]["stamp_SED"]) ###shoule be stamp-SED
# self._load_SED_lib_stamps() ###shoule be stamp-SED
self.tempSed_gal, self.tempRed_gal = seds(
"galaxy.list", seddir="/public/home/chengliang/CSSOSDataProductsSims/testCats/Templates/Galaxy/") # only for test
self._add_output_columns_header()
self._get_healpix_list()
self._load()
def _add_output_columns_header(self):
self.add_hdr = " model_tag teff logg feh"
self.add_hdr += " bulgemass diskmass detA e1 e2 kappa g1 g2 size galType veldisp "
self.add_fmt = " %10s %8.4f %8.4f %8.4f"
self.add_fmt += " %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f %4d %8.4f "
self.chip_output.update_output_header(additional_column_names=self.add_hdr)
self.chip_output.update_output_header(
additional_column_names=self.add_hdr)
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]))
......@@ -94,7 +98,7 @@ class Catalog(CatalogBase):
def load_norm_filt(self, obj):
if obj.type == "stamp":
return self.normF_galaxy ###normalize_filter for stamp
return self.normF_galaxy # normalize_filter for stamp
else:
return None
......@@ -102,35 +106,38 @@ class Catalog(CatalogBase):
print("debug:: load_stamps")
nstamps = len(stamps['filename'])
self.rng_sedGal = random.Random()
self.rng_sedGal.seed(float(pix_id)) # Use healpix index as the random seed
# Use healpix index as the random seed
self.rng_sedGal.seed(float(pix_id))
self.ud = galsim.UniformDeviate(pix_id)
for istamp in range(nstamps):
print("debug::", istamp)
fitsfile = os.path.join(self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8'))
fitsfile = os.path.join(
self.cat_dir, "stampCats/"+stamps['filename'][istamp].decode('utf-8'))
print("debug::", istamp, fitsfile)
hdu=fitsio.open(fitsfile)
hdu = fitsio.open(fitsfile)
param = self.initialize_param()
param['id'] = hdu[0].header['index'] #istamp
param['id'] = hdu[0].header['index'] # istamp
param['star'] = 3 # Stamp type in .cat file
param['ra'] = hdu[0].header['ra']
param['dec']= hdu[0].header['dec']
param['pixScale']= hdu[0].header['pixScale']
#param['srcGalaxyID'] = hdu[0].header['srcGID']
#param['mu']= hdu[0].header['mu']
#param['PA']= hdu[0].header['PA']
#param['bfrac']= hdu[0].header['bfrac']
#param['z']= hdu[0].header['z']
param['mag_use_normal'] = hdu[0].header['mag_g'] #gals['mag_true_g_lsst']
param['dec'] = hdu[0].header['dec']
param['pixScale'] = hdu[0].header['pixScale']
# param['srcGalaxyID'] = hdu[0].header['srcGID']
# param['mu']= hdu[0].header['mu']
# param['PA']= hdu[0].header['PA']
# param['bfrac']= hdu[0].header['bfrac']
# param['z']= hdu[0].header['z']
# gals['mag_true_g_lsst']
param['mag_use_normal'] = hdu[0].header['mag_g']
# Apply astrometric modeling
# in C3 case only aberration
param['ra_orig'] = param['ra']
param['dec_orig']= param['dec']
param['dec_orig'] = param['dec']
if self.config["obs_setting"]["enable_astrometric_model"]:
ra_list = [param['ra']] #ra_arr.tolist()
dec_list= [param['dec']] #dec_arr.tolist()
ra_list = [param['ra']] # ra_arr.tolist()
dec_list = [param['dec']] # dec_arr.tolist()
pmra_list = np.zeros(1).tolist()
pmdec_list = np.zeros(1).tolist()
rv_list = np.zeros(1).tolist()
......@@ -157,19 +164,20 @@ class Catalog(CatalogBase):
input_time_str=time_str
)
param['ra'] = ra_arr[0]
param['dec']= dec_arr[0]
param['dec'] = dec_arr[0]
# Assign each galaxy a template SED
param['sed_type'] = sed_assign(phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
param['sed_type'] = sed_assign(
phz=param['z'], btt=param['bfrac'], rng=self.rng_sedGal)
param['redden'] = self.tempRed_gal[param['sed_type']]
param['av'] = 0.0
param['redden'] = 0
param['mu'] = 1
#param["CSSTmag"]= True
#param["mag_r"] = 20.
#param['']
###more keywords for stamp###
# param["CSSTmag"]= True
# param["mag_r"] = 20.
# param['']
### more keywords for stamp###
param['image'] = hdu[0].data
param['image'] = param['image']/(np.sum(param['image']))
obj = Stamp(param)
......@@ -181,12 +189,12 @@ class Catalog(CatalogBase):
if "stamp_cat" in self.config["catalog_options"]["input_path"] and self.config["catalog_options"]["input_path"]["stamp_cat"] and self.config["catalog_options"]["stamp_yes"]:
stamps_cat = h5.File(self.stamp_path, 'r')['Stamps']
print("debug::",stamps_cat.keys())
print("debug::", stamps_cat.keys())
for pix in self.pix_list:
try:
stamps = stamps_cat[str(pix)]
print("debug::",stamps.keys())
print("debug::", stamps.keys())
self._load_stamps(stamps, pix_id=pix)
del stamps
except Exception as e:
......@@ -194,12 +202,12 @@ class Catalog(CatalogBase):
print(e)
if self.logger is not None:
self.logger.info("maximum galaxy size: %.4f"%(self.max_size))
self.logger.info("number of objects in catalog: %d"%(len(self.objs)))
self.logger.info("maximum galaxy size: %.4f" % (self.max_size))
self.logger.info("number of objects in catalog: %d" %
(len(self.objs)))
else:
print("number of objects in catalog: ", len(self.objs))
def load_sed(self, obj, **kwargs):
if obj.type == 'stamp':
sed_data = getObservedSED(
......@@ -217,7 +225,7 @@ class Catalog(CatalogBase):
# erg/s/cm2/A --> photon/s/m2/A
all_sed = y * lamb / (cons.h.value * cons.c.value) * 1e-13
sed = Table(np.array([lamb, all_sed]).T, names=('WAVELENGTH', 'FLUX'))
del wave
del flux
return sed
ObservationSim/ObservationSim.py → observation_sim/ObservationSim.py
View file @ 36189a3e
......@@ -8,24 +8,26 @@ from datetime import datetime
import traceback
from ObservationSim.Config import ChipOutput
from ObservationSim.Instrument import Telescope, Filter, FilterParam, FocalPlane, Chip
from ObservationSim.Instrument.Chip import Effects
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from ObservationSim.Astrometry.Astrometry_util import on_orbit_obs_position
from ObservationSim.sim_steps import SimSteps, SIM_STEP_TYPES
from observation_sim.config import ChipOutput
from observation_sim.instruments import Telescope, Filter, FilterParam, FocalPlane, Chip
from observation_sim.instruments.chip import effects
from observation_sim.instruments.chip import chip_utils as chip_utils
from observation_sim.astrometry.Astrometry_util import on_orbit_obs_position
from observation_sim.sim_steps import SimSteps, SIM_STEP_TYPES
class Observation(object):
def __init__(self, config, Catalog, work_dir=None, data_dir=None):
self.config = config
self.tel = Telescope()
self.filter_param = FilterParam()
self.filter_param = FilterParam()
self.Catalog = Catalog
def prepare_chip_for_exposure(self, chip, ra_cen, dec_cen, pointing, wcs_fp=None):
# Get WCS for the focal plane
if wcs_fp == None:
wcs_fp = self.focal_plane.getTanWCS(ra_cen, dec_cen, pointing.img_pa, chip.pix_scale)
wcs_fp = self.focal_plane.getTanWCS(
ra_cen, dec_cen, pointing.img_pa, chip.pix_scale)
# Create chip Image
chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
......@@ -35,29 +37,36 @@ class Observation(object):
# Get random generators for this chip
chip.rng_poisson, chip.poisson_noise = chip_utils.get_poisson(
seed=int(self.config["random_seeds"]["seed_poisson"]) + pointing.id*30 + chip.chipID, sky_level=0.)
# Get flat, shutter, and PRNU images
chip.flat_img, _ = chip_utils.get_flat(img=chip.img, seed=int(self.config["random_seeds"]["seed_flat"]))
chip.flat_img, _ = chip_utils.get_flat(
img=chip.img, seed=int(self.config["random_seeds"]["seed_flat"]))
if chip.chipID > 30:
chip.shutter_img = np.ones_like(chip.img.array)
else:
chip.shutter_img = Effects.ShutterEffectArr(chip.img, t_shutter=1.3, dist_bearing=735, dt=1E-3)
chip.prnu_img = Effects.PRNU_Img(xsize=chip.npix_x, ysize=chip.npix_y, sigma=0.01,
chip.shutter_img = effects.ShutterEffectArr(
chip.img, t_shutter=1.3, dist_bearing=735, dt=1E-3)
chip.prnu_img = effects.PRNU_Img(xsize=chip.npix_x, ysize=chip.npix_y, sigma=0.01,
seed=int(self.config["random_seeds"]["seed_prnu"]+chip.chipID))
return chip
def run_one_chip(self, chip, filt, pointing, chip_output, wcs_fp=None, psf_model=None, cat_dir=None, sed_dir=None):
chip_output.Log_info(':::::::::::::::::::Current Pointing Information::::::::::::::::::')
chip_output.Log_info(
':::::::::::::::::::Current Pointing Information::::::::::::::::::')
chip_output.Log_info("RA: %f, DEC; %f" % (pointing.ra, pointing.dec))
chip_output.Log_info("Time: %s" % datetime.utcfromtimestamp(pointing.timestamp).isoformat())
chip_output.Log_info("Time: %s" % datetime.utcfromtimestamp(
pointing.timestamp).isoformat())
chip_output.Log_info("Exposure time: %f" % pointing.exp_time)
chip_output.Log_info("Satellite Position (x, y, z): (%f, %f, %f)" % (pointing.sat_x, pointing.sat_y, pointing.sat_z))
chip_output.Log_info("Satellite Velocity (x, y, z): (%f, %f, %f)" % (pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
chip_output.Log_info("Satellite Position (x, y, z): (%f, %f, %f)" % (
pointing.sat_x, pointing.sat_y, pointing.sat_z))
chip_output.Log_info("Satellite Velocity (x, y, z): (%f, %f, %f)" % (
pointing.sat_vx, pointing.sat_vy, pointing.sat_vz))
chip_output.Log_info("Position Angle: %f" % pointing.img_pa.deg)
chip_output.Log_info('Chip : %d' % chip.chipID)
chip_output.Log_info(':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
chip_output.Log_info(
':::::::::::::::::::::::::::END:::::::::::::::::::::::::::::::::::')
# Apply astrometric simulation for pointing
if self.config["obs_setting"]["enable_astrometric_model"]:
......@@ -91,32 +100,35 @@ class Observation(object):
chip = self.prepare_chip_for_exposure(chip, ra_cen, dec_cen, pointing)
# Initialize SimSteps
sim_steps = SimSteps(overall_config=self.config, chip_output=chip_output, all_filters=self.all_filters)
sim_steps = SimSteps(overall_config=self.config,
chip_output=chip_output, all_filters=self.all_filters)
for step in pointing.obs_param["call_sequence"]:
if self.config["run_option"]["out_cat_only"]:
if step != "scie_obs":
continue
chip_output.Log_info("Starting simulation step: %s, calling function: %s"%(step, SIM_STEP_TYPES[step]))
chip_output.Log_info("Starting simulation step: %s, calling function: %s" % (
step, SIM_STEP_TYPES[step]))
obs_param = pointing.obs_param["call_sequence"][step]
step_name = SIM_STEP_TYPES[step]
try:
step_func = getattr(sim_steps, step_name)
chip, filt, tel, pointing = step_func(
chip=chip,
filt=filt,
chip=chip,
filt=filt,
tel=self.tel,
pointing=pointing,
catalog=self.Catalog,
pointing=pointing,
catalog=self.Catalog,
obs_param=obs_param)
chip_output.Log_info("Finished simulation step: %s"%(step))
chip_output.Log_info("Finished simulation step: %s" % (step))
except Exception as e:
traceback.print_exc()
chip_output.Log_error(e)
chip_output.Log_error("Failed simulation on step: %s"%(step))
chip_output.Log_error("Failed simulation on step: %s" % (step))
break
chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB"%(pointing.id, chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024) ))
chip_output.Log_info("check running:1: pointing-%d chip-%d pid-%d memory-%6.2fGB" % (pointing.id,
chip.chipID, os.getpid(), (psutil.Process(os.getpid()).memory_info().rss / 1024 / 1024 / 1024)))
del chip.img
def runExposure_MPI_PointingList(self, pointing_list, chips=None):
......@@ -131,9 +143,11 @@ class Observation(object):
# pointing_ID = pointing.id
pointing_ID = pointing.obs_id
pointing.make_output_pointing_dir(overall_config=self.config, copy_obs_config=True)
pointing.make_output_pointing_dir(
overall_config=self.config, copy_obs_config=True)
self.focal_plane = FocalPlane(chip_list=pointing.obs_param["run_chips"])
self.focal_plane = FocalPlane(
chip_list=pointing.obs_param["run_chips"])
# Make Chip & Filter lists
self.chip_list = []
self.filter_list = []
......@@ -167,7 +181,7 @@ class Observation(object):
run_chips.append(chip)
run_filts.append(filt)
nchips_per_fp = len(chips)
for ichip in range(nchips_per_fp):
i_process = process_counter + ichip
if i_process % num_thread != ind_thread:
......@@ -176,20 +190,22 @@ class Observation(object):
chip = run_chips[ichip]
filt = run_filts[ichip]
chip_output = ChipOutput(
config = self.config,
chip = chip,
filt = filt,
pointing = pointing
config=self.config,
chip=chip,
filt=filt,
pointing=pointing
)
chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..."%(int(pointing_ID), chip.chipID, pid))
chip_output.Log_info("running pointing#%d, chip#%d, at PID#%d..." % (
int(pointing_ID), chip.chipID, pid))
self.run_one_chip(
chip=chip,
filt=filt,
chip_output=chip_output,
pointing=pointing)
chip_output.Log_info("finished running chip#%d..."%(chip.chipID))
chip_output.Log_info(
"finished running chip#%d..." % (chip.chipID))
for handler in chip_output.logger.handlers[:]:
chip_output.logger.removeHandler(handler)
gc.collect()
......
ObservationSim/PSF/FieldDistortion.py → observation_sim/PSF/FieldDistortion.py
View file @ 36189a3e
......@@ -2,6 +2,7 @@ import galsim
import numpy as np
import cmath
class FieldDistortion(object):
def __init__(self, chip, fdModel=None, fdModel_path=None, img_rot=0.):
......@@ -13,7 +14,8 @@ class FieldDistortion(object):
with open(fdModel_path, "rb") as f:
self.fdModel = pickle.load(f)
else:
raise ValueError("Error: no field distortion model has been specified!")
raise ValueError(
"Error: no field distortion model has been specified!")
else:
self.fdModel = fdModel
self.img_rot = img_rot
......@@ -21,19 +23,20 @@ class FieldDistortion(object):
self.ixfdModel = self.ifdModel["xImagePos"]
self.iyfdModel = self.ifdModel["yImagePos"]
# first-order derivatives of the global field distortion model
self.ifx_dx = self.ixfdModel.partial_derivative(1,0)
self.ifx_dy = self.ixfdModel.partial_derivative(0,1)
self.ify_dx = self.iyfdModel.partial_derivative(1,0)
self.ify_dy = self.iyfdModel.partial_derivative(0,1)
self.ifx_dx = self.ixfdModel.partial_derivative(1, 0)
self.ifx_dy = self.ixfdModel.partial_derivative(0, 1)
self.ify_dx = self.iyfdModel.partial_derivative(1, 0)
self.ify_dy = self.iyfdModel.partial_derivative(0, 1)
if "residual" in self.fdModel["wave1"]:
self.irsModel = self.fdModel["wave1"]["residual"]["ccd" + chip.getChipLabel(chipID=chip.chipID)]
self.irsModel = self.fdModel["wave1"]["residual"]["ccd" +
chip.getChipLabel(chipID=chip.chipID)]
self.ixrsModel = self.irsModel["xResidual"]
self.iyrsModel = self.irsModel["yResidual"]
# first-order derivatives of the residual field distortion model
self.irx_dx = self.ixrsModel.partial_derivative(1,0)
self.irx_dy = self.ixrsModel.partial_derivative(0,1)
self.iry_dx = self.iyrsModel.partial_derivative(1,0)
self.iry_dy = self.iyrsModel.partial_derivative(0,1)
self.irx_dx = self.ixrsModel.partial_derivative(1, 0)
self.irx_dy = self.ixrsModel.partial_derivative(0, 1)
self.iry_dx = self.iyrsModel.partial_derivative(1, 0)
self.iry_dy = self.iyrsModel.partial_derivative(0, 1)
else:
self.irsModel = None
......@@ -95,8 +98,8 @@ class FieldDistortion(object):
ix_dy = self.ifx_dy(x, y)
iy_dx = self.ify_dx(x, y)
iy_dy = self.ify_dy(x, y)
g1k_fd = 0.0 + (iy_dy - ix_dx) / (iy_dy + ix_dx)
g2k_fd = 0.0 - (iy_dx + ix_dy) / (iy_dy + ix_dx)
g1k_fd = 0.0 + (iy_dy - ix_dx) / (iy_dy + ix_dx)
g2k_fd = 0.0 - (iy_dx + ix_dy) / (iy_dy + ix_dx)
# [TODO] [TESTING] Rotate the shear:
g_abs = np.sqrt(g1k_fd**2 + g2k_fd**2)
......@@ -107,4 +110,3 @@ class FieldDistortion(object):
fd_shear = galsim.Shear(g1=g1k_fd, g2=g2k_fd)
return galsim.PositionD(x, y), fd_shear
observation_sim/PSF/PSFGauss.py 0 → 100644
View file @ 36189a3e
import galsim
import sep
import numpy as np
from scipy.interpolate import interp1d
from observation_sim.PSF.PSFModel import PSFModel
class PSFGauss(PSFModel):
def __init__(self, chip, fwhm=0.187, sigSpin=0., psfRa=None):
self.pix_size = chip.pix_scale
self.chip = chip
if psfRa is None:
self.fwhm = fwhm
self.sigGauss = 0.15
else:
self.fwhm = self.fwhmGauss(r=psfRa)
self.sigGauss = psfRa # 80% light radius
self.sigSpin = sigSpin
self.psf = galsim.Gaussian(flux=1.0, fwhm=fwhm)
def perfGauss(self, r, sig):
"""
pseudo-error function, i.e. Cumulative distribution function of Gaussian distribution
Parameter:
r: radius
sig: sigma of the Gaussian distribution
Return:
the value of the pseudo CDF
"""
def gaussFun(sigma, r): return 1.0/(np.sqrt(2.0*np.pi)
* sigma) * np.exp(-r**2/(2.0*sigma**2))
nxx = 1000
rArr = np.linspace(0.0, r, nxx)
gauss = gaussFun(sig, rArr)
erf = 2.0*np.trapz(gauss, rArr)
return erf
def fracGauss(self, sig, r=0.15, pscale=None):
"""
For a given Gaussian PSF with sigma=sig,
derive the flux ratio ar the given radius r
Parameters:
sig: sigma of the Gauss PSF Function in arcsec
r: radius in arcsec
pscale: pixel scale
Return: the flux ratio
"""
if pscale == None:
pscale = self.pix_size
gaussx = galsim.Gaussian(flux=1.0, sigma=sig)
gaussImg = gaussx.drawImage(scale=pscale, method='no_pixel')
gaussImg = gaussImg.array
size = np.size(gaussImg, axis=0)
cxy = 0.5*(size-1)
flux, ferr, flag = sep.sum_circle(
gaussImg, [cxy], [cxy], [r/pscale], subpix=0)
return flux
def fwhmGauss(self, r=0.15, fr=0.8, pscale=None):
"""
Given a total flux ratio 'fr' within a fixed radius 'r',
estimate the fwhm of the Gaussian function
return the fwhm in arcsec
"""
if pscale == None:
pscale = self.pix_size
err = 1.0e-3
nxx = 100
sig = np.linspace(0.5*pscale, 1.0, nxx)
frA = np.zeros(nxx)
for i in range(nxx):
frA[i] = self.fracGauss(sig[i], r=r, pscale=pscale)
index = [i for i in range(nxx-1) if (fr-frA[i])
* (fr-frA[i+1]) <= 0.0][0]
while abs(frA[index]-fr) > 1.0e-3:
sig = np.linspace(sig[index], sig[index+1], nxx)
for i in range(nxx):
frA[i] = self.fracGauss(sig[i], r=r, pscale=pscale)
index = [i for i in range(
nxx-1) if (fr-frA[i])*(fr-frA[i+1]) <= 0.0][0]
fwhm = 2.35482*sig[index]
return fwhm
def get_PSF(self, pos_img, chip=None, bandpass=None, folding_threshold=5.e-3):
dx = pos_img.x - self.chip.cen_pix_x
dy = pos_img.y - self.chip.cen_pix_y
return self.PSFspin(dx, dy)
def PSFspin(self, x, y):
"""
The PSF profile at a given image position relative to the axis center
Parameters:
theta : spin angles in a given exposure in unit of [arcsecond]
dx, dy: relative position to the axis center in unit of [pixels]
Return:
Spinned PSF: g1, g2 and axis ratio 'a/b'
"""
a2Rad = np.pi/(60.0*60.0*180.0)
ff = self.sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
rc = np.sqrt(x*x + y*y)
cpix = rc*(self.sigSpin*a2Rad)
beta = (np.arctan2(y, x) + np.pi/2)
ell = cpix**2/(2.0*ff**2+cpix**2)
# ell *= 10.0
qr = np.sqrt((1.0+ell)/(1.0-ell))
# psfShape = galsim.Shear(e=ell, beta=beta)
# g1, g2 = psfShape.g1, psfShape.g2
# qr = np.sqrt((1.0+ell)/(1.0-ell))
# return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return self.psf.shear(PSFshear), PSFshear
ObservationSim/PSF/PSFInterp.py → observation_sim/PSF/PSFInterp.py
View file @ 36189a3e
......@@ -12,14 +12,14 @@ import scipy.spatial as spatial
import galsim
import h5py
from ObservationSim.PSF.PSFModel import PSFModel
from observation_sim.PSF.PSFModel import PSFModel
NPSF = 900 #***# 30*30
PixSizeInMicrons = 5. #***# in microns
NPSF = 900 # ***# 30*30
PixSizeInMicrons = 5. # ***# in microns
###find neighbors-KDtree###
### find neighbors-KDtree###
def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
"""
find nearest neighbors by 2D-KDTree
......@@ -38,7 +38,7 @@ def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
datay = py
tree = spatial.KDTree(list(zip(datax.ravel(), datay.ravel())))
dataq=[]
dataq = []
rr = dr
if OnlyDistance == True:
dataq = tree.query_ball_point([tx, ty], rr)
......@@ -51,7 +51,9 @@ def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
dataq = np.array(dataq)[ddSortindx[0:dn]]
return dataq
###find neighbors-hoclist###
### find neighbors-hoclist###
def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy):
if np.max(partx) > nhocx*dhocx:
print('ERROR')
......@@ -60,8 +62,8 @@ def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy):
print('ERROR')
sys.exit()
npart = partx.size
hoclist= np.zeros(npart, dtype=np.int32)-1
npart = partx.size
hoclist = np.zeros(npart, dtype=np.int32)-1
hoc = np.zeros([nhocy, nhocx], dtype=np.int32)-1
for ipart in range(npart):
ix = int(partx[ipart]/dhocx)
......@@ -70,18 +72,20 @@ def hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy):
hoc[iy, ix] = ipart
return hoc, hoclist
def hocFind(px, py, dhocx, dhocy, hoc, hoclist):
ix = int(px/dhocx)
iy = int(py/dhocy)
neigh=[]
neigh = []
it = hoc[iy, ix]
while it != -1:
neigh.append(it)
it = hoclist[it]
return neigh
def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None):
def findNeighbors_hoclist(px, py, tx=None, ty=None, dn=4, hoc=None, hoclist=None):
nhocy = nhocx = 20
pxMin = np.min(px)
......@@ -91,21 +95,21 @@ def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None)
dhocx = (pxMax - pxMin)/(nhocx-1)
dhocy = (pyMax - pyMin)/(nhocy-1)
partx = px - pxMin +dhocx/2
party = py - pyMin +dhocy/2
partx = px - pxMin + dhocx/2
party = py - pyMin + dhocy/2
if hoc is None:
hoc, hoclist = hocBuild(partx, party, nhocx, nhocy, dhocx, dhocy)
return hoc, hoclist
if hoc is not None:
tx = tx - pxMin +dhocx/2
ty = ty - pyMin +dhocy/2
tx = tx - pxMin + dhocx/2
ty = ty - pyMin + dhocy/2
itx = int(tx/dhocx)
ity = int(ty/dhocy)
ps = [-1, 0, 1]
neigh=[]
neigh = []
for ii in range(3):
for jj in range(3):
ix = itx + ps[ii]
......@@ -119,23 +123,23 @@ def findNeighbors_hoclist(px, py, tx=None,ty=None, dn=4, hoc=None, hoclist=None)
if iy > nhocy-1:
continue
#neightt = myUtil.hocFind(ppx, ppy, dhocx, dhocy, hoc, hoclist)
# neightt = myUtil.hocFind(ppx, ppy, dhocx, dhocy, hoc, hoclist)
it = hoc[iy, ix]
while it != -1:
neigh.append(it)
it = hoclist[it]
#neigh.append(neightt)
#ll = [i for k in neigh for i in k]
# neigh.append(neightt)
# ll = [i for k in neigh for i in k]
if dn != -1:
ptx = np.array(partx[neigh])
pty = np.array(party[neigh])
dd = np.hypot(ptx-tx, pty-ty)
dd = np.hypot(ptx-tx, pty-ty)
idx = np.argsort(dd)
neigh= np.array(neigh)[idx[0:dn]]
neigh = np.array(neigh)[idx[0:dn]]
return neigh
###PSF-IDW###
### PSF-IDW###
def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=None, hoclist=None, PSFCentroidWgt=False):
"""
psf interpolation by IDW
......@@ -161,9 +165,11 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=Tru
if OnlyNeighbors == True:
if hoc is None:
neigh = findNeighbors(px, py, cen_col, cen_row, dr=5., dn=4, OnlyDistance=False)
neigh = findNeighbors(px, py, cen_col, cen_row,
dr=5., dn=4, OnlyDistance=False)
if hoc is not None:
neigh = findNeighbors_hoclist(cen_col, cen_row, tx=px,ty=py, dn=4, hoc=hoc, hoclist=hoclist)
neigh = findNeighbors_hoclist(
cen_col, cen_row, tx=px, ty=py, dn=4, hoc=hoc, hoclist=hoclist)
neighFlag = np.zeros(npsf)
neighFlag[neigh] = 1
......@@ -173,7 +179,8 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=Tru
if neighFlag[ipsf] != 1:
continue
dist = np.sqrt((ref_col - cen_col[ipsf])**2 + (ref_row - cen_row[ipsf])**2)
dist = np.sqrt((ref_col - cen_col[ipsf])
** 2 + (ref_row - cen_row[ipsf])**2)
if IDWindex == 1:
psfWeight[ipsf] = dist
if IDWindex == 2:
......@@ -186,7 +193,7 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=Tru
psfWeight[ipsf] = 1./psfWeight[ipsf]
psfWeight /= np.sum(psfWeight)
psfMaker = np.zeros([ngy, ngx], dtype=np.float32)
psfMaker = np.zeros([ngy, ngx], dtype=np.float32)
for ipsf in range(npsf):
if OnlyNeighbors == True:
if neighFlag[ipsf] != 1:
......@@ -201,15 +208,14 @@ def psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=Tru
return psfMaker
###define PSFInterp###
### define PSFInterp###
class PSFInterp(PSFModel):
def __init__(self, chip, npsf=NPSF, PSF_data=None, PSF_data_file=None, PSF_data_prefix="", sigSpin=0, psfRa=0.15, HocBuild=False, LOG_DEBUG=False):
self.LOG_DEBUG = LOG_DEBUG
if self.LOG_DEBUG:
print('===================================================')
print('DEBUG: psf module for csstSim ' \
+time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
print('DEBUG: psf module for csstSim '
+ time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
print('===================================================')
self.sigSpin = sigSpin
......@@ -221,53 +227,62 @@ class PSFInterp(PSFModel):
print('Error - PSF_data_file is None')
sys.exit()
self.nwave= self._getPSFwave(self.iccd, PSF_data_file, PSF_data_prefix)
self.nwave = self._getPSFwave(
self.iccd, PSF_data_file, PSF_data_prefix)
self.npsf = npsf
self.PSF_data = self._loadPSF(self.iccd, PSF_data_file, PSF_data_prefix)
self.PSF_data = self._loadPSF(
self.iccd, PSF_data_file, PSF_data_prefix)
if self.LOG_DEBUG:
print('nwave-{:} on ccd-{:}::'.format(self.nwave, self.iccd), flush=True)
print('nwave-{:} on ccd-{:}::'.format(self.nwave,
self.iccd), flush=True)
print('self.PSF_data ... ok', flush=True)
print('Preparing self.[psfMat,cen_col,cen_row] for psfMaker ... ', end='', flush=True)
ngy, ngx = self.PSF_data[0][0]['psfMat'].shape
self.psfMat = np.zeros([self.nwave, self.npsf, ngy, ngx], dtype=np.float32)
self.cen_col= np.zeros([self.nwave, self.npsf], dtype=np.float32)
self.cen_row= np.zeros([self.nwave, self.npsf], dtype=np.float32)
self.hoc =[]
self.hoclist=[]
print(
'Preparing self.[psfMat,cen_col,cen_row] for psfMaker ... ', end='', flush=True)
ngy, ngx = self.PSF_data[0][0]['psfMat'].shape
self.psfMat = np.zeros(
[self.nwave, self.npsf, ngy, ngx], dtype=np.float32)
self.cen_col = np.zeros([self.nwave, self.npsf], dtype=np.float32)
self.cen_row = np.zeros([self.nwave, self.npsf], dtype=np.float32)
self.hoc = []
self.hoclist = []
for twave in range(self.nwave):
for tpsf in range(self.npsf):
self.psfMat[twave, tpsf, :, :] = self.PSF_data[twave][tpsf]['psfMat']
self.PSF_data[twave][tpsf]['psfMat'] = 0 ###free psfMat
self.psfMat[twave, tpsf, :,
:] = self.PSF_data[twave][tpsf]['psfMat']
self.PSF_data[twave][tpsf]['psfMat'] = 0 # free psfMat
self.pixsize = self.PSF_data[twave][tpsf]['pixsize']*1e-3 ##mm
self.cen_col[twave, tpsf] = self.PSF_data[twave][tpsf]['image_x'] + self.PSF_data[twave][tpsf]['centroid_x']
self.cen_row[twave, tpsf] = self.PSF_data[twave][tpsf]['image_y'] + self.PSF_data[twave][tpsf]['centroid_y']
self.pixsize = self.PSF_data[twave][tpsf]['pixsize']*1e-3 # mm
self.cen_col[twave, tpsf] = self.PSF_data[twave][tpsf]['image_x'] + \
self.PSF_data[twave][tpsf]['centroid_x']
self.cen_row[twave, tpsf] = self.PSF_data[twave][tpsf]['image_y'] + \
self.PSF_data[twave][tpsf]['centroid_y']
if HocBuild:
#hoclist on twave for neighborsFinding
hoc,hoclist = findNeighbors_hoclist(self.cen_col[twave], self.cen_row[twave])
# hoclist on twave for neighborsFinding
hoc, hoclist = findNeighbors_hoclist(
self.cen_col[twave], self.cen_row[twave])
self.hoc.append(hoc)
self.hoclist.append(hoclist)
if self.LOG_DEBUG:
print('ok', flush=True)
def _getPSFwave(self, iccd, PSF_data_file, PSF_data_prefix):
# fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r')
fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r')
fq = h5py.File(PSF_data_file+'/' + PSF_data_prefix +
'psfCube_{:}.h5'.format(iccd), 'r')
nwave = len(fq.keys())
fq.close()
return nwave
def _loadPSF(self, iccd, PSF_data_file, PSF_data_prefix):
psfSet = []
# fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_ccd{:}.h5'.format(iccd), 'r')
fq = h5py.File(PSF_data_file+'/' +PSF_data_prefix +'psfCube_{:}.h5'.format(iccd), 'r')
fq = h5py.File(PSF_data_file+'/' + PSF_data_prefix +
'psfCube_{:}.h5'.format(iccd), 'r')
for ii in range(self.nwave):
iwave = ii+1
psfWave = []
......@@ -276,30 +291,30 @@ class PSFInterp(PSFModel):
for jj in range(self.npsf):
ipsf = jj+1
psfInfo = {}
psfInfo['wavelength']= fq_iwave['wavelength'][()]
psfInfo['wavelength'] = fq_iwave['wavelength'][()]
fq_iwave_ipsf = fq_iwave['psf_{:}'.format(ipsf)]
psfInfo['pixsize'] = PixSizeInMicrons
psfInfo['field_x'] = fq_iwave_ipsf['field_x'][()]
psfInfo['field_y'] = fq_iwave_ipsf['field_y'][()]
psfInfo['image_x'] = fq_iwave_ipsf['image_x'][()]
psfInfo['image_y'] = fq_iwave_ipsf['image_y'][()]
psfInfo['centroid_x']= fq_iwave_ipsf['cx'][()]
psfInfo['centroid_y']= fq_iwave_ipsf['cy'][()]
psfInfo['psfMat'] = fq_iwave_ipsf['psfMat'][()]
psfInfo['pixsize'] = PixSizeInMicrons
psfInfo['field_x'] = fq_iwave_ipsf['field_x'][()]
psfInfo['field_y'] = fq_iwave_ipsf['field_y'][()]
psfInfo['image_x'] = fq_iwave_ipsf['image_x'][()]
psfInfo['image_y'] = fq_iwave_ipsf['image_y'][()]
psfInfo['centroid_x'] = fq_iwave_ipsf['cx'][()]
psfInfo['centroid_y'] = fq_iwave_ipsf['cy'][()]
psfInfo['psfMat'] = fq_iwave_ipsf['psfMat'][()]
psfWave.append(psfInfo)
psfSet.append(psfWave)
fq.close()
if self.LOG_DEBUG:
print('psfSet has been loaded:', flush=True)
print('psfSet[iwave][ipsf][keys]:', psfSet[0][0].keys(), flush=True)
print('psfSet[iwave][ipsf][keys]:',
psfSet[0][0].keys(), flush=True)
return psfSet
def _findWave(self, bandpass):
if isinstance(bandpass,int):
if isinstance(bandpass, int):
twave = bandpass
return twave
......@@ -308,7 +323,6 @@ class PSFInterp(PSFModel):
if bandpass.blue_limit < bandwave and bandwave < bandpass.red_limit:
return twave
return -1
def get_PSF(self, chip, pos_img, bandpass, galsimGSObject=True, findNeighMode='treeFind', folding_threshold=5.e-3, pointing_pa=0.0):
"""
......@@ -323,7 +337,7 @@ class PSFInterp(PSFModel):
Returns:
PSF: A 'galsim.GSObject'.
"""
pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize
pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 # set psf pixsize
# assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID'
twave = self._findWave(bandpass)
......@@ -331,16 +345,18 @@ class PSFInterp(PSFModel):
print("!!!PSF bandpass does not match.")
exit()
PSFMat = self.psfMat[twave]
cen_col= self.cen_col[twave]
cen_row= self.cen_row[twave]
cen_col = self.cen_col[twave]
cen_row = self.cen_row[twave]
px = (pos_img.x - chip.cen_pix_x)*0.01
py = (pos_img.y - chip.cen_pix_y)*0.01
if findNeighMode == 'treeFind':
imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, PSFCentroidWgt=True)
imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row,
IDWindex=2, OnlyNeighbors=True, PSFCentroidWgt=True)
if findNeighMode == 'hoclistFind':
assert(self.hoc != 0), 'hoclist should be built correctly!'
imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True, hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True)
assert (self.hoc != 0), 'hoclist should be built correctly!'
imPSF = psfMaker_IDW(px, py, PSFMat, cen_col, cen_row, IDWindex=2, OnlyNeighbors=True,
hoc=self.hoc[twave], hoclist=self.hoclist[twave], PSFCentroidWgt=True)
'''
############TEST: START
......@@ -353,20 +369,21 @@ class PSFInterp(PSFModel):
'''
if galsimGSObject:
imPSFt = np.zeros([257,257])
imPSFt = np.zeros([257, 257])
imPSFt[0:256, 0:256] = imPSF
# imPSFt[120:130, 0:256] = 1.
img = galsim.ImageF(imPSFt, scale=pixSize)
gsp = galsim.GSParams(folding_threshold=folding_threshold)
############TEST: START
# TEST: START
# Use sheared PSF to test the PSF orientation
# self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.)
############TEST: END
# TEST: END
self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
wcs = chip.img.wcs.local(pos_img)
scale = galsim.PixelScale(0.074)
self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img))
self.psf = wcs.toWorld(scale.toImage(
self.psf), image_pos=(pos_img))
# return self.PSFspin(x=px/0.01, y=py/0.01)
return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
......
ObservationSim/PSF/PSFInterpSLS.py → observation_sim/PSF/PSFInterpSLS.py
View file @ 36189a3e
......@@ -4,6 +4,7 @@ PSF interpolation for CSST-Sim
NOTE: [iccd, iwave, ipsf] are counted from 1 to n, but [tccd, twave, tpsf] are counted from 0 to n-1
'''
import yaml
import sys
import time
import copy
......@@ -12,8 +13,9 @@ import scipy.spatial as spatial
import galsim
import h5py
from ObservationSim.PSF.PSFModel import PSFModel
from ObservationSim.Instrument.Chip import ChipUtils as chip_utils
from observation_sim.instruments import Filter, FilterParam, Chip
from observation_sim.PSF.PSFModel import PSFModel
from observation_sim.instruments.chip import chip_utils
import os
from astropy.io import fits
......@@ -21,12 +23,12 @@ from astropy.modeling.models import Gaussian2D
from scipy import signal
LOG_DEBUG = False #***#
NPSF = 900 #***# 30*30
PixSizeInMicrons = 5. #***# in microns
LOG_DEBUG = False # ***#
NPSF = 900 # ***# 30*30
PIX_SIZE_MICRON = 5. # ***# in microns
###find neighbors-KDtree###
### find neighbors-KDtree###
# def findNeighbors(tx, ty, px, py, dr=0.1, dn=1, OnlyDistance=True):
# """
# find nearest neighbors by 2D-KDTree
......@@ -208,36 +210,35 @@ PixSizeInMicrons = 5. #***# in microns
# return psfMaker
###define PSFInterp###
### define PSFInterp###
class PSFInterpSLS(PSFModel):
def __init__(self, chip, filt,PSF_data_prefix="", sigSpin=0, psfRa=0.15, pix_size = 0.005):
def __init__(self, chip, filt, PSF_data_prefix="", sigSpin=0, psfRa=0.15, pix_size=0.005):
if LOG_DEBUG:
print('===================================================')
print('DEBUG: psf module for csstSim ' \
+time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
print('DEBUG: psf module for csstSim '
+ time.strftime("(%Y-%m-%d %H:%M:%S)", time.localtime()), flush=True)
print('===================================================')
self.sigSpin = sigSpin
self.sigGauss = psfRa
self.grating_ids = chip_utils.getChipSLSGratingID(chip.chipID)
_,self.grating_type = chip.getChipFilter(chipID=chip.chipID)
_, self.grating_type = chip.getChipFilter(chipID=chip.chipID)
self.data_folder = PSF_data_prefix
self.getPSFDataFromFile(filt)
self.pixsize = pix_size # um
self.pixsize = pix_size # um
def getPSFDataFromFile(self, filt):
gratingInwavelist = {'GU':0,'GV':1,'GI':2}
grating_orders = ['0','1']
gratingInwavelist = {'GU': 0, 'GV': 1, 'GI': 2}
grating_orders = ['0', '1']
waveListFn = self.data_folder + '/wavelist.dat'
wavelists = np.loadtxt(waveListFn)
self.waveList = wavelists[:,gratingInwavelist[self.grating_type]]
bandranges = np.zeros([4,2])
self.waveList = wavelists[:, gratingInwavelist[self.grating_type]]
bandranges = np.zeros([4, 2])
midBand = (self.waveList[0:3] + self.waveList[1:4])/2.*10000.
bandranges[0,0] = filt.blue_limit
bandranges[1:4,0] = midBand
bandranges[0, 0] = filt.blue_limit
bandranges[1:4, 0] = midBand
bandranges[0:3, 1] = midBand
bandranges[3,1] = filt.red_limit
bandranges[3, 1] = filt.red_limit
self.bandranges = bandranges
......@@ -246,7 +247,7 @@ class PSFInterpSLS(PSFModel):
for g_order in grating_orders:
g_folder_order = g_folder + 'PSF_Order_' + g_order + '/'
grating_order_data = {}
for bandi in [1,2,3,4]:
for bandi in [1, 2, 3, 4]:
subBand_data = {}
subBand_data['bandrange'] = bandranges[bandi-1]
final_folder = g_folder_order + str(bandi) + '/'
......@@ -305,7 +306,7 @@ class PSFInterpSLS(PSFModel):
# psfInfo['wavelength']= fq_iwave['wavelength'][()]
#
# fq_iwave_ipsf = fq_iwave['psf_{:}'.format(ipsf)]
# psfInfo['pixsize'] = PixSizeInMicrons
# psfInfo['pixsize'] = PIX_SIZE_MICRON
# psfInfo['field_x'] = fq_iwave_ipsf['field_x'][()]
# psfInfo['field_y'] = fq_iwave_ipsf['field_y'][()]
# psfInfo['image_x'] = fq_iwave_ipsf['image_x'][()]
......@@ -342,7 +343,7 @@ class PSFInterpSLS(PSFModel):
offset = int(np.ceil(sigma * 3))
g_size = 2 * offset + 1
m_cen = int(g_size / 2)
print('-----',g_size)
print('-----', g_size)
g_PSF_ = Gaussian2D(1, m_cen, m_cen, sigma, sigma)
yp, xp = np.mgrid[0:g_size, 0:g_size]
g_PSF = g_PSF_(xp, yp)
......@@ -351,8 +352,7 @@ class PSFInterpSLS(PSFModel):
convImg = convImg/np.sum(convImg)
return convImg
def get_PSF(self, chip, pos_img_local = [1000,1000], bandNo = 1, galsimGSObject=True, folding_threshold=5.e-3, g_order = 'A', grating_split_pos=3685):
def get_PSF(self, chip, pos_img_local=[1000, 1000], bandNo=1, galsimGSObject=True, folding_threshold=5.e-3, g_order='A', grating_split_pos=3685):
"""
Get the PSF at a given image position
......@@ -365,8 +365,9 @@ class PSFInterpSLS(PSFModel):
Returns:
PSF: A 'galsim.GSObject'.
"""
order_IDs = {'A': '1', 'B': '0' ,'C': '0', 'D': '0', 'E': '0'}
contam_order_sigma = {'C':0.28032344707964174,'D':0.39900182912061344,'E':1.1988309797685412} #arcsec
order_IDs = {'A': '1', 'B': '0', 'C': '0', 'D': '0', 'E': '0'}
contam_order_sigma = {'C': 0.28032344707964174,
'D': 0.39900182912061344, 'E': 1.1988309797685412} # arcsec
x_start = chip.x_cen/chip.pix_size - chip.npix_x / 2.
y_start = chip.y_cen/chip.pix_size - chip.npix_y / 2.
# print(pos_img.x - x_start)
......@@ -382,7 +383,6 @@ class PSFInterpSLS(PSFModel):
# if grating_order in ['-2','-1','2']:
# grating_order = '1'
# if grating_order in ['0', '1']:
psf_order = psf_data['order'+grating_order]
psf_order_b = psf_order['band'+str(bandNo)]
......@@ -396,31 +396,34 @@ class PSFInterpSLS(PSFModel):
px = pos_img.x*chip.pix_size
py = pos_img.y*chip.pix_size
dist2=(pos_p[:,1] - px)*(pos_p[:,1] - px) + (pos_p[:,0] - py)*(pos_p[:,0] - py)
temp_sort_dist = np.zeros([dist2.shape[0],2])
temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0],1)
dist2 = (pos_p[:, 1] - px)*(pos_p[:, 1] - px) + \
(pos_p[:, 0] - py)*(pos_p[:, 0] - py)
temp_sort_dist = np.zeros([dist2.shape[0], 2])
temp_sort_dist[:, 0] = np.arange(0, dist2.shape[0], 1)
temp_sort_dist[:, 1] = dist2
# print(temp_sort_dist)
dits2_sortlist = sorted(temp_sort_dist, key=lambda x:x[1])
dits2_sortlist = sorted(temp_sort_dist, key=lambda x: x[1])
# print(dits2_sortlist)
nearest4p = np.zeros([4,2])
pc_coeff_4p = np.zeros([pc_coeff.data.shape[0],4])
nearest4p = np.zeros([4, 2])
pc_coeff_4p = np.zeros([pc_coeff.data.shape[0], 4])
for i in np.arange(4):
smaller_ids = int(dits2_sortlist[i][0])
nearest4p[i, 0] = pos_p[smaller_ids, 1]
nearest4p[i, 1] = pos_p[smaller_ids, 0]
pc_coeff_4p[:,i] = pc_coeff[:,smaller_ids]
idw_dist = 1/(np.sqrt((px-nearest4p[:,0]) * (px-nearest4p[:,0]) + (py-nearest4p[:,1]) * (py-nearest4p[:,1])))
pc_coeff_4p[:, i] = pc_coeff[:, smaller_ids]
idw_dist = 1/(np.sqrt((px-nearest4p[:, 0]) * (px-nearest4p[:, 0]) + (
py-nearest4p[:, 1]) * (py-nearest4p[:, 1])))
coeff_int = np.zeros(pc_coeff.data.shape[0])
for i in np.arange(4):
coeff_int = coeff_int + pc_coeff_4p[:,i]*idw_dist[i]
coeff_int = coeff_int + pc_coeff_4p[:, i]*idw_dist[i]
coeff_int = coeff_int / np.sum(coeff_int)
npc = 10
m_size = int(pcs.shape[0]**0.5)
PSF_int = np.dot(pcs[:,0:npc],coeff_int[0:npc]).reshape(m_size,m_size)
PSF_int = np.dot(pcs[:, 0:npc], coeff_int[0:npc]
).reshape(m_size, m_size)
# PSF_int = PSF_int/np.sum(PSF_int)
PSF_int_trans = np.flipud(np.fliplr(PSF_int))
......@@ -434,7 +437,6 @@ class PSFInterpSLS(PSFModel):
# from astropy.io import fits
# fits.writeto(str(bandNo) + '_' + g_order+ '_psf_o.fits', PSF_int_trans)
# if g_order in ['C','D','E']:
# g_simgma = contam_order_sigma[g_order]/pixel_size_arc
# PSF_int_trans = self.convolveWithGauss(PSF_int_trans,g_simgma)
......@@ -459,24 +461,24 @@ class PSFInterpSLS(PSFModel):
pixel_size_arc = np.rad2deg(self.pixsize * 1e-3 / 28) * 3600
img = galsim.ImageF(PSF_int_trans, scale=pixel_size_arc)
gsp = galsim.GSParams(folding_threshold=folding_threshold)
############TEST: START
# TEST: START
# Use sheared PSF to test the PSF orientation
# self.psf = galsim.InterpolatedImage(img, gsparams=gsp).shear(g1=0.8, g2=0.)
############TEST: END
# TEST: END
self.psf = galsim.InterpolatedImage(img, gsparams=gsp)
# if g_order in ['C','D','E']:
# add_psf = galsim.Gaussian(sigma=contam_order_sigma[g_order], flux=1.0)
# self.psf = galsim.Convolve(self.psf, add_psf)
wcs = chip.img.wcs.local(pos_img)
scale = galsim.PixelScale(0.074)
self.psf = wcs.toWorld(scale.toImage(self.psf), image_pos=(pos_img))
self.psf = wcs.toWorld(scale.toImage(
self.psf), image_pos=(pos_img))
# return self.PSFspin(x=px/0.01, y=py/0.01)
return self.psf, galsim.Shear(e=0., beta=(np.pi/2)*galsim.radians)
return PSF_int_trans, PSF_int
# pixSize = np.rad2deg(self.pixsize*1e-3/28)*3600 #set psf pixsize
#
# # assert self.iccd == int(chip.getChipLabel(chipID=chip.chipID)), 'ERROR: self.iccd != chip.chipID'
......@@ -547,24 +549,23 @@ class PSFInterpSLS(PSFModel):
# PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
# return self.psf.shear(PSFshear), PSFshear
from ObservationSim.Instrument import Filter, FilterParam, Chip
import yaml
if __name__ == '__main__':
configfn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_new_sim/csst-simulation/config/config_C6_dev.yaml'
with open(configfn, "r") as stream:
try:
config = yaml.safe_load(stream)
for key, value in config.items():
print (key + " : " + str(value))
print(key + " : " + str(value))
except yaml.YAMLError as exc:
print(exc)
chip = Chip(chipID=1,config=config)
chip = Chip(chipID=1, config=config)
filter_id, filter_type = chip.getChipFilter()
filt = Filter(filter_id=filter_id,
filter_type=filter_type,
filter_param=FilterParam())
psf_i = PSFInterpSLS(chip, filt,PSF_data_prefix="/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/")
psf_i = PSFInterpSLS(
chip, filt, PSF_data_prefix="/Volumes/EAGET/CSST_PSF_data/SLS_PSF_PCA_fp/")
pos_img = galsim.PositionD(x=25155, y=-22060)
psf_im = psf_i.get_PSF(chip, pos_img = pos_img, g_order = '1')
psf_im = psf_i.get_PSF(chip, pos_img=pos_img, g_order='1')
observation_sim/PSF/PSFModel.py 0 → 100755
View file @ 36189a3e
import galsim
import sep
import numpy as np
from scipy.interpolate import interp1d
import pylab as pl
import os
import sys
class PSFModel(object):
def __init__(self, sigSpin=0., psfRa=0.15):
# TODO: what are the nesseary fields in PSFModel class?
pass
def PSFspin(self, psf, sigSpin, sigGauss, dx, dy):
"""
The PSF profile at a given image position relative to the axis center
Parameters:
theta : spin angles in a given exposure in unit of [arcsecond]
dx, dy: relative position to the axis center in unit of [pixels]
Return:
Spinned PSF: g1, g2 and axis ratio 'a/b'
"""
a2Rad = np.pi/(60.0*60.0*180.0)
ff = sigGauss * 0.107 * (1000.0/10.0) # in unit of [pixels]
rc = np.sqrt(dx*dx + dy*dy)
cpix = rc*(sigSpin*a2Rad)
beta = (np.arctan2(dy, dx) + np.pi/2)
ell = cpix**2/(2.0*ff**2+cpix**2)
# ell *= 10.0
qr = np.sqrt((1.0+ell)/(1.0-ell))
# psfShape = galsim.Shear(e=ell, beta=beta)
# g1, g2 = psfShape.g1, psfShape.g2
# qr = np.sqrt((1.0+ell)/(1.0-ell))
# return ell, beta, qr
PSFshear = galsim.Shear(e=ell, beta=beta*galsim.radians)
return psf.shear(PSFshear), PSFshear
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