Tibeuleu/FOC_Reduction
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

correct wcs handling for celestial wcs

Browse Source
This commit is contained in:
Thibault Barnouin
2024-04-12 14:41:48 +02:00
parent ce927c3409
commit 42ca956b25
3 changed files with 19 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
src/lib/fits.py
View File

@@ -50,13 +50,13 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
# force WCS to convention PCi_ja unitary, cdelt in deg # force WCS to convention PCi_ja unitary, cdelt in deg
for header in headers: 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(): if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt[:2] == np.array([1., 1.])).all():
# Update WCS with relevant information # Update WCS with relevant information
if new_wcs.wcs.has_cd(): 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 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: for key in keys:
header.remove(key, ignore_missing=True) header.remove(key, ignore_missing=True)
new_cdelt = np.linalg.eig(old_cd)[0] 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 # 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) 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: if np.unique(cdelt[np.logical_not(is_pol60)], axis=0).size != 2:
print(np.unique(cdelt[np.logical_not(is_pol60)], axis=0)) print(np.unique(cdelt[np.logical_not(is_pol60)], axis=0))
raise ValueError("Not all images have same pixel size") raise ValueError("Not all images have same pixel size")

14
src/lib/plots.py
View File

@@ -471,13 +471,13 @@ class align_maps(object):
self.map_data = fits.getdata(self.map_path) self.map_data = fits.getdata(self.map_path)
self.other_data = fits.getdata(self.other_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: if len(self.map_data.shape) == 4:
self.map_data = self.map_data[0, 0] self.map_data = self.map_data[0, 0]
elif len(self.map_data.shape) == 3: elif len(self.map_data.shape) == 3:
self.map_data = self.map_data[0] 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: if len(self.other_data.shape) == 4:
self.other_data = self.other_data[0, 0] self.other_data = self.other_data[0, 0]
elif len(self.other_data.shape) == 3: elif len(self.other_data.shape) == 3:
@@ -797,7 +797,7 @@ class overplot_chandra(align_maps):
pang = self.Stokes_UV['POL_ANG'].data pang = self.Stokes_UV['POL_ANG'].data
other_data = deepcopy(self.other_data) other_data = deepcopy(self.other_data)
other_wcs = deepcopy(self.other_wcs) other_wcs = self.other_wcs.deepcopy()
if zoom != 1: if zoom != 1:
other_data = sc_zoom(other_data, zoom) other_data = sc_zoom(other_data, zoom)
other_wcs.wcs.crpix *= zoom other_wcs.wcs.crpix *= zoom
@@ -1044,8 +1044,8 @@ class align_pol(object):
maps = np.array(maps)[order] maps = np.array(maps)[order]
self.ref_map, self.other_maps = maps[0], maps[1:] self.ref_map, self.other_maps = maps[0], maps[1:]
self.wcs = WCS(self.ref_map[0].header) self.wcs = WCS(self.ref_map[0].header).celestial.deepcopy()
self.wcs_other = np.array([WCS(map[0].header) for map in self.other_maps]) 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) self.aligned = np.zeros(self.other_maps.shape[0], dtype=bool)
@@ -1175,7 +1175,7 @@ class crop_map(object):
self.cropped = False self.cropped = False
self.hdul = hdul self.hdul = hdul
self.header = deepcopy(self.hdul[0].header) 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) self.data = deepcopy(self.hdul[0].data)
try: try:
@@ -2179,7 +2179,7 @@ class pol_map(object):
@property @property
def wcs(self): def wcs(self):
return WCS(self.Stokes[0].header).celestial return WCS(self.Stokes[0].header).celestial.deepcopy()
@property @property
def I(self): def I(self):

16
src/lib/reduction.py
View File

@@ -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_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]] 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 # 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]]) 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].update(curr_wcs.to_header())
crop_headers[i]['naxis1'], crop_headers[i]['naxis2'] = crop_array[i].shape 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)) pxsize = np.zeros((data_array.shape[0], 2))
for i, header in enumerate(headers): for i, header in enumerate(headers):
# Get current pixel size # Get current pixel size
w = WCS(header).deepcopy() w = WCS(header).celestial.deepcopy()
pxsize[i] = np.round(w.wcs.cdelt/3600., 15) pxsize[i] = np.round(w.wcs.cdelt/3600., 15)
if (pxsize != pxsize[0]).any(): if (pxsize != pxsize[0]).any():
raise ValueError("Not all images in array have same pixel size") 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)) Dxy_arr = np.ones((data_array.shape[0], 2))
for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))):
# Get current pixel size # Get current pixel size
w = WCS(header).deepcopy() w = WCS(header).celestial.deepcopy()
new_header = deepcopy(header) new_header = deepcopy(header)
# Compute binning ratio # 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))): for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))):
# Get current pixel size # Get current pixel size
w = WCS(header).deepcopy() w = WCS(header).celestial.deepcopy()
new_header = deepcopy(header) new_header = deepcopy(header)
Dxy = image.shape/new_shape 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) errors = np.array(errors)
# Update headers CRPIX value # 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] new_crpix = np.array([wcs.wcs.crpix for wcs in headers_wcs]) + shifts[:, ::-1] + res_shift[::-1]
for i in range(len(headers_wcs)): for i in range(len(headers_wcs)):
headers_wcs[i].wcs.crpix = new_crpix[0] 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)) pxsize = np.zeros((data_array.shape[0], 2))
for i, header in enumerate(headers): for i, header in enumerate(headers):
# Get current pixel size # Get current pixel size
w = WCS(header).deepcopy() w = WCS(header).celestial.deepcopy()
pxsize[i] = np.round(w.wcs.cdelt*3600., 4) pxsize[i] = np.round(w.wcs.cdelt*3600., 4)
if (pxsize != pxsize[0]).any(): if (pxsize != pxsize[0]).any():
raise ValueError("Not all images in array have same pixel size") 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: for header in headers:
new_header = deepcopy(header) new_header = deepcopy(header)
new_header['orientat'] = header['orientat'] + ang 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.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] 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 = deepcopy(header)
new_header['orientat'] = header['orientat'] + ang 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.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] new_wcs.wcs.crpix[:2] = np.dot(mrot, new_wcs.wcs.crpix[:2] - old_center[::-1]) + new_center[::-1]