diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png index 0e361fd..01e9a11 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_IQU.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_IQU.png index afd23cf..358bc0f 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_IQU.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_IQU.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_I_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_I_err.png index 156ac05..3a692fe 100644 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_I_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_I_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png index af21fbb..bd50888 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png index 117f340..5bbf717 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_err.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png index c75c255..72b3538 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_P_flux.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png index 8407d69..96ee04f 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRi.png differ diff --git a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png index 8b925f8..1039a71 100755 Binary files a/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png and b/plots/NGC1068_x274020/NGC1068_FOC_combine_FWHM020_rot_SNRp.png differ diff --git a/src/lib/fits.py b/src/lib/fits.py index 6d0b6ec..25364ec 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -39,7 +39,7 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): with fits.open(data_folder+infiles[i]) as f: headers.append(f[0].header) data_array.append(f[0].data) - data_array = np.array(data_array) + data_array = np.array(data_array,dtype=np.double) # Prevent negative count value in imported data for i in range(len(data_array)): diff --git a/src/lib/plots.py b/src/lib/plots.py index ccb077d..659dfd4 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -327,8 +327,10 @@ def polarization_map(Stokes, data_mask, rectangle=None, SNRp_cut=3., SNRi_cut=30 SNRp = pol.data/pol_err.data SNRp[np.isnan(SNRp)] = 0. pol.data[SNRp < SNRp_cut] = np.nan - SNRi = stkI.data/np.sqrt(stk_cov.data[0,0]) - SNRi[np.isnan(SNRi)] = 0. + + maskI = stk_cov.data[0,0] > 0 + SNRi = np.zeros(stkI.data.shape) + SNRi[maskI] = stkI.data[maskI]/np.sqrt(stk_cov.data[0,0][maskI]) pol.data[SNRi < SNRi_cut] = np.nan data_mask = (1.-data_mask).astype(bool) @@ -373,7 +375,7 @@ def polarization_map(Stokes, data_mask, rectangle=None, SNRp_cut=3., SNRi_cut=30 # Display polarization degree error map vmin, vmax = 0., 10. p_err = pol_err.data.copy() - p_err[p_err*100. > vmax] = np.nan + p_err[p_err > vmax/100.] = np.nan im = ax.imshow(p_err*100.,extent=[-pol_err.data.shape[1]/2.,pol_err.data.shape[1]/2.,-pol_err.data.shape[0]/2.,pol_err.data.shape[0]/2.], vmin=vmin, vmax=vmax, aspect='auto', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]") elif display.lower() in ['s_i','i_err']: diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 35892a5..a09d0c1 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -46,6 +46,9 @@ from datetime import datetime from scipy.ndimage import rotate as sc_rotate from scipy.ndimage import shift as sc_shift from astropy.wcs import WCS +from astropy import log +log.setLevel('ERROR') +import warnings from lib.deconvolve import deconvolve_im, gaussian_psf from lib.convex_hull import image_hull from lib.plots import plot_obs @@ -826,7 +829,7 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., # Define gaussian stdev stdev = FWHM/(2.*np.sqrt(2.*np.log(2))) - fmax = np.finfo(np.float64).max + fmax = np.finfo(np.double).max if smoothing.lower() in ['combine','combining']: # Smooth using N images combination algorithm @@ -843,7 +846,9 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., for c in range(smoothed.shape[1]): # Compute distance from current pixel dist_rc = np.where(data_mask, fmax, np.sqrt((r-xx)**2+(c-yy)**2)) - g_rc = np.array([np.exp(-0.5*(dist_rc/stdev)**2),]*data_array.shape[0]) + # Catch expected "OverflowWarning" as we overflow values that are not in the image + with warnings.catch_warnings(record=True) as w: + g_rc = np.array([np.exp(-0.5*(dist_rc/stdev)**2),]*data_array.shape[0]) # Apply weighted combination smoothed[r,c] = (1.-data_mask[r,c])*np.sum(data_array*weight*g_rc)/np.sum(weight*g_rc) error[r,c] = np.sqrt(np.sum(weight*g_rc**2))/np.sum(weight*g_rc) @@ -862,7 +867,9 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1., for r in range(image.shape[0]): for c in range(image.shape[1]): dist_rc = np.where(data_mask, fmax, np.sqrt((r-xx)**2+(c-yy)**2)) - g_rc = np.exp(-0.5*(dist_rc/stdev)**2)/(2.*np.pi*stdev**2) + # Catch expected "OverflowWarning" as we overflow values that are not in the image + with warnings.catch_warnings(record=True) as w: + g_rc = np.exp(-0.5*(dist_rc/stdev)**2)/(2.*np.pi*stdev**2) smoothed[i][r,c] = (1.-data_mask[r,c])*np.sum(image*g_rc) error[i][r,c] = np.sqrt(np.sum(error_array[i]*g_rc**2)) @@ -1137,7 +1144,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers, mask = (Q_stokes**2 + U_stokes**2) > I_stokes**2 if mask.any(): - print("WARNING : I_pol > I_stokes : ", I_stokes[mask].size) + print("WARNING : found {0:d} pixels for which I_pol > I_stokes".format(I_stokes[mask].size)) #Stokes covariance matrix Stokes_cov = np.zeros((3,3,I_stokes.shape[0],I_stokes.shape[1])) @@ -1265,22 +1272,34 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): for the new orientation angle. """ #Polarization degree and angle computation - I_pol = np.sqrt(Q_stokes**2 + U_stokes**2) - P = I_pol/I_stokes - P[I_stokes <= 0.] = 0. - PA = (90./np.pi)*np.arctan2(U_stokes,Q_stokes) + mask = I_stokes>0. + I_pol = np.zeros(I_stokes.shape) + I_pol[mask] = np.sqrt(Q_stokes[mask]**2 + U_stokes[mask]**2) + P = np.zeros(I_stokes.shape) + P[mask] = I_pol[mask]/I_stokes[mask] + PA = np.zeros(I_stokes.shape) + PA[mask] = (90./np.pi)*np.arctan2(U_stokes[mask],Q_stokes[mask]) if (P>1).any(): - print("WARNING : found pixels for which P > 1", P[P>1.].size) - + print("WARNING : found {0:d} pixels for which P > 1".format(P[P>1.].size)) + #Associated errors - s_P = (1/I_stokes)*np.sqrt((Q_stokes**2*Stokes_cov[1,1] + U_stokes**2*Stokes_cov[2,2] + 2.*Q_stokes*U_stokes*Stokes_cov[1,2])/(Q_stokes**2 + U_stokes**2) + ((Q_stokes/I_stokes)**2 + (U_stokes/I_stokes)**2)*Stokes_cov[0,0] - 2.*(Q_stokes/I_stokes)*Stokes_cov[0,1] - 2.*(U_stokes/I_stokes)*Stokes_cov[0,2]) - s_PA = (90./(np.pi*(Q_stokes**2 + U_stokes**2)))*np.sqrt(U_stokes**2*Stokes_cov[1,1] + Q_stokes**2*Stokes_cov[2,2] - 2.*Q_stokes*U_stokes*Stokes_cov[1,2]) + fmax = np.finfo(np.float64).max + + s_P = np.ones(I_stokes.shape)*fmax + s_P[mask] = (1/I_stokes[mask])*np.sqrt((Q_stokes[mask]**2*Stokes_cov[1,1][mask] + U_stokes[mask]**2*Stokes_cov[2,2][mask] + 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask])/(Q_stokes[mask]**2 + U_stokes[mask]**2) + ((Q_stokes[mask]/I_stokes[mask])**2 + (U_stokes[mask]/I_stokes[mask])**2)*Stokes_cov[0,0][mask] - 2.*(Q_stokes[mask]/I_stokes[mask])*Stokes_cov[0,1][mask] - 2.*(U_stokes[mask]/I_stokes[mask])*Stokes_cov[0,2][mask]) + s_P[np.isnan(s_P)] = fmax + s_PA = np.ones(I_stokes.shape)*fmax + s_PA[mask] = (90./(np.pi*(Q_stokes[mask]**2 + U_stokes[mask]**2)))*np.sqrt(U_stokes[mask]**2*Stokes_cov[1,1][mask] + Q_stokes[mask]**2*Stokes_cov[2,2][mask] - 2.*Q_stokes[mask]*U_stokes[mask]*Stokes_cov[1,2][mask]) - debiased_P = np.sqrt(P**2 - s_P**2) + # Catch expected "OverflowWarning" as wrong pixel have an overflowing error + with warnings.catch_warnings(record=True) as w: + mask2 = P**2 >= s_P**2 + debiased_P = np.zeros(I_stokes.shape) + debiased_P[mask2] = np.sqrt(P[mask2]**2 - s_P[mask2]**2) if (debiased_P>1.).any(): - print("WARNING : found pixels for which debiased_P > 100%", debiased_P[debiased_P>1.].size) + print("WARNING : found {0:d} pixels for which debiased_P > 100%".format(debiased_P[debiased_P>1.].size)) #Compute the total exposure time so that #I_stokes*exp_tot = N_tot the total number of events @@ -1289,18 +1308,18 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers): N_obs = I_stokes*exp_tot #Errors on P, PA supposing Poisson noise - s_P_P = np.sqrt(2.)/np.sqrt(N_obs)*100. - s_PA_P = s_P_P/(2.*P)*180./np.pi + s_P_P = np.ones(I_stokes.shape)*fmax + s_P_P[mask] = np.sqrt(2.)/np.sqrt(N_obs[mask])*100. + s_PA_P = np.ones(I_stokes.shape)*fmax + s_PA_P[mask2] = s_P_P[mask2]/(2.*P[mask2])*180./np.pi - # Nan handling : - fmax = np.finfo(np.float64).max - - P[np.isnan(P)] = 0. - s_P[np.isnan(s_P)] = fmax - s_PA[np.isnan(s_PA)] = fmax - debiased_P[np.isnan(debiased_P)] = 0. - s_P_P[np.isnan(s_P_P)] = fmax - s_PA_P[np.isnan(s_PA_P)] = fmax +# # Nan handling : +# P[np.isnan(P)] = 0. +# s_P[np.isnan(s_P)] = fmax +# s_PA[np.isnan(s_PA)] = fmax +# debiased_P[np.isnan(debiased_P)] = 0. +# s_P_P[np.isnan(s_P_P)] = fmax +# s_PA_P[np.isnan(s_PA_P)] = fmax return P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P @@ -1385,6 +1404,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, for j in range(3): 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]) #Update headers to new angle new_headers = []