diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png index 19354e4..3c35b76 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png index af19c59..bc91f9d 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_IQU.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png index c71cef6..d1adfb7 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_I_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png index c82ba62..b589d76 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png index 0c60036..3cf814b 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_err.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png index 9f89e15..5995159 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_P_flux.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png index a2d9ebe..fe79fb0 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRi.png differ diff --git a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png index 1bb629d..fe3b8a5 100644 Binary files a/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png and b/plots/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_SNRp.png differ diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 451956a..fddb39d 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -130,7 +130,7 @@ def main(): if (data < 0.).any(): print("ETAPE 1 : ", data) # Crop data to remove outside blank margins. - data_array, error_array = proj_red.crop_array(data_array, headers, step=5, null_val=0., inside=True, display=display_crop, savename=figname, plots_folder=plots_folder) + data_array, error_array, headers = proj_red.crop_array(data_array, headers, step=5, null_val=0., inside=True, display=display_crop, savename=figname, plots_folder=plots_folder) for data in data_array: if (data < 0.).any(): print("ETAPE 2 : ", data) @@ -138,7 +138,7 @@ def main(): if deconvolve: data_array = proj_red.deconvolve_array(data_array, headers, psf=psf, FWHM=psf_FWHM, scale=psf_scale, shape=psf_shape, iterations=iterations) # Estimate error from data background, estimated from sub-image of desired sub_shape. - data_array, error_array = proj_red.get_error(data_array, sub_shape=error_sub_shape, display=display_error, headers=headers, savename=figname+"_errors", plots_folder=plots_folder) + data_array, error_array, headers = proj_red.get_error(data_array, headers, sub_shape=error_sub_shape, display=display_error, savename=figname+"_errors", plots_folder=plots_folder) for data in data_array: if (data < 0.).any(): print("ETAPE 3 : ", data) @@ -150,7 +150,7 @@ def main(): if (data < 0.).any(): print("ETAPE 4 : ", data) # Align and rescale images with oversampling. - data_array, error_array, data_mask = proj_red.align_data(data_array, headers, error_array, upsample_factor=int(Dxy.min()), ref_center=align_center, return_shifts=False) + data_array, error_array, headers, data_mask = proj_red.align_data(data_array, headers, error_array, upsample_factor=int(Dxy.min()), ref_center=align_center, return_shifts=False) for data in data_array: if (data < 0.).any(): print("ETAPE 5 : ", data) @@ -166,7 +166,7 @@ def main(): mrot = np.array([[np.cos(-alpha), -np.sin(-alpha)], [np.sin(-alpha), np.cos(-alpha)]]) rectangle[0:2] = np.dot(mrot, np.asarray(rectangle[0:2]))+np.array(data_array.shape[1:])/2 rectangle[4] = alpha - data_array, error_array, data_mask, headers = proj_red.rotate_data(data_array, error_array, data_mask, headers, -ref_header['orientat']) + data_array, error_array, headers, data_mask = proj_red.rotate_data(data_array, error_array, data_mask, headers, -ref_header['orientat']) for data in data_array: if (data < 0.).any(): print("ETAPE 6 : ", data) @@ -191,7 +191,7 @@ def main(): [np.sin(-alpha), np.cos(-alpha)]]) rectangle[0:2] = np.dot(mrot, np.asarray(rectangle[0:2]))+np.array(data_array.shape[1:])/2 rectangle[4] = alpha - I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, -ref_header['orientat'], SNRi_cut=None) + I_stokes, Q_stokes, U_stokes, Stokes_cov, headers, data_mask = proj_red.rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, -ref_header['orientat'], SNRi_cut=None) # Compute polarimetric parameters (polarization degree and angle). P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P = proj_red.compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 352812f..97229d7 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -216,6 +216,10 @@ def crop_array(data_array, headers, error_array=None, step=5, null_val=None, Returns: cropped_array : numpy.ndarray Array containing the observationnal data homogeneously cropped. + headers : header list + Updated headers associated with the images in data_array. + cropped_error : numpy.ndarray + Array containing the error on the observationnal data homogeneously cropped. """ if error_array is None: error_array = np.zeros(data_array.shape) @@ -282,11 +286,16 @@ def crop_array(data_array, headers, error_array=None, step=5, null_val=None, crop_array = np.zeros((data_array.shape[0],new_shape[0],new_shape[1])) crop_error_array = np.zeros((data_array.shape[0],new_shape[0],new_shape[1])) - for i,image in enumerate(data_array): #Put the image data in the cropped array + for i,image in enumerate(data_array): + #Put the image data in the cropped array 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(headers[i])) + curr_wcs.wcs.crpix = curr_wcs.wcs.crpix - np.array([v_array[0], v_array[2]]) + headers[i].update(curr_wcs.to_header()) - return crop_array, crop_error_array + return crop_array, crop_error_array, headers def deconvolve_array(data_array, headers, psf='gaussian', FWHM=1., scale='px', @@ -353,7 +362,7 @@ def deconvolve_array(data_array, headers, psf='gaussian', FWHM=1., scale='px', return deconv_array -def get_error(data_array, sub_shape=(15,15), display=False, headers=None, +def get_error(data_array, headers, sub_shape=(15,15), display=False, savename=None, plots_folder="", return_background=False): """ Look for sub-image of shape sub_shape that have the smallest integrated @@ -363,6 +372,8 @@ def get_error(data_array, sub_shape=(15,15), display=False, headers=None, Inputs: data_array : numpy.ndarray Array containing the data to study (2D float arrays). + headers : header list + Headers associated with the images in data_array. sub_shape : tuple, optional Shape of the sub-image to look for. Must be odd. Defaults to (15,15). @@ -370,10 +381,6 @@ def get_error(data_array, sub_shape=(15,15), display=False, headers=None, If True, data_array will be displayed with a rectangle around the sub-image selected for background computation. Defaults to False. - headers : header list, optional - Headers associated with the images in data_array. Will only be used if - display is True. - Defaults to None. savename : str, optional Name of the figure the map should be saved to. If None, the map won't be saved (only displayed). Only used if display is True. @@ -390,6 +397,8 @@ def get_error(data_array, sub_shape=(15,15), display=False, headers=None, Returns: data_array : numpy.ndarray Array containing the data to study minus the background. + headers : header list + Updated headers associated with the images in data_array. error_array : numpy.ndarray Array containing the background values associated to the images in data_array. @@ -399,7 +408,7 @@ def get_error(data_array, sub_shape=(15,15), display=False, headers=None, Only returned if return_background is True. """ # Crop out any null edges - data, error = crop_array(data_array, headers, step=5, null_val=0., inside=False) + data, error, headers = crop_array(data_array, headers, step=5, null_val=0., inside=False) sub_shape = np.array(sub_shape) # Make sub_shape of odd values @@ -449,7 +458,6 @@ def get_error(data_array, sub_shape=(15,15), display=False, headers=None, print(data_array[i]) if display: - convert_flux = headers[0]['photflam'] date_time = np.array([headers[i]['date-obs']+';'+headers[i]['time-obs'] for i in range(len(headers))]) @@ -497,9 +505,9 @@ def get_error(data_array, sub_shape=(15,15), display=False, headers=None, plt.show() if return_background: - return data_array, error_array, np.array([error_array[i][0,0] for i in range(error_array.shape[0])]) + return data_array, error_array, headers, np.array([error_array[i][0,0] for i in range(error_array.shape[0])]) else: - return data_array, error_array + return data_array, error_array, headers def rebin_array(data_array, error_array, headers, pxsize, scale, @@ -656,19 +664,21 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1., raise ValueError("All images in data_array must have same shape as\ ref_data") if error_array is None: - _, error_array, background = get_error(data_array, return_background=True) + _, error_array, headers, background = get_error(data_array, headers, return_background=True) else: - _, _, background = get_error(data_array, return_background=True) + _, _, headers, background = get_error(data_array, headers, return_background=True) # Crop out any null edges #(ref_data must be cropped as well) full_array = np.concatenate((data_array,[ref_data]),axis=0) + full_headers = deepcopy(headers) + full_headers.append(headers[0]) err_array = np.concatenate((error_array,[np.zeros(ref_data.shape)]),axis=0) - full_array, err_array = crop_array(full_array, headers, err_array, step=5, + full_array, err_array, full_headers = crop_array(full_array, full_headers, err_array, step=5, inside=False) - data_array, ref_data = full_array[:-1], full_array[-1] + data_array, ref_data, headers = full_array[:-1], full_array[-1], full_headers[:-1] error_array = err_array[:-1] if ref_center is None: @@ -723,14 +733,22 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1., shifts.append(shift) errors.append(error) - + shifts = np.array(shifts) errors = np.array(errors) + # Update headers CRPIX value + added_pix = (np.array(rescaled_image.shape[1:]) - np.array(data_array.shape[1:]))/2 + headers_wcs = [deepcopy(WCS(header)) for header in headers] + new_crpix = np.array([wcs.wcs.crpix for wcs in headers_wcs]) + shifts + added_pix + for i in range(len(headers_wcs)): + headers_wcs[i].wcs.crpix = new_crpix.mean(axis=0) + headers[i].update(headers_wcs[i].to_header()) + if return_shifts: - return rescaled_image, rescaled_error, rescaled_mask.any(axis=0), shifts, errors + return rescaled_image, rescaled_error, headers, rescaled_mask.any(axis=0), shifts, errors else: - return rescaled_image, rescaled_error, rescaled_mask.any(axis=0) + return rescaled_image, rescaled_error, headers, rescaled_mask.any(axis=0) def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., @@ -1316,6 +1334,8 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, new_Stokes_cov[i,j] = sc_rotate(new_Stokes_cov[i,j], ang, reshape=False, cval=0.) new_Stokes_cov[i,i] = np.abs(new_Stokes_cov[i,i]) + center = np.array(new_I_stokes.shape)/2 + print('c',center) #Update headers to new angle new_headers = [] @@ -1324,7 +1344,6 @@ 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() if new_wcs.wcs.has_cd(): # CD matrix # Update WCS with relevant information @@ -1344,6 +1363,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, elif new_wcs.wcs.has_pc(): # PC matrix + CDELT newpc = np.dot(mrot, new_wcs.wcs.get_pc()) new_wcs.wcs.pc = newpc + new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - center) + center new_wcs.wcs.set() new_header.update(new_wcs.to_header()) @@ -1359,9 +1379,9 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, new_U_stokes[np.isnan(new_U_stokes)] = 0. new_Stokes_cov[np.isnan(new_Stokes_cov)] = fmax - return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers + return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_headers, new_data_mask -def rotate_data(data_array, error_array, data_mask, headers, ang): +def rotate_data(data_array, error_array, headers, data_mask, ang): """ Use scipy.ndimage.rotate to rotate I_stokes to an angle, and a rotation matrix to rotate Q, U of a given angle in degrees and update header @@ -1437,4 +1457,4 @@ def rotate_data(data_array, error_array, data_mask, headers, ang): new_headers.append(new_header) - return new_data_array, new_error_array, new_data_mask, new_headers + return new_data_array, new_error_array, new_headers, new_data_mask