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

plots clean-up, correct crop, prepare analysis

Browse Source
This commit is contained in:
Thibault Barnouin
2022-04-04 11:15:25 +02:00
parent 1ddff076b8
commit c4311bbb4b
399 changed files with 256 additions and 90 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
src/lib/convex_hull.py
View File

@@ -7,6 +7,33 @@ from copy import deepcopy
import numpy as np
def clean_ROI(image):
H,J = [],[]
shape = np.array(image.shape)
row, col = np.indices(shape)
for i in range(0,shape[0]):
r = row[i,:][image[i,:]>0.]
c = col[i,:][image[i,:]>0.]
if len(r)>1 and len(c)>1:
H.append((r[0],c[0]))
H.append((r[-1],c[-1]))
H = np.array(H)
for j in range(0,shape[1]):
r = row[:,j][image[:,j]>0.]
c = col[:,j][image[:,j]>0.]
if len(r)>1 and len(c)>1:
J.append((r[0],c[0]))
J.append((r[-1],c[-1]))
J = np.array(J)
xmin = np.min([H[:,1].min(),J[:,1].min()])
xmax = np.max([H[:,1].max(),J[:,1].max()])+1
ymin = np.min([H[:,0].min(),J[:,0].min()])
ymax = np.max([H[:,0].max(),J[:,0].max()])+1
return np.array([xmin,xmax,ymin,ymax])
# Define angle and vectors operations
def vector(A, B):
"""

35
src/lib/fits.py
View File

@@ -12,7 +12,7 @@ prototypes :
import numpy as np
from astropy.io import fits
from astropy import wcs
from lib.convex_hull import image_hull
from lib.convex_hull import image_hull, clean_ROI
def get_obs_data(infiles, data_folder="", compute_flux=False):
@@ -78,33 +78,6 @@ def get_obs_data(infiles, data_folder="", compute_flux=False):
return data_array, headers
def clean_ROI(image):
H,J = [],[]
shape = np.array(image.shape)
row, col = np.indices(shape)
for i in range(0,shape[0]):
r = row[i,:][image[i,:]>0.]
c = col[i,:][image[i,:]>0.]
if len(r)>1 and len(c)>1:
H.append((r[0],c[0]))
H.append((r[-1],c[-1]))
H = np.array(H)
for j in range(0,shape[1]):
r = row[:,j][image[:,j]>0.]
c = col[:,j][image[:,j]>0.]
if len(r)>1 and len(c)>1:
J.append((r[0],c[0]))
J.append((r[-1],c[-1]))
J = np.array(J)
xmin = np.min([H[:,1].min(),J[:,1].min()])-1
xmax = np.max([H[:,1].max(),J[:,1].max()])+1
ymin = np.min([H[:,0].min(),J[:,0].min()])-1
ymax = np.max([H[:,0].max(),J[:,0].max()])+1
return np.array([xmin,xmax,ymin,ymax])
def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
s_P_P, PA, s_PA, s_PA_P, headers, data_mask, filename, data_folder="",
return_hdul=False):
@@ -145,10 +118,10 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
exp_tot = np.array([header['exptime'] for header in headers]).sum()
new_wcs = wcs.WCS(ref_header).deepcopy()
vertex = clean_ROI(1.-data_mask)
vertex = clean_ROI(data_mask)
shape = vertex[1::2]-vertex[0::2]
new_wcs.array_shape = shape
new_wcs.wcs.crpix -= vertex[0::2]
new_wcs.wcs.crpix = np.array(new_wcs.wcs.crpix) - vertex[0::-2]
header = new_wcs.to_header()
header['instrume'] = (ref_header['instrume'], 'Instrument from which data is reduced')
@@ -179,7 +152,7 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
new_Stokes_cov = np.zeros((3,3,shape[0],shape[1]))
for i in range(3):
for j in range(3):
Stokes_cov[i,j][data_mask] = 0.
Stokes_cov[i,j][(1-data_mask).astype(bool)] = 0.
new_Stokes_cov[i,j] = Stokes_cov[i,j][vertex[0]:vertex[1],vertex[2]:vertex[3]]
Stokes_cov = new_Stokes_cov

191
src/lib/plots.py
View File

@@ -13,7 +13,7 @@ from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from matplotlib.widgets import RectangleSelector, Button
from matplotlib.widgets import RectangleSelector, Button, Slider
from matplotlib.transforms import Bbox
import matplotlib.font_manager as fm
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar, AnchoredDirectionArrows
@@ -219,23 +219,19 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
stkQ = Stokes[np.argmax([Stokes[i].header['datatype']=='Q_stokes' for i in range(len(Stokes))])]
stkU = Stokes[np.argmax([Stokes[i].header['datatype']=='U_stokes' for i in range(len(Stokes))])]
stk_cov = Stokes[np.argmax([Stokes[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes))])]
#pol = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg' for i in range(len(Stokes))])]
pol = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes))])]
pol_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes))])]
#pol_err_Poisson = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_deg_err_Poisson_noise' for i in range(len(Stokes))])]
pang = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_ang' for i in range(len(Stokes))])]
pang_err = Stokes[np.argmax([Stokes[i].header['datatype']=='Pol_ang_err' for i in range(len(Stokes))])]
if data_mask is None:
data_mask = Stokes[np.argmax([Stokes[i].header['datatype']=='Data_mask' for i in range(len(Stokes))])].data.astype(bool)
try:
if data_mask is None:
data_mask = Stokes[np.argmax([Stokes[i].header['datatype']=='Data_mask' for i in range(len(Stokes))])].data.astype(bool)
except KeyError:
data_mask = np.ones(stkI.shape).astype(bool)
pivot_wav = Stokes[0].header['photplam']
convert_flux = Stokes[0].header['photflam']
wcs = WCS(Stokes[0]).deepcopy()
#Get image mask
if data_mask is None:
data_mask = np.zeros(stkI.shape).astype(bool)
#Plot Stokes parameters map
if display is None or display.lower() == 'default':
plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder)
@@ -252,7 +248,6 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
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)
mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut)
# Look for pixel of max polarization
@@ -339,13 +334,13 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
if step_vec == 0:
pol.data[np.isfinite(pol.data)] = 1./2.
step_vec = 1
X, Y = np.meshgrid(np.linspace(0,stkI.data.shape[1],stkI.data.shape[1])-0.5, np.linspace(0,stkI.data.shape[0],stkI.data.shape[0])-0.5)
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.)
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=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w')
pol_sc = AnchoredSizeBar(ax.transData, 2., r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops)
ax.add_artist(pol_sc)
north_dir = AnchoredDirectionArrows(ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-Stokes[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
north_dir = AnchoredDirectionArrows(ax.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., angle=-Stokes[0].header['orientat'], color='white', text_props={'ec': None, 'fc': 'w', 'alpha': 1, 'lw': 0.4}, arrow_props={'ec': None,'fc':'w','alpha': 1,'lw': 1})
ax.add_artist(north_dir)
# Display instrument FOV
@@ -383,6 +378,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
plt.show()
return fig, ax
class align_maps(object):
"""
Class to interactively align maps with different WCS.
@@ -430,7 +426,7 @@ class align_maps(object):
self.ax1.add_artist(px_sc)
try:
north_dir1 = AnchoredDirectionArrows(self.ax1.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.map[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
north_dir1 = AnchoredDirectionArrows(self.ax1.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., angle=-self.map[0].header['orientat'], color='white', text_props={'ec': None, 'fc': 'w', 'alpha': 1, 'lw': 0.4}, arrow_props={'ec': None,'fc':'w','alpha': 1,'lw': 1})
self.ax1.add_artist(north_dir1)
except KeyError:
pass
@@ -452,7 +448,7 @@ class align_maps(object):
self.ax2.add_artist(px_sc)
try:
north_dir2 = AnchoredDirectionArrows(self.ax2.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.other_map[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
north_dir2 = AnchoredDirectionArrows(self.ax2.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.other_map[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
self.ax2.add_artist(north_dir2)
except KeyError:
pass
@@ -527,6 +523,8 @@ class overplot_radio(align_maps):
Inherit from class align_maps in order to get the same WCS on both maps.
"""
def overplot(self, other_levels, SNRp_cut=3., SNRi_cut=30., savename=None):
self.Stokes_UV = self.map
self.wcs_UV = self.wcs_map
#Get Data
obj = self.Stokes_UV[0].header['targname']
stkI = self.Stokes_UV[np.argmax([self.Stokes_UV[i].header['datatype']=='I_stokes' for i in range(len(self.Stokes_UV))])]
@@ -568,7 +566,7 @@ class overplot_radio(align_maps):
pol.data[np.isfinite(pol.data)] = 1./2.
step_vec = 1
X, Y = np.meshgrid(np.linspace(0,stkI.data.shape[0],stkI.data.shape[0]), np.linspace(0,stkI.data.shape[1],stkI.data.shape[1]))
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.)
Q = self.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=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w')
self.ax.autoscale(False)
@@ -584,7 +582,7 @@ class overplot_radio(align_maps):
px_size = self.wcs_UV.wcs.get_cdelt()[0]*3600.
px_sc = AnchoredSizeBar(self.ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops)
self.ax.add_artist(px_sc)
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
self.ax.add_artist(north_dir)
@@ -599,6 +597,7 @@ class overplot_radio(align_maps):
self.overplot(other_levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename)
plt.show(block=True)
class overplot_pol(align_maps):
"""
Class to overplot maps from different observations.
@@ -646,7 +645,7 @@ class overplot_pol(align_maps):
#Display full size polarization vectors
pol.data[np.isfinite(pol.data)] = 1./2.
step_vec = 1
X, Y = np.meshgrid(np.linspace(0,stkI.data.shape[0],stkI.data.shape[0]), np.linspace(0,stkI.data.shape[1],stkI.data.shape[1]))
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.)
Q = self.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=0.5,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,color='w')
@@ -661,7 +660,7 @@ class overplot_pol(align_maps):
px_size = self.wcs_other.wcs.get_cdelt()[0]*3600.
px_sc = AnchoredSizeBar(self.ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w', fontproperties=fontprops)
self.ax.add_artist(px_sc)
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': 'w', 'fc': 'w', 'alpha': 1,'lw': 2})
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.03, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, angle=-self.Stokes_UV[0].header['orientat'], color='w', arrow_props={'ec': None, 'fc': 'w', 'alpha': 1,'lw': 2})
self.ax.add_artist(north_dir)
@@ -750,7 +749,7 @@ class crop_map(object):
self.wcs_crop = self.wcs.deepcopy()
self.wcs_crop.array_shape = shape
if self.crpix_in_RS():
self.wcs_crop.wcs.crpix -= self.RSextent[::2]
self.wcs_crop.wcs.crpix = np.array(self.wcs_crop.wcs.crpix) - self.RSextent[::2]
else:
self.wcs_crop.wcs.crval = self.wcs.wcs_pix2world([self.RScenter],1)[0]
self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2]
@@ -793,6 +792,7 @@ class crop_map(object):
def writeto(self, filename):
self.hdul_crop.writeto(filename,overwrite=True)
class crop_Stokes(crop_map):
"""
Class to interactively crop a polarization map to desired Region of Interest.
@@ -810,7 +810,7 @@ class crop_Stokes(crop_map):
self.wcs_crop = self.wcs.deepcopy()
self.wcs_crop.array_shape = shape
if self.crpix_in_RS():
self.wcs_crop.wcs.crpix -= self.RSextent[::2]
self.wcs_crop.wcs.crpix = np.array(self.wcs_crop.wcs.crpix) - self.RSextent[::2]
else:
self.wcs_crop.wcs.crval = self.wcs.wcs_pix2world([self.RScenter],1)[0]
self.wcs_crop.wcs.crpix = self.RScenter-self.RSextent[::2]
@@ -849,4 +849,151 @@ class crop_Stokes(crop_map):
@property
def data_mask(self):
return self.hdul_crop[-1].data
return self.hdul_crop[-1].data
class pol_map(object):
"""
Class to interactively study polarization maps.
"""
def __init__(self,Stokes, SNRp_cut=3., SNRi_cut=30.):
self.Stokes = Stokes
self.wcs = deepcopy(WCS(Stokes[0].header))
self.SNRp_cut = SNRp_cut
self.SNRi_cut = SNRi_cut
self.SNRi = deepcopy(self.SNRi_cut)
self.SNRp = deepcopy(self.SNRp_cut)
self.region = None
#Get data
self.pivot_wav = self.Stokes[0].header['photplam']
self.convert_flux = self.Stokes[0].header['photflam']
self.I = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='I_stokes' for i in range(len(Stokes))])].data
self.Q = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Q_stokes' for i in range(len(Stokes))])].data
self.U = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='U_stokes' for i in range(len(Stokes))])].data
self.IQU_cov = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='IQU_cov_matrix' for i in range(len(Stokes))])].data
self.P = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Pol_deg_debiased' for i in range(len(Stokes))])].data
self.s_P = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Pol_deg_err' for i in range(len(Stokes))])].data
self.PA = self.Stokes[np.argmax([self.Stokes[i].header['datatype']=='Pol_ang' for i in range(len(Stokes))])].data
#Create figure
fontprops = fm.FontProperties(size=16)
self.fig = plt.figure(figsize=(15,15))
self.fig.subplots_adjust(hspace=0, wspace=0, right=0.9)
self.ax = self.fig.add_subplot(111,projection=self.wcs)
self.ax.set_facecolor('black')
self.ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
self.ax.coords[0].set_axislabel('Right Ascension (J2000)')
self.ax.coords[0].set_axislabel_position('t')
self.ax.coords[0].set_ticklabel_position('t')
self.ax.coords[1].set_axislabel('Declination (J2000)')
self.ax.coords[1].set_axislabel_position('l')
self.ax.coords[1].set_ticklabel_position('l')
self.ax.axis('equal')
#Display total flux
im = self.ax.imshow(self.I*self.convert_flux, vmin=0., vmax=self.I[self.I > 0.].max()*self.convert_flux, aspect='auto', cmap='inferno')
cbar_ax = self.fig.add_axes([0.95, 0.12, 0.01, 0.75])
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
#Display polarization vectors in SNR_cut
self.pol_vector()
#Display scales and orientation
px_size = self.wcs.wcs.cdelt[0]*3600.
px_sc = AnchoredSizeBar(self.ax.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='white', fontproperties=fontprops)
self.ax.add_artist(px_sc)
pol_sc = AnchoredSizeBar(self.ax.transData, 2., r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='white', fontproperties=fontprops)
self.ax.add_artist(pol_sc)
north_dir = AnchoredDirectionArrows(self.ax.transAxes, "E", "N", length=-0.08, fontsize=0.025, loc=1, aspect_ratio=-1, sep_y=0.01, sep_x=0.01, back_length=0., head_length=10., head_width=10., angle=-self.Stokes[0].header['orientat'], color='white', text_props={'ec': None, 'fc': 'w', 'alpha': 1, 'lw': 0.4}, arrow_props={'ec': None,'fc':'w','alpha': 1,'lw': 1})
self.ax.add_artist(north_dir)
#Display integrated values in ROI
self.pol_int
#Set axes for sliders (SNRp_cut, SNRi_cut)
ax_I_cut = self.fig.add_axes([0.125, 0.080, 0.35, 0.01])
ax_P_cut = self.fig.add_axes([0.125, 0.055, 0.35, 0.01])
ax_reset = self.fig.add_axes([0.125, 0.020, 0.05, 0.02])
SNRi_max = np.max(self.I[self.IQU_cov[0,0]>0.]/np.sqrt(self.IQU_cov[0,0][self.IQU_cov[0,0]>0.]))
SNRp_max = np.max(self.P[self.s_P>0.]/self.s_P[self.s_P > 0.])
s_I_cut = Slider(ax_I_cut,r"$SNR^{I}_{cut}$",1.,SNRi_max,valstep=1,valinit=self.SNRi_cut)
s_P_cut = Slider(ax_P_cut,r"$SNR^{P}_{cut}$",1.,SNRp_max,valstep=1,valinit=self.SNRp_cut)
b_reset = Button(ax_reset,"Reset")
def update_snri(val):
self.SNRi = val
self.quiver.remove()
self.pol_vector()
self.fig.canvas.draw_idle()
def update_snrp(val):
self.SNRp = val
self.quiver.remove()
self.pol_vector()
self.fig.canvas.draw_idle()
def reset(event):
s_I_cut.reset()
s_P_cut.reset()
s_I_cut.on_changed(update_snri)
s_P_cut.on_changed(update_snrp)
b_reset.on_clicked(reset)
plt.show(block=True)
@property
def cut(self):
s_I = np.sqrt(self.IQU_cov[0,0])
SNRp_mask, SNRi_mask = np.zeros(self.P.shape).astype(bool), np.zeros(self.I.shape).astype(bool)
SNRp_mask[self.s_P > 0.] = self.P[self.s_P > 0.] / self.s_P[self.s_P > 0.] > self.SNRp
SNRi_mask[s_I > 0.] = self.I[s_I > 0.] / s_I[s_I > 0.] > self.SNRi
return np.logical_and(SNRi_mask,SNRp_mask)
def pol_vector(self):
P_cut = np.ones(self.P.shape)*np.nan
P_cut[self.cut] = self.P[self.cut]
X, Y = np.meshgrid(np.arange(self.I.shape[1]),np.arange(self.I.shape[0]))
XY_U, XY_V = P_cut*np.cos(np.pi/2. + self.PA*np.pi/180.), P_cut*np.sin(np.pi/2. + self.PA*np.pi/180.)
self.quiver = self.ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='white')
return self.quiver
@property
def pol_int(self):
if self.region is None:
n_pix = self.I.size
s_I = np.sqrt(self.IQU_cov[0,0])
I_reg = self.I.sum()
I_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_I**2))
P_reg = self.Stokes[0].header['P_int']
P_reg_err = self.Stokes[0].header['P_int_err']
PA_reg = self.Stokes[0].header['PA_int']
PA_reg_err = self.Stokes[0].header['PA_int_err']
else:
n_pix = self.I[self.region].size
s_I = np.sqrt(self.IQU_cov[0,0])
s_Q = np.sqrt(self.IQU_cov[1,1])
s_U = np.sqrt(self.IQU_cov[2,2])
s_IQ = self.IQU_cov[0,1]
s_IU = self.IQU_cov[0,2]
s_QU = self.IQU_cov[1,2]
I_reg = self.I[self.region].sum()
Q_reg = self.Q[self.region].sum()
U_reg = self.U[self.region].sum()
I_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_I[self.region]**2))
Q_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_Q[self.region]**2))
U_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_U[self.region]**2))
IQ_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_IQ[self.region]**2))
IU_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_IU[self.region]**2))
QU_reg_err = np.sqrt(n_pix)*np.sqrt(np.sum(s_QU[self.region]**2))
P_reg = np.sqrt(Q_reg**2+U_reg**2)/I_reg
P_reg_err = np.sqrt((Q_reg**2*Q_reg_err**2 + U_reg**2*U_reg_err**2 + 2.*Q_reg*U_reg*QU_reg_err)/(Q_reg**2 + U_reg**2) + ((Q_reg/I_reg)**2 + (U_reg/I_reg)**2)*I_reg_err**2 - 2.*(Q_reg/I_reg)*IQ_reg_err - 2.*(U_reg/I_reg)*IU_reg_err)/I_reg
PA_reg = princ_angle((90./np.pi)*np.arctan2(U_reg,Q_reg))
PA_reg_err = (90./(np.pi*(Q_reg**2+U_reg**2)))*np.sqrt(U_reg**2*Q_reg_err**2 + Q_reg**2*U_reg_err**2 - 2.*Q_reg*U_reg*QU_reg_err)
return self.ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(self.pivot_wav,sci_not(I_reg*self.convert_flux,I_reg_err*self.convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,P_reg_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,PA_reg_err), color='white', fontsize=12, xy=(0.01, 0.90), xycoords='axes fraction')

41
src/lib/reduction.py
View File

@@ -50,7 +50,7 @@ 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.convex_hull import image_hull, clean_ROI
from lib.plots import plot_obs
from lib.cross_correlation import phase_cross_correlation
@@ -723,8 +723,11 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1.,
res_shape = int(np.ceil(np.sqrt(2.5)*np.max(shape[1:])))
rescaled_image = np.zeros((shape[0],res_shape,res_shape))
rescaled_error = np.ones((shape[0],res_shape,res_shape))
rescaled_mask = np.ones((shape[0],res_shape,res_shape),dtype=bool)
rescaled_mask = np.zeros((shape[0],res_shape,res_shape),dtype=bool)
res_center = (np.array(rescaled_image.shape[1:])/2).astype(int)
res_shift = res_center-ref_center
res_mask = np.zeros((res_shape,res_shape),dtype=bool)
res_mask[res_shift[0]:res_shift[0]+shape[1], res_shift[1]:res_shift[1]+shape[2]] = True
shifts, errors = [], []
for i,image in enumerate(data_array):
@@ -734,21 +737,19 @@ def align_data(data_array, headers, error_array=None, upsample_factor=1.,
shift, error, phase_diff = phase_cross_correlation(ref_data, image,
upsample_factor=upsample_factor)
# Rescale image to requested output
center = np.fix(ref_center-shift).astype(int)
res_shift = res_center-ref_center
rescaled_image[i,res_shift[0]:res_shift[0]+shape[1],
res_shift[1]:res_shift[1]+shape[2]] = deepcopy(image)
rescaled_error[i,res_shift[0]:res_shift[0]+shape[1],
res_shift[1]:res_shift[1]+shape[2]] = deepcopy(error_array[i])
rescaled_mask[i,res_shift[0]:res_shift[0]+shape[1],
res_shift[1]:res_shift[1]+shape[2]] = False
# Shift images to align
rescaled_image[i] = sc_shift(rescaled_image[i], shift, cval=0.)
rescaled_error[i] = sc_shift(rescaled_error[i], shift, cval=background[i])
rescaled_mask[i] = sc_shift(rescaled_mask[i], shift, cval=True)
curr_mask = sc_shift(res_mask, shift, cval=False)
mask_vertex = clean_ROI(curr_mask)
rescaled_mask[i,mask_vertex[2]:mask_vertex[3],mask_vertex[0]:mask_vertex[1]] = True
rescaled_image[i][rescaled_image[i] < 0.] = 0.
rescaled_image[i][1-rescaled_mask[i]] = 0.
rescaled_image[i][(1-rescaled_mask[i]).astype(bool)] = 0.
# Uncertainties from shifting
prec_shift = np.array([1.,1.])/upsample_factor
@@ -840,13 +841,13 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1.,
for r in range(smoothed.shape[0]):
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))
dist_rc = np.where(data_mask, np.sqrt((r-xx)**2+(c-yy)**2), fmax)
# 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] = (1.-data_mask[r,c])*np.sqrt(np.sum(weight*g_rc**2))/np.sum(weight*g_rc)
smoothed[r,c] = np.where(data_mask[r,c], np.sum(data_array*weight*g_rc)/np.sum(weight*g_rc), 0.)
error[r,c] = np.where(data_mask[r,c], np.sqrt(np.sum(weight*g_rc**2))/np.sum(weight*g_rc), 0.)
# Nan handling
error[np.isnan(smoothed)] = 0.
@@ -861,12 +862,12 @@ def smooth_data(data_array, error_array, data_mask, headers, FWHM=1.,
xx, yy = np.indices(image.shape)
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))
dist_rc = np.where(data_mask, np.sqrt((r-xx)**2+(c-yy)**2), fmax)
# 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] = (1.-data_mask[r,c])*np.sqrt(np.sum(error_array[i]**2*g_rc**2))
smoothed[i][r,c] = np.where(data_mask[r,c], np.sum(image*g_rc), 0.)
error[i][r,c] = np.where(data_mask[r,c], np.sqrt(np.sum(error_array[i]**2*g_rc**2)), 0.)
# Nan handling
error[i][np.isnan(smoothed[i])] = 0.
@@ -1185,7 +1186,7 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
Stokes_cov[0,0], Stokes_cov[1,1], Stokes_cov[2,2] = Stokes_error**2
#Compute integrated values for P, PA before any rotation
mask = (1-data_mask).astype(bool)
mask = deepcopy(data_mask).astype(bool)
n_pix = I_stokes[mask].size
I_diluted = I_stokes[mask].sum()
Q_diluted = Q_stokes[mask].sum()
@@ -1383,7 +1384,9 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
new_I_stokes = sc_rotate(new_I_stokes, ang, order=5, reshape=False, cval=0.)
new_Q_stokes = sc_rotate(new_Q_stokes, ang, order=5, reshape=False, cval=0.)
new_U_stokes = sc_rotate(new_U_stokes, ang, order=5, reshape=False, cval=0.)
new_data_mask = (sc_rotate(data_mask.astype(float)*10., ang, order=5, reshape=False, cval=10.).astype(int)).astype(bool)
new_data_mask = sc_rotate(data_mask.astype(float)*10., ang, order=5, reshape=False, cval=0.)
new_data_mask[new_data_mask < 2] = 0.
new_data_mask = new_data_mask.astype(bool)
for i in range(3):
for j in range(3):
new_Stokes_cov[i,j] = sc_rotate(new_Stokes_cov[i,j], ang,
@@ -1435,7 +1438,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
new_Stokes_cov[np.isnan(new_Stokes_cov)] = fmax
#Compute updated integrated values for P, PA
mask = (1-new_data_mask).astype(bool)
mask = deepcopy(new_data_mask).astype(bool)
n_pix = new_I_stokes[mask].size
I_diluted = new_I_stokes[mask].sum()
Q_diluted = new_Q_stokes[mask].sum()
@@ -1500,7 +1503,9 @@ def rotate_data(data_array, error_array, data_mask, headers, ang):
new_error_array.append(sc_rotate(error_array[i], ang, order=5, reshape=False,
cval=error_array.mean()))
new_data_array = np.array(new_data_array)
new_data_mask = sc_rotate(data_mask, ang, order=5, reshape=False, cval=True)
new_data_mask = sc_rotate(data_mask*10., ang, order=5, reshape=False, cval=0.)
new_data_mask[new_data_mask < 2] = 0.
new_data_mask = new_data_mask.astype(bool)
new_error_array = np.array(new_error_array)
for i in range(new_data_array.shape[0]):

18
src/lib/test_overplot.py
View File

@@ -1,11 +1,15 @@
#!/usr/bin/python3
from astropy.io import fits
import numpy as np
from copy import deepcopy
from plots import overplot_radio
Stokes_UV = fits.open("../../data/IC5063_x3nl030/IC5063_FOC_combine_FWHM020_pol.fits")
Stokes_UV = fits.open("../../data/IC5063_x3nl030/IC5063_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")
Stokes_229GHz = fits.open("../../data/IC5063_x3nl030/radio/I5063_229GHz.fits")
Stokes_357GHz = fits.open("../../data/IC5063_x3nl030/radio/I5063_357GHz.fits")
levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
@@ -18,3 +22,15 @@ A.plot(levels=levels18GHz, SNRp_cut=10.0, SNRi_cut=100.0, savename='../../plots/
levels24GHz = levelsMorganti*0.46*1e-3
B = overplot_radio(Stokes_UV, Stokes_24GHz)
B.plot(levels=levels24GHz, SNRp_cut=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/24GHz_overplot_forced.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=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/103GHz_overplot_forced.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=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/229GHz_overplot_forced.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=10.0, SNRi_cut=100.0, savename='../../plots/IC5063_x3nl030/357GHz_overplot_forced.png')