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

do some reduction over M87 POS1 & POS3 observations from 95 to 99

Browse Source
This commit is contained in:
Thibault Barnouin
2022-08-16 17:24:13 +02:00
parent e61984ae3d
commit d3c11eb4af
78 changed files with 57 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files

18
src/FOC_reduction.py
View File

@@ -90,22 +90,22 @@ from astropy.wcs import WCS
#globals()['plots_folder'] = "../plots/3C273_x0u20/"
#BEWARE: 5 observations separated by 1 year each (1995, 1996, 1997, 1998, 1999)
#globals()['data_folder'] = "../data/M87/POS1/"
globals()['data_folder'] = "../data/M87/POS1/"
#globals()['infiles'] = ['x2py010ct_c0f.fits','x2py010dt_c0f.fits','x2py010et_c0f.fits','x2py010ft_c0f.fits'] #1995
#globals()['infiles'] = ['x3be010ct_c0f.fits','x3be010dt_c0f.fits','x3be010et_c0f.fits','x3be010ft_c0f.fits'] #1996
#globals()['infiles'] = ['x43r010km_c0f.fits','x43r010mm_c0f.fits','x43r010om_c0f.fits','x43r010rm_c0f.fits'] #1997
#globals()['infiles'] = ['x43r110kr_c0f.fits','x43r110mr_c0f.fits','x43r110or_c0f.fits','x43r110rr_c0f.fits'] #1998
#globals()['infiles'] = ['x43r210kr_c0f.fits','x43r210mr_c0f.fits','x43r210or_c0f.fits','x43r210rr_c0f.fits'] #1999
#globals()['plots_folder'] = "../plots/M87/POS1/"
globals()['infiles'] = ['x43r210kr_c0f.fits','x43r210mr_c0f.fits','x43r210or_c0f.fits','x43r210rr_c0f.fits'] #1999
globals()['plots_folder'] = "../plots/M87/POS1/"
#BEWARE: 5 observations separated by 1 year each (1995, 1996, 1997, 1998, 1999)
globals()['data_folder'] = "../data/M87/POS3/"
globals()['infiles'] = ['x2py030at_c0f.fits','x2py030bt_c0f.fits','x2py030ct_c0f.fits','x2py0309t_c0f.fits'] #1995
#globals()['data_folder'] = "../data/M87/POS3/"
#globals()['infiles'] = ['x2py030at_c0f.fits','x2py030bt_c0f.fits','x2py030ct_c0f.fits','x2py0309t_c0f.fits'] #1995
#globals()['infiles'] = ['x3be030at_c0f.fits','x3be030bt_c0f.fits','x3be030ct_c0f.fits','x3be0309t_c0f.fits'] #1996
#globals()['infiles'] = ['x43r030em_c0f.fits','x43r030gm_c0f.fits','x43r030im_c0f.fits','x43r030lm_c0f.fits'] #1997
#globals()['infiles'] = ['x43r130er_c0f.fits','x43r130fr_c0f.fits','x43r130ir_c0f.fits','x43r130lr_c0f.fits'] #1998
#globals()['infiles'] = ['x43r230er_c0f.fits','x43r230fr_c0f.fits','x43r230ir_c0f.fits','x43r230lr_c0f.fits'] #1999
globals()['plots_folder'] = "../plots/M87/POS3/"
#globals()['plots_folder'] = "../plots/M87/POS3/"
def main():
@@ -136,7 +136,7 @@ def main():
display_data = False
# Smoothing
smoothing_function = 'combine' #gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine
smoothing_FWHM = 0.07 #If None, no smoothing is done
smoothing_FWHM = 0.10 #If None, no smoothing is done
smoothing_scale = 'arcsec' #pixel or arcsec
# Rotation
rotate_stokes = True #rotation to North convention can give erroneous results
@@ -145,8 +145,8 @@ def main():
crop = False #Crop to desired ROI
final_display = True
# Polarization map output
figname = 'M87_POS3_1995_FOC' #target/intrument name
figtype = '_combine_FWHM005' #additionnal informations
figname = 'M87_POS1_1999_FOC' #target/intrument name
figtype = '_combine_FWHM010' #additionnal informations
SNRp_cut = 3. #P measurments with SNR>3
SNRi_cut = 30. #I measurments with SNR>30, which implies an uncertainty in P of 4.7%.
step_vec = 0 #plot all vectors in the array. if step_vec = 2, then every other vector will be plotted

37
src/lib/plots.py
View File

