<jupyter_start><jupyter_text>Stitching PipelineIn this notebook the input WiFeS spectra (that has already been processed by the PyWiFeS pipeline) is modified so that the source is chosen, background corrected, the red and blue sides are aligned and stitched together, and the variance and bad pixels are also identified. The final result is made into a fits file ready for input into the calibSpec pipeline. If the observations are of the same object they can be combined into the same fits file, and if you would like you can also get plots of the raw data/selected source spaxels.Pipeline created by <NAME>. The Spectrumv18 object class is from the OzDES RM CalibSpec code written by <NAME>.For this code to work you must input the file names of the pyWiFES processed fits files into two arrays, for the red and blue, ensuring the order is the same so that observations of the same dates are put together.You should also define an output name, and choose the flags to be true or false if you want plots or to combine the fits files. You can also choose the number of pixels to use as the source and the background (these can be viewed in an output plot too)<jupyter_code>from astropy.io import fits
import numpy as np
from scipy.interpolate import interp1d
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.interpolate import UnivariateSpline

from scipy.spatial.distance import pdist, cdist, squareform
from sklearn.gaussian_process import GaussianProcessRegressor, kernels
import matplotlib.pyplot as plt
import sys
from astropy.time import Time
plt.rcParams['text.usetex'] = True

class Spectrumv18(object):
    def __init__(self, filepath=None):
        assert filepath is not None
        self.filepath = filepath
        try:
            self.data = fits.open(filepath)
        except IOError:
            print("Error: file {0} could not be found".format(filepath))
            exit()
        data = fits.open(filepath)
        self.combinedFlux = data[0]
        self.combinedVariance = data[1]
        self.combinedPixels = data[2]
        self.numEpochs = int((np.size(data)-2)/2)+1 
        self.cdelt1 = self.combinedFlux.header['cdelt3']  # Wavelength interval between subsequent pixels
        self.crpix1 = 1 #self.combinedFlux.header['crpix3']
        self.crval1 = self.combinedFlux.header['crval3']
        self.n_pix = self.combinedFlux.header['NAXIS3']
        self.RA = self.combinedFlux.header['RA']
        self.DEC = self.combinedFlux.header['DEC']

        self.fluxCoadd = self.combinedFlux.data
        self.varianceCoadd = self.combinedVariance.data
        self.badpixCoadd = self.combinedPixels.data

        self._wavelength = None
        self._flux = None
        self._variance = None
        self._badpix = None
        self._dates = None
        self._run = None
        self._ext = None
        self._qc = None
        self._exposed = None

    @property
    def wavelength(self):
        """Define wavelength solution."""
        if getattr(self, '_wavelength', None) is None:
            wave = ((np.arange(self.n_pix) - self.crpix1) * self.cdelt1) + self.crval1
            self._wavelength = wave
        return self._wavelength

    @property
    def flux(self):
        if getattr(self, '_flux', None) is None:
            self._flux = np.zeros((len(self.data[0].data), self.numEpochs), dtype=float) #2848 or 5000
            for i in range(self.numEpochs):
                self._flux[:, i] = self.data [i * 2].data     
        return self._flux

    @property
    def variance(self):
        if getattr(self, '_variance', None) is None:
            self._variance = np.zeros((len(self.data[0].data), self.numEpochs), dtype=float)
            for i in range(self.numEpochs):
                self._variance[:, i] = self.data[i * 2 + 1].data
        return self._variance

    @property
    def badpix(self):
        if getattr(self, '_badpix', None) is None:
            self._badpix = np.zeros((len(self.data[0].data), self.numEpochs), dtype=float)
            for i in range(self.numEpochs):
                self._badpix[:, i] = self.data[i * 2 + 2].data
        return self._badpix

    @property
    def dates(self):
        if getattr(self, '_dates', None) is None:
            self._dates = np.zeros(self.numEpochs, dtype=float)
            for i in range(self.numEpochs):
                self._dates[i] = round(Time(self.data[i*3].header['DATE-OBS'], format='isot', scale='utc').mjd,3)
                print('date:')
                print(self._dates)
                # this give Modified Julian Date (UTC) that observation was taken
        return self._dates
    
    @property
    def exposed(self):
        if getattr(self, '_exposed', None) is None:
            self._exposed = []
            for i in range(self.numEpochs):
                self._exposed.append(self.data[i * 2].header['EXPTIME']) 
                # this will give you the exposure time of each observation
        return self._exposed
