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

fix background estimation in get_error

Browse Source
This commit is contained in:
Tibeuleu
2022-11-25 16:40:36 +01:00
parent 1052784286
commit 93f43394e2
118 changed files with 43 additions and 38 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
src/FOC_reduction.py
View File

@@ -123,8 +123,8 @@ def main():
# Initial crop
display_crop = False
# Error estimation
error_sub_shape = (80,80)
display_error = False
error_sub_shape = (15,15)
display_error = True
# Data binning
rebin = True
if rebin:
@@ -143,7 +143,7 @@ def main():
rotate_data = False #rotation to North convention can give erroneous results
# Final crop
crop = False #Crop to desired ROI
final_display = False
final_display = True
# Polarization map output
figname = 'NGC1068_FOC' #target/intrument name
figtype = '_combine_FWHM010' #additionnal informations
@@ -163,10 +163,12 @@ def main():
# Deconvolve data using Richardson-Lucy iterative algorithm with a gaussian PSF of given FWHM.
if deconvolve:
data_array = proj_red.deconvolve_array(data_array, headers, psf=psf, FWHM=psf_FWHM, scale=psf_scale, shape=psf_shape, iterations=iterations, algo=algo)
# Estimate error from data background, estimated from sub-image of desired sub_shape.
background = None
if px_scale.lower() not in ['full','integrate']:
#data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_shape=error_sub_shape, display=display_error, savename=figname+"_errors", plots_folder=plots_folder, return_background=True)
data_array, error_array, headers, background = proj_red.get_error2(data_array, headers, error_array, display=display_error, savename=figname+"_errors", plots_folder=plots_folder, return_background=True)
data_array, error_array, headers, background = proj_red.get_error_hist(data_array, headers, error_array, display=display_error, savename=figname+"_errors", plots_folder=plots_folder, return_background=True)
# data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, display=display_error, savename=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)

2
src/comparison_Kishimoto.py
View File

@@ -20,7 +20,7 @@ data_K = {}
data_S = {}
for d,i in zip(['I','Q','U','P','PA','sI','sQ','sU','sP','sPA'],[0,1,2,5,8,(3,0,0),(3,1,1),(3,2,2),6,9]):
data_K[d] = np.loadtxt(path_join(root_dir_K,d+'.txt'))
with fits.open(path_join(root_dir_data_S,'NGC1068_K_FOC_bin10px.fits')) as f:
with fits.open(path_join(root_dir_data_S,'NGC1068_old_FOC_bin10px.fits')) as f:
if not type(i) is int:
data_S[d] = np.sqrt(f[i[0]].data[i[1],i[2]])
else:

14
src/lib/plots.py
View File

@@ -309,8 +309,8 @@ 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 = 0., 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.)
vmin, vmax = np.min(stkI.data[mask]*convert_flux)/5., 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)
print("Total flux contour levels : ", levelsI)
@@ -320,8 +320,8 @@ 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 = 0., np.max(stkI.data[pf_mask]*convert_flux*pol.data[pf_mask])
im = ax.imshow(stkI.data*convert_flux*pol.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
vmin, vmax = np.min(stkI.data[mask]*convert_flux)/5., 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)
print("Polarized flux contour levels : ", levelsPf)
@@ -375,7 +375,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
#ax.clabel(cont,inline=True,fontsize=6)
else:
# Defaults to intensity map
vmin, vmax = np.min(stkI.data[SNRi > SNRi_cut]*convert_flux)/10., np.max(stkI.data[SNRi > SNRi_cut]*convert_flux)
vmin, vmax = np.min(stkI.data[SNRi > SNRi_cut]*convert_flux)/5., np.max(stkI.data[SNRi > SNRi_cut]*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.)
@@ -1730,12 +1730,12 @@ class pol_map(object):
self.display_selection = "total_flux"
if self.display_selection.lower() in ['total_flux']:
self.data = self.I*self.convert_flux
vmin, vmax = np.min(self.data[self.cut])/10., np.max(self.data[self.data > 0.])
vmin, vmax = np.min(self.data[self.cut])/5., np.max(self.data[self.data > 0.])
norm = LogNorm(vmin, vmax)
label = r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]"
elif self.display_selection.lower() in ['pol_flux']:
self.data = self.I*self.convert_flux*self.P
vmin, vmax = np.min(self.I[self.cut]*self.convert_flux)/10., np.max(self.I[self.data > 0.]*self.convert_flux)
vmin, vmax = np.min(self.I[self.cut]*self.convert_flux)/5., np.max(self.I[self.data > 0.]*self.convert_flux)
norm = LogNorm(vmin, vmax)
label = r"$F_{\lambda} \cdot P$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]"
elif self.display_selection.lower() in ['pol_deg']:

37
src/lib/reduction.py
View File