@@ -526,6 +526,11 @@ class align_maps(object):
self.axreset = self.fig.add_axes([0.60, 0.01, 0.1, 0.04])
self.breset = Button(self.axreset, 'Leave as is')
self.breset.label.set_fontsize(8)
self.enter = self.fig.canvas.mpl_connect('key_press_event', self.on_key)
def on_key(self, event):
if event.key.lower() == "enter":
self.on_close_align(event)
def get_aligned_wcs(self):
return self.wcs_map, self.wcs_other
@@ -556,10 +561,16 @@ class align_maps(object):
self.aligned = True
def apply_align(self, event):
self.wcs_map.wcs.crpix = np.array(self.cr_map.get_data())
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]
else:
self.wcs_map.wcs.crpix = np.array(self.cr_map.get_data())
if np.array(self.cr_other.get_data()).shape == (2,1):
self.wcs_other.wcs.crpix = np.array(self.cr_other.get_data())[:,0]
else:
self.wcs_other.wcs.crpix = np.array(self.cr_other.get_data())
self.wcs_map.wcs.crval = np.array(self.wcs_map.pixel_to_world_values(*self.wcs_map.wcs.crpix))
self.wcs_other.wcs.crpix = np.array(self.cr_other.get_data())
self.wcs_other.wcs.crval = self.wcs_map.wcs.crval
self.fig.canvas.draw_idle()
@@ -569,8 +580,9 @@ class align_maps(object):
self.aligned = True
def on_close_align(self, event):
self.aligned = True
#print(self.get_aligned_wcs())
if not self.aligned:
self.aligned = True
self.apply_align()
def align(self):
self.fig.canvas.draw()
@@ -808,8 +820,6 @@ class align_pol(object):
if not ax_lim is None:
lim = np.concatenate([wcs.world_to_pixel(ax_lim[i]) for i in range(len(ax_lim))])
x_lim, y_lim = lim[0::2], lim[1::2]
print(x_lim[0], y_lim[0], wcs.pixel_to_world(x_lim[0], y_lim[0]))
print(x_lim[1], y_lim[1], wcs.pixel_to_world(x_lim[1], y_lim[1]))
ax.set(xlim=x_lim,ylim=y_lim)
if v_lim is None:
@@ -820,6 +830,8 @@ class align_pol(object):
for key, value in [["cmap",[["cmap","inferno"]]], ["norm",[["vmin",vmin],["vmax",vmax]]]]:
try:
test = kwargs[key]
if str(type(test)) == "<class 'matplotlib.colors.LogNorm'>":
kwargs[key] = LogNorm(vmin, vmax)
except KeyError:
for key_i, val_i in value:
kwargs[key_i] = val_i
@@ -856,17 +868,16 @@ class align_pol(object):
def plot(self, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs):
while not self.aligned.all():
self.align()
vmin = np.min([np.min(curr_map[0].data[curr_map[0].data > 0.]) for curr_map in self.other_maps])
vmax = np.max([np.max(curr_map[0].data[curr_map[0].data > 0.]) for curr_map in self.other_maps])
vmin, vmax = np.min([vmin, np.min(self.ref_map[0].data[self.ref_map[0].data > 0.])]), np.max([vmax, np.max(self.ref_map[0].data[self.ref_map[0].data > 0.])])
eps = 1e-35
vmin = np.min([np.min(curr_map[0].data[curr_map[0].data > SNRi_cut*np.max([eps*np.ones(curr_map[0].data.shape),np.sqrt(curr_map[3].data[0,0])],axis=0)]) for curr_map in self.other_maps])/2.5
vmax = np.max([np.max(curr_map[0].data[curr_map[0].data > SNRi_cut*np.max([eps*np.ones(curr_map[0].data.shape),np.sqrt(curr_map[3].data[0,0])],axis=0)]) for curr_map in self.other_maps])
vmin = np.min([vmin, np.min(self.ref_map[0].data[self.ref_map[0].data > SNRi_cut*np.max([eps*np.ones(self.ref_map[0].data.shape),np.sqrt(self.ref_map[3].data[0,0])],axis=0)])])/2.5
vmax = np.max([vmax, np.max(self.ref_map[0].data[self.ref_map[0].data > SNRi_cut*np.max([eps*np.ones(self.ref_map[0].data.shape),np.sqrt(self.ref_map[3].data[0,0])],axis=0)])])
v_lim = np.array([vmin, vmax])
fig, ax = self.single_plot(self.ref_map, self.wcs, v_lim = v_lim, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename+'_0', **kwargs)
x_lim, y_lim = ax.get_xlim(), ax.get_ylim()
ax_lim = np.array([self.wcs.pixel_to_world(x_lim[i], y_lim[i]) for i in range(len(x_lim))])
print(x_lim[0], y_lim[0], ax_lim[0])
print(x_lim[1], y_lim[1], ax_lim[1])
for i, curr_map in enumerate(self.other_maps):
self.single_plot(curr_map, self.wcs_other[i], v_lim=v_lim, ax_lim=ax_lim, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename+'_'+str(i+1), **kwargs)

