From 762b85772089e906b5e5728456d0cac738b4c5c6 Mon Sep 17 00:00:00 2001 From: Thibault Barnouin Date: Tue, 12 Mar 2024 17:06:28 +0100 Subject: [PATCH] bug fix for crop_stokes.write_to and background subtraction --- src/FOC_reduction.py | 6 +++--- src/lib/background.py | 10 +++------- src/lib/plots.py | 3 --- src/lib/reduction.py | 2 +- 4 files changed, 7 insertions(+), 14 deletions(-) diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index 3afaf6e..2683df8 100755 --- a/src/FOC_reduction.py +++ b/src/FOC_reduction.py @@ -33,7 +33,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 = 1.00 + subtract_error = 0.50 display_bkg = False # Data binning @@ -181,13 +181,13 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data", crop= figtype += "_crop" stokescrop = proj_plots.crop_Stokes(deepcopy(Stokes_test), norm=LogNorm()) stokescrop.crop() - stokescrop.writeto("/".join([data_folder, "_".join([figname, figtype+".fits"])])) + stokescrop.write_to("/".join([data_folder, "_".join([figname, figtype+".fits"])])) Stokes_test, data_mask, headers = stokescrop.hdul_crop, stokescrop.data_mask, [dataset.header for dataset in stokescrop.hdul_crop] print("F_int({0:.0f} Angs) = ({1} ± {2})e{3} ergs.cm^-2.s^-1.Angs^-1".format(headers[0]['photplam'], *proj_plots.sci_not( Stokes_test[0].data[data_mask].sum()*headers[0]['photflam'], np.sqrt(Stokes_test[3].data[0, 0][data_mask].sum())*headers[0]['photflam'], 2, out=int))) print("P_int = {0:.1f} ± {1:.1f} %".format(headers[0]['p_int']*100., np.ceil(headers[0]['p_int_err']*1000.)/10.)) - print("PA_int = {0:.1f} ±t {1:.1f} °".format(headers[0]['pa_int'], np.ceil(headers[0]['pa_int_err']*10.)/10.)) + print("PA_int = {0:.1f} ± {1:.1f} °".format(headers[0]['pa_int'], np.ceil(headers[0]['pa_int_err']*10.)/10.)) # Background values 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))) diff --git a/src/lib/background.py b/src/lib/background.py index 0316c9c..8d07c25 100755 --- a/src/lib/background.py +++ b/src/lib/background.py @@ -240,7 +240,7 @@ def bkg_fit(data, error, mask, headers, subtract_error=True, display=False, save # Substract background if subtract_error > 0: n_data_array[i][mask] = n_data_array[i][mask] - bkg - n_data_array[i][np.logical_and(mask, n_data_array[i] <= 0.01*bkg)] = 0.01*bkg + n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3*bkg)] = 1e-3*bkg std_bkg[i] = image[np.abs(image-bkg)/bkg < 1.].std() background[i] = bkg @@ -327,10 +327,6 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis histograms.append(hist) binning.append(np.exp(bin_centers(bin_edges))) - # Take the background as the count-rate with the maximum number of pixels - # hist_max = binning[-1][np.argmax(hist)] - # bkg = np.sqrt(np.sum(image[np.abs(image-hist_max)/hist_max<0.5]**2)/image[np.abs(image-hist_max)/hist_max<0.5].size) - # Fit a gaussian to the log-intensity histogram bins_stdev = binning[-1][hist > hist.max()/2.] stdev = bins_stdev[-1]-bins_stdev[0] @@ -346,7 +342,7 @@ def bkg_hist(data, error, mask, headers, sub_type=None, subtract_error=True, dis # Substract background if subtract_error > 0: n_data_array[i][mask] = n_data_array[i][mask] - bkg - n_data_array[i][np.logical_and(mask, n_data_array[i] < 0.)] = 0. + n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3*bkg)] = 1e-3*bkg std_bkg[i] = image[np.abs(image-bkg)/bkg < 1.].std() background[i] = bkg @@ -443,7 +439,7 @@ def bkg_mini(data, error, mask, headers, sub_shape=(15, 15), subtract_error=True # Substract background if subtract_error > 0.: n_data_array[i][mask] = n_data_array[i][mask] - bkg - n_data_array[i][np.logical_and(mask, n_data_array[i] <= 0.01*bkg)] = 0.01*bkg + n_data_array[i][np.logical_and(mask, n_data_array[i] <= 1e-3*bkg)] = 1e-3*bkg std_bkg[i] = image[np.abs(image-bkg)/bkg < 1.].std() background[i] = bkg diff --git a/src/lib/plots.py b/src/lib/plots.py index f149235..5a141ec 100755 --- a/src/lib/plots.py +++ b/src/lib/plots.py @@ -1438,9 +1438,6 @@ class crop_Stokes(crop_map): def data_mask(self): return self.hdul_crop[-1].data.astype(int) - def write_to(self, filename): - save_Stokes(self.hdul_crop['I_stokes'], self.hdul_crop['Q_stokes'], self.hdul_crop['U_stokes'], self.hdul_crop['IQU_cov_matrix'], self.hdul_crop['Pol_deg'], self.hdul_crop['Pol_deg_debiased'], self.hdul_crop['Pol_deg_err'], self.hdul_crop['Pol_deg_err_Poisson_noise'], self.hdul_crop['Pol_ang'], self.hdul_crop['Pol_ang_err'], self.hdul_crop['Pol_ang_err_Poisson_noise'], [hdu.header for hdu in self.hdul_crop], self.hdul_crop['data_mask'], filename, data_folder="", return_hdul=False) - class image_lasso_selector(object): def __init__(self, img, fig=None, ax=None): diff --git a/src/lib/reduction.py b/src/lib/reduction.py index 4ea48da..6f504bf 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -164,7 +164,7 @@ def bin_ndarray(ndarray, new_shape, operation='sum'): [342 350 358 366 374]] """ - if not operation.lower() in ['sum', 'mean', 'average', 'avg']: + if operation.lower() not in ['sum', 'mean', 'average', 'avg']: raise ValueError("Operation not supported.") if ndarray.ndim != len(new_shape): raise ValueError("Shape mismatch: {} -> {}".format(ndarray.shape,