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

rework variable names on align class

Browse Source
This commit is contained in:
Thibault Barnouin
2024-02-22 18:05:08 +01:00
parent 98a414e32d
commit 118ba25122
6 changed files with 320 additions and 310 deletions
Download Patch File Download Diff File Expand all files Collapse all files

26
src/FOC_reduction.py
View File

@@ -50,7 +50,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
# 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_FWHM = 0.10 #If None, no smoothing is done
smoothing_scale = 'arcsec' #pixel or arcsec
# Rotation
@@ -157,7 +157,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
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, SNRi_cut=None)
I_bkg, Q_bkg, U_bkg, S_cov_bkg, _, _ = proj_red.rotate_Stokes(I_bkg, Q_bkg, U_bkg, S_cov_bkg, np.array(True).reshape(1,1), headers, SNRi_cut=None)
# Compute polarimetric parameters (polarization degree and angle).
# Compute polarimetric parameters (polarisation 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)
P_bkg, debiased_P_bkg, s_P_bkg, s_P_P_bkg, PA_bkg, s_PA_bkg, s_PA_P_bkg = proj_red.compute_pol(I_bkg, Q_bkg, U_bkg, S_cov_bkg, headers)
@@ -182,19 +182,19 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=
print("F_bkg({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'],*proj_plots.sci_not(I_bkg[0,0]*headers[0]['photflam'],np.sqrt(S_cov_bkg[0,0][0,0])*headers[0]['photflam'],2,out=int)))
print("P_bkg = {0:.1f} ± {1:.1f} %".format(debiased_P_bkg[0,0]*100.,np.ceil(s_P_bkg[0,0]*1000.)/10.))
print("PA_bkg = {0:.1f} ± {1:.1f} °".format(PA_bkg[0,0],np.ceil(s_PA_bkg[0,0]*10.)/10.))
# Plot polarization map (Background is either total Flux, Polarization degree or Polarization degree error).
# Plot polarisation map (Background is either total Flux, Polarization degree or Polarization degree error).
if px_scale.lower() not in ['full','integrate'] and not interactive:
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder)
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P"]), plots_folder=plots_folder, display='Pol_deg')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"PA"]), plots_folder=plots_folder, display='Pol_ang')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I_err"]), plots_folder=plots_folder, display='I_err')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi')
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype]), plots_folder=plots_folder)
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I"]), plots_folder=plots_folder, display='Intensity')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_flux"]), plots_folder=plots_folder, display='Pol_Flux')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P"]), plots_folder=plots_folder, display='Pol_deg')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"PA"]), plots_folder=plots_folder, display='Pol_ang')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"I_err"]), plots_folder=plots_folder, display='I_err')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"P_err"]), plots_folder=plots_folder, display='Pol_deg_err')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRi"]), plots_folder=plots_folder, display='SNRi')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim, step_vec=step_vec, vec_scale=vec_scale, savename="_".join([figname,figtype,"SNRp"]), plots_folder=plots_folder, display='SNRp')
elif not interactive:
proj_plots.polarization_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename="_".join([figname,figtype]), plots_folder=plots_folder, display='integrate')
proj_plots.polarisation_map(deepcopy(Stokes_test), data_mask, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename="_".join([figname,figtype]), plots_folder=plots_folder, display='integrate')
elif px_scale.lower() not in ['full', 'integrate']:
pol_map = proj_plots.pol_map(Stokes_test, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim)

15
src/lib/background.py
View File