1
src/lib/reduction.py
View File

@@ -1503,7 +1503,6 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers,
new_wcs = WCS(header).deepcopy()
new_wcs.wcs.pc = np.dot(mrot, new_wcs.wcs.pc)
print(new_wcs.wcs.pc)
new_wcs.wcs.crpix = np.dot(mrot, new_wcs.wcs.crpix - old_center[::-1]) + new_center[::-1]
new_wcs.wcs.set()
for key, val in new_wcs.to_header().items():

48
src/overplot.py
View File

@@ -45,28 +45,28 @@ from matplotlib.colors import LogNorm
#G = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
#G.plot(SNRp_cut=3.0, SNRi_cut=60.0, savename='../plots/IC5063_x3nl030/IR_overplot_forced.png', norm=LogNorm(vmin=1e-17,vmax=5e-15), cmap='inferno_r')
#data_folder1 = "../data/M87/POS1/"
#plots_folder1 = "../plots/M87/POS1/"
#basename1 = "test"
#M87_1_95 = fits.open(data_folder1+"M87_POS1_1995_FOC_combine_FWHM005.fits")
#M87_1_96 = fits.open(data_folder1+"M87_POS1_1996_FOC_combine_FWHM005.fits")
#M87_1_97 = fits.open(data_folder1+"M87_POS1_1997_FOC_combine_FWHM005.fits")
#M87_1_98 = fits.open(data_folder1+"M87_POS1_1998_FOC_combine_FWHM005.fits")
#M87_1_99 = fits.open(data_folder1+"M87_POS1_1999_FOC_combine_FWHM005.fits")
data_folder1 = "../data/M87/POS1/"
plots_folder1 = "../plots/M87/POS1/"
basename1 = "M87_015_log_core"
M87_1_95 = fits.open(data_folder1+"M87_POS1_1995_FOC_combine_FWHM010.fits")
M87_1_96 = fits.open(data_folder1+"M87_POS1_1996_FOC_combine_FWHM010.fits")
M87_1_97 = fits.open(data_folder1+"M87_POS1_1997_FOC_combine_FWHM010.fits")
M87_1_98 = fits.open(data_folder1+"M87_POS1_1998_FOC_combine_FWHM010.fits")
M87_1_99 = fits.open(data_folder1+"M87_POS1_1999_FOC_combine_FWHM010.fits")
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=3.0, SNRi_cut=30.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))
#data_folder3 = "../data/M87/POS3/"
#plots_folder3 = "../plots/M87/POS3/"
#basename3 = "M87_005_log_star"
#M87_3_95 = fits.open(data_folder3+"M87_POS3_1995_FOC_combine_FWHM005.fits")
#M87_3_96 = fits.open(data_folder3+"M87_POS3_1996_FOC_combine_FWHM005.fits")
#M87_3_97 = fits.open(data_folder3+"M87_POS3_1997_FOC_combine_FWHM005.fits")
#M87_3_98 = fits.open(data_folder3+"M87_POS3_1998_FOC_combine_FWHM005.fits")
#M87_3_99 = fits.open(data_folder3+"M87_POS3_1999_FOC_combine_FWHM005.fits")
#
#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=3.0, SNRi_cut=30.0, savename=plots_folder1+'animated_loop/'+basename1, norm=LogNorm())
#command("convert -delay 20 -loop 0 {0:s}animated_loop/{1:s}*.png {0:s}animated_loop/{1:s}.gif && rm {0:s}animated_loop/{1:s}*.png".format(plots_folder1, basename1))
data_folder3 = "../data/M87/POS3/"
plots_folder3 = "../plots/M87/POS3/"
basename3 = "test"
M87_3_95 = fits.open(data_folder3+"M87_POS3_1995_FOC_combine_FWHM005.fits")
M87_3_96 = fits.open(data_folder3+"M87_POS3_1996_FOC_combine_FWHM005.fits")
M87_3_97 = fits.open(data_folder3+"M87_POS3_1997_FOC_combine_FWHM005.fits")
M87_3_98 = fits.open(data_folder3+"M87_POS3_1998_FOC_combine_FWHM005.fits")
M87_3_99 = fits.open(data_folder3+"M87_POS3_1999_FOC_combine_FWHM005.fits")
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=3.0, SNRi_cut=30.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 && rm {0:s}animated_loop/{1:s}*.png".format(plots_folder3, basename3))
#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=3.0, SNRi_cut=30.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))