Initial masking

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/). 

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)
+
+

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')

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)
+

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)