diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index e765147..99b44eb 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -46,8 +46,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): display_data = False # Smoothing - smoothing_function = 'combine' #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 0.20 #If None, no smoothing is done + smoothing_function = 'gaussian' #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine + smoothing_FWHM = None #If None, no smoothing is done smoothing_scale = 'arcsec' #pixel or arcsec # Rotation @@ -77,8 +77,8 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): target = input("Target name:\n>") else: target, products = retrieve_products(target,proposal_id,output_dir=output_dir) - prod = products.pop() - for prods in products: + prod = products[0] + for prods in products[1:]: main(target=target,infiles=["/".join(pr) for pr in prods],output_dir=output_dir) data_folder = prod[0,0] try: @@ -89,7 +89,15 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): data_array, headers = proj_fits.get_obs_data(infiles, data_folder=data_folder, compute_flux=True) figname = "_".join([target,"FOC"]) - figtype = "_".join(["".join([s[0] for s in smoothing_function.split("_")]),"{0:.2f}".format(smoothing_FWHM).replace(".","")]) #additionnal informations + if smoothing_FWHM is None: + if px_scale in ['px','pixel','pixels']: + figtype = "".join(["b_",str(pxsize),'px']) + elif px_scale in ['arcsec','arcseconds','arcs']: + figtype = "".join(["b_","{0:.2f}".format(pxsize).replace(".",""),'arcsec']) + else: + figtype = "full" + else: + figtype = "_".join(["".join([s[0] for s in smoothing_function.split("_")]),"".join(["{0:.2f}".format(smoothing_FWHM).replace(".",""),smoothing_scale])]) #additionnal informations # Crop data to remove outside blank margins. 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) @@ -100,7 +108,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): # Estimate error from data background, estimated from sub-image of desired sub_shape. background = None - data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_error, savename="_".join([figname,"_errors"]), plots_folder=plots_folder, return_background=True) + data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_error, savename="_".join([figname,"errors"]), plots_folder=plots_folder, return_background=True) # Align and rescale images with oversampling. data_array, error_array, headers, data_mask = proj_red.align_data(data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=False) diff --git a/src/comparison_Kishimoto.py b/src/comparison_Kishimoto.py index b81edb5..a81982f 100755 --- a/src/comparison_Kishimoto.py +++ b/src/comparison_Kishimoto.py @@ -16,8 +16,8 @@ import matplotlib.pyplot as plt root_dir = path_join('/home/t.barnouin/Documents/Thesis/HST') root_dir_K = path_join(root_dir,'Kishimoto','output') root_dir_S = path_join(root_dir,'FOC_Reduction','output') -root_dir_data_S = path_join(root_dir,'FOC_Reduction','data','NGC1068_x274020') -root_dir_plot_S = path_join(root_dir,'FOC_Reduction','plots','NGC1068_x274020') +root_dir_data_S = path_join(root_dir,'FOC_Reduction','data','NGC1068','5144') +root_dir_plot_S = path_join(root_dir,'FOC_Reduction','plots','NGC1068','5144') filename_S = "NGC1068_FOC_b_10px.fits" plt.rcParams.update({'font.size': 15}) diff --git a/src/lib/background.py b/src/lib/background.py old mode 100644 new mode 100755 index 1c175b9..c94be17 --- a/src/lib/background.py +++ b/src/lib/background.py @@ -119,7 +119,7 @@ def sky_part(img): # Intensity range sky_med = np.median(rand_pix) sig = np.min([img[imgsky_med].std()]) - sky_range = [sky_med-2.*sig, sky_med+sig] + sky_range = [sky_med-2.*sig, np.max([sky_med+sig,7e-4])] #Detector background average FOC Data Handbook Sec. 7.6 sky = img[np.logical_and(img>=sky_range[0],img<=sky_range[1])] return sky, sky_range @@ -170,7 +170,7 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save 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. - Defaults to None. + Defaults to None.CNRS-Unistra Labo ObsAstroS plots_folder : str, optional Relative (or absolute) filepath to the folder in wich the map will be saved. Not used if savename is None. @@ -212,18 +212,7 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save error_bkg[i] *= bkg - # Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999) - #wavelength dependence of the polariser filters - #estimated to less than 1% - err_wav = data[i]*0.01 - #difference in PSFs through each polarizers - #estimated to less than 3% - err_psf = data[i]*0.03 - #flatfielding uncertainties - #estimated to less than 3% - err_flat = data[i]*0.03 - - n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2) + n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2) #Substract background if subtract_error>0: @@ -327,18 +316,7 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis error_bkg[i] *= bkg - # Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999) - #wavelength dependence of the polariser filters - #estimated to less than 1% - err_wav = data[i]*0.01 - #difference in PSFs through each polarizers - #estimated to less than 3% - err_psf = data[i]*0.03 - #flatfielding uncertainties - #estimated to less than 3% - err_flat = data[i]*0.03 - - n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2) + n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2) #Substract background if subtract_error > 0: @@ -435,18 +413,7 @@ def bkg_mini(data, error, mask, headers, sub_shape=(15,15), subtract_error=True, bkg = np.sqrt(np.sum(sub_image**2)/sub_image.size)*subtract_error if subtract_error>0 else np.sqrt(np.sum(sub_image**2)/sub_image.size) error_bkg[i] *= bkg - # Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999) - #wavelength dependence of the polariser filters - #estimated to less than 1% - err_wav = data[i]*0.01 - #difference in PSFs through each polarizers - #estimated to less than 3% - err_psf = data[i]*0.03 - #flatfielding uncertainties - #estimated to less than 3% - err_flat = data[i]*0.03 - - n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2) + n_error_array[i] = np.sqrt(n_error_array[i]**2 + error_bkg[i]**2) #Substract background if subtract_error>0.: diff --git a/src/lib/fits.py b/src/lib/fits.py index 97fc9ec..5b96ba9 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -1,3 +1,5 @@ +#!/usr/bin/python3 +#-*- coding:utf-8 -*- """ Library function for simplified fits handling. @@ -12,7 +14,7 @@ prototypes : import numpy as np from os.path import join as path_join from astropy.io import fits -from astropy import wcs +from astropy.wcs import WCS from lib.convex_hull import image_hull, clean_ROI from lib.plots import princ_angle import matplotlib.pyplot as plt @@ -51,7 +53,7 @@ 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.WCS(header).deepcopy() + new_wcs = WCS(header).deepcopy() if new_wcs.wcs.has_cd() or (new_wcs.wcs.cdelt == np.array([1., 1.])).all(): # Update WCS with relevant information if new_wcs.wcs.has_cd(): @@ -70,6 +72,16 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): header[key] = val header['orientat'] = princ_angle(float(header['orientat'])) + # 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.array([WCS(head).wcs.cdelt for head in headers]) + if np.unique(cdelt[np.logical_not(is_pol60)],axis=0).size!=2: + print(np.unique(cdelt[np.logical_not(is_pol60)],axis=0).size) + raise ValueError("Not all images have same pixel size") + else: + for head in np.array(headers,dtype=object)[is_pol60]: + head['cdelt1'],head['cdelt2'] = np.unique(cdelt[np.logical_not(is_pol60)],axis=0)[0] + if compute_flux: for i in range(len(infiles)): # Compute the flux in counts/sec @@ -118,7 +130,7 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P, #Create new WCS object given the modified images ref_header = headers[0] exp_tot = np.array([header['exptime'] for header in headers]).sum() - new_wcs = wcs.WCS(ref_header).deepcopy() + new_wcs = WCS(ref_header).deepcopy() if data_mask.shape != (1,1): vertex = clean_ROI(data_mask) diff --git a/src/lib/plots.py b/src/lib/plots.py index 89621f9..4ecd430 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -287,20 +287,18 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder) #Compute SNR and apply cuts - pol.data[pol.data == 0.] = np.nan - pol_err.data[pol_err.data == 0.] = np.nan - SNRp = pol.data/pol_err.data - SNRp[np.isnan(SNRp)] = 0. - pol.data[SNRp < SNRp_cut] = np.nan - pang.data[SNRp < SNRp_cut] = np.nan + poldata, pangdata = pol.data.copy(), pang.data.copy() + maskP = pol_err.data > 0 + SNRp = np.zeros(pol.data.shape) + SNRp[maskP] = pol.data[maskP]/pol_err.data[maskP] 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 - pang.data[SNRi < SNRi_cut] = np.nan mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) + poldata[np.logical_not(mask)] = np.nan + pangdata[np.logical_not(mask)] = np.nan # Look for pixel of max polarization if np.isfinite(pol.data).any(): @@ -321,7 +319,10 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c if display.lower() in ['intensity']: # If no display selected, show intensity map display='i' - vmin, vmax = 1/5.0*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + if mask.sum() > 0.: + vmin, vmax = 1/5.0*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + else: + vmin, vmax = 1/5.0*np.mean(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsI = np.linspace(vmax*0.01, vmax*0.99, 10) @@ -332,7 +333,10 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c # Display polarisation flux display='pf' pf_mask = (stkI.data > 0.) * (pol.data > 0.) - vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + if mask.sum() > 0.: + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + else: + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) im = ax.imshow(stkI.data*convert_flux*pol.data, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") levelsPf = np.linspace(vmax*0.01, vmax*0.99, 10) @@ -354,40 +358,55 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c elif display.lower() in ['s_p','pol_err','pol_deg_err']: # Display polarization degree error map display='s_p' - vmin, vmax = 0., np.max(pol_err.data[SNRp > SNRp_cut])*100. - p_err = deepcopy(pol_err.data) - p_err[p_err > vmax/100.] = np.nan - im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + if (SNRp>SNRp_cut).any(): + vmin, vmax = 0., np.max(pol_err.data[SNRp > SNRp_cut])*100. + p_err = deepcopy(pol_err.data) + p_err[p_err > vmax/100.] = np.nan + im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + else: + im = ax.imshow(pol_err.data*100., aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_P$ [%]") elif display.lower() in ['s_i','i_err']: # Display intensity error map display='s_i' - vmin, vmax = np.min(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux), np.max(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux) - im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + if (SNRi>SNRi_cut).any(): + vmin, vmax = np.min(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux), np.max(np.sqrt(stk_cov.data[0,0][stk_cov.data[0,0] > 0.])*convert_flux) + im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + else: + im = ax.imshow(np.sqrt(stk_cov.data[0,0])*convert_flux, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\sigma_I$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") elif display.lower() in ['snr','snri']: # Display I_stokes signal-to-noise map display='snri' vmin, vmax = 0., np.max(SNRi[np.isfinite(SNRi)]) - im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + if vmax*0.99 > SNRi_cut: + im = ax.imshow(SNRi, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + levelsSNRi = np.linspace(SNRi_cut, vmax*0.99, 10) + print("SNRi contour levels : ", levelsSNRi) + cont = ax.contour(SNRi, levels=levelsSNRi, colors='grey', linewidths=0.5) + #ax.clabel(cont,inline=True,fontsize=6) + else: + im = ax.imshow(SNRi, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$I_{Stokes}/\sigma_{I}$") - levelsSNRi = np.linspace(SNRi_cut, vmax*0.99, 10) - print("SNRi contour levels : ", levelsSNRi) - cont = ax.contour(SNRi, levels=levelsSNRi, colors='grey', linewidths=0.5) - #ax.clabel(cont,inline=True,fontsize=6) elif display.lower() in ['snrp']: # Display polarization degree signal-to-noise map display='snrp' vmin, vmax = 0., np.max(SNRp[np.isfinite(SNRp)]) - im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + if vmax*0.99 > SNRp_cut: + im = ax.imshow(SNRp, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + levelsSNRp = np.linspace(SNRp_cut, vmax*0.99, 10) + print("SNRp contour levels : ", levelsSNRp) + cont = ax.contour(SNRp, levels=levelsSNRp, colors='grey', linewidths=0.5) + #ax.clabel(cont,inline=True,fontsize=6) + else: + im = ax.imshow(SNRp, aspect='equal', cmap='inferno', alpha=1.) cbar = plt.colorbar(im, cax=cbar_ax, label=r"$P/\sigma_{P}$") - levelsSNRp = np.linspace(SNRp_cut, vmax*0.99, 10) - print("SNRp contour levels : ", levelsSNRp) - cont = ax.contour(SNRp, levels=levelsSNRp, colors='grey', linewidths=0.5) - #ax.clabel(cont,inline=True,fontsize=6) else: # Defaults to intensity map - vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + if mask.sum() > 0.: + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0,0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + else: + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov.data[0,0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) #im = ax.imshow(stkI.data*convert_flux, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) #cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]") im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin,vmax), aspect='equal', cmap='inferno', alpha=1.) @@ -409,11 +428,11 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c if display.lower() in ['i','s_i','snri','pf','p','pa','s_p','snrp']: if step_vec == 0: - pol.data[np.isfinite(pol.data)] = 1./2. + poldata[np.isfinite(poldata)] = 1./2. step_vec = 1 vec_scale = 2. X, Y = np.meshgrid(np.arange(stkI.data.shape[1]), np.arange(stkI.data.shape[0])) - U, V = pol.data*np.cos(np.pi/2.+pang.data*np.pi/180.), pol.data*np.sin(np.pi/2.+pang.data*np.pi/180.) + U, V = poldata*np.cos(np.pi/2.+pangdata*np.pi/180.), poldata*np.sin(np.pi/2.+pangdata*np.pi/180.) Q = ax.quiver(X[::step_vec,::step_vec],Y[::step_vec,::step_vec],U[::step_vec,::step_vec],V[::step_vec,::step_vec],units='xy',angles='uv',scale=1./vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,linewidth=0.5,color='w',edgecolor='k') pol_sc = AnchoredSizeBar(ax.transData, vec_scale, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') @@ -1418,7 +1437,7 @@ class pol_map(object): b_aper.label.set_fontsize(8) b_aper_reset = Button(ax_aper_reset,"Reset") b_aper_reset.label.set_fontsize(8) - s_aper_radius = Slider(ax_aper_radius, r"$R_{aper}$", 0.5, 3.5, valstep=0.1, valinit=1) + s_aper_radius = Slider(ax_aper_radius, r"$R_{aper}$", np.ceil(self.wcs.wcs.cdelt.max()/1.33*3.6e5)/1e2, 3.5, valstep=1e-2, valinit=1.) def select_aperture(event): if self.data is None: @@ -1456,6 +1475,7 @@ class pol_map(object): def reset_aperture(event): self.region = None + s_aper_radius.reset() self.pol_int() self.fig.canvas.draw_idle() diff --git a/src/lib/query.py b/src/lib/query.py index ef44269..8ef2326 100755 --- a/src/lib/query.py +++ b/src/lib/query.py @@ -132,7 +132,7 @@ def get_product_list(target=None, proposal_id=None): products['target_name'] = Column(observations['target_name']) for prod in products: - products['proposal_id'] = results['Proposal ID'][results['Dataset']==prod['productFilename'][:len(results['Dataset'][0])].upper()] + prod['proposal_id'] = results['Proposal ID'][results['Dataset']==prod['productFilename'][:len(results['Dataset'][0])].upper()][0] #for prod in products: # prod['target_name'] = observations['target_name'][observation['obsid']==prod['obsID']] @@ -169,7 +169,7 @@ def retrieve_products(target=None, proposal_id=None, output_dir='./data'): filepaths.append([obs_dir,file]) prodpaths.append(np.array(filepaths,dtype=str)) - return target, prodpaths + return target, np.array(prodpaths) if __name__ == "__main__": diff --git a/src/lib/reduction.py b/src/lib/reduction.py index fe9043c..83c4096 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -485,17 +485,30 @@ def get_error(data_array, headers, error_array=None, data_mask=None, mask = zeropad(mask_c, data[0].shape).astype(bool) background = np.zeros((data.shape[0])) + #wavelength dependence of the polariser filters + #estimated to less than 1% + err_wav = data*0.01 + #difference in PSFs through each polarizers + #estimated to less than 3% + err_psf = data*0.03 + #flatfielding uncertainties + #estimated to less than 3% + err_flat = data*0.03 + if (sub_type is None): - n_data_array, n_error_array, headers, background = bkg_hist(data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) + n_data_array, c_error_bkg, headers, background = bkg_hist(data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) elif type(sub_type)==str: if sub_type.lower() in ['auto']: - n_data_array, n_error_array, headers, background = bkg_fit(data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) + n_data_array, c_error_bkg, headers, background = bkg_fit(data, error, mask, headers, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) else: - n_data_array, n_error_array, headers, background = bkg_hist(data, error, mask, headers, sub_type=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) + n_data_array, c_error_bkg, headers, background = bkg_hist(data, error, mask, headers, sub_type=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) elif type(sub_type)==tuple: - n_data_array, n_error_array, headers, background = bkg_mini(data, error, mask, headers, sub_shape=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) + n_data_array, c_error_bkg, headers, background = bkg_mini(data, error, mask, headers, sub_shape=sub_type, subtract_error=subtract_error, display=display, savename=savename, plots_folder=plots_folder) else: print("Warning: Invalid subtype.") + + # Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999) + n_error_array = np.sqrt(err_wav**2+err_psf**2+err_flat**2+c_error_bkg**2) if return_background: return n_data_array, n_error_array, headers, background @@ -553,29 +566,39 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, (yet)".format(instr)) rebinned_data, rebinned_error, rebinned_headers = [], [], [] - Dxy = np.array([1, 1],dtype=int) + Dxy = np.array([1., 1.]) # Routine for the FOC instrument if instr == 'FOC': HST_aper = 2400. # HST aperture in mm + Dxy_arr = np.ones((data_array.shape[0],2)) for i, enum in enumerate(list(zip(data_array, error_array, headers))): image, error, header = enum # Get current pixel size w = WCS(header).deepcopy() + new_header = deepcopy(header) # Compute binning ratio if scale.lower() in ['px', 'pixel']: - Dxy = np.array([pxsize,]*2) + Dxy_arr[i] = np.array([pxsize,]*2) elif scale.lower() in ['arcsec','arcseconds']: - Dxy = np.floor(pxsize/np.abs(w.wcs.cdelt)/3600.).astype(int) + Dxy_arr[i] = np.array(pxsize/np.abs(w.wcs.cdelt)/3600.) elif scale.lower() in ['full','integrate']: - Dxy = np.floor(image.shape).astype(int) + Dxy_arr[i] = image.shape else: raise ValueError("'{0:s}' invalid scale for binning.".format(scale)) + new_shape = np.ceil(min(image.shape/Dxy_arr,key=lambda x:x[0]+x[1])).astype(int) + for i, enum in enumerate(list(zip(data_array, error_array, headers))): + image, error, header = enum + # Get current pixel size + w = WCS(header).deepcopy() + new_header = deepcopy(header) + + Dxy = image.shape/new_shape if (Dxy < 1.).any(): raise ValueError("Requested pixel size is below resolution.") - new_shape = (image.shape//Dxy).astype(int) + new_shape = np.ceil(image.shape/Dxy).astype(int) # Rebin data rebin_data = bin_ndarray(image, new_shape=new_shape, @@ -606,10 +629,15 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, rebinned_error.append(np.sqrt(rms_image**2 + new_error**2)) # Update header - w = w.slice((np.s_[::Dxy[0]], np.s_[::Dxy[1]])) - header['NAXIS1'],header['NAXIS2'] = w.array_shape - header.update(w.to_header()) - rebinned_headers.append(header) + #nw = w.slice((np.s_[::Dxy[0]], np.s_[::Dxy[1]])) + nw = w.deepcopy() + nw.wcs.cdelt *= Dxy + nw.wcs.crpix /= Dxy + nw.array_shape = new_shape + new_header['NAXIS1'],new_header['NAXIS2'] = nw.array_shape + for key, val in nw.to_header().items(): + new_header.set(key,val) + rebinned_headers.append(new_header) if not data_mask is None: data_mask = bin_ndarray(data_mask,new_shape=new_shape,operation='average') > 0.80 @@ -1180,6 +1208,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers, if mask.any(): print("WARNING : found {0:d} pixels for which I_pol > I_stokes".format(I_stokes[mask].size)) + # Statistical error: Poisson noise is assumed + sigma_flux = np.array([np.sqrt(flux/head['exptime']) for flux,head in zip(pol_flux,pol_headers)]) + s_I2_stat = np.sum([coeff_stokes[0,i]**2*sigma_flux[i]**2 for i in range(len(sigma_flux))],axis=0) + s_Q2_stat = np.sum([coeff_stokes[1,i]**2*sigma_flux[i]**2 for i in range(len(sigma_flux))],axis=0) + s_U2_stat = np.sum([coeff_stokes[2,i]**2*sigma_flux[i]**2 for i in range(len(sigma_flux))],axis=0) + # Compute the derivative of each Stokes parameter with respect to the polarizer orientation dI_dtheta1 = 2.*pol_eff[0]/A*(pol_eff[2]*np.cos(-2.*theta[2]+2.*theta[0])*(pol_flux[1]-I_stokes) - pol_eff[1]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes)) dI_dtheta2 = 2.*pol_eff[1]/A*(pol_eff[0]*np.cos(-2.*theta[0]+2.*theta[1])*(pol_flux[2]-I_stokes) - pol_eff[2]*np.cos(-2.*theta[1]+2.*theta[2])*(pol_flux[0]-I_stokes)) @@ -1205,9 +1239,9 @@ def compute_Stokes(data_array, error_array, data_mask, headers, #np.savetxt("output/sU_dir.txt", np.sqrt(s_U2_axis)) # Add quadratically the uncertainty to the Stokes covariance matrix - Stokes_cov[0,0] += s_I2_axis - Stokes_cov[1,1] += s_Q2_axis - Stokes_cov[2,2] += s_U2_axis + Stokes_cov[0,0] += s_I2_axis + s_I2_stat + Stokes_cov[1,1] += s_Q2_axis + s_Q2_stat + Stokes_cov[2,2] += s_U2_axis + s_U2_stat #Compute integrated values for P, PA before any rotation mask = np.logical_and(data_mask.astype(bool), (I_stokes > 0.))