@@ -407,7 +407,7 @@ def deconvolve_array(data_array, headers, psf='gaussian', FWHM=1., scale='px',
return deconv_array
def get_error2(data_array, headers, error_array=None, data_mask=None,
def get_error_hist(data_array, headers, error_array=None, data_mask=None,
display=False, savename=None, plots_folder="",
return_background=False):
"""
@@ -472,12 +472,15 @@ def get_error2(data_array, headers, error_array=None, data_mask=None,
background = np.zeros((data.shape[0]))
if display:
plt.rcParams.update({'font.size': 15})
filt_obs = {"POL0":0, "POL60":0, "POL120":0}
fig_h, ax_h = plt.subplots(figsize=(10,6), constrained_layout=True)
date_time = np.array([headers[i]['date-obs']+';'+headers[i]['time-obs']
for i in range(len(headers))])
date_time = np.array([datetime.strptime(d,'%Y-%m-%d;%H:%M:%S')
for d in date_time])
for i, image in enumerate(data):
filt_obs[headers[i]['filtnam1']] += 1
#Compute the Count-rate histogram for the image
n_mask = np.logical_and(mask,image>0.)
@@ -496,8 +499,9 @@ def get_error2(data_array, headers, error_array=None, data_mask=None,
#bkg = np.percentile(image[image<hist_max],25.)
#bkg = 0.95*hist_max
if display:
ax_h.plot(bin_centers,hist,'+',color="C{0:d}".format(i),alpha=0.8,label=headers[i]['filtnam1']+' ('+str(date_time[i])+") with n_bins = {0:d}".format(n_bins))
ax_h.plot(bin_centers,hist,'+',color="C{0:d}".format(i),alpha=0.8,label=headers[i]['filtnam1']+' (Obs '+str(filt_obs[headers[i]['filtnam1']])+')')
ax_h.plot([bkg,bkg],[hist.min(), hist.max()],'x--',color="C{0:d}".format(i),alpha=0.8)
print(headers[i]['filtnam1']+' ('+str(date_time[i])+') : n_bins =',n_bins,'; bkg = {0:.2e}'.format(bkg))
error_bkg[i] *= bkg
# Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999)
@@ -515,9 +519,9 @@ def get_error2(data_array, headers, error_array=None, data_mask=None,
#Substract background
n_data_array[i][data_mask] = n_data_array[i][data_mask] - bkg
n_data_array[i][np.logical_and(data_mask,n_data_array[i] <= 0.01*bkg)] = 0.01*bkg#n_data_array[i][np.logical_and(data_mask,n_data_array[i] > 0.)].min()
n_data_array[i][np.logical_and(data_mask,n_data_array[i] <= 0.01*bkg)] = 0.01*bkg
std_bkg[i] = image[image<2*bkg].std()
std_bkg[i] = image[np.abs(image-bkg)/bkg<1.].std()
background[i] = bkg
if (data_array[i] < 0.).any():
@@ -536,7 +540,6 @@ def get_error2(data_array, headers, error_array=None, data_mask=None,
ax_h.set_title("Histogram for each observation")
plt.legend()
plt.rcParams.update({'font.size': 15})
convert_flux = np.array([head['photflam'] for head in headers])
filt = np.array([headers[i]['filtnam1'] for i in range(len(headers))])
dict_filt = {"POL0":'r', "POL60":'g', "POL120":'b'}
@@ -596,7 +599,7 @@ def get_error2(data_array, headers, error_array=None, data_mask=None,
plt.show()
if return_background:
return n_data_array, n_error_array, headers, background #np.array([n_error_array[i][0,0] for i in range(n_error_array.shape[0])])
return n_data_array, n_error_array, headers, background
else:
return n_data_array, n_error_array, headers
@@ -698,7 +701,7 @@ def get_error(data_array, headers, error_array=None, data_mask=None,
# Compute error : root mean square of the background
sub_image = image[minima[0]:minima[0]+sub_shape[0],minima[1]:minima[1]+sub_shape[1]]
#bkg = np.std(sub_image) # Previously computed using standard deviation over the background
bkg = np.sqrt(np.sum((sub_image-sub_image.mean())**2)/sub_image.size)
bkg = np.sqrt(np.sum(sub_image**2)/sub_image.size)
error_bkg[i] *= bkg
# Quadratically add uncertainties in the "correction factors" (see Kishimoto 1999)
@@ -712,13 +715,13 @@ def get_error(data_array, headers, error_array=None, data_mask=None,
#estimated to less than 3%
err_flat = data_array[i]*0.03
error_array[i] = np.sqrt(error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2)
n_error_array[i] = np.sqrt(error_array[i]**2 + error_bkg[i]**2 + err_wav**2 + err_psf**2 + err_flat**2)
#Substract background
n_data_array[i][data_mask] = n_data_array[i][data_mask] - bkg
n_data_array[i][np.logical_and(data_mask,n_data_array[i] <= 0.)] = n_data_array[i][np.logical_and(data_mask,n_data_array[i] > 0.)].min()
n_data_array[i][np.logical_and(data_mask,n_data_array[i] <= 0.01*bkg)] = 0.01*bkg
std_bkg[i] = image[image<2*bkg].std()
std_bkg[i] = image[np.abs(image-bkg)/bkg<1.].std()
background[i] = bkg
if (data_array[i] < 0.).any():
@@ -800,9 +803,9 @@ def get_error(data_array, headers, error_array=None, data_mask=None,
plt.show()
if return_background:
return data_array, error_array, headers, np.array([error_array[i][0,0] for i in range(error_array.shape[0])])
return n_data_array, n_error_array, headers, background
else:
return data_array, error_array, headers
return n_data_array, n_error_array, headers
def rebin_array(data_array, error_array, headers, pxsize, scale,
@@ -990,7 +993,7 @@ def align_data(data_array, headers, error_array=None, background=None,
raise ValueError("All images in data_array must have same shape as\
ref_data")
if (error_array is None) or (background is None):
_, error_array, headers, background = get_error2(data_array, headers, return_background=True)
_, error_array, headers, background = get_error(data_array, headers, return_background=True)
# Crop out any null edges
#(ref_data must be cropped as well)
@@ -1519,9 +1522,9 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
for header in headers:
header['P_int'] = (P_diluted, 'Integrated polarization degree')
header['P_int_err'] = (P_diluted_err, 'Integrated polarization degree error')
header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error')
header['PA_int'] = (PA_diluted, 'Integrated polarization angle')
header['PA_int_err'] = (PA_diluted_err, 'Integrated polarization angle error')
header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
return I_stokes, Q_stokes, U_stokes, Stokes_cov
@@ -1778,9 +1781,9 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
for header in new_headers:
header['P_int'] = (P_diluted, 'Integrated polarization degree')
header['P_int_err'] = (P_diluted_err, 'Integrated polarization degree error')
header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error')
header['PA_int'] = (PA_diluted, 'Integrated polarization angle')
header['PA_int_err'] = (PA_diluted_err, 'Integrated polarization angle error')
header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers

16
src/overplot.py
View File

@@ -6,7 +6,7 @@ from copy import deepcopy
from lib.plots import overplot_radio, overplot_pol, align_pol
from matplotlib.colors import LogNorm
Stokes_UV = fits.open("../data/IC5063_x3nl030/IC5063_FOC_combine_FWHM020.fits")
Stokes_UV = fits.open("../data/IC5063_x3nl030/IC5063_old_FOC_combine_FWHM020.fits")
Stokes_18GHz = fits.open("../data/IC5063_x3nl030/radio/IC5063.18GHz.fits")
Stokes_24GHz = fits.open("../data/IC5063_x3nl030/radio/IC5063.24GHz.fits")
Stokes_103GHz = fits.open("../data/IC5063_x3nl030/radio/I5063_103GHz.fits")
@@ -20,30 +20,30 @@ levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
#levels18GHz = np.array([0.6, 1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_18GHz[0].data.max()
levels18GHz = levelsMorganti*0.28*1e-3
A = overplot_radio(Stokes_UV, Stokes_18GHz)
A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/18GHz_overplot.png')
A.plot(levels=levels18GHz, SNRp_cut=3.0, SNRi_cut=30.0, savename='../plots/IC5063_x3nl030/18GHz_overplot_old.png')
#levels24GHz = np.array([1.,1.5, 3, 6, 12, 24, 48, 96])/100.*Stokes_24GHz[0].data.max()
levels24GHz = levelsMorganti*0.46*1e-3
B = overplot_radio(Stokes_UV, Stokes_24GHz)
B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/24GHz_overplot.png')
B.plot(levels=levels24GHz, SNRp_cut=3.0, SNRi_cut=30.0, savename='../plots/IC5063_x3nl030/24GHz_overplot_old.png')
levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.]))
C = overplot_radio(Stokes_UV, Stokes_103GHz)
C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/103GHz_overplot.png')
C.plot(levels=levels103GHz, SNRp_cut=3.0, SNRi_cut=30.0, savename='../plots/IC5063_x3nl030/103GHz_overplot_old.png')
levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.]))
D = overplot_radio(Stokes_UV, Stokes_229GHz)
D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/229GHz_overplot.png')
D.plot(levels=levels229GHz, SNRp_cut=3.0, SNRi_cut=30.0, savename='../plots/IC5063_x3nl030/229GHz_overplot_old.png')
levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.]))
E = overplot_radio(Stokes_UV, Stokes_357GHz)
E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/357GHz_overplot.png')
E.plot(levels=levels357GHz, SNRp_cut=3.0, SNRi_cut=30.0, savename='../plots/IC5063_x3nl030/357GHz_overplot_old.png')
#F = overplot_pol(Stokes_UV, Stokes_S2)
#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='../plots/IC5063_x3nl030/S2_overplot_old.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
G = overplot_pol(Stokes_UV, Stokes_IR, cmap='inferno')
G.plot(SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/IR_overplot.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
G.plot(SNRp_cut=3.0, SNRi_cut=30.0, savename='../plots/IC5063_x3nl030/IR_overplot_old.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
#data_folder1 = "../data/M87/POS1/"
#plots_folder1 = "../plots/M87/POS1/"