diff --git a/decomp_checks/README b/decomp_checks/README new file mode 100644 index 0000000..cbd1b34 --- /dev/null +++ b/decomp_checks/README @@ -0,0 +1,24 @@ +class MaskClass + +This class creates masks for each individual parameter for each scousepy decomposition. It then creates super masks that are the combinations of the parameter-specific masks + +The class assumes the following directory structure: + +the MaskClass script (mask_class.py) is located in a sub-directory of the directory pulled from github. The path to the directory containing the stage_4 data files from the directory containing the class is then: '../HC3N_TP_7m_12m_feather/stage_4/' + +The class outputs the masks to the directory 'mask_dir/' + +The class can be run using the test_class.py script. When prompted, enter True for AIC and PB and if all files are in the expected directories a mask should be created for AIC and PB as well as a super mask combining the two parameter-specific masks. + +Below is a description of the functions within the class: + +FUNCTIONS: + +__init__(): Initialized the MaskClass instance. As part of the initialization, a dictionary of key-value pairs is initialized which is used to determine which parameters to use for the mask generation. The possible mask options are based on the scousepy output parameters (AIC, rms, chisq, etc) plus the primary beam response (PB). By default, all of the key value pairs are false (e.g. {'AIC' : False}), but as part of the initialization the user will be prompted whether they wish to mask any of the parameters. If the user types True for any of the parameters, the corresponding dictionary entry will be updated to True, and that parameter will be used to make a mask. + +make_mask(): For each scousepy decomposition generate a mask for each desired parameter. These parameter-specific masks are saved in the mask_dir/ folder. Thresholds used to create the masks are initialized using the __init_() function. + +super_mask(v0_mlist,v1_mlist,v2_mlist): For each decomposition, create a super mask incorporating all of the parameter-specific masks generated from the make_mask() function. + +TBD: +Add a function which applies the masks to the stage 4 .dat files. This function will then save the new .dat files after masking to the data directory (HC3N_TP_7m_12m_feather/stage_4/). diff --git a/decomp_checks/archive/AIC_comp.py b/decomp_checks/archive/AIC_comp.py new file mode 100644 index 0000000..ed9b45a --- /dev/null +++ b/decomp_checks/archive/AIC_comp.py @@ -0,0 +1,64 @@ +""" +@author: Dylan Pare +@date: 03/27/2023 +@name: AIC_comp.py +@description: This code compares the AIC values obtained from the different scousepy decompositions to determine +the most likely of the three. It also determines the likelihood of the other models to minimize information loss. +""" + +# Import needed packages +import numpy as np +import matplotlib.pyplot as plt +from astropy.io import fits +from aplpy import FITSFigure +import os +import glob + +# Now read in the AIC files from the decomposistions after masking +cwd = os.getcwd() +parent = os.path.dirname(cwd) +dat_dir = parent + '/HC3N_TP_7m_12m_feather/stage_4/' +aic_files = glob.glob(dat_dir + '*AIC*_mask.fits') +aic_names = [os.path.basename(x) for x in glob.glob(dat_dir + '*AIC*_mask.fits')] +decomp_v = ['v0','v1','v2'] + +# Now iterate over the non-masked lines of sight and determine which AIC value is the minimum +test_im = fits.open(aic_files[0]) +test_dat = test_im[0].data +header = test_im[0].header + +v0_v1_diff = np.empty(test_dat.shape) +v0_v2_diff = np.empty(test_dat.shape) +v1_v2_diff = np.empty(test_dat.shape) + +v0_v1_diff[:] = np.nan +v1_v2_diff[:] = np.nan +v0_v2_diff[:] = np.nan + +v0_head = fits.open(aic_files[0])[0].header +v1_head = fits.open(aic_files[1])[0].header +v2_head = fits.open(aic_files[2])[0].header + +name_v0_v1 = 'v0_v1_AIC_diff.fits' +name_v0_v2 = 'v0_v2_AIC_diff.fits' +name_v1_v2 = 'v1_v2_AIC_diff.fits' + +im_v0 = fits.open(aic_files[0]) +im_v1 = fits.open(aic_files[1]) +im_v2 = fits.open(aic_files[2]) +dat_v0 = im_v0[0].data +dat_v1 = im_v1[0].data +dat_v2 = im_v2[0].data +for x in range(len(test_dat)): + for y in range(len(test_dat[0])): + if (np.isnan(dat_v0[x][y]) != True) and (np.isnan(dat_v1[x][y]) != True) and (np.isnan(dat_v2[x][y]) != True): + v0_v1_diff[x][y] = dat_v0[x][y] - dat_v1[x][y] + v0_v2_diff[x][y] = dat_v0[x][y] - dat_v2[x][y] + v1_v2_diff[x][y] = dat_v1[x][y] - dat_v2[x][y] + +# Save the minimum AIC maps for future use +fits.writeto(name_v0_v1,data=v0_v1_diff,header=v0_head,overwrite=True) +fits.writeto(name_v0_v2,data=v0_v2_diff,header=v1_head,overwrite=True) +fits.writeto(name_v1_v2,data=v1_v2_diff,header=v2_head,overwrite=True) + + diff --git a/decomp_checks/archive/apply_mask.py b/decomp_checks/archive/apply_mask.py new file mode 100644 index 0000000..fcb0ba3 --- /dev/null +++ b/decomp_checks/archive/apply_mask.py @@ -0,0 +1,86 @@ +""" +@author: Dylan Pare +@date: 03/14/2023 +@name: apply_mask.py +@description: apply the masks to the decomposition data files, removing the entries for the pixels that are marked by the mask file. +""" + +# Import needed packages +import numpy as np +from astropy.io import ascii +from astropy.table import Table +from astropy.io import fits +import glob +import os + +# Get paths to needed files and file basenames +path = os.getcwd() +parent = os.path.dirname(path) +dat_dir = parent + '/HC3N_TP_7m_12m_feather/stage_4/' +dat_files = glob.glob(dat_dir + '*refit.dat') +mask_files = glob.glob(dat_dir + '*mask.fits') +dat_names = [os.path.basename(x) for x in glob.glob(dat_dir+'*refit.dat')] + +# Read in the mask and the corresponding data file and remove the data entries that have been marked by the mask +for file in range(len(dat_files)): + tab_dat = Table.read(dat_files[file],format='ascii') + ncomps = tab_dat['ncomps'] + x = tab_dat['x'] + y = tab_dat['y'] + amp = tab_dat['amplitude'] + err_amp = tab_dat['err amplitude'] + shift = tab_dat['shift'] + err_shift = tab_dat['err shift'] + width = tab_dat['width'] + err_width = tab_dat['err width'] + rms = tab_dat['rms'] + residual = tab_dat['residual'] + chisq = tab_dat['chisq'] + dof = tab_dat['dof'] + redchisq = tab_dat['redchisq'] + aic = tab_dat['aic'] + + mask_name = mask_files[file] + prefix = dat_names[file][:-4] + mask_im = fits.open(mask_name) + mask_dat = mask_im[0].data + mask_head = mask_im[0].header + + ncomps_thresh = [] + x_thresh = [] + y_thresh = [] + amp_thresh = [] + err_amp_thresh = [] + shift_thresh = [] + err_shift_thresh = [] + width_thresh = [] + err_width_thresh = [] + rms_thresh = [] + residual_thresh = [] + chisq_thresh = [] + dof_thresh = [] + redchisq_thresh = [] + aic_thresh = [] + for xi in range(len(x)): + if mask_dat[int(y[xi])][int(x[xi])] == 0: + ncomps_thresh.append(ncomps[xi]) + x_thresh.append(x[xi]) + y_thresh.append(y[xi]) + amp_thresh.append(amp[xi]) + err_amp_thresh.append(err_amp[xi]) + shift_thresh.append(shift[xi]) + err_shift_thresh.append(err_shift[xi]) + width_thresh.append(width[xi]) + err_width_thresh.append(err_width[xi]) + rms_thresh.append(rms[xi]) + residual_thresh.append(residual[xi]) + chisq_thresh.append(chisq[xi]) + dof_thresh.append(dof[xi]) + redchisq_thresh.append(redchisq[xi]) + aic_thresh.append(aic[xi]) + + print(len(x)) + print(len(x_thresh)) + + tab_dat_new = Table([ncomps_thresh,x_thresh,y_thresh,amp_thresh,err_amp_thresh,shift_thresh,err_shift_thresh,width_thresh,err_width_thresh,rms_thresh,residual_thresh,chisq_thresh,dof_thresh,redchisq_thresh,aic_thresh],names=('ncomps','x','y','amplitude','err amplitude','shift','err shift','width','err width','rms','residual','chisq','dof','redchisq','aic')) + tab_dat_new.write(prefix+'_mask.dat',format='ascii') diff --git a/decomp_checks/archive/make_mask.py b/decomp_checks/archive/make_mask.py new file mode 100644 index 0000000..e9f53ec --- /dev/null +++ b/decomp_checks/archive/make_mask.py @@ -0,0 +1,44 @@ +""" +@author: Dylan Pare +@name: make_mask.py +@date: 03/14/2023 +@description: Create a mask to remove artifact pixels caused by e.g. edge pixels. +""" + +# Import needed packages +import numpy as np +from astropy.io import fits +import os +import glob + +# Get set of residstd files generated by scousepy +path = os.getcwd() +parent = os.path.dirname(path) +print(parent) +fits_dir = parent + '/HC3N_TP_7m_12m_feather/stage_4/' +rms_files = glob.glob(fits_dir + '*rms*refit.fits') +rms_names = [os.path.basename(x) for x in glob.glob(fits_dir+'*rms*refit.fits')] +pb_file = parent + '/SgrB2_PrimaryBeam_new.fits' +pb_im = fits.open(pb_file) +pb_data = pb_im[0].data +pb_thresh = 0.5 + +mask_prefix = ['v0','v1','v2'] + +# Read in PB file and create a mask file based on PB response threshold +for file in range(len(rms_files)): + + print(rms_names[file]) + rms_file = rms_files[file] + rms_image = fits.open(rms_file) + rms_data = rms_image[0].data + rms_header = rms_image[0].header + + mask = np.zeros(rms_data.shape) + for x in range(len(rms_data)): + for y in range(len(rms_data[0])): + if (pb_data[x][y] <= pb_thresh) or (np.isnan(pb_data[x][y]) == True): + mask[x][y] = 1.0 + + fits.writeto(mask_prefix[file]+'_mask.fits',data=mask,header=rms_header,overwrite=True) + diff --git a/decomp_checks/archive/plot_maskdist.py b/decomp_checks/archive/plot_maskdist.py new file mode 100644 index 0000000..382e5a3 --- /dev/null +++ b/decomp_checks/archive/plot_maskdist.py @@ -0,0 +1,92 @@ +""" +@author: Dylan Pare +@date: 03/14/2023 +@name: plot_maskdist.py +@description: Plot the distributions obtained from the masked scouse data set. +""" + +# Import needed packages +import numpy as np +from astropy.io import fits +from astropy.io import ascii +from astropy.table import Table +import glob +import os + +# Get paths to needed files and file basenames +path = os.getcwd() +parent = os.path.dirname(path) +dat_dir = parent + '/HC3N_TP_7m_12m_feather/stage_4/' +dat_files = glob.glob(dat_dir + '*_mask.dat') +mask_files = [dat_dir + 'v0_mask.fits',dat_dir+'v1_mask.fits',dat_dir+'v2_mask.fits'] +dat_names = [os.path.basename(x) for x in glob.glob(dat_dir+'*_mask.dat')] +aic_files = glob.glob(dat_dir + '*AIC*.refit.fits') +chisq_files = glob.glob(dat_dir + '*chisq*.refit.fits') +redchisq_files = glob.glob(dat_dir + '*redchisq*.refit.fits') +ncomps_files = glob.glob(dat_dir + '*ncomps*.refit.fits') +rms_files = glob.glob(dat_dir + '*rms*.refit.fits') + +# Read in the data file and plot new distributions based on the masked data +for file in range(len(dat_files)): + mask_name = mask_files[file] + mask_im = fits.open(mask_name) + mask_data = mask_im[0].data + + aic_name = aic_files[file] + chisq_name = chisq_files[file] + redchisq_name = redchisq_files[file] + ncomps_name = ncomps_files[file] + rms_name = rms_files[file] + AIC_head = fits.open(aic_name)[0].header + chisq_head = fits.open(chisq_name)[0].header + redchisq_head = fits.open(redchisq_name)[0].header + ncomps_head = fits.open(ncomps_name)[0].header + rms_head = fits.open(rms_name)[0].header + + aic_prefix = aic_name[:-5] + chisq_prefix = chisq_name[:-5] + redchisq_prefix = redchisq_name[:-5] + ncomps_prefix = ncomps_name[:-5] + rms_prefix = rms_name[:-5] + + AIC = np.empty(mask_data.shape) + AIC[:] = np.nan + chisq = np.empty(mask_data.shape) + chisq[:] = np.nan + ncomps = np.empty(mask_data.shape) + ncomps[:] = np.nan + redchisq = np.empty(mask_data.shape) + redchisq[:] = np.nan + rms = np.empty(mask_data.shape) + rms[:] = np.nan + + tab_dat = Table.read(dat_files[file],format='ascii') + ncomps_ind = tab_dat['ncomps'] + x = tab_dat['x'] + y = tab_dat['y'] + amp_ind = tab_dat['amplitude'] + err_amp_ind = tab_dat['err amplitude'] + shift_ind = tab_dat['shift'] + err_shift_ind = tab_dat['err shift'] + width_ind = tab_dat['width'] + err_width_ind = tab_dat['err width'] + rms_ind = tab_dat['rms'] + residual_ind = tab_dat['residual'] + chisq_ind = tab_dat['chisq'] + dof_ind = tab_dat['dof'] + redchisq_ind = tab_dat['redchisq'] + aic_ind = tab_dat['aic'] + + for xi in range(len(x)): + AIC[int(y[xi])][int(x[xi])] = aic_ind[xi] + chisq[int(y[xi])][int(x[xi])] = chisq_ind[xi] + ncomps[int(y[xi])][int(x[xi])] = ncomps_ind[xi] + redchisq[int(y[xi])][int(x[xi])] = redchisq_ind[xi] + rms[int(y[xi])][int(x[xi])] = rms_ind[xi] + + # Now save the new distributions using the masked data + fits.writeto(aic_prefix+'_mask.fits',AIC,AIC_head,overwrite=True) + fits.writeto(chisq_prefix+'_mask.fits',chisq,chisq_head,overwrite=True) + fits.writeto(ncomps_prefix+'_mask.fits',ncomps,ncomps_head,overwrite=True) + fits.writeto(redchisq_prefix+'_mask.fits',redchisq,redchisq_head,overwrite=True) + fits.writeto(rms_prefix+'_mask.fits',rms,rms_head,overwrite=True) diff --git a/decomp_checks/mask_class.py b/decomp_checks/mask_class.py new file mode 100644 index 0000000..d45aba5 --- /dev/null +++ b/decomp_checks/mask_class.py @@ -0,0 +1,164 @@ +""" +@author: Dylan Pare +@date: 03/30/2023 +@name: mask_class.py +@description: Make a generate class allowing for the construction of different masks based on +different parameters. Also have a routine to combine all generated masks to produce a final mask +for distribution. +""" +# Import needed packages +from astropy.io import fits +import numpy as np +import glob +import os + +class MaskClass: + """Class for making and combining various masks.""" + + def __init__(self): + """Initialization function for MaskClass. + Inputs: + None + Outputs: + None + """ + print("Please specify which masks to use [Mark True to use, else leave blank]: ") + #prompt_thresh = input("Do you want to update thresholds? [defaults: ncomps = 8, AIC = -360.0, PB = 0.5, rms = 3.0, chisq = 10.0, redchisq = 1.0, residstd = 1.0]: ") + + params = {'ncomps': False,\ + 'AIC': False,\ + 'PB': False,\ + 'rms' : False,\ + 'chisq' : False,\ + 'redchisq' : False,\ + 'residstd' : False} + + thresh = {'ncomps' : 8,\ + 'AIC' : -360.0,\ + 'PB' : 0.5,\ + 'rms' : 3.0,\ + 'chisq' : 10.0,\ + 'redchisq' : 1.0,\ + 'residstd' : 1.0} + + ncomps = input("ncomps: ") + AIC = input("AIC: ") + PB = input("PB: ") + rms = input("rms: ") + chisq = input('chisq: ') + redchisq = input('redchisq: ') + residstd = input('residstd: ') + + if ncomps == "True": + params.update({"ncomps" : True}) + if AIC == "True": + params.update({"AIC" : True}) + if PB == "True": + params.update({"PB" : True}) + if rms == "True": + params.update({"rms" : True}) + if chisq == "True": + params.update({"chisq" : True}) + if redchisq == "True": + params.update({"redchisq" : True}) + if residstd == "True": + params.update({"residstd" : True}) + + print("Default threshold levels for the different parameters: ncomps = 8, AIC = -360.0, PB = 0.5, rms = 3.0, chisq = 10.0, residstd = 1.0 [all except for AIC and PB are placeholders]") + + parent = os.path.dirname(os.getcwd()) + data_dir = parent + '/HC3N_TP_7m_12m_feather/stage_4/' + mask_dir = os.getcwd() + '/mask_dir/' + decomps = ['v0','v1','v2'] + + self.params = params + self.thresh = thresh + self.data_dir = data_dir + self.mask_dir = mask_dir + self.decomps = decomps + return + + def make_mask(self): + """Make a mask for each parameter flagged as True. + Inputs: + None + Outputs: + mask_list = aarray of the names for the individual masks generated + """ + v0_mlist = [] + v1_mlist = [] + v2_mlist = [] + value_list = list(self.thresh.values()) + do_mask = list(self.params.values()) + key_list = list(self.params.keys()) + for decomp in range(len(sorted(self.decomps))): + for el in range(len(value_list)): + if do_mask[el] == True: + if key_list[el] != 'PB': + print("Making mask based on " + key_list[el] + ' ' + self.decomps[decomp] + " threshold.") + print('') + data_name = self.data_dir + 'stage_4_' + key_list[el] + '.' + self.decomps[decomp] + '.refit.fits' + data_im = fits.open(data_name)[0].data + cutoff = value_list[el] + mask_name = key_list[el] + '_' + self.decomps[decomp] + '_mask.fits' + mask_dat = np.zeros(data_im.shape) + mask_head = fits.open(data_name)[0].header + for ra in range(len(data_im)): + for dec in range(len(data_im[0])): + if data_im[ra][dec] <= cutoff: + mask_dat[ra][dec] = 1 + fits.writeto(self.mask_dir + mask_name,data=mask_dat,header=mask_head,overwrite=True) + else: + print("Making mask based on " + key_list[el] + ' ' + self.decomps[decomp] + " threshold.") + print('') + data_name = self.data_dir + 'SgrB2_PB.fits' + data_im = fits.open(data_name)[0].data + cutoff = value_list[el] + mask_name = key_list[el] + '_' + self.decomps[decomp] + '_mask.fits' + mask_dat = np.zeros(data_im.shape) + mask_head = fits.open(data_name)[0].header + for ra in range(len(data_im)): + for dec in range(len(data_im[0])): + if (data_im[ra][dec] <= cutoff) or (np.isnan(data_im[ra][dec]) == True): + mask_dat[ra][dec] = 1 + fits.writeto(self.mask_dir + mask_name,data=mask_dat,header=mask_head,overwrite=True) + if decomp == 0: + v0_mlist.append(mask_name) + elif decomp == 1: + v1_mlist.append(mask_name) + else: + v2_mlist.append(mask_name) + return v0_mlist, v1_mlist, v2_mlist + + def super_mask(self,v0_mlist,v1_mlist,v2_mlist): + """Create a super mask incorporating all of the sub-masks generated by the make_mask routine. + Inputs: + mask_list = array of string names referencing the sub-masks from make_mask. + """ + print("Making super masks based on mask parameters used.") + mask_list = np.array([]) + for decomp in range(len(sorted(self.decomps))): + if decomp == 0: + mask_list = v0_mlist + elif decomp == 1: + mask_list = v1_mlist + else: + mask_list = v2_mlist + + # Iterate over the masks to create a super mask + m_name = self.mask_dir + mask_list[0] + m_shape = fits.open(m_name)[0].data + super_mask = np.zeros(m_shape.shape) + for mask in range(len(mask_list)): + mask_name = self.mask_dir + mask_list[mask] + mask_dat = fits.open(mask_name)[0].data + mask_hdr = fits.open(mask_name)[0].header + for ra in range(len(mask_dat)): + for dec in range(len(mask_dat[0])): + if super_mask[ra][dec] == 1: + continue + elif mask_dat[ra][dec] == 1: + super_mask[ra][dec] = 1 + + fits.writeto(self.mask_dir + 'super_' + self.decomps[decomp] + '_mask.fits',data=super_mask,header=mask_hdr,overwrite=True) + return diff --git a/decomp_checks/reproject_pb.py b/decomp_checks/reproject_pb.py new file mode 100644 index 0000000..9cc6003 --- /dev/null +++ b/decomp_checks/reproject_pb.py @@ -0,0 +1,30 @@ +""" +@author: Dylan Pare +@date: 03/24/2023 +@name: reproject_pb.py +@description: Rescale the primary beam obtained from Ginsburg (on Globus) to +match WCS with the SgrB2 pilot data. +""" + +# Import needed packages +import numpy as np +from reproject import reproject_interp +from astropy.io import fits +from astropy import wcs +import glob +import os + +# Reproject the PB distribution to match WCS construction of the Sgr B2 data +path = os.getcwd() +parent = os.path.dirname(path) +pb_name = parent + '/SgrB2_PrimaryBeam.fits' +WCS_name = parent + '/HC3N_TP_7m_12m_feather/stage_4/stage_4_AIC.fits' # Any image from stage 4 to get WCS +pb_im = fits.open(pb_name) +WCS_im = fits.open(WCS_name) +WCS_head = WCS_im[0].header + +pb_im_new, footprint = reproject_interp(pb_im,WCS_head) + +# Write out the new primary beam image with the appropriate pixel / WCS scaling +fits.writeto('Sgr_PrimaryBeam_new.fits',data=pb_im_new,header=WCS_head) + diff --git a/decomp_checks/test_class.py b/decomp_checks/test_class.py new file mode 100644 index 0000000..a5ed3f7 --- /dev/null +++ b/decomp_checks/test_class.py @@ -0,0 +1,16 @@ +""" +@author: Dylan Pare +@date: 03/30/2023 +@name: test_class.py +@description: Test the make_mask class. +""" + +# Import needed packages +import mask_class + +test = mask_class.MaskClass() +print(test.params) + +v0_mnames, v1_mnames, v2_mnames = test.make_mask() + +test.super_mask(v0_mnames, v1_mnames, v2_mnames)