From 845c6c30da144bd529553add5d9767440511f0fc Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Thu, 2 May 2024 12:36:17 +0200 Subject: [PATCH 01/11] Fix wrong orientation in some cases --- src/lib/fits.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib/fits.py b/src/lib/fits.py index 9116c83..ec90801 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -67,7 +67,7 @@ def get_obs_data(infiles, data_folder="", compute_flux=False): new_wcs.wcs.cdelt = new_cdelt for key, val in new_wcs.to_header().items(): header[key] = val - header['orientat'] = princ_angle(float(header['orientat'])) + # header['orientat'] = princ_angle(float(header['orientat'])) # force WCS for POL60 to have same pixel size as POL0 and POL120 is_pol60 = np.array([head['filtnam1'].lower() == 'pol60' for head in headers], dtype=bool) From a931e20f5f9ea4a3e852eb2c9fa4dd800a1a800f Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Thu, 2 May 2024 12:36:53 +0200 Subject: [PATCH 02/11] Fix MAST target query --- src/lib/query.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/lib/query.py b/src/lib/query.py index 421e55a..5c244c5 100755 --- a/src/lib/query.py +++ b/src/lib/query.py @@ -16,14 +16,6 @@ def divide_proposal(products): """ Divide observation in proposals by time or filter """ - for pid in np.unique(products['Proposal ID']): - obs = products[products['Proposal ID'] == pid].copy() - close_date = np.unique(np.array([TimeDelta(np.abs(obs['Start'][0].unix-date.unix), format='sec') - < 7.*u.d for date in obs['Start']], dtype=bool), axis=0) - if len(close_date) > 1: - for date in close_date: - products['Proposal ID'][np.any([products['Dataset'] == dataset for dataset in obs['Dataset'][date]], axis=0) - ] = "_".join([obs['Proposal ID'][date][0], str(obs['Start'][date][0])[:10]]) for pid in np.unique(products['Proposal ID']): obs = products[products['Proposal ID'] == pid].copy() same_filt = np.unique(np.array(np.sum([obs['Filters'][:, 1:] == filt[1:] for filt in obs['Filters']], axis=2) < 3, dtype=bool), axis=0) @@ -31,6 +23,14 @@ def divide_proposal(products): for filt in same_filt: products['Proposal ID'][np.any([products['Dataset'] == dataset for dataset in obs['Dataset'][filt]], axis=0)] = "_".join( [obs['Proposal ID'][filt][0], "_".join([fi for fi in obs['Filters'][filt][0][1:] if fi[:-1] != "CLEAR"])]) + for pid in np.unique(products['Proposal ID']): + obs = products[products['Proposal ID'] == pid].copy() + close_date = np.unique([[np.abs(TimeDelta(obs['Start'][i].unix-date.unix, format='sec')) + < 7.*u.d for i in range(len(obs))] for date in obs['Start']], axis=0) + if len(close_date) > 1: + for date in close_date: + products['Proposal ID'][np.any([products['Dataset'] == dataset for dataset in obs['Dataset'][date]], axis=0) + ] = "_".join([obs['Proposal ID'][date][0], str(obs['Start'][date][0])[:10]]) return products @@ -147,8 +147,6 @@ def get_product_list(target=None, proposal_id=None): for prod in products: prod['target_name'] = observations['target_name'][observations['obsid'] == prod['obsID']][0] tab = unique(products, ['target_name', 'proposal_id']) - if len(tab) > 1 and np.all(tab['target_name'] == tab['target_name'][0]): - target = tab['target_name'][0] products["Obs"] = [np.argmax(np.logical_and(tab['proposal_id'] == data['proposal_id'], tab['target_name'] == data['target_name']))+1 for data in products] return target, products From cbe6aeda4524f93daeb02318d990471df5155f4e Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Thu, 2 May 2024 12:37:40 +0200 Subject: [PATCH 03/11] Uniform plots and better handled maps for plots --- src/FOC_reduction.py | 12 ++--- src/lib/plots.py | 109 ++++++++++++++++++++++++------------------- src/lib/reduction.py | 15 +++--- 3 files changed, 73 insertions(+), 63 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 2e60d98..0b029d9 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -18,7 +18,7 @@ from matplotlib.colors import LogNorm def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop=False, interactive=False): # Reduction parameters # Deconvolution - deconvolve = True + deconvolve = False if deconvolve: # from lib.deconvolve import from_file_psf psf = 'gaussian' # Can be user-defined as well @@ -33,9 +33,9 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= display_crop = False # Background estimation - error_sub_type = 'sturges' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.25 - display_bkg = True + error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) + subtract_error = 0.50 + display_bkg = False # Data binning rebin = True @@ -45,12 +45,12 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Alignement align_center = 'center' # If None will not align the images - display_align = True + display_align = False display_data = False # Smoothing smoothing_function = 'combine' # gaussian_after, weighted_gaussian_after, gaussian, weighted_gaussian or combine - smoothing_FWHM = 0.10 # If None, no smoothing is done + smoothing_FWHM = 0.20 # If None, no smoothing is done smoothing_scale = 'arcsec' # pixel or arcsec # Rotation diff --git a/src/lib/plots.py b/src/lib/plots.py index fa12277..cca65ee 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -90,7 +90,7 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" plt.rcParams.update({'font.size': 10}) nb_obs = np.max([np.sum([head['filtnam1'] == curr_pol for head in headers]) for curr_pol in ['POL0', 'POL60', 'POL120']]) shape = np.array((3, nb_obs)) - fig, ax = plt.subplots(shape[0], shape[1], figsize=(10, 10), dpi=200, + fig, ax = plt.subplots(shape[0], shape[1], figsize=(3*shape[1], 3*shape[0]), dpi=200, layout='constrained', sharex=True, sharey=True) r_pol = dict(pol0=0, pol60=1, pol120=2) c_pol = dict(pol0=0, pol60=0, pol120=0) @@ -107,7 +107,11 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" else: ax_curr = ax[r_ax] # plots - vmin, vmax = convert*data[data > 0.].min()/10., convert*data[data > 0.].max() + if ('vmin' in kwargs.keys() or 'vmax' in kwargs.keys()): + vmin, vmax = kwargs['vmin'], kwargs['vmax'] + del kwargs['vmin'], kwargs['vmax'] + else: + vmin, vmax = convert*data[data > 0.].min()/10., convert*data[data > 0.].max() for key, value in [["cmap", [["cmap", "gray"]]], ["norm", [["norm", LogNorm(vmin, vmax)]]]]: try: _ = kwargs[key] @@ -127,7 +131,7 @@ def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="" ax_curr.annotate(filt, color='white', fontsize=10, xy=(0.01, 0.01), xycoords='axes fraction', verticalalignment='bottom', horizontalalignment='left') ax_curr.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=1.02) + # 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): @@ -160,6 +164,10 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): stkI = Stokes['I_stokes'].data.copy() stkQ = Stokes['Q_stokes'].data.copy() stkU = Stokes['U_stokes'].data.copy() + data_mask = Stokes['Data_mask'].data.astype(bool) + + for dataset in [stkI, stkQ, stkU]: + dataset[np.logical_not(data_mask)] = np.nan wcs = WCS(Stokes[0]).deepcopy() @@ -244,17 +252,23 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c The figure and ax created for interactive contour maps. """ # Get data - stkI = Stokes[np.argmax([Stokes[i].header['datatype'] == 'I_stokes' for i in range(len(Stokes))])] - stk_cov = Stokes[np.argmax([Stokes[i].header['datatype'] == 'IQU_cov_matrix' for i in range(len(Stokes))])] - pol = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Pol_deg_debiased' for i in range(len(Stokes))])] - pol_err = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Pol_deg_err' for i in range(len(Stokes))])] - pang = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Pol_ang' for i in range(len(Stokes))])] + stkI = Stokes['I_stokes'].data.copy() + stk_cov = Stokes['IQU_cov_matrix'].data.copy() + pol = Stokes['Pol_deg_debiased'].data.copy() + pol_err = Stokes['Pol_deg_err'].data.copy() + pang = Stokes['Pol_ang'].data.copy() try: if data_mask is None: - data_mask = Stokes[np.argmax([Stokes[i].header['datatype'] == 'Data_mask' for i in range(len(Stokes))])].data.astype(bool) + data_mask = Stokes['Data_mask'].data.astype(bool).copy() except KeyError: data_mask = np.ones(stkI.shape).astype(bool) + for dataset in [stkI, pol, pol_err, pang]: + dataset[np.logical_not(data_mask)] = np.nan + for i in range(3): + for j in range(3): + stk_cov[i][j][np.logical_not(data_mask)] = np.nan + pivot_wav = Stokes[0].header['photplam'] convert_flux = Stokes[0].header['photflam'] wcs = WCS(Stokes[0]).deepcopy() @@ -264,95 +278,94 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c plot_Stokes(Stokes, savename=savename, plots_folder=plots_folder) # Compute SNR and apply cuts - poldata, pangdata = pol.data.copy(), pang.data.copy() - maskP = pol_err.data > 0 - SNRp = np.zeros(pol.data.shape) - SNRp[maskP] = pol.data[maskP]/pol_err.data[maskP] + poldata, pangdata = pol.copy(), pang.copy() + maskP = pol_err > 0 + SNRp = np.ones(pol.shape)*np.nan + SNRp[maskP] = pol[maskP]/pol_err[maskP] - maskI = stk_cov.data[0, 0] > 0 - SNRi = np.zeros(stkI.data.shape) - SNRi[maskI] = stkI.data[maskI]/np.sqrt(stk_cov.data[0, 0][maskI]) + maskI = stk_cov[0, 0] > 0 + SNRi = np.ones(stkI.shape)*np.nan + SNRi[maskI] = stkI[maskI]/np.sqrt(stk_cov[0, 0][maskI]) mask = (SNRp > SNRp_cut) * (SNRi > SNRi_cut) poldata[np.logical_not(mask)] = np.nan pangdata[np.logical_not(mask)] = np.nan # Look for pixel of max 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) + if np.isfinite(pol).any(): + p_max = np.max(pol[np.isfinite(pol)]) + x_max, y_max = np.unravel_index(np.argmax(pol == p_max), pol.shape) else: print("No pixel with polarisation information above requested SNR.") # Plot the map plt.rcParams.update({'font.size': 10}) plt.rcdefaults() - fig, ax = plt.subplots(figsize=(10, 10), subplot_kw=dict(projection=wcs)) + fig, ax = plt.subplots(figsize=(10, 10), layout='constrained', subplot_kw=dict(projection=wcs)) ax.set(aspect='equal', fc='k') - fig.subplots_adjust(hspace=0, wspace=0, left=0.102, right=1.02) + # 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 display = 'i' if flux_lim is None: if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: vmin, vmax = flux_lim - im = ax.imshow(stkI.data*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(stkI*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) 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.array([0.8, 2., 5., 10., 20., 50.])/100.*vmax print("Total flux contour levels : ", levelsI) - ax.contour(stkI.data*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) + ax.contour(stkI*convert_flux, levels=levelsI, colors='grey', linewidths=0.5) elif display.lower() in ['pol_flux']: # Display polarisation flux display = 'pf' if flux_lim is None: if mask.sum() > 0.: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][mask])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][stkI.data > 0.])*convert_flux), np.max(stkI.data[stkI.data > 0.]*convert_flux) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) else: vmin, vmax = flux_lim - im = ax.imshow(stkI.data*convert_flux*pol.data, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(stkI*convert_flux*pol, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) 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) - ax.contour(stkI.data*convert_flux*pol.data, levels=levelsPf, colors='grey', linewidths=0.5) + ax.contour(stkI*convert_flux*pol, levels=levelsPf, colors='grey', linewidths=0.5) elif display.lower() in ['p', 'pol', 'pol_deg']: # 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.) + im = ax.imshow(pol*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) 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 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.) + im = ax.imshow(princ_angle(pang), vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) 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 polarisation degree error map display = 's_p' if (SNRp > SNRp_cut).any(): - vmin, vmax = 0., np.max(pol_err.data[SNRp > SNRp_cut])*100. - p_err = deepcopy(pol_err.data) - p_err[p_err > vmax/100.] = np.nan - im = ax.imshow(p_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno', alpha=1.) + vmin, vmax = 0., np.max([pol_err[SNRp > SNRp_cut].max(), 1.])*100. + im = ax.imshow(pol_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno_r', alpha=1.) else: - im = ax.imshow(pol_err.data*100., aspect='equal', cmap='inferno', alpha=1.) + vmin, vmax = 0., 100. + im = ax.imshow(pol_err*100., vmin=vmin, vmax=vmax, aspect='equal', cmap='inferno_r', alpha=1.) 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' if (SNRi > SNRi_cut).any(): - vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov.data[0, 0][stk_cov.data[0, 0] > 0.]) * - convert_flux), np.max(np.sqrt(stk_cov.data[0, 0][stk_cov.data[0, 0] > 0.])*convert_flux) - im = ax.imshow(np.sqrt(stk_cov.data[0, 0])*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + vmin, vmax = 1./2.*np.median(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.]) * + convert_flux), np.max(np.sqrt(stk_cov[0, 0][stk_cov[0, 0] > 0.])*convert_flux) + im = ax.imshow(np.sqrt(stk_cov[0, 0])*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno_r', alpha=1.) else: - im = ax.imshow(np.sqrt(stk_cov.data[0, 0])*convert_flux, aspect='equal', cmap='inferno', alpha=1.) + im = ax.imshow(np.sqrt(stk_cov[0, 0])*convert_flux, aspect='equal', cmap='inferno', alpha=1.) 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 @@ -381,15 +394,15 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c 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) + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov[0, 0][mask])*convert_flux), np.max(stkI[stkI > 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, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) + vmin, vmax = 1.*np.mean(np.sqrt(stk_cov[0, 0][stkI > 0.])*convert_flux), np.max(stkI[stkI > 0.]*convert_flux) + im = ax.imshow(stkI*convert_flux, norm=LogNorm(vmin, vmax), aspect='equal', cmap='inferno', alpha=1.) 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 - I_diluted = stkI.data[data_mask].sum() - I_diluted_err = np.sqrt(np.sum(stk_cov.data[0, 0][data_mask])) + I_diluted = stkI[data_mask].sum() + I_diluted_err = np.sqrt(np.sum(stk_cov[0, 0][data_mask])) P_diluted = Stokes[0].header['P_int'] P_diluted_err = Stokes[0].header['P_int_err'] @@ -407,7 +420,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c poldata[np.isfinite(poldata)] = 1./2. step_vec = 1 vec_scale = 2. - X, Y = np.meshgrid(np.arange(stkI.data.shape[1]), np.arange(stkI.data.shape[0])) + X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) U, V = poldata*np.cos(np.pi/2.+pangdata*np.pi/180.), poldata*np.sin(np.pi/2.+pangdata*np.pi/180.) ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', angles='uv', scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='w', edgecolor='k') @@ -429,7 +442,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c # Display instrument FOV if not (rectangle is None): x, y, width, height, angle, color = rectangle - x, y = np.array([x, y]) - np.array(stkI.data.shape)/2. + x, y = np.array([x, y]) - np.array(stkI.shape)/2. ax.add_patch(Rectangle((x, y), width, height, angle=angle, edgecolor=color, fill=False)) @@ -2053,7 +2066,7 @@ class pol_map(object): def submit_save(expression): ax_text_save.set(visible=False) if expression != '': - save_fig, save_ax = plt.subplots(figsize=(12, 10), layout='tight', subplot_kw=dict(projection=self.wcs)) + save_fig, save_ax = plt.subplots(figsize=(12, 10), layout='constrained', subplot_kw=dict(projection=self.wcs)) self.ax_cosmetics(ax=save_ax) self.display(fig=save_fig, ax=save_ax) self.pol_vector(fig=save_fig, ax=save_ax) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index ba7df4a..23551b5 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -276,7 +276,7 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu if display: plt.rcParams.update({'font.size': 15}) - fig, ax = plt.subplots(figsize=(10, 10)) + fig, ax = plt.subplots(figsize=(10, 10), layout='constrained') convert_flux = headers[0]['photflam'] data = deepcopy(data_array[0]*convert_flux) data[data <= data[data > 0.].min()] = data[data > 0.].min() @@ -301,16 +301,13 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu xycoords='axes fraction') ax.annotate(str(exptime)+" s", color='white', fontsize=10, xy=(0.80, 0.02), xycoords='axes fraction') - ax.set( # title="Location of cropped image.", - xlabel='pixel offset', - ylabel='pixel offset') + ax.set(title="Location of cropped image.", xlabel='pixel offset', ylabel='pixel offset') - fig.subplots_adjust(hspace=0, wspace=0, 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, wspace=0, right=0.85) + # cbar_ax = fig.add_axes([0.9, 0.12, 0.02, 0.75]) + fig.colorbar(im, ax=ax, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if savename is not None: - # fig.suptitle(savename+'_'+filt+'_crop_region') fig.savefig("/".join([plots_folder, savename+'_'+filt+'_crop_region.png']), bbox_inches='tight') plot_obs(data_array, headers, vmin=convert_flux*data_array[data_array > 0.].mean()/5., @@ -1406,7 +1403,7 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, for i, head in enumerate(headers): ang[i] = -head['orientat'] ang = ang.mean() - alpha = np.radians(ang) + alpha = np.pi/180.*ang mrot = np.array([[1., 0., 0.], [0., np.cos(2.*alpha), np.sin(2.*alpha)], [0, -np.sin(2.*alpha), np.cos(2.*alpha)]]) From a1767de88471528507d0e229474e7e262a3e2d95 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Thu, 2 May 2024 17:12:41 +0200 Subject: [PATCH 04/11] propagate options to multiple reductions, default SNRi_cut to 3 --- src/FOC_reduction.py | 10 +++++----- src/lib/plots.py | 20 ++++++++++---------- src/lib/reduction.py | 4 ++++ 3 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 0b029d9..4b3bc84 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -35,11 +35,11 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) subtract_error = 0.50 - display_bkg = False + display_bkg = True # Data binning rebin = True - pxsize = 0.10 + pxsize = 0.05 px_scale = 'arcsec' # pixel, arcsec or full rebin_operation = 'sum' # sum or average @@ -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.075 # If None, no smoothing is done smoothing_scale = 'arcsec' # pixel or arcsec # Rotation @@ -58,7 +58,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= rotate_stokes = True # Polarization map output - SNRp_cut = 3. # P measurments with SNR>3 + SNRp_cut = 5. # P measurments with SNR>3 SNRi_cut = 3. # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None vec_scale = 3 @@ -78,7 +78,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= target, products = retrieve_products(target, proposal_id, output_dir=output_dir) prod = products.pop() for prods in products: - main(target=target, infiles=["/".join(pr) for pr in prods], output_dir=output_dir) + main(target=target, infiles=["/".join(pr) for pr in prods], output_dir=output_dir, crop=crop, interactive=interactive) data_folder = prod[0][0] try: plots_folder = data_folder.replace("data", "plots") diff --git a/src/lib/plots.py b/src/lib/plots.py index cca65ee..e357348 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -200,7 +200,7 @@ def plot_Stokes(Stokes, savename=None, plots_folder=""): return 0 -def polarisation_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=3., flux_lim=None, step_vec=1, vec_scale=2., savename=None, plots_folder="", display="default"): """ Plots polarisation map from Stokes HDUList. @@ -673,7 +673,7 @@ class overplot_radio(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2, savename=None, **kwargs): + def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=3., vec_scale=2, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data @@ -780,7 +780,7 @@ class overplot_radio(align_maps): self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs) -> None: + def plot(self, levels=None, SNRp_cut=3., SNRi_cut=3., savename=None, **kwargs) -> None: while not self.aligned: self.align() self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, savename=savename, **kwargs) @@ -793,7 +793,7 @@ class overplot_chandra(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2, zoom=1, savename=None, **kwargs): + def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=3., vec_scale=2, zoom=1, savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data @@ -899,7 +899,7 @@ class overplot_chandra(align_maps): self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3., SNRi_cut=30., zoom=1, savename=None, **kwargs) -> None: + def plot(self, levels=None, SNRp_cut=3., SNRi_cut=3., zoom=1, savename=None, **kwargs) -> None: while not self.aligned: self.align() self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, zoom=zoom, savename=savename, **kwargs) @@ -912,7 +912,7 @@ class overplot_pol(align_maps): Inherit from class align_maps in order to get the same WCS on both maps. """ - def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2., savename=None, **kwargs): + def overplot(self, levels=None, SNRp_cut=3., SNRi_cut=3., vec_scale=2., savename=None, **kwargs): self.Stokes_UV = self.map self.wcs_UV = self.map_wcs # Get Data @@ -1019,7 +1019,7 @@ class overplot_pol(align_maps): self.fig_overplot.canvas.draw() - def plot(self, levels=None, SNRp_cut=3., SNRi_cut=30., vec_scale=2., savename=None, **kwargs) -> None: + def plot(self, levels=None, SNRp_cut=3., SNRi_cut=3., vec_scale=2., savename=None, **kwargs) -> None: while not self.aligned: self.align() self.overplot(levels=levels, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, vec_scale=vec_scale, savename=savename, **kwargs) @@ -1059,7 +1059,7 @@ class align_pol(object): self.kwargs = kwargs - def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs): + def single_plot(self, curr_map, wcs, v_lim=None, ax_lim=None, SNRp_cut=3., SNRi_cut=3., savename=None, **kwargs): # Get data stkI = curr_map['I_STOKES'].data stk_cov = curr_map['IQU_COV_MATRIX'].data @@ -1150,7 +1150,7 @@ class align_pol(object): self.wcs, self.wcs_other[i] = curr_align.align() self.aligned[i] = curr_align.aligned - def plot(self, SNRp_cut=3., SNRi_cut=30., savename=None, **kwargs): + def plot(self, SNRp_cut=3., SNRi_cut=3., savename=None, **kwargs): while not self.aligned.all(): self.align() eps = 1e-35 @@ -1737,7 +1737,7 @@ class pol_map(object): Class to interactively study polarisation maps. """ - def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=30., flux_lim=None, selection=None): + def __init__(self, Stokes, SNRp_cut=3., SNRi_cut=3., flux_lim=None, selection=None): if isinstance(Stokes, str): Stokes = fits.open(Stokes) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 23551b5..8d46219 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -952,6 +952,10 @@ def polarizer_avg(data_array, error_array, data_mask, headers, FWHM=None, scale= err60_array = error_array[is_pol60] err120_array = error_array[is_pol120] + # For a single observation, combination amount to a weighted gaussian + if np.max([is_pol0.sum(), is_pol60.sum(), is_pol120.sum()]) == 1 and smoothing.lower() in ['combine', 'combining']: + smoothing = 'weighted_gaussian' + if (FWHM is not None) and (smoothing.lower() in ['combine', 'combining']): # Smooth by combining each polarizer images pol0, err0 = smooth_data(pol0_array, err0_array, data_mask, headers0, FWHM=FWHM, scale=scale, smoothing=smoothing) From ea0ac013b730db82eff774a56e763aba7c4053ce Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 13 May 2024 16:20:53 +0200 Subject: [PATCH 05/11] revert background estimation to gaussian fit --- src/lib/background.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/lib/background.py b/src/lib/background.py index 8d07c25..74f41e6 100755 --- a/src/lib/background.py +++ b/src/lib/background.py @@ -80,7 +80,8 @@ def display_bkg(data, background, std_bkg, headers, histograms=None, binning=Non label=headers[i]['filtnam1']+' (Obs '+str(filt_obs[headers[i]['filtnam1']])+')') ax_h.plot([background[i]*convert_flux[i], background[i]*convert_flux[i]], [hist.min(), hist.max()], 'x--', color="C{0:d}".format(i), alpha=0.8) if not (coeff is None): - ax_h.plot(bins*convert_flux[i], gausspol(bins, *coeff[i]), '--', color="C{0:d}".format(i), alpha=0.8) + # ax_h.plot(bins*convert_flux[i], gausspol(bins, *coeff[i]), '--', color="C{0:d}".format(i), alpha=0.8) + ax_h.plot(bins*convert_flux[i], gauss(bins, *coeff[i]), '--', color="C{0:d}".format(i), alpha=0.8) ax_h.set_xscale('log') ax_h.set_ylim([0., np.max([hist.max() for hist in histograms])]) ax_h.set_xlim([np.min(background*convert_flux)*1e-2, np.max(background*convert_flux)*1e2]) @@ -159,12 +160,15 @@ def bkg_estimate(img, bins=None, chi2=None, coeff=None): peak = binning[np.argmax(hist)] bins_stdev = binning[hist > hist.max()/2.] stdev = bins_stdev[-1]-bins_stdev[0] - p0 = [hist.max(), peak, stdev, 1e-3, 1e-3, 1e-3, 1e-3] + # p0 = [hist.max(), peak, stdev, 1e-3, 1e-3, 1e-3, 1e-3] + p0 = [hist.max(), peak, stdev] try: - popt, pcov = curve_fit(gausspol, binning, hist, p0=p0) + # popt, pcov = curve_fit(gausspol, binning, hist, p0=p0) + popt, pcov = curve_fit(gauss, binning, hist, p0=p0) except RuntimeError: popt = p0 - chi2.append(np.sum((hist - gausspol(binning, *popt))**2)/hist.size) + # chi2.append(np.sum((hist - gausspol(binning, *popt))**2)/hist.size) + chi2.append(np.sum((hist - gauss(binning, *popt))**2)/hist.size) coeff.append(popt) return bins, chi2, coeff @@ -330,8 +334,10 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis # Fit a gaussian to the log-intensity histogram bins_stdev = binning[-1][hist > hist.max()/2.] stdev = bins_stdev[-1]-bins_stdev[0] - p0 = [hist.max(), binning[-1][np.argmax(hist)], stdev, 1e-3, 1e-3, 1e-3, 1e-3] - popt, pcov = curve_fit(gausspol, binning[-1], hist, p0=p0) + # p0 = [hist.max(), binning[-1][np.argmax(hist)], stdev, 1e-3, 1e-3, 1e-3, 1e-3] + p0 = [hist.max(), binning[-1][np.argmax(hist)], stdev] + # popt, pcov = curve_fit(gausspol, binning[-1], hist, p0=p0) + popt, pcov = curve_fit(gauss, binning[-1], hist, p0=p0) coeff.append(popt) bkg = popt[1]+np.abs(popt[2])*subtract_error From a515182b2cc78d3539b6d9c0565e73e96354dc43 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Mon, 13 May 2024 16:22:12 +0200 Subject: [PATCH 06/11] fix crop and clean up --- src/lib/convex_hull.py | 54 ++++++++++++++++++++++-------------------- src/lib/fits.py | 1 - src/lib/reduction.py | 4 ++-- 3 files changed, 30 insertions(+), 29 deletions(-) diff --git a/src/lib/convex_hull.py b/src/lib/convex_hull.py index 7cdfcd8..5e576fe 100755 --- a/src/lib/convex_hull.py +++ b/src/lib/convex_hull.py @@ -270,7 +270,7 @@ def del_aligned(s, Omega): # Implement Graham algorithm -def convex_hull(A): +def convex_hull(H): """ Implement Graham algorithm to find the convex hull of a given list of points, handling aligned points. @@ -283,17 +283,20 @@ def convex_hull(A): List of points defining the convex hull of the input list A. """ S = empty_stack() - Omega = min_lexico(A) - sort_angles_distances(Omega, A) - A = del_aligned(A, Omega) - stack(S, A[0]) - stack(S, A[1]) - stack(S, A[2]) - for i in range(3, len(A)): - while pseudo_angle(stack_sub_top(S), stack_top(S), A[i]) <= 0: - unstack(S) - stack(S, A[i]) - return S + Omega = min_lexico(H) + sort_angles_distances(Omega, H) + A = del_aligned(H, Omega) + if len(A) < 3: + return H + else: + stack(S, A[0]) + stack(S, A[1]) + stack(S, A[2]) + for i in range(3, len(A)): + while pseudo_angle(stack_sub_top(S), stack_top(S), A[i]) <= 0: + unstack(S) + stack(S, A[i]) + return S def image_hull(image, step=5, null_val=0., inside=True): @@ -328,25 +331,24 @@ def image_hull(image, step=5, null_val=0., inside=True): H = [] shape = np.array(image.shape) row, col = np.indices(shape) - for i in range(0, shape[0], step): - r = row[i, :][image[i, :] > null_val] - c = col[i, :][image[i, :] > null_val] - if len(r) > 1 and len(c) > 1: - H.append((r[0], c[0])) - H.append((r[-1], c[-1])) - for j in range(0, shape[1], step): - r = row[:, j][image[:, j] > null_val] - c = col[:, j][image[:, j] > null_val] - if len(r) > 1 and len(c) > 1: - if not ((r[0], c[0]) in H): - H.append((r[0], c[0])) - if not ((r[-1], c[-1]) in H): - H.append((r[-1], c[-1])) + for i in range(0, int(min(shape)/2), step): + r1, r2 = row[i, :][image[i, :] > null_val], row[-i, :][image[-i, :] > null_val] + c1, c2 = col[i, :][image[i, :] > null_val], col[-i, :][image[-i, :] > null_val] + if r1.shape[0] > 1: + H.append((r1[0], c1[0])) + H.append((r1[-1], c1[-1])) + if r2.shape[0] > 1: + H.append((r2[0], c2[0])) + H.append((r2[-1], c2[-1])) S = np.array(convex_hull(H)) x_min, y_min = S[:, 0] < S[:, 0].mean(), S[:, 1] < S[:, 1].mean() x_max, y_max = S[:, 0] > S[:, 0].mean(), S[:, 1] > S[:, 1].mean() # Get the 4 extrema + # S0 = S[x_min*y_min][np.argmin(S[x_min*y_min][:, 0])] + # S1 = S[x_min*y_max][np.argmax(S[x_min*y_max][:, 1])] + # S2 = S[x_max*y_min][np.argmin(S[x_max*y_min][:, 1])] + # S3 = S[x_max*y_max][np.argmax(S[x_max*y_max][:, 0])] S0 = S[x_min*y_min][np.abs(0-S[x_min*y_min].sum(axis=1)).min() == np.abs(0-S[x_min*y_min].sum(axis=1))][0] S1 = S[x_min*y_max][np.abs(shape[1]-S[x_min*y_max].sum(axis=1)).min() == np.abs(shape[1]-S[x_min*y_max].sum(axis=1))][0] S2 = S[x_max*y_min][np.abs(shape[0]-S[x_max*y_min].sum(axis=1)).min() == np.abs(shape[0]-S[x_max*y_min].sum(axis=1))][0] diff --git a/src/lib/fits.py b/src/lib/fits.py index ec90801..8f9b389 100755 --- a/src/lib/fits.py +++ b/src/lib/fits.py @@ -14,7 +14,6 @@ from os.path import join as path_join from astropy.io import fits from astropy.wcs import WCS from lib.convex_hull import clean_ROI -from lib.utils import princ_angle def get_obs_data(infiles, data_folder="", compute_flux=False): diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 8d46219..2c194b4 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -442,11 +442,11 @@ def get_error(data_array, headers, error_array=None, data_mask=None, sub_type=No Returns: data_array : numpy.ndarray Array containing the data to study minus the background. - headers : header list - Updated headers associated with the images in data_array. error_array : numpy.ndarray Array containing the background values associated to the images in data_array. + headers : header list + Updated headers associated with the images in data_array. background : numpy.ndarray Array containing the pixel background value for each image in data_array. From d3915b370672702fb55608a9d4529093cca6e9b6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 21 May 2024 17:11:09 +0200 Subject: [PATCH 07/11] replace src/analysis with execution from src/lib/plots and move unused scripts to utils --- src/lib/plots.py | 37 +++++++++++++++----- {src => utils}/analysis.py | 4 +-- {src => utils}/comparison_Kishimoto.py | 10 +++--- utils/get_cdelt.py | 47 ++++++++++++++++++++++++++ 4 files changed, 82 insertions(+), 16 deletions(-) rename {src => utils}/analysis.py (95%) rename {src => utils}/comparison_Kishimoto.py (98%) create mode 100755 utils/get_cdelt.py diff --git a/src/lib/plots.py b/src/lib/plots.py index e357348..3386c41 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -1,3 +1,4 @@ +#!/usr/bin/python """ Library functions for displaying informations using matplotlib @@ -54,7 +55,7 @@ from astropy.wcs import WCS from astropy.io import fits from astropy.coordinates import SkyCoord from scipy.ndimage import zoom as sc_zoom -from lib.utils import rot2D, princ_angle, sci_not +from utils import rot2D, princ_angle, sci_not def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): @@ -423,7 +424,7 @@ def polarisation_map(Stokes, data_mask=None, rectangle=None, SNRp_cut=3., SNRi_c X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) U, V = poldata*np.cos(np.pi/2.+pangdata*np.pi/180.), poldata*np.sin(np.pi/2.+pangdata*np.pi/180.) ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', angles='uv', - scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='w', edgecolor='k') + scale=1./vec_scale, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='w', edgecolor='k') pol_sc = AnchoredSizeBar(ax.transData, vec_scale, r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') ax.add_artist(pol_sc) @@ -734,7 +735,7 @@ class overplot_radio(align_maps): 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="{0:s} polarisation map".format(self.map_observer)) + scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) self.ax_overplot.autoscale(False) # Display other map as contours @@ -853,7 +854,7 @@ class overplot_chandra(align_maps): 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="{0:s} polarisation map".format(self.map_observer)) + scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black', label="{0:s} polarisation map".format(self.map_observer)) self.ax_overplot.autoscale(False) # Display other map as contours @@ -971,8 +972,8 @@ class overplot_pol(align_maps): 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="{0:s} polarisation map".format(self.map_observer)) + 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.5, linewidth=0.8, 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: @@ -1128,7 +1129,7 @@ class align_pol(object): X, Y = np.meshgrid(np.arange(stkI.shape[1]), np.arange(stkI.shape[0])) U, V = pol*np.cos(np.pi/2.+pang*np.pi/180.), pol*np.sin(np.pi/2.+pang*np.pi/180.) ax.quiver(X[::step_vec, ::step_vec], Y[::step_vec, ::step_vec], U[::step_vec, ::step_vec], V[::step_vec, ::step_vec], units='xy', - angles='uv', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.1, color='w') + angles='uv', scale=0.5, scale_units='xy', pivot='mid', headwidth=0., headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='w') pol_sc = AnchoredSizeBar(ax.transData, 2., r"$P$= 100 %", 4, pad=0.5, sep=5, borderpad=0.5, frameon=False, size_vertical=0.005, color='w') ax.add_artist(pol_sc) @@ -2349,12 +2350,12 @@ class pol_map(object): if hasattr(self, 'quiver'): self.quiver.remove() self.quiver = ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='white', edgecolor='black') + headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black') fig.canvas.draw_idle() return self.quiver else: ax.quiver(X, Y, XY_U, XY_V, units='xy', scale=1./self.vec_scale, scale_units='xy', pivot='mid', headwidth=0., - headlength=0., headaxislength=0., width=0.2, linewidth=0.3, color='white', edgecolor='black') + headlength=0., headaxislength=0., width=0.5, linewidth=0.8, color='white', edgecolor='black') fig.canvas.draw_idle() def pol_int(self, fig=None, ax=None): @@ -2463,3 +2464,21 @@ class pol_map(object): if self.region is not None: ax.contour(self.region.astype(float), levels=[0.5], colors='white', linewidths=0.8) fig.canvas.draw_idle() + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description='Interactively plot the pipeline products') + parser.add_argument('-f', '--file', metavar='path', required=False, help='the full or relative path to the data product', type=str, default=None) + parser.add_argument('-p', '--snrp', metavar='snrp_cut', required=False, help='the cut in signal-to-noise for the polarisation degree', type=float, default=3.) + parser.add_argument('-i', '--snri', metavar='snri_cut', required=False, help='the cut in signal-to-noise for the intensity', type=float, default=3.) + parser.add_argument('-l', '--lim', metavar='flux_lim', nargs=2, required=False, help='limits for the intensity map', default=None) + args = parser.parse_args() + + if args.file is not None: + Stokes_UV = fits.open(args.file, mode='readonly') + p = pol_map(Stokes_UV, SNRp_cut=args.snrp, SNRi_cut=args.snri, flux_lim=args.lim) + + else: + print("python3 plots.py -f -p -i -l ") diff --git a/src/analysis.py b/utils/analysis.py similarity index 95% rename from src/analysis.py rename to utils/analysis.py index e2d2fea..6daecb8 100755 --- a/src/analysis.py +++ b/utils/analysis.py @@ -1,4 +1,4 @@ -#!/usr/bin/python3 +#!/usr/bin/python from getopt import getopt, error as get_error from sys import argv @@ -30,7 +30,7 @@ except get_error as err: if fits_path is not None: from astropy.io import fits - from lib.plots import pol_map + from src.lib.plots import pol_map Stokes_UV = fits.open(fits_path) p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) diff --git a/src/comparison_Kishimoto.py b/utils/comparison_Kishimoto.py similarity index 98% rename from src/comparison_Kishimoto.py rename to utils/comparison_Kishimoto.py index 56b67ae..75b7073 100755 --- a/src/comparison_Kishimoto.py +++ b/utils/comparison_Kishimoto.py @@ -1,8 +1,8 @@ -# !/usr/bin/env python -from lib.background import gauss, bin_centers -from lib.deconvolve import zeropad -from lib.reduction import align_data -from lib.plots import princ_angle +#!/usr/bin/python +from src.lib.background import gauss, bin_centers +from src.lib.deconvolve import zeropad +from src.lib.reduction import align_data +from src.lib.plots import princ_angle from matplotlib.colors import LogNorm from os.path import join as path_join from astropy.io import fits diff --git a/utils/get_cdelt.py b/utils/get_cdelt.py new file mode 100755 index 0000000..45e526b --- /dev/null +++ b/utils/get_cdelt.py @@ -0,0 +1,47 @@ +#!/usr/bin/python + +def main(infiles=None): + """ + Retrieve native spatial resolution from given observation. + """ + from os.path import join as path_join + from warnings import catch_warnings, filterwarnings + from astropy.io.fits import getheader + from astropy.wcs import WCS, FITSFixedWarning + from numpy.linalg import eig + + if infiles is None: + print("Usage: \"python get_cdelt.py -f infiles\"") + return 1 + prod = [["/".join(filepath.split('/')[:-1]), filepath.split('/')[-1]] for filepath in infiles] + data_folder = prod[0][0] + infiles = [p[1] for p in prod] + + cdelt = {} + size = {} + for currfile in infiles: + with catch_warnings(): + filterwarnings('ignore', message="'datfix' made the change", category=FITSFixedWarning) + wcs = WCS(getheader(path_join(data_folder, currfile))).celestial + key = currfile[:-5] + size[key] = wcs.array_shape + if wcs.wcs.has_cd(): + cdelt[key] = eig(wcs.wcs.cd)[0]*3600. + else: + cdelt[key] = wcs.wcs.cdelt*3600. + + print("Image name, native resolution in arcsec and shape") + for currfile in infiles: + key = currfile[:-5] + print(key, cdelt[key], size[key]) + + return 0 + + +if __name__ == "__main__": + import argparse + + parser = argparse.ArgumentParser(description='Query MAST for target products') + parser.add_argument('-f', '--files', metavar='path', required=False, nargs='*', help='the full or relative path to the data products', default=None) + args = parser.parse_args() + exitcode = main(infiles=args.files) From 1f43a3194d319bad7346bfe172a77361423ebe29 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 21 May 2024 17:12:59 +0200 Subject: [PATCH 08/11] default snrp to 3 and image output to pdf everywhere --- src/lib/reduction.py | 16 ++++++++++------ src/overplot_MRK463E.py | 7 ++++--- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 2c194b4..0e74d5a 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -308,8 +308,8 @@ def crop_array(data_array, headers, error_array=None, data_mask=None, step=5, nu fig.colorbar(im, ax=ax, label=r"Flux [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") if savename is not None: - fig.savefig("/".join([plots_folder, savename+'_'+filt+'_crop_region.png']), - bbox_inches='tight') + fig.savefig("/".join([plots_folder, savename+'_'+filt+'_crop_region.pdf']), + bbox_inches='tight', dpi=200) plot_obs(data_array, headers, vmin=convert_flux*data_array[data_array > 0.].mean()/5., vmax=convert_flux*data_array[data_array > 0.].max(), rectangle=[rectangle,]*len(headers), savename=savename+'_crop_region', plots_folder=plots_folder) @@ -1228,10 +1228,12 @@ def compute_Stokes(data_array, error_array, data_mask, headers, FWHM=None, scale QU_diluted_err = np.sqrt(np.sum(Stokes_cov[1, 2][mask]**2)) P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted - P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted ** + 2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted, Q_diluted)) - 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) + 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 polarisation degree') @@ -1487,10 +1489,12 @@ def rotate_Stokes(I_stokes, Q_stokes, U_stokes, Stokes_cov, data_mask, headers, QU_diluted_err = np.sqrt(np.sum(new_Stokes_cov[1, 2][mask]**2)) P_diluted = np.sqrt(Q_diluted**2+U_diluted**2)/I_diluted - P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted**2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) + P_diluted_err = (1./I_diluted)*np.sqrt((Q_diluted**2*Q_diluted_err**2 + U_diluted**2*U_diluted_err**2 + 2.*Q_diluted*U_diluted*QU_diluted_err)/(Q_diluted**2 + U_diluted ** + 2) + ((Q_diluted/I_diluted)**2 + (U_diluted/I_diluted)**2)*I_diluted_err**2 - 2.*(Q_diluted/I_diluted)*IQ_diluted_err - 2.*(U_diluted/I_diluted)*IU_diluted_err) PA_diluted = princ_angle((90./np.pi)*np.arctan2(U_diluted, Q_diluted)) - 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) + 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 polarisation degree') diff --git a/src/overplot_MRK463E.py b/src/overplot_MRK463E.py index e2070df..79de0ca 100755 --- a/src/overplot_MRK463E.py +++ b/src/overplot_MRK463E.py @@ -6,7 +6,7 @@ from matplotlib.colors import LogNorm Stokes_UV = fits.open("./data/MRK463E/5960/MRK463E_FOC_b0.05arcsec_c0.10arcsec.fits") Stokes_IR = fits.open("./data/MRK463E/WFPC2/IR_rot_crop.fits") -Stokes_Xr = fits.open("./data/MRK463E/Chandra/4913/primary/acisf04913N004_cntr_img2.fits") +Stokes_Xr = fits.open("./data/MRK463E/Chandra/X_ray_crop.fits") levels = np.geomspace(1., 99., 7) @@ -14,12 +14,13 @@ levels = np.geomspace(1., 99., 7) # A.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=20.0, zoom=1, savename='./plots/MRK463E/Chandra_overplot.pdf') B = overplot_chandra(Stokes_UV, Stokes_Xr, norm=LogNorm()) -B.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=30.0, vec_scale=3, zoom=1, savename='./plots/MRK463E/Chandra_overplot_forced.pdf') +B.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, zoom=1, savename='./plots/MRK463E/Chandra_overplot_forced.pdf') B.write_to(path1="./data/MRK463E/FOC_data_Chandra.fits", path2="./data/MRK463E/Chandra_data.fits", suffix="aligned") # C = overplot_pol(Stokes_UV, Stokes_IR) # C.plot(SNRp_cut=3.0, SNRi_cut=20.0, savename='./plots/MRK463E/IR_overplot.pdf') +levels = np.array([0.8, 2, 5, 10, 20, 50])/100.*Stokes_UV[0].header['photflam'] D = overplot_pol(Stokes_UV, Stokes_IR, norm=LogNorm()) -D.plot(SNRp_cut=3.0, SNRi_cut=30.0, vec_scale=3, norm=LogNorm(1e-18, 1e-15), savename='./plots/MRK463E/IR_overplot_forced.pdf') +D.plot(levels=levels, SNRp_cut=3.0, SNRi_cut=3.0, vec_scale=5, norm=LogNorm(8.5e-18, 2.5e-15), savename='./plots/MRK463E/IR_overplot_forced.pdf') D.write_to(path1="./data/MRK463E/FOC_data_WFPC.fits", path2="./data/MRK463E/WFPC_data.fits", suffix="aligned") From 371ec8c1ea82bef5251b1d0f1bd85ce2123f7b35 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 21 May 2024 17:14:19 +0200 Subject: [PATCH 09/11] better handling of image border while background computation, clean imports --- src/FOC_reduction.py | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 4b3bc84..6b586e0 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -1,4 +1,4 @@ -#!/usr/bin/python3 +#!/usr/bin/python # -*- coding:utf-8 -*- """ Main script where are progressively added the steps for the FOC pipeline reduction. @@ -7,10 +7,11 @@ Main script where are progressively added the steps for the FOC pipeline reducti # Project libraries import numpy as np from copy import deepcopy +from os import system +from os.path import exists as path_exists import lib.fits as proj_fits # Functions to handle fits files import lib.reduction as proj_red # Functions used in reduction pipeline import lib.plots as proj_plots # Functions for plotting data -from lib.query import retrieve_products, path_exists, system from lib.utils import sci_not, princ_angle from matplotlib.colors import LogNorm @@ -26,20 +27,20 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= psf_FWHM = 3.1 psf_scale = 'px' psf_shape = None # (151, 151) - iterations = 3 - algo = "richardson" + iterations = 1 + algo = "conjgrad" # Initial crop display_crop = False # Background estimation error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.50 + subtract_error = 0.01 display_bkg = True # Data binning rebin = True - pxsize = 0.05 + pxsize = 0.10 px_scale = 'arcsec' # pixel, arcsec or full rebin_operation = 'sum' # sum or average @@ -50,7 +51,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.075 # If None, no smoothing is done + smoothing_FWHM = 0.200 # If None, no smoothing is done smoothing_scale = 'arcsec' # pixel or arcsec # Rotation @@ -58,10 +59,10 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= rotate_stokes = True # Polarization map output - SNRp_cut = 5. # P measurments with SNR>3 + SNRp_cut = 3. # P measurments with SNR>3 SNRi_cut = 3. # I measurments with SNR>30, which implies an uncertainty in P of 4.7%. flux_lim = None # lowest and highest flux displayed on plot, defaults to bkg and maximum in cut if None - vec_scale = 3 + vec_scale = 5 step_vec = 1 # plot all vectors in the array. if step_vec = 2, then every other vector will be plotted if step_vec = 0 then all vectors are displayed at full length # Pipeline start @@ -75,6 +76,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= if target is None: target = input("Target name:\n>") else: + from lib.query import retrieve_products target, products = retrieve_products(target, proposal_id, output_dir=output_dir) prod = products.pop() for prods in products: @@ -107,6 +109,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Crop data to remove outside blank margins. data_array, error_array, headers = proj_red.crop_array(data_array, headers, step=5, null_val=0., inside=True, display=display_crop, savename=figname, plots_folder=plots_folder) + data_mask = np.ones(data_array[0].shape, dtype=bool) # Deconvolve data using Richardson-Lucy iterative algorithm with a gaussian PSF of given FWHM. if deconvolve: @@ -114,8 +117,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Estimate error from data background, estimated from sub-image of desired sub_shape. background = None - data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, sub_type=error_sub_type, subtract_error=subtract_error, display=display_bkg, savename="_".join([ - figname, "errors"]), plots_folder=plots_folder, return_background=True) + data_array, error_array, headers, background = proj_red.get_error(data_array, headers, error_array, data_mask=data_mask, sub_type=error_sub_type, subtract_error=subtract_error, display=display_bkg, savename="_".join([figname, "errors"]), plots_folder=plots_folder, return_background=True) # Align and rescale images with oversampling. data_array, error_array, headers, data_mask = proj_red.align_data( From 4e17f4053486f002e75806348cab6154aec68ad0 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Wed, 22 May 2024 11:00:26 +0200 Subject: [PATCH 10/11] temp fix for ModuleNotFound error in lib.plots --- src/lib/plots.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/lib/plots.py b/src/lib/plots.py index 3386c41..dcde2be 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -55,7 +55,10 @@ from astropy.wcs import WCS from astropy.io import fits from astropy.coordinates import SkyCoord from scipy.ndimage import zoom as sc_zoom -from utils import rot2D, princ_angle, sci_not +try: + from utils import rot2D, princ_angle, sci_not +except ModuleNotFoundError: + from lib.utils import rot2D, princ_angle, sci_not def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): From fb4849df25732edaa975ea70c523596666640df4 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Fri, 24 May 2024 12:15:06 +0200 Subject: [PATCH 11/11] package-style architecture --- {src => package}/FOC_reduction.py | 4 ++-- package/__init__.py | 2 ++ package/lib/__init__.py | 9 +++++++++ {src => package}/lib/background.py | 2 +- {src => package}/lib/convex_hull.py | 0 {src => package}/lib/cross_correlation.py | 0 {src => package}/lib/deconvolve.py | 0 {src => package}/lib/fits.py | 2 +- {src => package}/lib/plots.py | 8 ++++---- {src => package}/lib/query.py | 0 {src => package}/lib/reduction.py | 12 ++++++------ {src => package}/lib/utils.py | 0 {src => package}/overplot_IC5063.py | 0 {src => package}/overplot_MRK463E.py | 0 package/src/__init__.py | 0 {utils => package/src}/analysis.py | 2 +- {utils => package/src}/comparison_Kishimoto.py | 0 {utils => package/src}/get_cdelt.py | 0 18 files changed, 26 insertions(+), 15 deletions(-) rename {src => package}/FOC_reduction.py (99%) create mode 100644 package/__init__.py create mode 100644 package/lib/__init__.py rename {src => package}/lib/background.py (99%) rename {src => package}/lib/convex_hull.py (100%) rename {src => package}/lib/cross_correlation.py (100%) rename {src => package}/lib/deconvolve.py (100%) rename {src => package}/lib/fits.py (99%) rename {src => package}/lib/plots.py (99%) rename {src => package}/lib/query.py (100%) rename {src => package}/lib/reduction.py (99%) rename {src => package}/lib/utils.py (100%) rename {src => package}/overplot_IC5063.py (100%) rename {src => package}/overplot_MRK463E.py (100%) create mode 100644 package/src/__init__.py rename {utils => package/src}/analysis.py (96%) rename {utils => package/src}/comparison_Kishimoto.py (100%) rename {utils => package/src}/get_cdelt.py (100%) diff --git a/src/FOC_reduction.py b/package/FOC_reduction.py similarity index 99% rename from src/FOC_reduction.py rename to package/FOC_reduction.py index 6b586e0..565ff7a 100755 --- a/src/FOC_reduction.py +++ b/package/FOC_reduction.py @@ -35,7 +35,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= # Background estimation error_sub_type = 'freedman-diaconis' # sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51, 51)) - subtract_error = 0.01 + subtract_error = 0.50 display_bkg = True # Data binning @@ -51,7 +51,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.200 # If None, no smoothing is done + smoothing_FWHM = 0.150 # If None, no smoothing is done smoothing_scale = 'arcsec' # pixel or arcsec # Rotation diff --git a/package/__init__.py b/package/__init__.py new file mode 100644 index 0000000..80d0eb7 --- /dev/null +++ b/package/__init__.py @@ -0,0 +1,2 @@ +from . import lib +from . import src diff --git a/package/lib/__init__.py b/package/lib/__init__.py new file mode 100644 index 0000000..b866f6e --- /dev/null +++ b/package/lib/__init__.py @@ -0,0 +1,9 @@ +from . import background +from . import convex_hull +from . import cross_correlation +from . import deconvolve +from . import fits +from . import plots +from . import query +from . import reduction +from . import utils diff --git a/src/lib/background.py b/package/lib/background.py similarity index 99% rename from src/lib/background.py rename to package/lib/background.py index 74f41e6..fb0ae83 100755 --- a/src/lib/background.py +++ b/package/lib/background.py @@ -17,7 +17,7 @@ import matplotlib.dates as mdates from matplotlib.colors import LogNorm from matplotlib.patches import Rectangle from datetime import datetime -from lib.plots import plot_obs +from .plots import plot_obs from scipy.optimize import curve_fit diff --git a/src/lib/convex_hull.py b/package/lib/convex_hull.py similarity index 100% rename from src/lib/convex_hull.py rename to package/lib/convex_hull.py diff --git a/src/lib/cross_correlation.py b/package/lib/cross_correlation.py similarity index 100% rename from src/lib/cross_correlation.py rename to package/lib/cross_correlation.py diff --git a/src/lib/deconvolve.py b/package/lib/deconvolve.py similarity index 100% rename from src/lib/deconvolve.py rename to package/lib/deconvolve.py diff --git a/src/lib/fits.py b/package/lib/fits.py similarity index 99% rename from src/lib/fits.py rename to package/lib/fits.py index 8f9b389..e966f43 100755 --- a/src/lib/fits.py +++ b/package/lib/fits.py @@ -13,7 +13,7 @@ import numpy as np from os.path import join as path_join from astropy.io import fits from astropy.wcs import WCS -from lib.convex_hull import clean_ROI +from .convex_hull import clean_ROI def get_obs_data(infiles, data_folder="", compute_flux=False): diff --git a/src/lib/plots.py b/package/lib/plots.py similarity index 99% rename from src/lib/plots.py rename to package/lib/plots.py index dcde2be..0ea1354 100755 --- a/src/lib/plots.py +++ b/package/lib/plots.py @@ -56,9 +56,9 @@ from astropy.io import fits from astropy.coordinates import SkyCoord from scipy.ndimage import zoom as sc_zoom try: + from .utils import rot2D, princ_angle, sci_not +except ImportError: from utils import rot2D, princ_angle, sci_not -except ModuleNotFoundError: - from lib.utils import rot2D, princ_angle, sci_not def plot_obs(data_array, headers, rectangle=None, savename=None, plots_folder="", **kwargs): @@ -975,8 +975,8 @@ class overplot_pol(align_maps): 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.5, linewidth=0.8, color='white', edgecolor='black', transform=self.ax_overplot.get_transform(self.wcs_UV), label="{0:s} polarisation map".format(self.map_observer)) + 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=2.0, linewidth=1.0, 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: diff --git a/src/lib/query.py b/package/lib/query.py similarity index 100% rename from src/lib/query.py rename to package/lib/query.py diff --git a/src/lib/reduction.py b/package/lib/reduction.py similarity index 99% rename from src/lib/reduction.py rename to package/lib/reduction.py index 0e74d5a..2553f13 100755 --- a/src/lib/reduction.py +++ b/package/lib/reduction.py @@ -49,12 +49,12 @@ from scipy.signal import fftconvolve from astropy.wcs import WCS from astropy import log import warnings -from lib.deconvolve import deconvolve_im, gaussian_psf, gaussian2d, zeropad -from lib.convex_hull import image_hull, clean_ROI -from lib.background import bkg_fit, bkg_hist, bkg_mini -from lib.plots import plot_obs -from lib.utils import princ_angle -from lib.cross_correlation import phase_cross_correlation +from .deconvolve import deconvolve_im, gaussian_psf, gaussian2d, zeropad +from .convex_hull import image_hull, clean_ROI +from .background import bkg_fit, bkg_hist, bkg_mini +from .plots import plot_obs +from .utils import princ_angle +from .cross_correlation import phase_cross_correlation log.setLevel('ERROR') diff --git a/src/lib/utils.py b/package/lib/utils.py similarity index 100% rename from src/lib/utils.py rename to package/lib/utils.py diff --git a/src/overplot_IC5063.py b/package/overplot_IC5063.py similarity index 100% rename from src/overplot_IC5063.py rename to package/overplot_IC5063.py diff --git a/src/overplot_MRK463E.py b/package/overplot_MRK463E.py similarity index 100% rename from src/overplot_MRK463E.py rename to package/overplot_MRK463E.py diff --git a/package/src/__init__.py b/package/src/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/utils/analysis.py b/package/src/analysis.py similarity index 96% rename from utils/analysis.py rename to package/src/analysis.py index 6daecb8..1cb3bc9 100755 --- a/utils/analysis.py +++ b/package/src/analysis.py @@ -30,7 +30,7 @@ except get_error as err: if fits_path is not None: from astropy.io import fits - from src.lib.plots import pol_map + from lib.plots import pol_map Stokes_UV = fits.open(fits_path) p = pol_map(Stokes_UV, SNRp_cut=SNRp_cut, SNRi_cut=SNRi_cut, flux_lim=flux_lim) diff --git a/utils/comparison_Kishimoto.py b/package/src/comparison_Kishimoto.py similarity index 100% rename from utils/comparison_Kishimoto.py rename to package/src/comparison_Kishimoto.py diff --git a/utils/get_cdelt.py b/package/src/get_cdelt.py similarity index 100% rename from utils/get_cdelt.py rename to package/src/get_cdelt.py