@@ -100,19 +100,18 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non
exptime = headers[0]['exptime']
filt = headers[0]['filtnam1']
#plots
im = ax2.imshow(data0, norm=LogNorm(data0[data0>0.].mean()/10.,data0.max()), origin='lower', cmap='gray')
im2 = ax2.imshow(data0, norm=LogNorm(data0[data0>0.].mean()/10.,data0.max()), origin='lower', cmap='gray')
bkg_im = ax2.imshow(bkg_data0, origin='lower', cmap='Reds', alpha=0.5)
if not(rectangle is None):
x, y, width, height, angle, color = rectangle[0]
ax2.add_patch(Rectangle((x, y),width,height,edgecolor=color,fill=False,lw=2))
ax2.annotate(instr+":"+rootname, color='white', fontsize=10, xy=(0.02, 0.95), xycoords='axes fraction')
ax2.annotate(filt, color='white', fontsize=14, xy=(0.02, 0.02), xycoords='axes fraction')
ax2.annotate(str(exptime)+" s", color='white', fontsize=10, xy=(0.80, 0.02), xycoords='axes fraction')
ax2.set(xlabel='pixel offset', ylabel='pixel offset')
ax2.annotate(instr+":"+rootname, color='white', fontsize=10, xy=(0.01, 1.00), xycoords='axes fraction',verticalalignment='top', horizontalalignment='left')
ax2.annotate(filt, color='white', fontsize=14, xy=(0.01, 0.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='left')
ax2.annotate(str(exptime)+" s", color='white', fontsize=10, xy=(1.00, 0.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='right')
ax2.set(xlabel='pixel offset', ylabel='pixel offset',aspect='equal')
fig2.subplots_adjust(hspace=0, wspace=0, right=0.85)
cbar_ax = fig2.add_axes([0.9, 0.12, 0.02, 0.75])
fig2.colorbar(im, cax=cbar_ax, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
fig2.subplots_adjust(hspace=0, wspace=0, right=1.0)
fig2.colorbar(im2, ax=ax2, location='right', aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
if not(savename is None):
this_savename = deepcopy(savename)

11
src/lib/fits.py
View File

@@ -139,7 +139,8 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
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')
header['telescop'] = (ref_header['telescop'] if 'TELESCOP' in list(ref_header.keys()) else 'HST', 'telescope used to acquire data')
header['instrume'] = (ref_header['instrume'] if 'INSTRUME' in list(ref_header.keys()) else 'FOC', 'identifier for instrument used to acuire data')
header['photplam'] = (ref_header['photplam'], 'Pivot Wavelength')
header['photflam'] = (ref_header['photflam'], 'Inverse Sensitivity in DN/sec/cm**2/Angst')
header['exptot'] = (exp_tot, 'Total exposure time in sec')
@@ -147,10 +148,10 @@ def save_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, P, debiased_P, s_P,
header['targname'] = (ref_header['targname'], 'Target name')
header['orientat'] = (ref_header['orientat'], 'Angle between North and the y-axis of the image')
header['filename'] = (filename, 'Original filename')
header['P_int'] = (ref_header['P_int'], 'Integrated polarization degree')
header['P_int_err'] = (ref_header['P_int_err'], 'Integrated polarization degree error')
header['PA_int'] = (ref_header['PA_int'], 'Integrated polarization angle')
header['PA_int_err'] = (ref_header['PA_int_err'], 'Integrated polarization angle error')
header['P_int'] = (ref_header['P_int'], 'Integrated polarisation degree')
header['P_int_err'] = (ref_header['P_int_err'], 'Integrated polarisation degree error')
header['PA_int'] = (ref_header['PA_int'], 'Integrated polarisation angle')
header['PA_int_err'] = (ref_header['PA_int_err'], 'Integrated polarisation angle error')
#Crop Data to mask
if data_mask.shape != (1,1):

481
src/lib/plots.py
View File

@@ -8,8 +8,8 @@ prototypes :
- plot_Stokes(Stokes, savename, plots_folder)
Plot the I/Q/U maps from the Stokes HDUList.
- polarization_map(Stokes, data_mask, rectangle, SNRp_cut, SNRi_cut, step_vec, savename, plots_folder, display) -> fig, ax
Plots polarization map of polarimetric parameters saved in an HDUList.
- polarisation_map(Stokes, data_mask, rectangle, SNRp_cut, SNRi_cut, step_vec, savename, plots_folder, display) -> fig, ax
Plots polarisation map of polarimetric parameters saved in an HDUList.
class align_maps(map, other_map, **kwargs)
Class to interactively align maps with different WCS.
@@ -21,13 +21,13 @@ prototypes :
Class inherited from align_maps to overplot chandra data as contours.
class overplot_pol(align_maps)
Class inherited from align_maps to overplot UV polarization vectors on other maps.
Class inherited from align_maps to overplot UV polarisation vectors on other maps.
class crop_map(hdul, fig, ax)
Class to interactively crop a region of interest of a HDUList.
class crop_Stokes(crop_map)
Class inherited from crop_map to work on polarization maps.
Class inherited from crop_map to work on polarisation maps.
class image_lasso_selector(img, fig, ax)
Class to interactively select part of a map to work on.
@@ -36,7 +36,7 @@ prototypes :
Class to interactively simulate aperture integration.
class pol_map(Stokes, SNRp_cut, SNRi_cut, selection)
Class to interactively study polarization maps making use of the cropping and selecting tools.
Class to interactively study polarisation maps making use of the cropping and selecting tools.
"""
from copy import deepcopy
@@ -155,16 +155,18 @@ def plot_obs(data_array, headers, shape=None, vmin=None, vmax=None, rectangle=No
color='grey', alpha=0.5)
axe.plot([0,data.shape[1]-1], [data.shape[1]/2, data.shape[1]/2], '--', lw=1,
color='grey', alpha=0.5)
axe.annotate(instr+":"+rootname,color='white',fontsize=5,xy=(0.02, 0.95),
xycoords='axes fraction')
axe.annotate(filt,color='white',fontsize=10,xy=(0.02, 0.02),
xycoords='axes fraction')
axe.annotate(exptime,color='white',fontsize=5,xy=(0.80, 0.02),
xycoords='axes fraction')
axe.annotate(instr+":"+rootname,color='white',fontsize=5,xy=(0.01, 1.00),
xycoords='axes fraction', verticalalignment='top',
horizontalalignment='left')
axe.annotate(filt,color='white',fontsize=10,xy=(0.01, 0.01),
xycoords='axes fraction', verticalalignment='bottom',
horizontalalignment='left')
axe.annotate(exptime,color='white',fontsize=5,xy=(1.00, 0.01),
xycoords='axes fraction', verticalalignment='bottom',
horizontalalignment='right')
fig.subplots_adjust(hspace=0.01, wspace=0.01, right=0.85)
cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75])
fig.colorbar(im, cax=cbar_ax, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
fig.subplots_adjust(hspace=0.01, wspace=0.01, right=1.02)
fig.colorbar(im, ax=ax[:,:], location='right', shrink=0.75, aspect=50, pad=0.025, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
if not (savename is None):
#fig.suptitle(savename)
@@ -193,30 +195,29 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
Defaults to current folder.
"""
# Get data
stkI = Stokes[np.argmax([Stokes[i].header['datatype']=='I_stokes' for i in range(len(Stokes))])].data
stkQ = Stokes[np.argmax([Stokes[i].header['datatype']=='Q_stokes' for i in range(len(Stokes))])].data
stkU = Stokes[np.argmax([Stokes[i].header['datatype']=='U_stokes' for i in range(len(Stokes))])].data
stkI = Stokes['I_stokes'].data.copy()
stkQ = Stokes['Q_stokes'].data.copy()
stkU = Stokes['U_stokes'].data.copy()
wcs = WCS(Stokes[0]).deepcopy()
# Plot figure
plt.rcParams.update({'font.size': 10})
fig = plt.figure(figsize=(15,5))
fig, (axI,axQ,axU) = plt.subplots(ncols=3,figsize=(20,6),subplot_kw=dict(projection=wcs))
fig.subplots_adjust(hspace=0,wspace=0.75,bottom=0.01,top=0.99,left=0.08,right=0.95)
fig.suptitle("I, Q, U Stokes parameters")
ax = fig.add_subplot(131, projection=wcs)
im = ax.imshow(stkI, origin='lower', cmap='inferno')
plt.colorbar(im)
ax.set(xlabel="RA", ylabel="DEC", title=r"$I_{stokes}$")
imI = axI.imshow(stkI, origin='lower', cmap='inferno')
fig.colorbar(imI,ax=axI, aspect=50, shrink=0.50, pad=0.025,label='counts/sec')
axI.set(xlabel="RA", ylabel='DEC', title=r"$I_{stokes}$")
ax = fig.add_subplot(132, projection=wcs)
im = ax.imshow(stkQ, origin='lower', cmap='inferno')
plt.colorbar(im)
ax.set(xlabel="RA", ylabel="DEC", title=r"$Q_{stokes}$")
imQ = axQ.imshow(stkQ, origin='lower', cmap='inferno')
fig.colorbar(imQ,ax=axQ, aspect=50, shrink=0.50, pad=0.025,label='counts/sec')
axQ.set(xlabel="RA", ylabel='DEC', title=r"$Q_{stokes}$")
ax = fig.add_subplot(133, projection=wcs)
im = ax.imshow(stkU, origin='lower', cmap='inferno')
plt.colorbar(im)
ax.set(xlabel="RA", ylabel="DEC", title=r"$U_{stokes}$")
imU = axU.imshow(stkU, origin='lower', cmap='inferno')
fig.colorbar(imU,ax=axU, aspect=50, shrink=0.50, pad=0.025,label='counts/sec')
axU.set(xlabel="RA", ylabel='DEC', title=r"$U_{stokes}$")
if not (savename is None):
#fig.suptitle(savename+"_IQU")
@@ -229,10 +230,10 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""):
return 0
def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=30.,
def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_cut=30.,
flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"):
"""
Plots polarization map from Stokes HDUList.
Plots polarisation map from Stokes HDUList.
----------
Inputs:
Stokes : astropy.io.fits.hdu.hdulist.HDUList
@@ -255,12 +256,12 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
Defaults to None, limits are computed on the background value and the
maximum value in the cut.
step_vec : int, optional
Number of steps between each displayed polarization vector.
Number of steps between each displayed polarisation vector.
If step_vec = 2, every other vector will be displayed.
Defaults to 1
vec_scale : float, optional
Pixel length of displayed 100% polarization vector.
If vec_scale = 2, a vector of 50% polarization will be 1 pixel wide.
Pixel length of displayed 100% polarisation vector.
If vec_scale = 2, a vector of 50% polarisation will be 1 pixel wide.
Defaults to 2.
savename : str, optional
Name of the figure the map should be saved to. If None, the map won't
@@ -271,8 +272,8 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
be saved. Not used if savename is None.
Defaults to current folder.
display : str, optional
Choose the map to display between intensity (default), polarization
degree ('p','pol','pol_deg') or polarization degree error ('s_p',
Choose the map to display between intensity (default), polarisation
degree ('p','pol','pol_deg') or polarisation degree error ('s_p',
'pol_err','pol_deg_err').
Defaults to None (intensity).
----------
@@ -316,21 +317,19 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
poldata[np.logical_not(mask)] = np.nan
pangdata[np.logical_not(mask)] = np.nan
# Look for pixel of max polarization
# Look for pixel of max polarisation
if np.isfinite(pol.data).any():
p_max = np.max(pol.data[np.isfinite(pol.data)])
x_max, y_max = np.unravel_index(np.argmax(pol.data==p_max),pol.data.shape)
else:
print("No pixel with polarization information above requested SNR.")
print("No pixel with polarisation information above requested SNR.")
#Plot the map
plt.rcParams.update({'font.size': 10})
plt.rcdefaults()
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection=wcs)
ax.set_facecolor('k')
fig.subplots_adjust(hspace=0, wspace=0, right=0.9)
cbar_ax = fig.add_axes([0.95, 0.12, 0.01, 0.75])
fig, ax = plt.subplots(figsize=(10,10),subplot_kw=dict(projection=wcs))
ax.set(aspect='equal',fc='k')
fig.subplots_adjust(hspace=0,wspace=0,left=0.102,right=1.02)
if display.lower() in ['intensity']:
# If no display selected, show intensity map
@@ -343,8 +342,8 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
else:
vmin, vmax = flux_lim
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)
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]")
levelsI = np.logspace(0.31,1.955,6)/100.*vmax
print("Total flux contour levels : ", levelsI)
cont = ax.contour(stkI.data*convert_flux, levels=levelsI, colors='grey', linewidths=0.5)
#ax.clabel(cont,inline=True,fontsize=6)
@@ -360,25 +359,25 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
else:
vmin, vmax = flux_lim
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}$]")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, 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)
cont = ax.contour(stkI.data*convert_flux*pol.data, levels=levelsPf, colors='grey', linewidths=0.5)
#ax.clabel(cont,inline=True,fontsize=6)
elif display.lower() in ['p','pol','pol_deg']:
# Display polarization degree map
# Display polarisation degree map
display='p'
vmin, vmax = 0., 100.
im = ax.imshow(pol.data*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$P$ [%]")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$P$ [%]")
elif display.lower() in ['pa','pang','pol_ang']:
# Display polarization degree map
# Display polarisation degree map
display='pa'
vmin, vmax = 0., 180.
im = ax.imshow(princ_angle(pang.data), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.)
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$\theta_P$ [°]")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\theta_P$ [°]")
elif display.lower() in ['s_p','pol_err','pol_deg_err']:
# Display polarization degree error map
# Display polarisation degree error map
display='s_p'
if (SNRp>SNRp_cut).any():
vmin, vmax = 0., np.max(pol_err.data[SNRp > SNRp_cut])*100.
@@ -387,7 +386,7 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
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$ [%]")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=r"$\sigma_P$ [%]")
elif display.lower() in ['s_i','i_err']:
# Display intensity error map
display='s_i'
@@ -396,43 +395,41 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
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}$]")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025 , 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)])
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)
levelsSNRi = np.linspace(SNRi_cut, vmax*0.99, 5)
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}$")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025 , label=r"$I_{Stokes}/\sigma_{I}$")
elif display.lower() in ['snrp']:
# Display polarization degree signal-to-noise map
# Display polarisation degree signal-to-noise map
display='snrp'
vmin, vmax = 0., np.max(SNRp[np.isfinite(SNRp)])
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)
levelsSNRp = np.linspace(SNRp_cut, vmax*0.99, 5)
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}$")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025 , label=r"$P/\sigma_{P}$")
else:
# Defaults to intensity map
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.)
cbar = plt.colorbar(im, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]")
cbar = fig.colorbar(im, ax=ax, aspect=50, shrink=0.75, pad=0.025 , label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA$]")
#Get integrated values from header
n_pix = stkI.data[data_mask].size
@@ -462,12 +459,12 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
ax.add_artist(px_sc)
ax.add_artist(north_dir)
ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted*100.,P_diluted_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted,PA_diluted_err), color='white', xy=(0.01, 0.92), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_diluted*100.,P_diluted_err*100.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_diluted,PA_diluted_err), color='white', xy=(0.01, 1.00), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')],verticalalignment='top', horizontalalignment='left')
else:
if display.lower() == 'default':
ax.add_artist(px_sc)
ax.add_artist(north_dir)
ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2)), color='white', xy=(0.01, 0.97), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
ax.annotate(r"$F_{{\lambda}}^{{int}}$({0:.0f} $\AA$) = {1} $ergs \cdot cm^{{-2}} \cdot s^{{-1}} \cdot \AA^{{-1}}$".format(pivot_wav,sci_not(I_diluted*convert_flux,I_diluted_err*convert_flux,2)), color='white', xy=(0.01, 1.00), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')],verticalalignment='top', horizontalalignment='left')
# Display instrument FOV
if not(rectangle is None):
@@ -478,13 +475,8 @@ def polarization_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c
#ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
ax.coords[0].set_axislabel('Right Ascension (J2000)')
ax.coords[0].set_axislabel_position('t')
ax.coords[0].set_ticklabel_position('t')
ax.coords[1].set_axislabel('Declination (J2000)')
ax.coords[1].set_axislabel_position('l')
ax.coords[1].set_ticklabel_position('l')
ax.axis('equal')
ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)',labelpad=-1)
if not savename is None:
if not savename[-4:] in ['.png', '.jpg', '.pdf']:
@@ -499,96 +491,102 @@ class align_maps(object):
"""
Class to interactively align maps with different WCS.
"""
def __init__(self, map1, other_map, **kwargs):
def __init__(self, map, other_map, **kwargs):
self.aligned = False
self.map = map1
self.other_map = other_map
self.map = map
self.other = other_map
self.map_path = self.map.fileinfo(0)['filename']
self.other_path = self.other.fileinfo(0)['filename']
self.map_header = fits.getheader(self.map_path)
self.other_header = fits.getheader(self.other_path)
self.map_data = fits.getdata(self.map_path)
self.other_data = fits.getdata(self.other_path)
self.wcs_map = deepcopy(WCS(self.map[0])).celestial
if len(self.map[0].data.shape) == 4:
self.map[0].data = self.map[0].data[0,0]
elif len(self.map[0].data.shape) == 3:
self.map[0].data = self.map[0].data[0]
self.map_wcs = deepcopy(WCS(self.map_header)).celestial
if len(self.map_data.shape) == 4:
self.map_data = self.map_data[0,0]
elif len(self.map_data.shape) == 3:
self.map_data = self.map_data[0]
self.wcs_other = deepcopy(WCS(self.other_map[0])).celestial
if len(self.other_map[0].data.shape) == 4:
self.other_map[0].data = self.other_map[0].data[0,0]
elif len(self.other_map[0].data.shape) == 3:
self.other_map[0].data = self.other_map[0].data[0]
self.other_wcs = deepcopy(WCS(self.other_header)).celestial
if len(self.other_data.shape) == 4:
self.other_data = self.other_data[0,0]
elif len(self.other_data.shape) == 3:
self.other_data = self.other_data[0]
self.convert_flux, self.map_unit = (float(self.map[0].header['photflam']), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$") if "PHOTFLAM" in list(self.map[0].header.keys()) else (1., self.map[0].header['bunit'] if 'BUNIT' in list(self.map[0].header.keys()) else "Arbitray Units")
self.other_convert, self.other_map_unit = (float(self.other_map[0].header['photflam']), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$") if "PHOTFLAM" in list(self.other_map[0].header.keys()) else (1., self.other_map[0].header['bunit'] if 'BUNIT' in list(self.other_map[0].header.keys()) else "Arbitray Units")
#Get data
data = self.map[0].data
other_data = self.other_map[0].data
self.map_convert, self.map_unit = (float(self.map_header['photflam']), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$") if "PHOTFLAM" in list(self.map_header.keys()) else (1., self.map_header['bunit'] if 'BUNIT' in list(self.map_header.keys()) else "Arbitray Units")
self.other_convert, self.other_unit = (float(self.other_map[0].header['photflam']), r"$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$") if "PHOTFLAM" in list(self.other_header.keys()) else (1., self.other_header['bunit'] if 'BUNIT' in list(self.other_header.keys()) else "Arbitray Units")
self.map_observer = "/".join([self.map_header['telescop'],self.map_header['instrume']]) if "INSTRUME" in list(self.map_header.keys()) else self.map_header['telescop']
self.other_observer = "/".join([self.other_header['telescop'],self.other_header['instrume']]) if "INSTRUME" in list(self.other_header.keys()) else self.other_header['telescop']
plt.rcParams.update({'font.size': 10})
fontprops = fm.FontProperties(size=16)
self.fig_align = plt.figure(figsize=(20,10))
self.ax_map = self.fig_align.add_subplot(121, projection=self.wcs_map)
self.ax_other = self.fig_align.add_subplot(122, projection=self.wcs_other)
self.map_ax = self.fig_align.add_subplot(121, projection=self.map_wcs)
self.other_ax = self.fig_align.add_subplot(122, projection=self.other_wcs)
#Plot the UV map
other_kwargs = deepcopy(kwargs)
vmin, vmax = data[data > 0.].max()/1e3*self.convert_flux, data[data > 0.].max()*self.convert_flux
vmin, vmax = self.map_data[self.map_data > 0.].max()/1e3*self.map_convert, self.map_data[self.map_data > 0.].max()*self.map_convert
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["norm",LogNorm(vmin,vmax)]]]]:
try:
test = kwargs[key]
except KeyError:
for key_i, val_i in value:
kwargs[key_i] = val_i
im1 = self.ax_map.imshow(data*self.convert_flux, aspect='equal', **kwargs)
im1 = self.map_ax.imshow(self.map_data*self.map_convert, aspect='equal', **kwargs)
if kwargs['cmap'] in ['inferno','magma','Greys_r','binary_r','gist_yarg_r','gist_gray','gray','bone','pink','hot','afmhot','gist_heat','copper','gist_earth','gist_stern','gnuplot','gnuplot2','CMRmap','cubehelix','nipy_spectral','gist_ncar','viridis']:
self.ax_map.set_facecolor('black')
self.ax_other.set_facecolor('black')
self.map_ax.set_facecolor('black')
self.other_ax.set_facecolor('black')
font_color="white"
else:
self.ax_map.set_facecolor('white')
self.ax_other.set_facecolor('white')
self.map_ax.set_facecolor('white')
self.other_ax.set_facecolor('white')
font_color="black"
px_size1 = self.wcs_map.wcs.get_cdelt()[0]*3600.
px_sc1 = AnchoredSizeBar(self.ax_map.transData, 1./px_size1, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.ax_map.add_artist(px_sc1)
px_size1 = self.map_wcs.wcs.get_cdelt()[0]*3600.
px_sc1 = AnchoredSizeBar(self.map_ax.transData, 1./px_size1, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.map_ax.add_artist(px_sc1)
if 'PHOTPLAM' in list(self.map[0].header.keys()):
annote1 = self.ax_map.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(self.map[0].header['photplam']), color=font_color, fontsize=12, xy=(0.01, 0.93), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
if 'ORIENTAT' in list(self.map[0].header.keys()):
north_dir1 = AnchoredDirectionArrows(self.ax_map.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=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1,'lw': 0.5})
self.ax_map.add_artist(north_dir1)
if 'PHOTPLAM' in list(self.map_header.keys()):
annote1 = self.map_ax.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(self.map_header['photplam']), color=font_color, fontsize=12, xy=(0.01, 0.93), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
if 'ORIENTAT' in list(self.map_header.keys()):
north_dir1 = AnchoredDirectionArrows(self.map_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.map_header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1,'lw': 0.5})
self.map_ax.add_artist(north_dir1)
self.cr_map, = self.ax_map.plot(*(self.wcs_map.wcs.crpix-(1.,1.)), 'r+')
self.cr_map, = self.map_ax.plot(*(self.map_wcs.wcs.crpix-(1.,1.)), 'r+')
self.ax_map.set_title("Click on selected point of reference.")
self.ax_map.set_xlabel(label="Right Ascension (J2000)")
self.ax_map.set_ylabel(label="Declination (J2000)",labelpad=-1)
self.map_ax.set_title("{0:s} observation\nClick on selected point of reference.".format(self.map_observer))
self.map_ax.set_xlabel(label="Right Ascension (J2000)")
self.map_ax.set_ylabel(label="Declination (J2000)",labelpad=-1)
#Plot the other map
vmin, vmax = other_data[other_data > 0.].max()/1e3*self.other_convert, other_data[other_data > 0.].max()*self.other_convert
vmin, vmax = self.other_data[self.other_data > 0.].max()/1e3*self.other_convert, self.other_data[self.other_data > 0.].max()*self.other_convert
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["norm",LogNorm(vmin,vmax)]]]]:
try:
test = other_kwargs[key]
except KeyError:
for key_i, val_i in value:
other_kwargs[key_i] = val_i
im2 = self.ax_other.imshow(other_data*self.other_convert, aspect='equal', **other_kwargs)
im2 = self.other_ax.imshow(self.other_data*self.other_convert, aspect='equal', **other_kwargs)
px_size2 = self.wcs_other.wcs.get_cdelt()[0]*3600.
px_sc2 = AnchoredSizeBar(self.ax_other.transData, 1./px_size2, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.ax_other.add_artist(px_sc2)
px_size2 = self.other_wcs.wcs.get_cdelt()[0]*3600.
px_sc2 = AnchoredSizeBar(self.other_ax.transData, 1./px_size2, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.other_ax.add_artist(px_sc2)
if 'PHOTPLAM' in list(self.other_map[0].header.keys()):
annote2 = self.ax_other.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(self.other_map[0].header['photplam']), color='white', fontsize=12, xy=(0.01, 0.93), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
if 'ORIENTAT' in list(self.other_map[0].header.keys()):
north_dir2 = AnchoredDirectionArrows(self.ax_map.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=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1,'lw': 0.5})
self.ax_other.add_artist(north_dir2)
if 'PHOTPLAM' in list(self.other_header.keys()):
annote2 = self.other_ax.annotate(r"$\lambda$ = {0:.0f} $\AA$".format(self.other_header['photplam']), color='white', fontsize=12, xy=(0.01, 0.93), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
if 'ORIENTAT' in list(self.other_header.keys()):
north_dir2 = AnchoredDirectionArrows(self.map_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.other_header['orientat'], color=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1,'lw': 0.5})
self.other_ax.add_artist(north_dir2)
self.cr_other, = self.ax_other.plot(*(self.wcs_other.wcs.crpix-(1.,1.)), 'r+')
self.cr_other, = self.other_ax.plot(*(self.other_wcs.wcs.crpix-(1.,1.)), 'r+')
self.ax_other.set_title("Click on selected point of reference.")
self.ax_other.set_xlabel(label="Right Ascension (J2000)")
self.ax_other.set_ylabel(label="Declination (J2000)",labelpad=-1)
self.other_ax.set_title("{0:s} observation\nClick on selected point of reference.".format(self.other_observer))
self.other_ax.set_xlabel(label="Right Ascension (J2000)")
self.other_ax.set_ylabel(label="Declination (J2000)",labelpad=-1)
#Selection button
self.axapply = self.fig_align.add_axes([0.80, 0.01, 0.1, 0.04])
@@ -604,18 +602,18 @@ class align_maps(object):
self.on_close_align(event)
def get_aligned_wcs(self):
return self.wcs_map, self.wcs_other
return self.map_wcs, self.other_wcs
def onclick_ref(self, event) -> None:
if self.fig_align.canvas.manager.toolbar.mode == '':
if (event.inaxes is not None) and (event.inaxes == self.ax_map):
if (event.inaxes is not None) and (event.inaxes == self.map_ax):
x = event.xdata
y = event.ydata
self.cr_map.set(data=[x,y])
self.fig_align.canvas.draw_idle()
if (event.inaxes is not None) and (event.inaxes == self.ax_other):
if (event.inaxes is not None) and (event.inaxes == self.other_ax):
x = event.xdata
y = event.ydata
@@ -623,8 +621,8 @@ class align_maps(object):
self.fig_align.canvas.draw_idle()
def reset_align(self, event):
self.wcs_map.wcs.crpix = WCS(self.map[0].header).wcs.crpix[:2]
self.wcs_other.wcs.crpix = WCS(self.other_map[0].header).wcs.crpix[:2]
self.map_wcs.wcs.crpix = WCS(self.map_header).wcs.crpix[:2]
self.other_wcs.wcs.crpix = WCS(self.other_header).wcs.crpix[:2]
self.fig_align.canvas.draw_idle()
if self.aligned:
@@ -634,15 +632,15 @@ class align_maps(object):
def apply_align(self, event=None):
if np.array(self.cr_map.get_data()).shape == (2,1):
self.wcs_map.wcs.crpix = np.array(self.cr_map.get_data())[:,0]+(1.,1.)
self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data())[:,0]+(1.,1.)
else:
self.wcs_map.wcs.crpix = np.array(self.cr_map.get_data())+(1.,1.)
self.map_wcs.wcs.crpix = np.array(self.cr_map.get_data())+(1.,1.)
if np.array(self.cr_other.get_data()).shape == (2,1):
self.wcs_other.wcs.crpix = np.array(self.cr_other.get_data())[:,0]+(1.,1.)
self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data())[:,0]+(1.,1.)
else:
self.wcs_other.wcs.crpix = np.array(self.cr_other.get_data())+(1.,1.)
self.wcs_map.wcs.crval = np.array(self.wcs_map.pixel_to_world_values(*self.wcs_map.wcs.crpix))
self.wcs_other.wcs.crval = self.wcs_map.wcs.crval
self.other_wcs.wcs.crpix = np.array(self.cr_other.get_data())+(1.,1.)
self.map_wcs.wcs.crval = np.array(self.map_wcs.pixel_to_world_values(*self.map_wcs.wcs.crpix))
self.other_wcs.wcs.crval = self.map_wcs.wcs.crval
self.fig_align.canvas.draw_idle()
if self.aligned:
@@ -664,6 +662,24 @@ class align_maps(object):
plt.show(block=True)
return self.get_aligned_wcs()
def write_map_to(self,path="map.fits",suffix="aligned",data_dir="."):
new_head = deepcopy(self.map_header)
new_head.update(self.map_wcs.to_header())
new_hdul = fits.HDUList(fits.PrimaryHDU(self.map_data,new_head))
new_hdul.writeto("_".join([path[:-5],suffix])+".fits")
return 0
def write_other_to(self,path="other_map.fits",suffix="aligned",data_dir="."):
new_head = deepcopy(self.other_header)
new_head.update(self.other_wcs.to_header())
new_hdul = fits.HDUList(fits.PrimaryHDU(self.other_data,new_head))
new_hdul.writeto("_".join([path[:-5],suffix])+".fits")
return 0
def write_to(self,path1="map.fits",path2="other_map.fits",suffix="aligned",data_dir="."):
self.write_map_to(path=path1,suffix=suffix,data_dir=data_dir)
self.write_other_to(path=path2,suffix=suffix,data_dir=data_dir)
return 0
class overplot_radio(align_maps):
"""
@@ -672,23 +688,23 @@ class overplot_radio(align_maps):
"""
def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2, savename=None, **kwargs):
self.Stokes_UV = self.map
self.wcs_UV = self.wcs_map
self.wcs_UV = self.map_wcs
#Get Data
obj = self.Stokes_UV[0].header['targname']
stkI = deepcopy(self.Stokes_UV['I_STOKES'].data)
stk_cov = deepcopy(self.Stokes_UV['IQU_COV_MATRIX'].data)
stkI = self.Stokes_UV['I_STOKES'].data
stk_cov = self.Stokes_UV['IQU_COV_MATRIX'].data
pol = deepcopy(self.Stokes_UV['POL_DEG_DEBIASED'].data)
pol_err = deepcopy(self.Stokes_UV['POL_DEG_ERR'].data)
pang = deepcopy(self.Stokes_UV['POL_ANG'].data)
pol_err = self.Stokes_UV['POL_DEG_ERR'].data
pang = self.Stokes_UV['POL_ANG'].data
other_data = self.other_map[0].data
other_data = self.other_data
self.other_convert = 1.
if self.other_map_unit.lower() == 'jy/beam':
self.other_map_unit = r"mJy/Beam"
if self.other_unit.lower() == 'jy/beam':
self.other_unit = r"mJy/Beam"
self.other_convert = 1e3
other_freq = self.other_map[0].header['crval3'] if 'CRVAL3' in list(self.other_map[0].header.keys()) else 1.
other_freq = self.other_header['crval3'] if 'CRVAL3' in list(self.other_header.keys()) else 1.
self.convert_flux = self.Stokes_UV[0].header['photflam']
self.map_convert = self.Stokes_UV[0].header['photflam']
#Compute SNR and apply cuts
pol[pol == 0.] = np.nan
@@ -703,8 +719,8 @@ class overplot_radio(align_maps):
self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(10,10), subplot_kw=dict(projection=self.wcs_UV))
self.fig_overplot.subplots_adjust(hspace=0,wspace=0,bottom=0.1,left=0.1,top=0.8,right=1)
#Display UV intensity map with polarization vectors
vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.convert_flux,stkI[np.isfinite(stkI)].max()*self.convert_flux
#Display UV intensity map with polarisation vectors
vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.map_convert,stkI[np.isfinite(stkI)].max()*self.map_convert
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["norm",LogNorm(vmin,vmax)]]]]:
try:
test = kwargs[key]
@@ -717,33 +733,33 @@ class overplot_radio(align_maps):
else:
self.ax_overplot.set_facecolor('white')
font_color="black"
self.im = self.ax_overplot.imshow(stkI*self.convert_flux, aspect='equal', **kwargs)
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{\lambda}$ [{0:s}]".format(self.map_unit))
self.im = self.ax_overplot.imshow(stkI*self.map_convert, aspect='equal',label="{0:s} observation".format(self.map_observer), **kwargs)
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit))
#Display full size polarization vectors
#Display full size polarisation vectors
if vec_scale is None:
self.vec_scale = 2.
pol[np.isfinite(pol)] = 1./2.
else:
self.vec_scale = vec_scale
step_vec = 1
px_scale = self.wcs_other.wcs.get_cdelt()[0]/self.wcs_UV.wcs.get_cdelt()[0]
px_scale = self.other_wcs.wcs.get_cdelt()[0]/self.wcs_UV.wcs.get_cdelt()[0]
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.)
self.Q = self.ax_overplot.quiver(self.X[::step_vec,::step_vec],self.Y[::step_vec,::step_vec],self.U[::step_vec,::step_vec],self.V[::step_vec,::step_vec],units='xy',angles='uv',scale=1./self.vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,linewidth=0.5,color='white',edgecolor='black',label="HST/FOC polarisation map")
self.Q = self.ax_overplot.quiver(self.X[::step_vec,::step_vec],self.Y[::step_vec,::step_vec],self.U[::step_vec,::step_vec],self.V[::step_vec,::step_vec],units='xy',angles='uv',scale=1./self.vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,linewidth=0.5,color='white',edgecolor='black',label="{0:s} polarisation map".format(self.map_observer))
self.ax_overplot.autoscale(False)
#Display other map as contours
if levels is None:
levels = np.logspace(np.log(3)/np.log(10),2.,5)/100.*other_data[other_data > 0.].max()*self.other_convert
other_cont = self.ax_overplot.contour(other_data*self.other_convert, transform=self.ax_overplot.get_transform(self.wcs_other.celestial), levels=levels*self.other_convert, colors='grey')
levels = np.logspace(np.log(3)/np.log(10),2.,5)/100.*other_data[other_data > 0.].max()
other_cont = self.ax_overplot.contour(other_data*self.other_convert, transform=self.ax_overplot.get_transform(self.other_wcs.celestial), levels=levels*self.other_convert, colors='grey')
self.ax_overplot.clabel(other_cont, inline=True, fontsize=5)
other_proxy = Rectangle((0,0),1,1,fc='w',ec=other_cont.collections[0].get_edgecolor()[0], label=r"{0:s} contour".format(self.other_map[0].header['telescop']))
other_proxy = Rectangle((0,0),1,1,fc='w',ec=other_cont.collections[0].get_edgecolor()[0], label=r"{0:s} contour".format(self.other_observer))
self.ax_overplot.add_patch(other_proxy)
self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
self.ax_overplot.set_ylabel(label="Declination (J2000)",labelpad=-1)
self.fig_overplot.suptitle("HST/FOC UV polarization map of {0:s} overplotted with {1:.2f}GHz map in {2:s}.".format(obj, other_freq*1e-9, self.other_map_unit),wrap=True)
self.fig_overplot.suptitle("{0:s} polarisation map of {1:s} overplotted with {2:s} {3:.2f}GHz map in {4:s}.".format(self.map_observer, obj, self.other_observer, other_freq*1e-9, self.other_unit),wrap=True)
#Display pixel scale and North direction
fontprops = fm.FontProperties(size=16)
@@ -755,10 +771,12 @@ class overplot_radio(align_maps):
pol_sc = AnchoredSizeBar(self.ax_overplot.transData, self.vec_scale, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.ax_overplot.add_artist(pol_sc)
self.cr_map, = self.ax_overplot.plot(*(self.wcs_map.celestial.wcs.crpix-(1.,1.)), 'r+')
self.cr_other, = self.ax_overplot.plot(*(self.wcs_other.celestial.wcs.crpix-(1.,1.)), 'g+', transform=self.ax_overplot.get_transform(self.wcs_other))
self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1.,1.)), 'r+')
self.cr_other, = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix-(1.,1.)), 'g+', transform=self.ax_overplot.get_transform(self.other_wcs))
self.legend = self.ax_overplot.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.)
h,l = self.ax_overplot.get_legend_handles_labels()
h[np.argmax([li=="{0:s} polarisation map".format(self.map_observer) for li in l])] = FancyArrowPatch((0,0),(0,1),arrowstyle='-',fc='w',ec='k',lw=2)
self.legend = self.ax_overplot.legend(handles=h,labels=l,bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.)
if not(savename is None):
if not savename[-4:] in ['.png', '.jpg', '.pdf']:
@@ -781,17 +799,17 @@ class overplot_chandra(align_maps):
"""
def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2, zoom=1, savename=None, **kwargs):
self.Stokes_UV = self.map
self.wcs_UV = self.wcs_map
self.wcs_UV = self.map_wcs
#Get Data
obj = self.Stokes_UV[0].header['targname']
stkI = deepcopy(self.Stokes_UV['I_STOKES'].data)
stk_cov = deepcopy(self.Stokes_UV['IQU_COV_MATRIX'].data)
stkI = self.Stokes_UV['I_STOKES'].data
stk_cov = self.Stokes_UV['IQU_COV_MATRIX'].data
pol = deepcopy(self.Stokes_UV['POL_DEG_DEBIASED'].data)
pol_err = deepcopy(self.Stokes_UV['POL_DEG_ERR'].data)
pang = deepcopy(self.Stokes_UV['POL_ANG'].data)
pol_err = self.Stokes_UV['POL_DEG_ERR'].data
pang = self.Stokes_UV['POL_ANG'].data
other_data = deepcopy(self.other_map[0].data)
other_wcs = deepcopy(self.wcs_other)
other_data = deepcopy(self.other_data)
other_wcs = deepcopy(self.other_wcs)
if zoom != 1:
other_data = sc_zoom(other_data,zoom)
other_wcs.wcs.crpix *= zoom
@@ -811,8 +829,8 @@ class overplot_chandra(align_maps):
self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11,10), subplot_kw=dict(projection=self.wcs_UV))
self.fig_overplot.subplots_adjust(hspace=0,wspace=0,bottom=0.1,left=0.1,top=0.8,right=1)
#Display UV intensity map with polarization vectors
vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.convert_flux,stkI[np.isfinite(stkI)].max()*self.convert_flux
#Display UV intensity map with polarisation vectors
vmin, vmax = stkI[np.isfinite(stkI)].max()/1e3*self.map_convert,stkI[np.isfinite(stkI)].max()*self.map_convert
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["norm",LogNorm(vmin,vmax)]]]]:
try:
test = kwargs[key]
@@ -825,10 +843,10 @@ class overplot_chandra(align_maps):
else:
self.ax_overplot.set_facecolor('white')
font_color="black"
self.im = self.ax_overplot.imshow(stkI*self.convert_flux, aspect='equal', **kwargs)
self.im = self.ax_overplot.imshow(stkI*self.map_convert, aspect='equal', **kwargs)
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=50, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.map_unit))
#Display full size polarization vectors
#Display full size polarisation vectors
if vec_scale is None:
self.vec_scale = 2.
pol[np.isfinite(pol)] = 1./2.
@@ -838,7 +856,7 @@ class overplot_chandra(align_maps):
px_scale = 1./self.wcs_UV.wcs.get_cdelt()[0]
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.)
self.Q = self.ax_overplot.quiver(self.X[::step_vec,::step_vec],self.Y[::step_vec,::step_vec],self.U[::step_vec,::step_vec],self.V[::step_vec,::step_vec],units='xy',angles='uv',scale=1./self.vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,linewidth=0.5,color='white',edgecolor='black',label="HST/FOC polarisation map")
self.Q = self.ax_overplot.quiver(self.X[::step_vec,::step_vec],self.Y[::step_vec,::step_vec],self.U[::step_vec,::step_vec],self.V[::step_vec,::step_vec],units='xy',angles='uv',scale=1./self.vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1,linewidth=0.5,color='white',edgecolor='black',label="{0:s} polarisation map".format(self.map_observer))
proxy_Q = FancyArrowPatch((0,0),(0,1),arrowstyle='-',fc='w',ec='k',lw=3)
self.ax_overplot.autoscale(False)
@@ -846,15 +864,15 @@ class overplot_chandra(align_maps):
if levels is None:
levels = np.logspace(np.log(3)/np.log(10),2.,5)/100.*other_data[other_data > 0.].max()*self.other_convert
elif zoom != 1:
levels *= other_data.max()/self.other_map[0].data.max()
levels *= other_data.max()/self.other_data.max()
other_cont = self.ax_overplot.contour(other_data*self.other_convert, transform=self.ax_overplot.get_transform(other_wcs), levels=levels, colors='grey')
self.ax_overplot.clabel(other_cont, inline=True, fontsize=8)
other_proxy = Rectangle((0,0),1.,1.,fc='w',ec=other_cont.collections[0].get_edgecolor()[0], lw=2, label=r"{0:s} contour in counts".format(self.other_map[0].header['telescop']))
other_proxy = Rectangle((0,0),1.,1.,fc='w',ec=other_cont.collections[0].get_edgecolor()[0], lw=2, label=r"{0:s} contour in counts".format(self.other_observer))
self.ax_overplot.add_patch(other_proxy)
self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
self.ax_overplot.set_ylabel(label="Declination (J2000)",labelpad=-1)
self.fig_overplot.suptitle("HST/FOC UV polarization map of {0:s} overplotted \nwith {1:s} contour in counts.".format(obj,self.other_map[0].header['telescop']),wrap=True)
self.fig_overplot.suptitle("{0:s} polarisation map of {1:s} overplotted\nwith {2:s} contour in counts.".format(self.map_observer,obj,self.other_observer),wrap=True)
#Display pixel scale and North direction
fontprops = fm.FontProperties(size=16)
@@ -866,10 +884,10 @@ class overplot_chandra(align_maps):
pol_sc = AnchoredSizeBar(self.ax_overplot.transData, self.vec_scale, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.ax_overplot.add_artist(pol_sc)
self.cr_map, = self.ax_overplot.plot(*(self.wcs_map.celestial.wcs.crpix-(1.,1.)), 'r+')
self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1.,1.)), 'r+')
self.cr_other, = self.ax_overplot.plot(*(other_wcs.celestial.wcs.crpix-(1.,1.)), 'g+', transform=self.ax_overplot.get_transform(other_wcs))
h,l = self.ax_overplot.get_legend_handles_labels()
h[np.argmax([li=='HST/FOC polarisation map' for li in l])] = FancyArrowPatch((0,0),(0,1),arrowstyle='-',fc='w',ec='k',lw=2)
h[np.argmax([li=="{0:s} polarisation map".format(self.map_observer) for li in l])] = FancyArrowPatch((0,0),(0,1),arrowstyle='-',fc='w',ec='k',lw=2)
self.legend = self.ax_overplot.legend(handles=h,labels=l,bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.)
if not(savename is None):
@@ -893,7 +911,7 @@ class overplot_pol(align_maps):
"""
def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2., savename=None, **kwargs):
self.Stokes_UV = self.map
self.wcs_UV = self.wcs_map
self.wcs_UV = self.map_wcs
#Get Data
obj = self.Stokes_UV[0].header['targname']
stkI = self.Stokes_UV['I_STOKES'].data
@@ -902,7 +920,7 @@ class overplot_pol(align_maps):
pol_err = self.Stokes_UV['POL_DEG_ERR'].data
pang = self.Stokes_UV['POL_ANG'].data
other_data = self.other_map[0].data
other_data = self.other_data
#Compute SNR and apply cuts
pol[pol == 0.] = np.nan
@@ -914,16 +932,16 @@ class overplot_pol(align_maps):
pol[SNRi < SNRi_cut] = np.nan
plt.rcParams.update({'font.size': 16})
self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11,10), subplot_kw=dict(projection=self.wcs_other))
self.fig_overplot.subplots_adjust(hspace=0,wspace=0,bottom=-0.02,left=0.12,top=0.90,right=1.02)
self.fig_overplot, self.ax_overplot = plt.subplots(figsize=(11,10), subplot_kw=dict(projection=self.other_wcs))
self.fig_overplot.subplots_adjust(hspace=0,wspace=0,bottom=0.1,left=0.1,top=0.80,right=1.02)
self.ax_overplot.set_xlabel(label="Right Ascension (J2000)")
self.ax_overplot.set_ylabel(label="Declination (J2000)",labelpad=-1)
self.fig_overplot.suptitle("{0:s} observation from {1:s} overplotted with polarization vectors and Stokes I contours from HST/FOC".format(obj,self.other_map[0].header['telescop']),wrap=True)
self.fig_overplot.suptitle("{0:s} observation from {1:s} overplotted with polarisation vectors and Stokes I contours from {2:s}".format(obj,self.other_observer,self.map_observer),wrap=True)
#Display "other" intensity map
vmin, vmax = other_data[other_data > 0.].max()/1e3*self.other_convert, other_data[other_data > 0.].max()*self.other_convert
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["norm",LogNorm(vmin,vmax)]]]]:
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["vmin",vmin],["vmax",vmax]]]]:
try:
test = kwargs[key]
except KeyError:
@@ -935,32 +953,32 @@ class overplot_pol(align_maps):
else:
self.ax_overplot.set_facecolor('white')
font_color="black"
self.im = self.ax_overplot.imshow(other_data*self.other_convert, alpha=1., label="{0:s} observation".format(self.other_map[0].header['telescop']), **kwargs)
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=80, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_map_unit))
self.im = self.ax_overplot.imshow(other_data*self.other_convert, alpha=1., label="{0:s} observation".format(self.other_observer), **kwargs)
self.cbar = self.fig_overplot.colorbar(self.im, ax=self.ax_overplot, aspect=80, shrink=0.75, pad=0.025, label=r"$F_{{\lambda}}$ [{0:s}]".format(self.other_unit))
#Display full size polarization vectors
#Display full size polarisation vectors
if vec_scale is None:
self.vec_scale = 2.
pol[np.isfinite(pol)] = 1./2.
else:
self.vec_scale = vec_scale
step_vec = 1
px_scale = self.wcs_other.wcs.get_cdelt()[0]/self.wcs_UV.wcs.get_cdelt()[0]
px_scale = self.other_wcs.wcs.get_cdelt()[0]/self.wcs_UV.wcs.get_cdelt()[0]
self.X, self.Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0]))
self.U, self.V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.)
self.Q = self.ax_overplot.quiver(self.X[::step_vec,::step_vec],self.Y[::step_vec,::step_vec],self.U[::step_vec,::step_vec],self.V[::step_vec,::step_vec],units='xy',angles='uv',scale=px_scale/self.vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1/px_scale,linewidth=0.5,color='white',edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV),label="HST/FOC polarisation map")
self.Q = self.ax_overplot.quiver(self.X[::step_vec,::step_vec],self.Y[::step_vec,::step_vec],self.U[::step_vec,::step_vec],self.V[::step_vec,::step_vec],units='xy',angles='uv',scale=px_scale/self.vec_scale,scale_units='xy',pivot='mid',headwidth=0.,headlength=0.,headaxislength=0.,width=0.1/px_scale,linewidth=0.5,color='white',edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV),label="{0:s} polarisation map".format(self.map_observer))
#Display Stokes I as contours
if levels is None:
levels = np.logspace(np.log(3)/np.log(10),2.,5)/100.*np.max(stkI[stkI > 0.])*self.convert_flux
cont_stkI = self.ax_overplot.contour(stkI*self.convert_flux, levels=levels, colors='grey', alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV))
self.ax_overplot.clabel(cont_stkI, inline=True, fontsize=5)
cont_proxy = Rectangle((0,0),1,1,fc='w',ec=cont_stkI.collections[0].get_edgecolor()[0], label="HST/FOC Stokes I contour")
levels = np.logspace(np.log(3)/np.log(10),2.,5)/100.*np.max(stkI[stkI > 0.])*self.map_convert
cont_stkI = self.ax_overplot.contour(stkI*self.map_convert, levels=levels, colors='grey', alpha=0.75, transform=self.ax_overplot.get_transform(self.wcs_UV))
#self.ax_overplot.clabel(cont_stkI, inline=True, fontsize=5)
cont_proxy = Rectangle((0,0),1,1,fc='w',ec=cont_stkI.collections[0].get_edgecolor()[0], label="{0:s} Stokes I contour".format(self.map_observer))
self.ax_overplot.add_patch(cont_proxy)
#Display pixel scale and North direction
fontprops = fm.FontProperties(size=16)
px_size = self.wcs_other.wcs.get_cdelt()[0]*3600.
px_size = self.other_wcs.wcs.get_cdelt()[0]*3600.
px_sc = AnchoredSizeBar(self.ax_overplot.transData, 1./px_size, '1 arcsec', 3, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.ax_overplot.add_artist(px_sc)
north_dir = AnchoredDirectionArrows(self.ax_overplot.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=font_color, arrow_props={'ec': 'k', 'fc': 'w', 'alpha': 1,'lw': 0.5})
@@ -968,18 +986,18 @@ class overplot_pol(align_maps):
pol_sc = AnchoredSizeBar(self.ax_overplot.transData, self.vec_scale/px_scale, r"$P$= 100%", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color=font_color, fontproperties=fontprops)
self.ax_overplot.add_artist(pol_sc)
self.cr_map, = self.ax_overplot.plot(*(self.wcs_map.celestial.wcs.crpix-(1.,1.)), 'r+', transform=self.ax_overplot.get_transform(self.wcs_UV))
self.cr_other, = self.ax_overplot.plot(*(self.wcs_other.celestial.wcs.crpix-(1.,1.)), 'g+')
self.cr_map, = self.ax_overplot.plot(*(self.map_wcs.celestial.wcs.crpix-(1.,1.)), 'r+', transform=self.ax_overplot.get_transform(self.wcs_UV))
self.cr_other, = self.ax_overplot.plot(*(self.other_wcs.celestial.wcs.crpix-(1.,1.)), 'g+')
if "PHOTPLAM" in list(self.other_map[0].header.keys()):
self.legend_title = r"{0:s} image at $\lambda$ = {1:.0f} $\AA$".format(self.other_map[0].header['telescop'],float(self.other_map[0].header['photplam']))
elif "CRVAL3" in list(self.other_map[0].header.keys()):
self.legend_title = "{0:s} image at {1:.2f} GHz".format(self.other_map[0].header['telescop'],float(self.other_map[0].header['crval3'])*1e-9)
if "PHOTPLAM" in list(self.other_header.keys()):
self.legend_title = r"{0:s} image at $\lambda$ = {1:.0f} $\AA$".format(self.other_map_observer,float(self.other_header['photplam']))
elif "CRVAL3" in list(self.other_header.keys()):
self.legend_title = "{0:s} image at {1:.2f} GHz".format(self.other_observer,float(self.other_header['crval3'])*1e-9)
else:
self.legend_title = r"{0:s} image".format(self.other_map[0].header['telescop'])
self.legend_title = r"{0:s} image".format(self.other_observer)
h,l = self.ax_overplot.get_legend_handles_labels()
h[np.argmax([li=='HST/FOC polarisation map' for li in l])] = FancyArrowPatch((0,0),(0,1),arrowstyle='-',fc='w',ec='k',lw=2)
h[np.argmax([li=="{0:s} polarisation map".format(self.map_observer) for li in l])] = FancyArrowPatch((0,0),(0,1),arrowstyle='-',fc='w',ec='k',lw=2)
self.legend = self.ax_overplot.legend(handles=h,labels=l,bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left', mode="expand", borderaxespad=0.)
if not(savename is None):
@@ -999,7 +1017,7 @@ class overplot_pol(align_maps):
if position == 'center':
position = np.array(self.X.shape)/2.
if type(position) == SkyCoord:
position = self.wcs_other.world_to_pixel(position)
position = self.other_wcs.world_to_pixel(position)
u, v = pol_deg*np.cos(np.radians(pol_ang)+np.pi/2.), pol_deg*np.sin(np.radians(pol_ang)+np.pi/2.)
for key, value in [["scale",[["scale",self.vec_scale]]], ["width",[["width",0.1]]], ["color",[["color",'k']]]]:
@@ -1149,9 +1167,9 @@ class crop_map(object):
self.data = deepcopy(self.hdul[0].data)
try:
self.convert_flux = self.header['photflam']
self.map_convert = self.header['photflam']
except KeyError:
self.convert_flux = 1.
self.map_convert = 1.
try:
self.kwargs = kwargs
except AttributeError:
@@ -1179,7 +1197,7 @@ class crop_map(object):
self.rect_selector = RectangleSelector(self.ax, self.onselect_crop,
button=[1])
self.embedded = True
self.display(self.data, self.wcs, self.convert_flux, **self.kwargs)
self.display(self.data, self.wcs, self.map_convert, **self.kwargs)
self.extent = np.array([0.,self.data.shape[0],0., self.data.shape[1]])
self.center = np.array(self.data.shape)/2
@@ -1193,7 +1211,7 @@ class crop_map(object):
if wcs is None:
wcs = self.wcs
if convert_flux is None:
convert_flux = self.convert_flux
convert_flux = self.map_convert
if kwargs is None:
kwargs = self.kwargs
else:
@@ -1317,7 +1335,7 @@ class crop_map(object):
class crop_Stokes(crop_map):
"""
Class to interactively crop a polarization map to desired Region of Interest.
Class to interactively crop a polarisation map to desired Region of Interest.
Inherit from crop_map.
"""
def apply_crop(self,event):
@@ -1396,10 +1414,10 @@ class crop_Stokes(crop_map):
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for dataset in self.hdul_crop:
dataset.header['P_int'] = (P_diluted, 'Integrated polarization degree')
dataset.header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarization degree error')
dataset.header['PA_int'] = (PA_diluted, 'Integrated polarization angle')
dataset.header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
dataset.header['P_int'] = (P_diluted, 'Integrated polarisation degree')
dataset.header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarisation degree error')
dataset.header['PA_int'] = (PA_diluted, 'Integrated polarisation angle')
dataset.header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarisation angle error')
self.fig.canvas.draw_idle()
@property
@@ -1578,7 +1596,7 @@ class aperture(object):
class pol_map(object):
"""
Class to interactively study polarization maps.
Class to interactively study polarisation maps.
"""
def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=30., flux_lim=None, selection=None):
@@ -1598,19 +1616,17 @@ class pol_map(object):
#Get data
self.targ = self.Stokes[0].header['targname']
self.pivot_wav = self.Stokes[0].header['photplam']
self.convert_flux = self.Stokes[0].header['photflam']
self.map_convert = self.Stokes[0].header['photflam']
#Create figure
plt.rcParams.update({'font.size': 10})
self.fig = plt.figure(figsize=(10,10))
self.fig.subplots_adjust(hspace=0, wspace=0, right=0.88)
self.ax = self.fig.add_subplot(111,projection=self.wcs)
self.fig, self.ax = plt.subplots(figsize=(10,10),subplot_kw=dict(projection=self.wcs))
self.fig.subplots_adjust(hspace=0, wspace=0, right=1.02)
self.ax_cosmetics()
self.cbar_ax = self.fig.add_axes([0.925, 0.13, 0.01, 0.74])
#Display selected data (Default to total flux)
self.display()
#Display polarization vectors in SNR_cut
#Display polarisation vectors in SNR_cut
self.pol_vector()
#Display integrated values in ROI
self.pol_int()
@@ -1869,7 +1885,7 @@ class pol_map(object):
dump_list = []
for i in range(shape[0]):
for j in range(shape[1]):
dump_list.append([x[i,j], y[i,j], self.I[i,j]*self.convert_flux, self.Q[i,j]*self.convert_flux, self.U[i,j]*self.convert_flux, P[i,j], PA[i,j]])
dump_list.append([x[i,j], y[i,j], self.I[i,j]*self.map_convert, self.Q[i,j]*self.map_convert, self.U[i,j]*self.map_convert, P[i,j], PA[i,j]])
self.data_dump = np.array(dump_list)
b_dump.on_clicked(dump)
@@ -1943,7 +1959,7 @@ class pol_map(object):
@property
def wcs(self):
return deepcopy(WCS(self.Stokes[0].header))
return WCS(self.Stokes[0].header).celestial
@property
def I(self):
return self.Stokes['I_STOKES'].data
@@ -1983,16 +1999,11 @@ class pol_map(object):
def ax_cosmetics(self, ax=None):
if ax is None:
ax = self.ax
ax.set_facecolor('black')
ax.set(aspect='equal',fc='black')
ax.coords.grid(True, color='white', ls='dotted', alpha=0.5)
ax.coords[0].set_axislabel('Right Ascension (J2000)')
ax.coords[0].set_axislabel_position('t')
ax.coords[0].set_ticklabel_position('t')
ax.coords[1].set_axislabel('Declination (J2000)')
ax.coords[1].set_axislabel_position('l')
ax.coords[1].set_ticklabel_position('l')
ax.axis('equal')
ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)',labelpad=-1)
#Display scales and orientation
fontprops = fm.FontProperties(size=14)
@@ -2017,7 +2028,7 @@ class pol_map(object):
if flux_lim is None:
flux_lim = self.flux_lim
if self.display_selection.lower() in ['total_flux']:
self.data = self.I*self.convert_flux
self.data = self.I*self.map_convert
if flux_lim is None:
vmin, vmax = 1./2.*np.median(self.data[self.data > 0.]), np.max(self.data[self.data > 0.])
else:
@@ -2025,9 +2036,9 @@ class pol_map(object):
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
self.data = self.I*self.map_convert*self.P
if flux_lim is None:
vmin, vmax = 1./2.*np.median(self.I[self.I > 0.]*self.convert_flux), np.max(self.I[self.I > 0.]*self.convert_flux)
vmin, vmax = 1./2.*np.median(self.I[self.I > 0.]*self.map_convert), np.max(self.I[self.I > 0.]*self.map_convert)
else:
vmin, vmax = flux_lim
norm = LogNorm(vmin, vmax)
@@ -2058,13 +2069,15 @@ class pol_map(object):
fig = self.fig
if ax is None:
ax = self.ax
if hasattr(self, 'cbar'):
self.cbar.remove()
if hasattr(self, 'im'):
self.im.remove()
if not norm is None:
self.im = ax.imshow(self.data, norm=norm, aspect='equal', cmap='inferno')
else:
self.im = ax.imshow(self.data, vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno')
self.cbar = plt.colorbar(self.im, cax=self.cbar_ax, label=label)
self.cbar = fig.colorbar(self.im, ax=ax, aspect=50, shrink=0.75, pad=0.025, label=label)
fig.canvas.draw_idle()
return self.im
else:
@@ -2185,15 +2198,13 @@ class pol_map(object):
ax = self.ax
if hasattr(self, 'an_int'):
self.an_int.remove()
self.an_int = 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.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.93), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
#self.an_int = 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.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.)+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100.,np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut,np.ceil(PA_cut_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.85), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
self.an_int = 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.map_convert,I_reg_err*self.map_convert,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')], verticalalignment='top', horizontalalignment='left')
if not self.region is None:
self.cont = ax.contour(self.region.astype(float),levels=[0.5], colors='white', linewidths=0.8)
fig.canvas.draw_idle()
return self.an_int
else:
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.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.94), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
#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.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.)+"\n"+r"$P^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_cut*100.,np.ceil(P_cut_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{cut}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_cut,np.ceil(PA_cut_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 0.90), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')])
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.map_convert,I_reg_err*self.map_convert,2))+"\n"+r"$P^{{int}}$ = {0:.1f} $\pm$ {1:.1f} %".format(P_reg*100.,np.ceil(P_reg_err*1000.)/10.)+"\n"+r"$\theta_{{P}}^{{int}}$ = {0:.1f} $\pm$ {1:.1f} °".format(PA_reg,np.ceil(PA_reg_err*10.)/10.), color='white', fontsize=12, xy=(0.01, 1.00), xycoords='axes fraction',path_effects=[pe.withStroke(linewidth=0.5,foreground='k')], verticalalignment='top', horizontalalignment='left')
if not self.region is None:
ax.contour(self.region.astype(float),levels=[0.5], colors='white', linewidths=0.8)
fig.canvas.draw_idle()

52
src/lib/reduction.py
View File

@@ -30,7 +30,7 @@ prototypes :
Compute Stokes parameters I, Q and U and their respective correlated errors from data_array.
- compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers) -> P, debiased_P, s_P, s_P_P, PA, s_PA, s_PA_P
Compute polarization degree (in %) and angle (in degree) and their respective errors.
Compute polarisation degree (in %) and angle (in degree) and their respective errors.
- rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, ang, SNRi_cut) -> I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers
Rotate I, Q, U given an angle in degrees using scipy functions.
@@ -992,7 +992,7 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None,
FWHM=FWHM, scale=scale, smoothing=smoothing)
else:
# Sum on each polarization filter.
# Sum on each polarisation filter.
pol0_t = np.sum([header['exptime'] for header in headers0])
pol60_t = np.sum([header['exptime'] for header in headers60])
pol120_t = np.sum([header['exptime'] for header in headers120])
@@ -1101,10 +1101,10 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
total intensity
Q_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for
vertical/horizontal linear polarization intensity
vertical/horizontal linear polarisation intensity
U_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for
+45/-45deg linear polarization intensity
+45/-45deg linear polarisation intensity
Stokes_cov : numpy.ndarray
Covariance matrix of the Stokes parameters I, Q, U.
"""
@@ -1257,17 +1257,17 @@ def compute_Stokes(data_array, error_array, data_mask, headers,
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in headers:
header['P_int'] = (P_diluted, 'Integrated polarization degree')
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'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
header['P_int'] = (P_diluted, 'Integrated polarisation degree')
header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarisation degree error')
header['PA_int'] = (PA_diluted, 'Integrated polarisation angle')
header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarisation angle error')
return I_stokes, Q_stokes, U_stokes, Stokes_cov
def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
"""
Compute the polarization degree (in %) and angle (in deg) and their
Compute the polarisation degree (in %) and angle (in deg) and their
respective errors from given Stokes parameters.
----------
Inputs:
@@ -1276,10 +1276,10 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
total intensity
Q_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for
vertical/horizontal linear polarization intensity
vertical/horizontal linear polarisation intensity
U_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for
+45/-45deg linear polarization intensity
+45/-45deg linear polarisation intensity
Stokes_cov : numpy.ndarray
Covariance matrix of the Stokes parameters I, Q, U.
headers : header list
@@ -1287,21 +1287,21 @@ def compute_pol(I_stokes, Q_stokes, U_stokes, Stokes_cov, headers):
----------
Returns:
P : numpy.ndarray
Image (2D floats) containing the polarization degree (in %).
Image (2D floats) containing the polarisation degree (in %).
debiased_P : numpy.ndarray
Image (2D floats) containing the debiased polarization degree (in %).
Image (2D floats) containing the debiased polarisation degree (in %).
s_P : numpy.ndarray
Image (2D floats) containing the error on the polarization degree.
Image (2D floats) containing the error on the polarisation degree.
s_P_P : numpy.ndarray
Image (2D floats) containing the Poisson noise error on the
polarization degree.
polarisation degree.
PA : numpy.ndarray
Image (2D floats) containing the polarization angle.
Image (2D floats) containing the polarisation angle.
s_PA : numpy.ndarray
Image (2D floats) containing the error on the polarization angle.
Image (2D floats) containing the error on the polarisation angle.
s_PA_P : numpy.ndarray
Image (2D floats) containing the Poisson noise error on the
polarization angle.
polarisation angle.
new_headers : header list
Updated list of headers corresponding to the reduced images accounting
for the new orientation angle.
@@ -1374,10 +1374,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
total intensity
Q_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for
vertical/horizontal linear polarization intensity
vertical/horizontal linear polarisation intensity
U_stokes : numpy.ndarray
Image (2D floats) containing the Stokes parameters accounting for
+45/-45deg linear polarization intensity
+45/-45deg linear polarisation intensity
Stokes_cov : numpy.ndarray
Covariance matrix of the Stokes parameters I, Q, U.
data_mask : numpy.ndarray
@@ -1399,10 +1399,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
accounting for total intensity
new_Q_stokes : numpy.ndarray
Rotated mage (2D floats) containing the rotated Stokes parameters
accounting for vertical/horizontal linear polarization intensity
accounting for vertical/horizontal linear polarisation intensity
new_U_stokes : numpy.ndarray
Rotated image (2D floats) containing the rotated Stokes parameters
accounting for +45/-45deg linear polarization intensity.
accounting for +45/-45deg linear polarisation intensity.
new_Stokes_cov : numpy.ndarray
Updated covariance matrix of the Stokes parameters I, Q, U.
new_headers : header list
@@ -1516,10 +1516,10 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
PA_diluted_err = (90./(np.pi*(Q_diluted**2 + U_diluted**2)))*np.sqrt(U_diluted**2*Q_diluted_err**2 + Q_diluted**2*U_diluted_err**2 - 2.*Q_diluted*U_diluted*QU_diluted_err)
for header in new_headers:
header['P_int'] = (P_diluted, 'Integrated polarization degree')
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'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarization angle error')
header['P_int'] = (P_diluted, 'Integrated polarisation degree')
header['P_int_err'] = (np.ceil(P_diluted_err*1000.)/1000., 'Integrated polarisation degree error')
header['PA_int'] = (PA_diluted, 'Integrated polarisation angle')
header['PA_int_err'] = (np.ceil(PA_diluted_err*10.)/10., 'Integrated polarisation angle error')
return new_I_stokes, new_Q_stokes, new_U_stokes, new_Stokes_cov, new_data_mask, new_headers

45
src/overplot_IC5063.py
View File

@@ -15,35 +15,34 @@ Stokes_UV = fits.open("./data/IC5063/5918/IC5063_FOC_b0.10arcsec_c0.20arcsec.fit
#Stokes_S2 = fits.open("./data/IC5063/POLARIZATION_COMPARISON/S2_rot_crop.fits")
Stokes_IR = fits.open("./data/IC5063/IR/u2e65g01t_c0f_rot.fits")
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
##levelsMorganti = np.array([1.,2.,3.,8.,16.,32.,64.,128.])
#levelsMorganti = np.logspace(0.,1.97,5)/100.
#
#levels18GHz = levelsMorganti*Stokes_18GHz[0].data.max()
#A = overplot_radio(Stokes_UV, Stokes_18GHz)
#A.plot(levels=levels18GHz, SNRp_cut=1.0, SNRi_cut=10.0, savename='./plots/IC5063/18GHz_overplot_forced.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
#A.plot(levels=levels18GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/18GHz_overplot_forced.pdf',vec_scale=None)
##
#levels24GHz = levelsMorganti*Stokes_24GHz[0].data.max()
#B = overplot_radio(Stokes_UV, Stokes_24GHz)
#B.plot(levels=levels24GHz, SNRp_cut=1.0, SNRi_cut=10.0, savename='./plots/IC5063/24GHz_overplot_forced.png')
#
#levels103GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_103GHz[0].data[Stokes_103GHz[0].data > 0.]))
#B.plot(levels=levels24GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/24GHz_overplot_forced.pdf',vec_scale=None)
##
#levels103GHz = levelsMorganti*Stokes_103GHz[0].data.max()
#C = overplot_radio(Stokes_UV, Stokes_103GHz)
#C.plot(levels=levels103GHz, SNRp_cut=1.0, SNRi_cut=10.0, savename='./plots/IC5063/103GHz_overplot_forced.png')
#
#levels229GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_229GHz[0].data[Stokes_229GHz[0].data > 0.]))
#C.plot(levels=levels103GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/103GHz_overplot_forced.pdf',vec_scale=None)
##
#levels229GHz = levelsMorganti*Stokes_229GHz[0].data.max()
#D = overplot_radio(Stokes_UV, Stokes_229GHz)
#D.plot(levels=levels229GHz, SNRp_cut=1.0, SNRi_cut=10.0, savename='./plots/IC5063/229GHz_overplot_forced.png')
#
#levels357GHz = np.linspace(1,99,11)/100.*np.max(deepcopy(Stokes_357GHz[0].data[Stokes_357GHz[0].data > 0.]))
#D.plot(levels=levels229GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/229GHz_overplot_forced.pdf',vec_scale=None)
##
#levels357GHz = levelsMorganti*Stokes_357GHz[0].data.max()
#E = overplot_radio(Stokes_UV, Stokes_357GHz)
#E.plot(levels=levels357GHz, SNRp_cut=1.0, SNRi_cut=10.0, savename='./plots/IC5063/357GHz_overplot_forced.png')
#
#E.plot(levels=levels357GHz, SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/357GHz_overplot_forced.pdf',vec_scale=None)
##
#F = overplot_pol(Stokes_UV, Stokes_S2)
#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot_forced.png', norm=LogNorm(vmin=5e-20,vmax=5e-18))
#F.plot(SNRp_cut=3.0, SNRi_cut=80.0, savename='./plots/IC5063/S2_overplot_forced.pdf', norm=LogNorm(vmin=5e-20,vmax=5e-18))
G = overplot_pol(Stokes_UV, Stokes_IR, cmap='inferno')
G.plot(SNRp_cut=1.0, SNRi_cut=10.0, vec_scale=3., savename='./plots/IC5063/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
G.plot(SNRp_cut=2.0, SNRi_cut=10.0, savename='./plots/IC5063/IR_overplot_forced.pdf',vec_scale=None,norm=LogNorm(Stokes_IR[0].data.max()*Stokes_IR[0].header['photflam']/1e3,Stokes_IR[0].data.max()*Stokes_IR[0].header['photflam']),cmap='inferno_r')
#data_folder1 = "./data/M87/POS1/"
#plots_folder1 = "./plots/M87/POS1/"
@@ -56,7 +55,7 @@ G.plot(SNRp_cut=1.0, SNRi_cut=10.0, vec_scale=3., savename='./plots/IC5063/IR_ov
#H = align_pol(np.array([M87_1_95,M87_1_96,M87_1_97,M87_1_98,M87_1_99]), norm=LogNorm())
#H.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder1+'animated_loop/'+basename1, norm=LogNorm())
#command("convert -delay 50 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif".format(plots_folder1, basename1))
#command("convert -delay 50 -loop 0 {0:s}animated_loop/{1:s}*.pdf {0:s}animated_loop/{1:s}.gif".format(plots_folder1, basename1))
#data_folder3 = "./data/M87/POS3/"
#plots_folder3 = "./plots/M87/POS3/"
@@ -69,4 +68,4 @@ G.plot(SNRp_cut=1.0, SNRi_cut=10.0, vec_scale=3., savename='./plots/IC5063/IR_ov
#I = align_pol(np.array([M87_3_95,M87_3_96,M87_3_97,M87_3_98,M87_3_99]), norm=LogNorm())
#I.plot(SNRp_cut=5.0, SNRi_cut=50.0, savename=plots_folder3+'animated_loop/'+basename3, norm=LogNorm())
#command("convert -delay 20 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif".format(plots_folder3, basename3))
#command("convert -delay 20 -loop 0 {0:s}animated_loop/{1:s}*.pdf {0:s}animated_loop/{1:s}.gif".format(plots_folder3, basename3))