diff --git a/tests/PSFInterpTest/readme b/tests/PSFInterpTest/readme deleted file mode 100644 index 4eca40ef16c9d63bc22a7f5b16ec11299999511b..0000000000000000000000000000000000000000 --- a/tests/PSFInterpTest/readme +++ /dev/null @@ -1,7 +0,0 @@ -unittest of PSF module: - --test_loadPSFSet.py - --testdata: S30x30_psfCube_1.h5 config_test.yaml - --test_PSFInterpModule_coverage.py - --testdata: S30x30_psfCube_1.h5 S20x20_psfCube_1.h5 config_test.yaml - --test_Convolve.py - --testdata: ccd1-w1 diff --git a/tests/PSFInterpTest/test_Convolve.py b/tests/PSFInterpTest/test_Convolve.py deleted file mode 100644 index 94aa93f0768f0c3795864af4a5cf3342abc0a147..0000000000000000000000000000000000000000 --- a/tests/PSFInterpTest/test_Convolve.py +++ /dev/null @@ -1,176 +0,0 @@ -''' -test galsim.interpolatedImage & galsim.convolve -''' - -import os -import unittest - -import numpy as np -import matplotlib.pyplot as plt -import galsim -import scipy.io -from scipy import ndimage - - -def psfEncircle(img, fraction=0.8, psfSampleSizeInMicrons=2.5, focalLengthInMeters=28, cenPix=None): - #imgMaxPix_x, imgMaxPix_y = findMaxPix(img) - y,x = ndimage.center_of_mass(img) #y-rows, x-cols - imgMaxPix_x = x #int(x) - imgMaxPix_y = y #int(y) - if cenPix != None: - imgMaxPix_x = cenPix[0] - imgMaxPix_y = cenPix[1] - - im1 = img.copy() - im1size = im1.shape - - dis = np.zeros_like(img) - for irow in range(im1size[0]): - for icol in range(im1size[1]): - dx = icol - imgMaxPix_x - dy = irow - imgMaxPix_y - dis[irow, icol] = np.hypot(dx, dy) - - nn = im1size[1]*im1size[0] - disX = dis.reshape(nn) - disXsortId = np.argsort(disX) - - imgX = img.reshape(nn) - imgY = imgX[disXsortId] - psfFrac = np.cumsum(imgY)/np.sum(imgY) - ind = np.where(psfFrac > fraction)[0][0] - - REE80 = np.rad2deg(dis[np.where(img == imgY[ind])]*psfSampleSizeInMicrons*1e-6/focalLengthInMeters)*3600 - return REE80 - - - -def check_galsimConvolve(path=None, plotImage=True): - #load psf data - data=scipy.io.loadmat(path) - imPSF = data['psf'] - pixSize = np.rad2deg(5.*1e-6/28)*3600 - imPSF = imPSF/np.sum(imPSF) - - #psf -> galsimInterpolatedImage - img = galsim.ImageF(imPSF, scale=pixSize) - imgt= galsim.InterpolatedImage(img) - - #imPSFt = imgt.drawImage(nx=256, ny=256, scale=pixSize, method='no_pixel') - imPSFt = imgt.drawImage(nx=256, ny=256, scale=pixSize) - - ree80 = psfEncircle(imPSF, fraction=0.8, psfSampleSizeInMicrons=5.) - ree80_pix = ree80/(np.rad2deg((5.*1e-6/28))*3600) - - sliceX = slice(128-int(np.round(ree80_pix[0])), 128+int(np.round(ree80_pix[0]))+1, 1) - - #set a point sorce - src = galsim.DeltaFunction(flux=1.0) - - result = galsim.Convolve(src, imgt) - - #drawImage with same pixSize - #tmp = result.drawImage(nx=256, ny=256, scale=pixSize, method='no_pixel') - tmp = result.drawImage(nx=256, ny=256, scale=pixSize) - - - if plotImage != True: - res = (imPSFt.array - imPSF)/imPSF - d0 = np.mean(res[sliceX, sliceX].flatten()) - - res = (tmp.array - imPSFt.array)/imPSFt.array - d1 = np.mean(res[sliceX, sliceX].flatten()) - return d0, d1 - - - #plot images - fig = plt.figure(figsize=(22, 5)) - ax=plt.subplot(1,3,1) - plt.imshow(imPSF[128-10:128+10, 128-10:128+10]) - plt.annotate("ORG", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w") - plt.colorbar() - - ax=plt.subplot(1,3,2) - plt.imshow(imPSFt.array[128-10:128+10, 128-10:128+10]) - plt.annotate("InterpolatedImage", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w") - plt.colorbar() - - ax=plt.subplot(1,3,3) - plt.imshow(tmp.array[128-10:128+10, 128-10:128+10]) - plt.annotate("ConvolvedImage", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w") - plt.colorbar() - plt.savefig(OUTPUTPATH+'/fig_check_galsimConvolve_1.pdf') - - fig = plt.figure(figsize=(13, 10)) - ax=plt.subplot(2,2,1) - res = (imPSFt.array - imPSF)/imPSF - plt.imshow(res[128-10:128+10, 128-10:128+10]) - plt.annotate("$\Delta_1$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w") - plt.colorbar() - - ax=plt.subplot(2,2,2) - plt.hist(res[sliceX, sliceX].flatten(), alpha=0.75, bins=4) - #plt.annotate("$\Delta_1^{\\rm REE80}$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="k") - plt.xlabel("$\Delta_1^{\\rm REE80}$", fontsize=16) - plt.ylabel("PDF", fontsize=16) - - ax=plt.subplot(2,2,3) - res = (tmp.array - imPSFt.array)/imPSFt.array - plt.imshow(res[128-10:128+10, 128-10:128+10]) - plt.annotate("$\Delta_2$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="w") - plt.colorbar() - - ax=plt.subplot(2,2,4) - plt.hist(res[sliceX, sliceX].flatten(), alpha=0.75, bins=4) - #plt.annotate("$\Delta_2^{\\rm REE80}$", [0.1, 0.9],xycoords="axes fraction", fontsize=16, color="k") - plt.xlabel("$\Delta_2^{\\rm REE80}$", fontsize=16) - plt.ylabel("PDF", fontsize=16) - - plt.savefig(OUTPUTPATH+'/fig_check_galsimConvolve_2.pdf') - - -def check_galsimConvolveALL(dataPath): - d0 = np.zeros(900) - d1 = np.zeros(900) - for ipsf in range(1,901,1): - print("ipsf={:}".format(ipsf), end="\r") - psfPath = dataPath+"/ccd1-w1/psf_{:}_centroidWgt_BC.mat".format(ipsf) - t0,t1=check_galsimConvolve(path = psfPath, plotImage=False) - d0[ipsf-1] = t0 - d1[ipsf-1] = t1 - - fig = plt.figure(figsize=(12,6)) - ax = plt.subplot(1,2,1) - #plt.scatter(np.linspace(1,900,900), d0) - plt.hist(d0, bins=8, alpha=0.75) - plt.xlabel("mean($\Delta_1^{\\rm REE80}$)", fontsize=16) - plt.ylabel("PDF", fontsize=16) - - ax = plt.subplot(1,2,2) - #plt.scatter(np.linspace(1,900,900), d1) - plt.hist(d1, bins=8, alpha=0.75) - plt.xlabel("mean($\Delta_2^{\\rm REE80}$)", fontsize=16) - plt.ylabel("PDF", fontsize=16) - - plt.savefig(OUTPUTPATH+'/fig_check_galsimConvolveALL.pdf') - - -class testConvolve(unittest.TestCase): - def __init__(self, methodName='runTest'): - super(testConvolve,self).__init__(methodName) - self.dataPath = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1') - - global OUTPUTPATH - OUTPUTPATH = os.path.join(self.dataPath, 'outputs') - - def test_galsimConvolve(self): - ipsf = 1 - psfPath = self.dataPath+"/ccd1-w1/psf_{:}_centroidWgt_BC.mat".format(ipsf) - check_galsimConvolve(path = psfPath) - - def test_galsimConvolveALL(self): - check_galsimConvolveALL(dataPath=self.dataPath) - - -if __name__ == "__main__": - unittest.main() diff --git a/tests/PSFInterpTest/test_PSFInterpModule_coverage.py b/tests/PSFInterpTest/test_PSFInterpModule_coverage.py deleted file mode 100644 index 6330c1efe1437e71717e9026d0420c3114718311..0000000000000000000000000000000000000000 --- a/tests/PSFInterpTest/test_PSFInterpModule_coverage.py +++ /dev/null @@ -1,403 +0,0 @@ -import unittest - -import sys,os,math -from itertools import islice - -import numpy as np -import matplotlib.pyplot as plt -import matplotlib as mpl -mpl.use('Agg') - -import yaml -from ObservationSim.Instrument import Chip -from ObservationSim.PSF.PSFInterp import PSFInterp - - -def defineCCD(iccd, config_file): - with open(config_file, "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) - chip = Chip(chipID=iccd, config=config) - #chip = Chip(chipID=iccd, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"], normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict['sls_dir'], config=config) - return chip - -def psfSecondMoments(psfMat, cenX, cenY, pixSize=1): - apr = 0.5 #arcsec, 0.5角秒内测量 - fl = 28. #meters - pxs = 2.5 #microns - apr = np.deg2rad(apr/3600.)*fl*1e6 - apr = apr/pxs - apr = int(np.ceil(apr)) - - I = psfMat - ncol = I.shape[1] - nrow = I.shape[0] - w = 0.0 - w11 = 0.0 - w12 = 0.0 - w22 = 0.0 - for icol in range(ncol): - for jrow in range(nrow): - x = icol*pixSize - cenX - y = jrow*pixSize - cenY - rr = np.sqrt(x*x + y*y) - wgt= 0.0 - if rr <= apr: - wgt = 1.0 - w += I[jrow, icol]*wgt - w11 += x*x*I[jrow, icol]*wgt - w12 += x*y*I[jrow, icol]*wgt - w22 += y*y*I[jrow, icol]*wgt - w11 /= w - w12 /= w - w22 /= w - sz = w11 + w22 - e1 = (w11 - w22)/sz - e2 = 2.0*w12/sz - - return sz, e1, e2 - - -def test_psfEll(iccd, iwave, psfMat): - psfMat_iwave = psfMat.psfMat[iwave-1, :,:,:] - npsf = np.shape(psfMat_iwave)[0] - - imx = np.zeros(npsf) - imy = np.zeros(npsf) - psf_e1 = np.zeros(npsf) - psf_e2 = np.zeros(npsf) - psf_sz = np.zeros(npsf) - - for ipsf in range(1, npsf+1): - print('ipsf-{:}'.format(ipsf), end='\r') - imx[ipsf-1] = psfMat.cen_col[iwave-1, ipsf-1] - imy[ipsf-1] = psfMat.cen_row[iwave-1, ipsf-1] - psfMat_iwave_ipsf = psfMat_iwave[ipsf-1, :, :] - cenX = 256 - cenY = 256 - sz, e1, e2 = psfSecondMoments(psfMat_iwave_ipsf, cenX, cenY, pixSize=1) - psf_e1[ipsf-1] = e1 - psf_e2[ipsf-1] = e2 - psf_sz[ipsf-1] = sz - #print('ell======', ipsf, np.sqrt(e1**2 + e2**2)) - ####### - arr = [imx, imy, psf_e1, psf_e2, psf_sz] - np.save(OUTPUTPATH+'/psfEll{:}_{:}_{:}'.format(int(np.sqrt(npsf)),iccd, iwave), arr) - - -def test_psfEllPlot(OVERPLOT=False): - #if ThisTask == 0: - if True: - prefix = 'psfEll30' - iccd = 1 - iwave= 1 - data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy') - imx= data[0] - imy= data[1] - psf_e1 = data[2] - psf_e2 = data[3] - print(np.shape(imx)) - npsf = np.shape(imx)[0] - - ####### - plt.cla() - plt.close("all") - - fig = plt.figure(figsize=(12, 12)) - plt.subplots_adjust(wspace=0.1, hspace=0.1) - ax = plt.subplot(1, 1, 1) - for ipsf in range(npsf): - plt.plot(imx[ipsf], imy[ipsf], 'r.') - - ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 - ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) - ell *= 15 - lcos = ell*np.cos(ang) - lsin = ell*np.sin(ang) - plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'r', lw=2) - ########### - ang = 0. - ell = 0.05 - ell*= 15 - lcos = ell*np.cos(ang) - lsin = ell*np.sin(ang) - plt.plot([imx[898]-lcos, imx[898]+lcos],[imy[898]+5.-lsin, imy[898]+5.+lsin],'k', lw=2) - plt.annotate('{:}'.format(ell/15), (imx[898]-2., imy[898]+6.), xycoords='data', fontsize=10) - plt.xlabel('CCD X (mm)') - plt.ylabel('CCD Y (mm)') - - if OVERPLOT == True: - prefix = 'psfEll20' - data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy') - imx= data[0] - imy= data[1] - psf_e1 = data[2] - psf_e2 = data[3] - - npsf = np.shape(imx)[0] - for ipsf in range(npsf): - plt.plot(imx[ipsf], imy[ipsf], 'b.') - - ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 - ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) - ell *= 15 - lcos = ell*np.cos(ang) - lsin = ell*np.sin(ang) - plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'b', lw=2) - - plt.gca().set_aspect(1) - if OVERPLOT == True: - prefix = 'psfEllOP' - plt.savefig(OUTPUTPATH+'/'+prefix+'_iccd{:}.pdf'.format(iccd)) - - -def test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB): - bandpass = iwave-1 - class pos_img(): - def __init__(self,x, y): - self.x = x*1e3/10. #in unit of pixels - self.y = y*1e3/10. - - psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:] - npsf = np.shape(psfMat_iwave)[0] - psf_e1 = np.zeros(npsf) - psf_e2 = np.zeros(npsf) - psf_sz = np.zeros(npsf) - - for ipsf in range(1, npsf+1): - print('ipsf:', ipsf, end='\r', flush=True) - tpos_img = pos_img(psfMatA.cen_col[iwave-1, ipsf-1], psfMatA.cen_row[iwave-1, ipsf-1]) - psfIDW = psfMatB.get_PSF(chip, tpos_img, bandpass, galsimGSObject=False, findNeighMode='treeFind') - np.save(OUTPUTPATH+'/psfIDW_{:}_{:}_{:}'.format(iccd, iwave, ipsf), psfIDW) - - cenX = 256 - cenY = 256 - sz, e1, e2 = psfSecondMoments(psfIDW, cenX, cenY, pixSize=1) - psf_e1[ipsf-1] = e1 - psf_e2[ipsf-1] = e2 - psf_sz[ipsf-1] = sz - arr = [psf_e1, psf_e2, psf_sz] - np.save(OUTPUTPATH+'/psfEll20IDW_{:}_{:}'.format(iccd, iwave), arr) - - -def test_psfResidualPlot(iccd, iwave, ipsf, psfMatA): - psfMat_iwave = psfMatA.psfMat[iwave-1, :,:,:] - - psfMatORG = psfMat_iwave[ipsf-1, :, :] - psfMatIDW = np.load(OUTPUTPATH+'/psfIDW_{:}_{:}_{:}.npy'.format(iccd, iwave, ipsf)) - - npix = psfMatORG.shape[0] - pixCutEdge= int(npix/2-15) - - img0 = psfMatORG[pixCutEdge:npix-pixCutEdge, pixCutEdge:npix-pixCutEdge] - img1 = psfMatIDW[pixCutEdge:npix-pixCutEdge, pixCutEdge:npix-pixCutEdge] - imgX = (img1 - img0)/img0 - - img0 = np.log10(img0) - img1 = np.log10(img1) - imgX = np.log10(np.abs(imgX)) - - fig = plt.figure(figsize=(18,4)) - ax = plt.subplot(1,3,1) - plt.imshow(img0, origin='lower', vmin=-7, vmax=-1.3) - plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--') - plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--') - plt.annotate('ORG', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15) - cticks=[-7, -6, -5, -4, -3, -2, -1] - cbar = plt.colorbar(ticks=cticks) - cbar.ax.set_yticklabels(['$10^{-7}$', '$10^{-6}$', '$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$']) - print(img0.min(), img0.max()) - - - ax = plt.subplot(1,3,2) - plt.imshow(img1, origin='lower', vmin=-7, vmax=-1.3) - plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--') - plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--') - plt.annotate('IDW', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15) - cticks=[-7, -6, -5, -4, -3, -2, -1] - cbar = plt.colorbar(ticks=cticks) - cbar.ax.set_yticklabels(['$10^{-7}$', '$10^{-6}$', '$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$']) - print(img1.min(), img1.max()) - - - ax = plt.subplot(1,3,3) - plt.imshow(imgX, origin='lower', vmin =-3, vmax =np.log10(3e-1)) - plt.plot([npix/2-pixCutEdge, npix/2-pixCutEdge],[0, (npix/2-pixCutEdge)*2-1],'w--') - plt.plot([0, (npix/2-pixCutEdge)*2-1],[npix/2-pixCutEdge, npix/2-pixCutEdge],'w--') - #plt.annotate('(IDW-ORG)/ORG', [0,(npix/2-pixCutEdge)*2-5], c='w', size=15) - cticks=[-5, -4, -3, -2, -1] - cbar = plt.colorbar(ticks=cticks) - cbar.ax.set_yticklabels(['$10^{-5}$','$10^{-4}$','$10^{-3}$','$10^{-2}$', '$10^{-1}$']) - - print(np.max((psfMatORG-psfMatIDW))) - plt.savefig(OUTPUTPATH+'/psfResidual_iccd{:}.pdf'.format(iccd)) - - -def test_psfEllIDWPlot(OVERPLOT=False): - #if ThisTask == 0: - if True: - prefix = 'psfEll20' - iccd = 1 - iwave= 1 - data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy') - imx= data[0] - imy= data[1] - psf_e1 = data[2] - psf_e2 = data[3] - print(np.shape(imx)) - npsf = np.shape(imx)[0] - - ####### - plt.cla() - plt.close("all") - - fig = plt.figure(figsize=(12, 12)) - plt.subplots_adjust(wspace=0.1, hspace=0.1) - ax = plt.subplot(1, 1, 1) - for ipsf in range(npsf): - plt.plot(imx[ipsf], imy[ipsf], 'b.') - - ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 - ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) - ell *= 15 - lcos = ell*np.cos(ang) - lsin = ell*np.sin(ang) - plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'b', lw=2) - ########### - ang = 0. - ell = 0.05 - ell*= 15 - lcos = ell*np.cos(ang) - lsin = ell*np.sin(ang) - #plt.plot([imx[898]-lcos, imx[898]+lcos],[imy[898]+5.-lsin, imy[898]+5.+lsin],'k', lw=2) - #plt.annotate('{:}'.format(ell/15), (imx[898]-2., imy[898]+6.), xycoords='data', fontsize=10) - plt.xlabel('CCD X (mm)') - plt.ylabel('CCD Y (mm)') - - if OVERPLOT == True: - prefix = 'psfEll20IDW' - data = np.load(OUTPUTPATH+'/'+prefix+'_1_1.npy') - #imx= data[0] - #imy= data[1] - psf_e1 = data[0] - psf_e2 = data[1] - - npsf = np.shape(imx)[0] - for ipsf in range(npsf): - #plt.plot(imx[ipsf], imy[ipsf], 'r.') - - ang = np.arctan2(psf_e2[ipsf], psf_e1[ipsf])/2 - ell = np.sqrt(psf_e1[ipsf]**2 + psf_e2[ipsf]**2) - ell *= 15 - lcos = ell*np.cos(ang) - lsin = ell*np.sin(ang) - plt.plot([imx[ipsf]-lcos, imx[ipsf]+lcos],[imy[ipsf]-lsin, imy[ipsf]+lsin],'r', lw=1) - - plt.gca().set_aspect(1) - if OVERPLOT == True: - prefix = 'psfEllOPIDW' - plt.savefig(OUTPUTPATH+'/'+prefix+'_iccd{:}.pdf'.format(iccd)) - - -def test_psfdEllabsPlot(iccd): - #if ThisTask == 0: - if True: - prefix = 'psfEll20' - #iccd = 1 - #iwave= 1 - data = np.load(OUTPUTPATH+'/'+prefix+'_{:}_1.npy'.format(iccd)) - imx= data[0] - imy= data[1] - psf_e1 = data[2] - psf_e2 = data[3] - psf_sz = data[4] - print(np.shape(imx)) - npsf = np.shape(imx)[0] - - ellX = np.sqrt(psf_e1**2 + psf_e2**2) - angX = np.arctan2(psf_e2, psf_e1)/2 - angX = np.rad2deg(angX) - szX = psf_sz - - ############################## - prefix = 'psfEll20IDW' - data = np.load(OUTPUTPATH+'/'+prefix+'_{:}_1.npy'.format(iccd)) - #imx= data[0] - #imy= data[1] - psf_e1 = data[0] - psf_e2 = data[1] - psf_sz = data[2] - - ellY = np.sqrt(psf_e1**2 + psf_e2**2) - angY = np.arctan2(psf_e2, psf_e1)/2 - angY = np.rad2deg(angY) - szY = psf_sz - - ############################## - fig=plt.figure(figsize=(6, 5)) - grid = plt.GridSpec(3,1,left=0.15, right=0.95, wspace=None, hspace=0.02) - #plt.subplots_adjust(left=None,bottom=None,right=None,top=None,wspace=None,hspace=0.02) - ax = plt.subplot(grid[0:2,0]) - plt.plot([0.01,0.1],[0.01,0.1], 'k--', lw=1. ) - plt.scatter(ellX, ellY, color='b', alpha=1., s=3., edgecolors='None') - #plt.xlim([0.015, 0.085]) - #plt.ylim([0.015, 0.085]) - plt.ylabel('$\epsilon_{\\rm IDW}$') - plt.gca().axes.get_xaxis().set_visible(False) - - ax = plt.subplot(grid[2,0]) - plt.plot([0.015,0.085],[0.,0.], 'k--', lw=1. ) - plt.scatter(ellX, (ellY-ellX), color='b', s=3., edgecolors='None') - #plt.xlim([0.015, 0.085]) - #plt.ylim([-0.0018, 0.0018]) - plt.xlabel('$\epsilon_{\\rm ORG}$') - plt.ylabel('$\Delta$') - plt.savefig(OUTPUTPATH+'/psfEllOPIDWPDF_{:}.pdf'.format(iccd)) - - fig=plt.figure(figsize=(6, 6)) - plt.hist((szY-szX)/szX, bins=20, color='r', alpha=0.5) - plt.xlabel('$(R_{\\rm IDW}-R_{\\rm ORG})/R_{\\rm ORG}$') - plt.ylabel('PDF') - plt.savefig(OUTPUTPATH+'/psfEllOPIDWPDF_dsz_{:}.pdf'.format(iccd)) - - -class PSFInterpModule_coverage(unittest.TestCase): - def __init__(self, methodName='runTest'): - super(PSFInterpModule_coverage,self).__init__(methodName) - self.dataPath = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1') - global OUTPUTPATH - OUTPUTPATH = os.path.join(self.dataPath, 'outputs') - - def test_psfEll_(self): - iccd = 1 - iwave= 1 - - config_file = os.path.join(self.dataPath, 'config_test.yaml') - chip = defineCCD(iccd, config_file) - print(chip.chipID) - print(chip.cen_pix_x, chip.cen_pix_y) - - psfMatA = PSFInterp(chip, npsf=400, PSF_data_file=self.dataPath, PSF_data_prefix="S20x20_") - psfMatB = PSFInterp(chip, npsf=900, PSF_data_file=self.dataPath, PSF_data_prefix="S30x30_") - - test_psfEll(iccd, iwave, psfMatA) - test_psfEll(iccd, iwave, psfMatB) - test_psfEllPlot(OVERPLOT=True) - - test_psfIDW(iccd, iwave, psfMatA, chip, psfMatB) - ipsf = 1 - test_psfResidualPlot(iccd, iwave, ipsf, psfMatA) - test_psfEllIDWPlot(OVERPLOT=True) - test_psfdEllabsPlot(iccd) - - -if __name__ == '__main__': - unittest.main() - print('#####haha#####') - - - diff --git a/tests/PSFInterpTest/test_loadPSFSet.py b/tests/PSFInterpTest/test_loadPSFSet.py deleted file mode 100644 index 6c65c79d458802f6f61fb90d0017dd2efe67f031..0000000000000000000000000000000000000000 --- a/tests/PSFInterpTest/test_loadPSFSet.py +++ /dev/null @@ -1,62 +0,0 @@ -import unittest - -import sys,os,math -from itertools import islice - -import numpy as np - -import yaml -from ObservationSim.Instrument import Chip -from ObservationSim.PSF.PSFInterp import PSFInterp - - -def defineCCD(iccd, config_file): - with open(config_file, "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) - chip = Chip(chipID=iccd, config=config) - #chip = Chip(chipID=iccd, ccdEffCurve_dir=path_dict["ccd_dir"], CRdata_dir=path_dict["CRdata_dir"], normalize_dir=path_dict["normalize_dir"], sls_dir=path_dict['sls_dir'], config=config) - return chip - -def loadPSFSet(iccd, dataPath): - config_file = os.path.join(dataPath, 'config_test.yaml') - chip = defineCCD(iccd, config_file) - print(chip.chipID) - print(chip.cen_pix_x, chip.cen_pix_y) - - psfMat= PSFInterp(chip, npsf=900, PSF_data_file=dataPath, PSF_data_prefix="S30x30_") - psfSet= psfMat._loadPSF(iccd, dataPath, PSF_data_prefix="S30x30_") - - twave = 0 #[0...3] - tpsf = 0 #[0...899] - field_x= psfSet[twave][tpsf]['field_x'] - field_y= psfSet[twave][tpsf]['field_y'] - image_x= psfSet[twave][tpsf]['image_x'] - image_y= psfSet[twave][tpsf]['image_y'] - centroid_x= psfSet[twave][tpsf]['centroid_x'] - centroid_y= psfSet[twave][tpsf]['centroid_y'] - print("pos_info:", field_x, field_y, image_x, image_y, centroid_x, centroid_y) - return psfSet - - - -class PSFInterpModule_coverage(unittest.TestCase): - def __init__(self, methodName='runTest'): - super(PSFInterpModule_coverage,self).__init__(methodName) - self.dataPath = os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1') - - def test_psfEll_(self): - iccd = 1 #[1...30] - psfSet = loadPSFSet(iccd, dataPath=self.dataPath) - - -if __name__ == '__main__': - unittest.main() - print('#####haha#####') - - - diff --git a/tests/UNIT_TEST_DATA/config_test.yaml b/tests/UNIT_TEST_DATA/config_test.yaml deleted file mode 100644 index 72ce1b23c3dc05e0c8df3d989ef2d15ce85babcd..0000000000000000000000000000000000000000 --- a/tests/UNIT_TEST_DATA/config_test.yaml +++ /dev/null @@ -1,153 +0,0 @@ ---- -############################################### -# -# Configuration file for CSST simulation -# Overall settings -# CSST-Sim Group, 2024/01/08 -# -############################################### - -# Base diretories and naming setup -# can add some of the command-line arguments here as well; -# ok to pass either way or both, as long as they are consistent -work_dir: "/public/home/chengliang/CSSOSDataProductsSims/csst-simulation/tests/UNIT_TEST_DATA/" -data_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/" -run_name: "testRun" - -# Project cycle and run counter are used to name the outputs -project_cycle: 9 -run_counter: 0 - -# Run options -run_option: - use_mpi: YES - # NOTE: "n_threads" paramters is currently not used in the backend - # simulation codes. It should be implemented later in the web frontend - # in order to config the number of threads to request from NAOC cluster - n_threads: 80 - - # Output catalog only? - # If yes, no imaging simulation will run - out_cat_only: NO - -############################################### -# Catalog setting -############################################### -# Configure the input catalog: options should be implemented -# in the corresponding (user defined) 'Catalog' class -catalog_options: - input_path: - # cat_dir: "Catalog_C6_20221212" - cat_dir: "" - star_cat: "starcat/" - galaxy_cat: "qsocat/cat2CSSTSim_bundle-50sqDeg/" - # AGN_cat: "AGN_C6_ross13_rand_pos_rmax-1.3.fits" - - SED_templates_path: - star_SED: "SpecLib.hdf5" - galaxy_SED: "sedlibs/" - AGN_SED: "qsocat/qsosed/" - # AGN_SED_WAVE: "wave_ross13.npy" - - # Only simulate stars? - star_only: NO - - # Only simulate galaxies? - galaxy_only: NO - - # rotate galaxy ellipticity - rotateEll: 0. # [degree] - - seed_Av: 121212 # Seed for generating random intrinsic extinction - -############################################### -# Observation setting -############################################### -obs_setting: - # (Optional) a file of point list - # if you just want to run default pointing: - # - pointing_dir: null - # - pointing_file: null - pointing_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/pointing_gir25/" - pointing_file: "pointing_50_1.dat" - - obs_config_file: "/public/home/chengliang/CSSOSDataProductsSims/csst-simulation/config/obs_config_SCI_WIDE_phot.yaml" - - # Run specific pointing(s): - # - give a list of indexes of pointings: [ip_1, ip_2...] - # - run all pointings: null - # Note: only valid when a pointing list is specified - run_pointings: [0] - - # Whether to enable astrometric modeling - enable_astrometric_model: True - - # Cut by saturation magnitude in which band? - cut_in_band: "z" - - # saturation magnitude margin - mag_sat_margin: -2.5 - # mag_sat_margin: -15. - - # limiting magnitude margin - mag_lim_margin: +1.0 - -############################################### -# PSF setting -############################################### -psf_setting: - - # Which PSF model to use: - # "Gauss": simple gaussian profile - # "Interp": Interpolated PSF from sampled ray-tracing data - psf_model: "Interp" - - # PSF size [arcseconds] - # radius of 80% energy encircled - # NOTE: only valid for "Gauss" PSF - # psf_rcont: 0.15 - - # path to PSF data - # NOTE: only valid for "Interp" PSF - # PSF models for photometry survey simulation - psf_pho_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/psfCube/set1_dynamic/" - # PSF models for slitless spectrum survey simulation - psf_sls_dir: "/public/share/yangxuliu/CSSOSDataProductsSims/data_50sqDeg/SLS_PSF_PCA_fp/" - -############################################### -# Shear setting -############################################### - -shear_setting: - # Options to generate mock shear field: - # "constant": all galaxies are assigned a constant reduced shear - # "catalog": get shear values from catalog - shear_type: "constant" - - # For constant shear field - reduced_g1: 0. - reduced_g2: 0. - -############################################### -# Output options -############################################### -output_setting: - output_format: "channels" # Whether to export as 16 channels (subimages) with pre- and over-scan ("image"/"channels") - shutter_output: OFF # Whether to export shutter effect 16-bit image - prnu_output: OFF # Whether to export the PRNU (pixel-to-pixel flat-fielding) files - -############################################### -# Random seeds -############################################### -random_seeds: - seed_poisson: 20210601 # Seed for Poisson noise - seed_CR: 20210317 # Seed for generating random cosmic ray maps - seed_flat: 20210101 # Seed for generating random flat fields - seed_prnu: 20210102 # Seed for photo-response non-uniformity - seed_gainNonUniform: 20210202 # Seed for gain nonuniformity - seed_biasNonUniform: 20210203 # Seed for bias nonuniformity - seed_rnNonUniform: 20210204 # Seed for readout-noise nonuniformity - seed_badcolumns: 20240309 # Seed for bad columns - seed_defective: 20210304 # Seed for defective (bad) pixels - seed_readout: 20210601 # Seed for read-out gaussian noise -... diff --git a/tests/UNIT_TEST_DATA/testCTE_image_before.fits b/tests/UNIT_TEST_DATA/testCTE_image_before.fits deleted file mode 100644 index b9defbff7fe01815c9f4f8594add481da9a86a7d..0000000000000000000000000000000000000000 Binary files a/tests/UNIT_TEST_DATA/testCTE_image_before.fits and /dev/null differ diff --git a/tests/test_prescan_overscan.py b/tests/test_prescan_overscan.py deleted file mode 100644 index 0e1291862fe6ca32e8809526dba2b29ff5cfa941..0000000000000000000000000000000000000000 --- a/tests/test_prescan_overscan.py +++ /dev/null @@ -1,63 +0,0 @@ -import unittest - -import sys,os,math -from itertools import islice -import numpy as np -import galsim -import yaml - -from ObservationSim.Instrument import Chip, Filter, FilterParam, FocalPlane -from ObservationSim.Instrument.Chip import ChipUtils as chip_utils - - -def defineCCD(iccd, config_file): - with open(config_file, "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) - chip = Chip(chipID=iccd, config=config) - chip.img = galsim.ImageF(chip.npix_x, chip.npix_y) - focal_plane = FocalPlane(chip_list=[iccd]) - chip.img.wcs= focal_plane.getTanWCS(192.8595, 27.1283, -113.4333*galsim.degrees, chip.pix_scale) - return chip - -def defineFilt(chip): - 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 - return bandpass_list - - -class detModule_coverage(unittest.TestCase): - def __init__(self, methodName='runTest'): - super(detModule_coverage, self).__init__(methodName) - self.dataPath = "/public/home/chengliang/CSSOSDataProductsSims/csst-simulation/tests/UNIT_TEST_DATA" ##os.path.join(os.getenv('UNIT_TEST_DATA_ROOT'), 'csst_fz_gc1') - self.iccd = 1 - - def test_add_prescan_overscan(self): - config_file = os.path.join(self.dataPath, 'config_test.yaml') - chip = defineCCD(self.iccd, config_file) - bandpass = defineFilt(chip) - print(chip.chipID) - print(chip.cen_pix_x, chip.cen_pix_y) - - chip.img = chip_utils.AddPreScan(GSImage=chip.img, - pre1=chip.prescan_x, - pre2=chip.prescan_y, - over1=chip.overscan_x, - over2=chip.overscan_y) - - self.assertTrue( (chip.prescan_x+chip.overscan_x)*8+chip.npix_x == np.shape(chip.img.array)[1] ) - self.assertTrue( (chip.prescan_y+chip.overscan_y)*2+chip.npix_y == np.shape(chip.img.array)[0] ) - - -if __name__ == '__main__': - unittest.main()