##CONFIGURATIONS##

PlotFlag = True #Choose whether to make plots
CombineFlag = True #choose whether it's the same object (combine files)
n_source = 1 #Number of pixels to use as the source - the brightest n_source pixels will be used. Unless you have a large source
            #avoid using more than 9 pixels.
n_back = 50 #the least bright n_back pixels will be used. Default 50. To minimise error don't use fewer than 20, or more than 80.
#Now input the file names of the fits files you would like to process. These should be the output of pyWiFeS. Make
#sure they are in the same order of dates.
sourceNamesblue = np.array(['AGNData/T2m3wb-20200530.101942-0029.p11.fits','AGNData/T2m3wb-20200531.103610-0028.p11.fits'])
sourceNamesred = np.array(['AGNData/T2m3wr-20200530.101942-0029.p11.fits','AGNData/T2m3wr-20200531.103610-0028.p11.fits'])
outName = 'AGNData/IRAS09149-6206-'
##PROCESSING##

#First we check that you've added the sourceNames correctly.
if len(sourceNamesblue) !=len(sourceNamesred):
    print('Error - check that your source names are the same length')

for j in np.arange(len(sourceNamesblue)):
    
    blue_spectra = Spectrumv18(sourceNamesblue[j])
    red_spectra = Spectrumv18(sourceNamesred[j])
    n_b = blue_spectra.n_pix
    n_r = red_spectra.n_pix
    ##Let's choose the source pixels (n_source brightest pixels) and load their data
    a =  np.mean(blue_spectra.fluxCoadd, axis = 0)[1:36,1:23]
    if n_source < 1:
        print("Error - n_source must be greater than or equal to 1")
    if n_source ==1:
        ind = np.unravel_index(np.argmax(a, axis=None), a.shape)
        blue = blue_spectra.fluxCoadd[np.arange(0,n_b,1), np.ones(n_b).astype(int)*(ind[0]+1), np.ones(n_b).astype(int)*(ind[1]+1)]
        red = red_spectra.fluxCoadd[np.arange(0,n_r,1), np.ones(n_r).astype(int)*(ind[0]+1), np.ones(n_r).astype(int)*(ind[1]+1)]
        if PlotFlag ==True:
            plt.imshow(a)
            plt.plot(np.array([ind[1]]),np.array([ind[0]]),color = 'r',marker = "s" )
            plt.title('Selected regions of the source'+str(j))
            plt.savefig(outName+'_selected_regions'+str(j))
            plt.clf()
        
    if n_source >1:    
        indices_source = np.unravel_index(np.argsort(a, axis=None)[-n_source:], a.shape)
        if PlotFlag ==True:
            image = np.zeros([36,23])
            image[indices_source]=1
            plt.contour(image, levels = [0.5])
            plt.imshow(a)
            plt.title('Selected regions of the source'+str(j))
            plt.savefig(outName+'_selected_regions'+str(j), dpi = 600)
            plt.clf()
    
        blue = [np.zeros(n_b)]
        for i in np.arange(0,n_source-1,1):
            arr_i = blue_spectra.fluxCoadd[np.arange(0,n_b,1), np.ones(n_b).astype(int)*indices_source[0][i]+1, np.ones(n_b).astype(int)*indices_source[1][i]+1]
            blue = np.append(blue,[arr_i], axis = 0)
        blue = np.mean(blue[1:], axis = 0)
        red = [np.zeros(n_r)]
        for i in np.arange(0,n_source-1,1):
            arr_i = red_spectra.fluxCoadd[np.arange(0,n_r,1), np.ones(n_r).astype(int)*indices_source[0][i]+1, np.ones(n_r).astype(int)*indices_source[1][i]+1]
            red = np.append(red,[arr_i], axis = 0)
        red = np.mean(red[1:], axis = 0)
    
    ##Now let's choose the background pixels by getting the coordinates of the n_back least intense spaxels
    def get_indices_of_k_smallest(arr, k):
        idx = np.argpartition(arr.ravel(), k)
        return tuple(np.array(np.unravel_index(idx, arr.shape))[:, range(min(k, 0), max(k, 0))])
    indices = get_indices_of_k_smallest(a, n_back)

    background_blue = [np.zeros(n_b)]
    for i in np.arange(0,n_back-1,1):
        arr_i = blue_spectra.fluxCoadd[np.arange(0,n_b,1), np.ones(n_b).astype(int)*indices[0][i]+1, np.ones(n_b).astype(int)*indices[1][i]+1]
        background_blue = np.append(background_blue,[arr_i], axis = 0)
    background_blue = np.mean(background_blue[1:], axis = 0)
    blue = blue - background_blue

    background_red = [np.zeros(n_r)]
    for i in np.arange(0,n_back-1,1):
        arr_i = red_spectra.fluxCoadd[np.arange(0,n_r,1), np.ones(n_r).astype(int)*indices[0][i]+1, np.ones(n_r).astype(int)*indices[1][i]+1]
        background_red = np.append(background_red,[arr_i], axis = 0)
    background_red = np.mean(background_red[1:], axis = 0)
    red = red - background_red

    blue[blue < 0] = 0
    red[red < 0]=0 
    #red and blue are the arrays of fluxes for the source pixels
    
    if PlotFlag == True:
        plt.plot(blue_spectra.wavelength, blue)
        plt.plot(red_spectra.wavelength, red,  color = 'r')
        plt.title('Background corrected spectrum of each side'+str(j))
        plt.xlabel('Wavelength ($\AA$)')
        plt.savefig(outName+'_background_corrected'+str(j), dpi = 600)
        plt.clf()
    
    #Now we stitch the red and blue sides together:
    overlap_blue = blue_spectra.wavelength[blue_spectra.wavelength > red_spectra.wavelength[0] ]
    overlap_red = red_spectra.wavelength[red_spectra.wavelength < blue_spectra.wavelength[-1] ]
    #The difference in the average of twice the overlap length on the long wavelength side of blue and short side of red
    difference = np.mean(red[len(overlap_red):2*len(overlap_red)]) - np.mean(blue[-2*len(overlap_blue):-len(overlap_blue)])
    red_shifted = red-difference
    #We then average the data in the overlap region
    interpolated_red_overlap = np.interp(blue_spectra.wavelength[-len(overlap_blue):],
                                     red_spectra.wavelength[:len(overlap_red)],red_shifted[:len(overlap_red)])
    overlap_average = np.mean([blue[-len(overlap_blue):],interpolated_red_overlap],axis =0)
    final_spectrum = np.append(np.append(blue[:-len(overlap_blue)],overlap_average),red_shifted[len(overlap_red):])
    final_wavelengths = np.append(blue_spectra.wavelength,red_spectra.wavelength[len(overlap_red):])
    
    ##Now we repeat for the source variance
    if n_source ==1:
        blue_variance = blue_spectra.varianceCoadd[np.arange(0,n_b,1), np.ones(n_b).astype(int)*(ind[0]+1), np.ones(n_b).astype(int)*(ind[1]+1)]
        red_variance = red_spectra.varianceCoadd[np.arange(0,n_r,1), np.ones(n_r).astype(int)*(ind[0]+1), np.ones(n_r).astype(int)*(ind[1]+1)]

    if n_source >1:
        blue_variance = [np.zeros(n_b)]
        for i in np.arange(0,n_source-1,1):
            arr_i = blue_spectra.varianceCoadd[np.arange(0,n_b,1), np.ones(n_b).astype(int)*indices_source[0][i]+1, np.ones(n_b).astype(int)*indices_source[1][i]+1]
            blue_variance = np.append(blue_variance,[arr_i], axis = 0)
        blue_variance = np.sqrt(np.sum((blue_variance[1:])**2, axis = 0))/np.sqrt(n_source)
        red_variance = [np.zeros(n_r)]
        for i in np.arange(0,n_source-1,1):
            arr_i = red_spectra.varianceCoadd[np.arange(0,n_r,1), np.ones(n_r).astype(int)*indices_source[0][i]+1, np.ones(n_r).astype(int)*indices_source[1][i]+1]
            red_variance = np.append(red_variance,[arr_i], axis = 0)
        red_variance = np.sqrt(np.sum((red_variance[1:])**2, axis = 0))/np.sqrt(n_source)

