diff --git a/src/FOC_reduction.py b/src/FOC_reduction.py index ecf559a..f27f145 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"): # Background estimation error_sub_type = 'freedman-diaconis' #sqrt, sturges, rice, scott, freedman-diaconis (default) or shape (example (51,51)) subtract_error = 1.00 - display_error = True + display_error = False # Data binning rebin = True @@ -93,15 +93,13 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): data_array, headers = proj_fits.get_obs_data(infiles, data_folder=data_folder, compute_flux=True) figname = "_".join([target,"FOC"]) - if smoothing_FWHM is None: - if px_scale in ['px','pixel','pixels']: - figtype = "".join(["b_",str(pxsize),'px']) - elif px_scale in ['arcsec','arcseconds','arcs']: - figtype = "".join(["b_","{0:.2f}".format(pxsize).replace(".",""),'arcsec']) + if rebin: + if not px_scale in ['full']: + figtype = "".join(["b","{0:.2f}".format(pxsize),px_scale]) #additionnal informations else: figtype = "full" - else: - figtype = "_".join(["".join([s[0] for s in smoothing_function.split("_")]),"".join(["{0:.2f}".format(smoothing_FWHM).replace(".",""),smoothing_scale])]) #additionnal informations + if not smoothing_FWHM is None: + figtype += "_"+"".join(["".join([s[0] for s in smoothing_function.split("_")]),"{0:.2f}".format(smoothing_FWHM),smoothing_scale]) #additionnal informations if align_center is None: figtype += "_not_aligned" @@ -123,7 +121,7 @@ def main(target=None, proposal_id=None, infiles=None, output_dir="./data"): data_array, error_array, headers, data_mask = proj_red.align_data(data_array, headers, error_array=error_array, background=background, upsample_factor=10, ref_center=align_center, return_shifts=False) if display_align: - proj_plots.plot_obs(data_array, headers, vmin=data_array[data_array>0.].min()*headers[0]['photflam'], vmax=data_array[data_array>0.].max()*headers[0]['photflam'], savename="_".join([figname,"center",str(align_center)]), plots_folder=plots_folder) + proj_plots.plot_obs(data_array, headers, vmin=data_array[data_array>0.].min()*headers[0]['photflam'], vmax=data_array[data_array>0.].max()*headers[0]['photflam'], savename="_".join([figname,str(align_center)]), plots_folder=plots_folder) # Rebin data to desired pixel size. if rebin: diff --git a/src/comparison_Kishimoto.py b/src/comparison_Kishimoto.py index ebdb05b..1b645bd 100755 --- a/src/comparison_Kishimoto.py +++ b/src/comparison_Kishimoto.py @@ -17,8 +17,8 @@ root_dir = path_join('/home/t.barnouin/Documents/Thesis/HST') root_dir_K = path_join(root_dir,'Kishimoto','output') root_dir_S = path_join(root_dir,'FOC_Reduction','output') root_dir_data_S = path_join(root_dir,'FOC_Reduction','data','NGC1068','5144') -root_dir_plot_S = path_join(root_dir,'FOC_Reduction','plots','NGC1068','5144') -filename_S = "NGC1068_FOC_b_10px.fits" +root_dir_plot_S = path_join(root_dir,'FOC_Reduction','plots','NGC1068','5144','notaligned') +filename_S = "NGC1068_FOC_b10.00pixel_not_aligned.fits" plt.rcParams.update({'font.size': 15}) SNRi_cut = 30. @@ -140,7 +140,8 @@ fig_dif_pa.savefig(path_join(root_dir_plot_S,"NGC1068_K_polang_diff.png"),bbox_i ##### ###display both polarization maps to check consistency ##### -fig = plt.figure(num="Polarization maps comparison") +#plt.rcParams.update({'font.size': 15}) +fig = plt.figure(num="Polarization maps comparison",figsize=(10,10)) ax = fig.add_subplot(111, projection=wcs) fig.subplots_adjust(right=0.85) cbar_ax = fig.add_axes([0.88, 0.12, 0.01, 0.75]) @@ -164,7 +165,6 @@ ax.coords[1].set_ticklabel_position('l') #ax.axis('equal') cbar = plt.colorbar(im0, cax=cbar_ax, label=r"$F_{\lambda}$ [$ergs \cdot cm^{-2} \cdot s^{-1} \cdot \AA^{-1}$]") -#plt.rcParams.update({'font.size': 8}) ax.legend(loc='upper right') fig.savefig(path_join(root_dir_plot_S,"NGC1068_K_comparison.png"),bbox_inches="tight",dpi=300) diff --git a/src/lib/query.py b/src/lib/query.py index 4a819a9..7424dcb 100755 --- a/src/lib/query.py +++ b/src/lib/query.py @@ -137,7 +137,7 @@ 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 np.all(tab['target_name']==tab['target_name'][0]): + 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] diff --git a/src/lib/reduction.py b/src/lib/reduction.py index aacfc87..c2c8563 100755 --- a/src/lib/reduction.py +++ b/src/lib/reduction.py @@ -575,8 +575,7 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, if instr == 'FOC': HST_aper = 2400. # HST aperture in mm Dxy_arr = np.ones((data_array.shape[0],2)) - for i, enum in enumerate(list(zip(data_array, error_array, headers))): - image, error, header = enum + for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): # Get current pixel size w = WCS(header).deepcopy() new_header = deepcopy(header) @@ -592,8 +591,7 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, raise ValueError("'{0:s}' invalid scale for binning.".format(scale)) new_shape = np.ceil(min(image.shape/Dxy_arr,key=lambda x:x[0]+x[1])).astype(int) - for i, enum in enumerate(list(zip(data_array, error_array, headers))): - image, error, header = enum + for i, (image, error, header) in enumerate(list(zip(data_array, error_array, headers))): # Get current pixel size w = WCS(header).deepcopy() new_header = deepcopy(header) @@ -617,21 +615,12 @@ def rebin_array(data_array, error_array, headers, pxsize, scale, if operation.lower() in ["mean", "average", "avg"]: new_error = np.sqrt(bin_ndarray(error**2, new_shape=new_shape, operation='average')) - #new_error[mask] = np.sqrt(bin_ndarray(error**2*image, - # new_shape=new_shape, operation='average')[mask]/sum_image[mask]) - #new_error[mask] = np.sqrt(bin_ndarray(error**2, - # new_shape=new_shape, operation='average')[mask]) else: new_error = np.sqrt(bin_ndarray(error**2, new_shape=new_shape, operation='sum')) - #new_error[mask] = np.sqrt(bin_ndarray(error**2*image, - # new_shape=new_shape, operation='sum')[mask]/sum_image[mask]) - #new_error[mask] = np.sqrt(bin_ndarray(error**2, - # new_shape=new_shape, operation='sum')[mask]) rebinned_error.append(np.sqrt(rms_image**2 + new_error**2)) # Update header - #nw = w.slice((np.s_[::Dxy[0]], np.s_[::Dxy[1]])) nw = w.deepcopy() nw.wcs.cdelt *= Dxy nw.wcs.crpix /= Dxy @@ -762,21 +751,20 @@ def align_data(data_array, headers, error_array=None, background=None, # Initialize rescaled images to background values rescaled_error[i] *= 0.01*background[i] # Get shifts and error by cross-correlation to ref_data - shift, error, phase_diff = phase_cross_correlation(ref_data/ref_data.max(), image/image.max(), - upsample_factor=upsample_factor) + if do_shift: + shift, error, _ = phase_cross_correlation(ref_data/ref_data.max(), image/image.max(), + upsample_factor=upsample_factor) + else: + shift = pol_shift[headers[i]['filtnam1'].lower()] + error = sigma_shift[headers[i]['filtnam1'].lower()] # Rescale image to requested output rescaled_image[i,res_shift[0]:res_shift[0]+shape[1], res_shift[1]:res_shift[1]+shape[2]] = deepcopy(image) rescaled_error[i,res_shift[0]:res_shift[0]+shape[1], res_shift[1]:res_shift[1]+shape[2]] = deepcopy(error_array[i]) # Shift images to align - if do_shift: - rescaled_image[i] = sc_shift(rescaled_image[i], shift, order=1, cval=0.) - rescaled_error[i] = sc_shift(rescaled_error[i], shift, order=1, cval=background[i]) - else: - shift = pol_shift[headers[i]['filtnam1'].lower()] - rescaled_image[i] = sc_shift(rescaled_image[i], shift, order=1, cval=0.) - rescaled_error[i] = sc_shift(rescaled_error[i], shift, order=1, cval=background[i]) + rescaled_image[i] = sc_shift(rescaled_image[i], shift, order=1, cval=0.) + rescaled_error[i] = sc_shift(rescaled_error[i], shift, order=1, cval=background[i]) curr_mask = sc_shift(res_mask, shift, order=1, cval=False) mask_vertex = clean_ROI(curr_mask) @@ -792,9 +780,6 @@ def align_data(data_array, headers, error_array=None, background=None, #sum quadratically the errors rescaled_error[i] = np.sqrt(rescaled_error[i]**2 + error_shift**2) - #if i==1: - #np.savetxt("output/s_shift.txt",error_shift) - shifts.append(shift) errors.append(error)