diff --git a/ObservationSim/Config/Header/ImageHeader.py b/ObservationSim/Config/Header/ImageHeader.py index 75229b718b932cc73623ddb0b50cc02ac3ec2ff9..0281dcddae182eafd299637289f9cf2c3b0ceead 100644 --- a/ObservationSim/Config/Header/ImageHeader.py +++ b/ObservationSim/Config/Header/ImageHeader.py @@ -87,6 +87,17 @@ def rotate_CD_matrix(cd, pa_aper): cd_rot = np.dot(mat, cd) return cd_rot +def calcaluteSLSRotSkyCoor(pix_xy = None,rot_angle = 1, xlen = 9216, ylen = 9232, w = None): + rad = np.deg2rad(rot_angle) + 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)]) + center = np.array([xlen/2, ylen/2]) + rot_pix = np.dot(mat, pix_xy-center) + center + skyCoor = w.wcs_pix2world(np.array([rot_pix]), 1) + + return skyCoor + # def Header_extention(xlen = 9216, ylen = 9232, gain = 1.0, readout = 5.0, dark = 0.02,saturation=90000, row_num = 1, col_num = 1): # @@ -153,7 +164,7 @@ def rotate_CD_matrix(cd, pa_aper): ##9232 9216 898 534 1309 60 -40 -23.4333 -def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1): +def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, ra = 60, dec = -40, pa = -23.433,psize = 0.074, row_num = 1, col_num = 1, filter = 'GI'): """ Creat a wcs frame for CCST with multiple extensions @@ -228,6 +239,63 @@ def WCS_def(xlen = 9216, ylen = 9232, gapy = 898.0, gapx1 = 534, gapx2 = 1309, r r_dat['CD1_2'] = cd_rot[0,1] r_dat['CD2_1'] = cd_rot[1,0] r_dat['CD2_2'] = cd_rot[1,1] + + if filter in ['GU', 'GV', 'GI']: + from astropy import wcs + + w = wcs.WCS(naxis=2) + w.wcs.crpix = [r_dat['CRPIX1'], r_dat['CRPIX2']] + w.wcs.cd = cd_rot + w.wcs.crval = [ra_ref, dec_ref] + w.wcs.ctype = [r_dat['CTYPE1'], r_dat['CTYPE2']] + + # test_center_o = w.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1) + + sls_rot = 1 + if i > 2: + sls_rot = -sls_rot + + sn_x = 30 + sn_y = 30 + x_pixs = np.zeros(sn_y * sn_x) + y_pixs = np.zeros(sn_y * sn_x) + xpixs_line = np.linspace(1, xlen, sn_x) + ypixs_line = np.linspace(1, ylen, sn_y) + + sky_coors = [] + + for n1, y in enumerate(ypixs_line): + for n2, x in enumerate(xpixs_line): + i_pix = n1 * sn_x + n2 + x_pixs[i_pix] = x + y_pixs[i_pix] = y + + pix_coor = np.array([x, y]) + sc1 = calcaluteSLSRotSkyCoor(pix_xy=pix_coor, rot_angle=sls_rot, w=w) + # print(sc1[0,0],sc1[0,1]) + sky_coors.append((sc1[0, 0], sc1[0, 1])) + + from astropy.coordinates import SkyCoord + from astropy.wcs.utils import fit_wcs_from_points + + wcs_new = fit_wcs_from_points(xy=np.array([x_pixs, y_pixs]), + world_coords=SkyCoord(sky_coors, frame="icrs", unit="deg"), projection='TAN') + + # print(wcs_new) + # test_center = wcs_new.wcs_pix2world(np.array([[xlen / 2, ylen / 2]]), 1) + # + # print(test_center - test_center_o) + + r_dat['CD1_1'] = wcs_new.wcs.cd[0, 0] + r_dat['CD1_2'] = wcs_new.wcs.cd[0, 1] + r_dat['CD2_1'] = wcs_new.wcs.cd[1, 0] + r_dat['CD2_2'] = wcs_new.wcs.cd[1, 1] + r_dat['CRPIX1'] = wcs_new.wcs.crpix[0] + r_dat['CRPIX2'] = wcs_new.wcs.crpix[1] + + r_dat['CRVAL1'] = wcs_new.wcs.crval[0] + r_dat['CRVAL2'] = wcs_new.wcs.crval[1] + return r_dat @@ -347,7 +415,7 @@ def generateExtensionHeader(xlen = 9216, ylen = 9232,ra = 60, dec = -40, pa = -2 h_ext['CCDCHIP'] = 'ccd' + CCDID[k-1].rjust(2,'0') h_ext['POS_ANG'] = pa header_wcs = WCS_def(xlen=xlen, ylen=ylen, gapy=898.0, gapx1=534, gapx2=1309, ra=ra, dec=dec, pa=pa, psize=psize, - row_num=row_num, col_num=col_num) + row_num=row_num, col_num=col_num, filter = h_ext['FILTER']) h_ext['CRPIX1'] = header_wcs['CRPIX1'] h_ext['CRPIX2'] = header_wcs['CRPIX2']