#background variance
    background_blue_variance = [np.zeros(n_b)]
    for i in np.arange(0,n_back-1,1):
        arr_i = blue_spectra.varianceCoadd[np.arange(0,n_b,1), np.ones(n_b).astype(int)*indices[0][i]+1, np.ones(n_b).astype(int)*indices[1][i]+1]
        background_blue_variance = np.append(background_blue_variance,[arr_i], axis = 0)
    background_blue_variance = np.sqrt(np.sum((background_blue_variance[1:])**2, axis = 0))/np.sqrt(n_back)

##I use standard error in the mean, which assumes samples come from same background population - this seems reasonable given would
## expect background spectrum to be uniform on the camera, except for any issues in the pixels themselves (which would probably 
##come up as a nonzero value in "badpix")

    blue_variance = np.sqrt(blue_variance**2 + background_blue_variance**2)
    background_red_variance = [np.zeros(n_r)]
    for i in np.arange(0,n_back-1,1):
        arr_i = red_spectra.varianceCoadd[np.arange(0,n_r,1), np.ones(n_r).astype(int)*indices[0][i]+1, np.ones(n_r).astype(int)*indices[1][i]+1]
        background_red_variance = np.append(background_red_variance,[arr_i], axis = 0)
    background_red_variance = np.sqrt(np.sum((background_red_variance[1:])**2, axis = 0))/np.sqrt(n_back)
    red_variance = np.sqrt(red_variance**2 + background_red_variance**2)
    blue_variance[blue_variance < 0] = 0
    red_variance[red_variance < 0]=0

    overlap_interpolated_variance = np.interp(blue_spectra.wavelength[-len(overlap_blue):],
                                     red_spectra.wavelength[:len(overlap_red)],red_variance[:len(overlap_red)])+blue_variance[-len(overlap_blue):]
    final_spectrum_variance = np.append(np.append(blue_variance[:-len(overlap_blue)],overlap_interpolated_variance),red_variance[len(overlap_red):])
    final_wavelengths_variance = np.append(blue_spectra.wavelength,red_spectra.wavelength[len(overlap_red):])
    
    ##Now we look at the bad pixels. Since WiFes has badpix for every wavelength, we average the badpix value (0 or 1) over all wavelengths
    #for each spaxel. If the spaxel has an average badpix value of above 0.1 it counts as a bad pixel
    blue_badpix = np.mean(np.mean(blue_spectra.badpixCoadd, axis = 1), axis = 1)
    red_badpix = np.mean(np.mean(red_spectra.badpixCoadd, axis = 1), axis = 1)
    spectrum_badpix = np.append(blue_badpix,red_badpix[len(overlap_red):])

    def filter_function(x):
        if x > 0.1:
            return 1
        else :
            return 0
    
    final_spectrum_badpix = list(map(lambda x: filter_function(x), spectrum_badpix))
    if PlotFlag ==True:
        
        fig, axs = plt.subplots(3, 1, sharex=True, figsize = (6,5))
        fig.subplots_adjust(hspace=0)
        fig.suptitle('Final stitched spectrum '+str(j))
        axs[0].plot(final_wavelengths, final_spectrum)
        axs[0].set_ylabel('Flux')

        axs[1].plot(final_wavelengths, final_spectrum_variance)
        axs[1].set_ylabel('Variance')
        axs[2].plot(final_wavelengths, final_spectrum_badpix)
        axs[2].set_xlabel('Wavelength ($\AA$)')
        axs[2].set_ylabel('Badpix')        
        
        fig.savefig(outName+'_spectrumvariancebadpix'+str(j), dpi = 600)
        plt.clf()
    
    #And now we make the data into a fits file, for input into the calibSpec pipeline.
    
    outputName = outName+str(j)+"_stitched.fits"
    print("Saving Data to " + outputName)
    hdulist = fits.HDUList(fits.PrimaryHDU())
    header = fits.Header()
    header['SOURCE'] = outputName
    header['RA'] = blue_spectra.RA
    header['DEC'] = blue_spectra.DEC
    header['CTYPE1'] = 'wavelength'
    header['CUNIT1'] = 'angstrom'
    header['DATE-OBS'] = blue_spectra.dates[0]
    header['EXPTIME'] = blue_spectra.exposed[0]
    hdulist[0].data = final_wavelengths
    hdulist.append(fits.ImageHDU(data=final_spectrum, header=header))
    hdulist.append(fits.ImageHDU(data=final_spectrum_variance, header=header))
    hdulist.append(fits.ImageHDU(data=final_spectrum_badpix, header=header))
    hdulist.writeto(outputName, overwrite=True)
    hdulist.close()

    #If desired we can put all the observations of the source together in one fits file:
