From 42ca956b25d58e2c0a994456295fe8b656c3efa6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 12 Apr 2024 14:41:48 +0200 Subject: [PATCH] correct wcs handling for celestial wcs --- src/lib/fits.py | 8 ++++---- src/lib/plots.py | 14 +++++++------- src/lib/reduction.py | 16 ++++++++-------- 3 files changed, 19 insertions(+), 19 deletions(-) diff --git a/src/lib/fits.py b/src/lib/fits.py index 2de034b..9116c83 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -50,13 +50,13 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): # force WCS to convention PCi_ja unitary, cdelt in deg for header in headers: - new_wcs = WCS(header).deepcopy() + new_wcs = WCS(header).celestial.deepcopy() if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt[:2] == np.array([1., 1.])).all(): # Update WCS with relevant information if new_wcs.wcs.has_cd(): - old_cd = new_wcs.wcs.cd[:2, :2] + old_cd = new_wcs.wcs.cd del new_wcs.wcs.cd - keys = list(new_wcs.to_header().keys())+['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2'] + keys = list(new_wcs.to_header().keys())+['CD1_1', 'CD1_2', 'CD1_3', 'CD2_1', 'CD2_2', 'CD2_3', 'CD3_1', 'CD3_2', 'CD3_3'] for key in keys: header.remove(key, ignore_missing=True) new_cdelt = np.linalg.eig(old_cd)[0] @@ -71,7 +71,7 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): # force WCS for POL60 to have same pixel size as POL0 and POL120 is_pol60 = np.array([head['filtnam1'].lower() == 'pol60' for head in headers], dtype=bool) - cdelt = np.round(np.array([WCS(head).wcs.cdelt for head in headers]), 14) + cdelt = np.round(np.array([WCS(head).wcs.cdelt[:2] for head in headers]), 14) if np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2: print(np.unique(cdelt[np.logical_not(is_pol60)], axis=0)) raise ValueError("Not all images have same pixel size") diff --git a/src/lib/plots.py b/src/lib/plots.py index d305784..93a931e 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -471,13 +471,13 @@ class align_maps(object): self.map_data = fits.getdata(self.map_path) self.other_data = fits.getdata(self.other_path) - self.map_wcs = deepcopy(WCS(self.map_header)).celestial + self.map_wcs = WCS(self.map_header).celestial.deepcopy() if len(self.map_data.shape) == 4: self.map_data = self.map_data[0, 0] elif len(self.map_data.shape) == 3: self.map_data = self.map_data[0] - self.other_wcs = deepcopy(WCS(self.other_header)).celestial + self.other_wcs = WCS(self.other_header).celestial.deepcopy() if len(self.other_data.shape) == 4: self.other_data = self.other_data[0, 0] elif len(self.other_data.shape) == 3: @@ -797,7 +797,7 @@ class overplot_chandra(align_maps): pang = self.Stokes_UV['POL_ANG'].data other_data = deepcopy(self.other_data) - other_wcs = deepcopy(self.other_wcs) + other_wcs = self.other_wcs.deepcopy() if zoom != 1: other_data = sc_zoom(other_data, zoom) other_wcs.wcs.crpix *= zoom @@ -1044,8 +1044,8 @@ class align_pol(object): maps = np.array(maps)[order] self.ref_map, self.other_maps = maps[0], maps[1:] - self.wcs = WCS(self.ref_map[0].header) - self.wcs_other = np.array([WCS(map[0].header) for map in self.other_maps]) + self.wcs = WCS(self.ref_map[0].header).celestial.deepcopy() + self.wcs_other = np.array([WCS(map[0].header).celestial.deepcopy() for map in self.other_maps]) self.aligned = np.zeros(self.other_maps.shape[0], dtype=bool) @@ -1175,7 +1175,7 @@ class crop_map(object): self.cropped = False self.hdul = hdul self.header = deepcopy(self.hdul[0].header) - self.wcs = WCS(self.header).deepcopy() + self.wcs = WCS(self.header).celestial.deepcopy() self.data = deepcopy(self.hdul[0].data) try: @@ -2179,7 +2179,7 @@ class pol_map(object): @property def wcs(self): - return WCS(self.Stokes[0].header).celestial + return WCS(self.Stokes[0].header).celestial.deepcopy() @property def I(self): diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 582b0a4..b91b3bd 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -269,7 +269,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu crop_array[i] = image[v_array[0]:v_array[1], v_array[2]:v_array[3]] crop_error_array[i] = error_array[i][v_array[0]:v_array[1], v_array[2]:v_array[3]] # Update CRPIX value in the associated header - curr_wcs = deepcopy(WCS(crop_headers[i])) + curr_wcs = WCS(crop_headers[i]).celestial.deepcopy() curr_wcs.wcs.crpix[:2] = curr_wcs.wcs.crpix[:2] - np.array([v_array[2], v_array[0]]) crop_headers[i].update(curr_wcs.to_header()) crop_headers[i]['naxis1'], crop_headers[i]['naxis2'] = crop_array[i].shape @@ -371,7 +371,7 @@ def deconvolve_array(data_array, headers, psf='gaussian', FWHM=1., scale='px', pxsize = np.zeros((data_array.shape[0], 2)) for i, header in enumerate(headers): # Get current pixel size - w = WCS(header).deepcopy() + w = WCS(header).celestial.deepcopy() pxsize[i] = np.round(w.wcs.cdelt/3600., 15) if (pxsize != pxsize[0]).any(): raise ValueError("Not all images in array have same pixel size") @@ -559,7 +559,7 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, operation='sum' Dxy_arr = np.ones((data_array.shape[0], 2)) for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): # Get current pixel size - w = WCS(header).deepcopy() + w = WCS(header).celestial.deepcopy() new_header = deepcopy(header) # Compute binning ratio @@ -575,7 +575,7 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, operation='sum' for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): # Get current pixel size - w = WCS(header).deepcopy() + w = WCS(header).celestial.deepcopy() new_header = deepcopy(header) Dxy = image.shape/new_shape @@ -759,7 +759,7 @@ def align_data(data_array, headers, error_array=None, background=None, upsample_ errors = np.array(errors) # Update headers CRPIX value - headers_wcs = [deepcopy(WCS(header)) for header in headers] + headers_wcs = [WCS(header).celestial.deepcopy() for header in headers] new_crpix = np.array([wcs.wcs.crpix for wcs in headers_wcs]) + shifts[:, ::-1] + res_shift[::-1] for i in range(len(headers_wcs)): headers_wcs[i].wcs.crpix = new_crpix[0] @@ -814,7 +814,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., scale='pix pxsize = np.zeros((data_array.shape[0], 2)) for i, header in enumerate(headers): # Get current pixel size - w = WCS(header).deepcopy() + w = WCS(header).celestial.deepcopy() pxsize[i] = np.round(w.wcs.cdelt*3600., 4) if (pxsize != pxsize[0]).any(): raise ValueError("Not all images in array have same pixel size") @@ -1447,7 +1447,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, for header in headers: new_header = deepcopy(header) new_header['orientat'] = header['orientat'] + ang - new_wcs = WCS(header).deepcopy() + new_wcs = WCS(header).celestial.deepcopy() new_wcs.wcs.pc = np.dot(mrot, new_wcs.wcs.pc) new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center[::-1]) + new_center[::-1] @@ -1560,7 +1560,7 @@ def rotate_data(data_array, error_array, data_mask, headers, ang): new_header = deepcopy(header) new_header['orientat'] = header['orientat'] + ang - new_wcs = WCS(header).deepcopy() + new_wcs = WCS(header).celestial.deepcopy() new_wcs.wcs.pc[:2, :2] = np.dot(mrot, new_wcs.wcs.pc[:2, :2]) new_wcs.wcs.crpix[:2] = np.dot(mrot, new_wcs.wcs.crpix[:2] - old_center[::-1]) + new_center[::-1]