if CombineFlag == True:
    hdulist = fits.HDUList(fits.PrimaryHDU())
    image_file = outName+str(0)+"_stitched.fits"
    hdul = fits.open(image_file)
    wavelengths = hdul[0]
    flux = hdul[1]
    variance = hdul[2]
    badpix = hdul[3]
    hdulist[0].data=wavelengths.data
    hdulist.append(fits.ImageHDU(data=flux.data, header=flux.header))
    hdulist.append(fits.ImageHDU(data=variance.data, header=variance.header))
    hdulist.append(fits.ImageHDU(data=badpix.data, header=badpix.header))   
    for j in np.arange(len(sourceNamesblue)):
        if j ==0:
            print(':)')
        else:
            image_file = outName+str(j)+"_stitched.fits"
            hdul = fits.open(image_file)
            wavelengths1 = hdul[0]
            flux1 = hdul[1]
            variance1 = hdul[2]
            badpix1 = hdul[3]
            if len(wavelengths1.data)==len(wavelengths.data): 
                hdulist.append(fits.ImageHDU(data=flux1.data, header=flux1.header))
                hdulist.append(fits.ImageHDU(data=variance1.data, header=variance1.header))
                hdulist.append(fits.ImageHDU(data=badpix1.data, header=badpix1.header))   
                
                ##Sometimes the observations have slightly different wavelength arrays. We choose to use the first observation
                ##wavelength since the difference will be only one or so wavelengths - negligible. If it is longer we remove from 
                ##the red side, if it is shorter we pad with the mean.
            elif len(wavelengths1.data) < len(wavelengths.data):
                flux2 = flux1.data
                variance2 = variance1.data
                badpix2 = badpix1.data
                flux2 = np.pad(flux2, (0, len(wavelengths.data)-len(wavelengths1.data)), 'mean')
                variance2 = np.pad(variance2, (0, len(wavelengths.data)-len(wavelengths1.data)), 'mean')
                badpix2 = np.pad(badpix2, (0, len(wavelengths.data)-len(wavelengths1.data)), 'constant',constant_values=0)
                hdulist.append(fits.ImageHDU(data=flux2, header=flux1.header))
                hdulist.append(fits.ImageHDU(data=variance2, header=variance1.header))
                hdulist.append(fits.ImageHDU(data=badpix2, header=badpix1.header))  
            elif len(wavelengths1.data) > len(wavelengths.data):
                flux2 = flux1.data
                variance2 = variance1.data
                badpix2 = badpix1.data
                flux2 = flux2[:len(wavelengths.data)]
                variance2 = variance2[:len(wavelengths.data)]
                badpix2 = badpix2[:len(wavelengths.data)]
                hdulist.append(fits.ImageHDU(data=flux2, header=flux1.header))
                hdulist.append(fits.ImageHDU(data=variance2, header=variance1.header))
                hdulist.append(fits.ImageHDU(data=badpix2, header=badpix1.header))  
                
    outputname = outName+"_stitched_combined.fits"
    print('Writing combined output to '+outputname)
    hdulist.writeto(outputname, overwrite=True)
    hdulist.close()<jupyter_output>Saving Data to AGNData/IRAS09149-6206-0_stitched.fits
date:
[58999.427]
Saving Data to AGNData/IRAS09149-6206-1_stitched.fits
date:
[59000.438]
:)
Writing combined output to AGNData/IRAS09149-6206-_stitched_combined.fits
{'fraction_non_alphanumeric': 0.12338989981598855, 'fraction_numerical': 0.02392148844817011, 'mean_word_length': 3.176093916755603, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '99649', 'license': 'Unlicense', 'n_tokens_mistral': 6998, 'n_tokens_neox': 6622, 'n_words': 1423, 'path': 'stitching pipeline.ipynb', 'repo_name': 'claireyung/WiFeS-RM', 'size': '19576'}