Applying barycentric corrections as a part of the DRP?

[themiyan]

Hi,

Following observations of same galaxies in different years/months we ran into an issue where the redshifts were off by ~50 km/s. We were gladly ignoring heliocentric/barycentric corrections to our data, forgetting the improved resolution of MOSFIRE + that values can range from positive to negative throughout the year. This meant that for some of our galaxies, the emission line centroid could be offset by 1-2 pixels. It took us a while to figure out why this happened, but once we did it was a quick fix.

One of the challenges is that, if you require to add data in 2D between different months, this offsets can produce additional smearing of emission lines, and apart from everything else would be disastrous for kinematic studies. The solution is to interpolate the 2D spectra to a common grid, which means that data will be interpolated at least twice (the first being during rectification). From what I understood from @csteidel in https://github.com/Keck-DataReductionPipelines/MosfireDRP/issues/42, interpolation itself might be problematic and could lead to differences between different algorithms.

So the best solution would be to apply these corrections at the end of the wavelength solution? The lambda_solution_wave_stack has the wavelength solution per pixel which will later be interpolated in to a common grid depending on the filter during the rectification step.

I've written a function within Wavelength.py to perform these corrections shown at the end. This requires some minor changes to other functions as as well: https://github.com/themiyan/MosfireDRP_Themiyan/tree/master/MOSFIRE

Anyways, just wanted to check if anyone would notice any obvious flaws in doing this using this method?

def bary_corr(files,wavename,options,maskname,extension=None):
    """This function will perfom barycentric corrections to the data. 
    Uses PyAstronomy module which may not be a part of the Ureka installation. Pandas can be replaced with
    Astropy ascii or numpy if required.

    Will calculate the median barycentric velocity for all the exposures in files and apply the 
    velocity correction to the wavename file.

    Args: 
        files: list of strings with file names
        wavename: path (relative or full) to the wavelength stack file, string
        options: passed across from Options file
        maskname: string with name of the mask
        extension: path to file that contains a well formated fits header
            this should be used only when the detector server fails
            to write the full FITS header
    Returns:
        1. modifies the existing wavename file to wavename_before_correction.fits
        2. saves the barycentric corrected wavelength solution as wavename
        3. writes a csv file with barycentric velocities for the exposures in files. 

    """
    try:
        from PyAstronomy import pyasl
    except ImportError:
        import pip
        pip.main(['install', 'PyAstronomy'])
        from PyAstronomy import pyasl
    import pandas as pd

    lam = IO.readfits(wavename, options)

    try:
        assert lam[0]['BARYCORR']==None, 'The lambda solutions file has already been corrected for barycentric velocity. Comment this step from the Driver file.'
    except KeyError:
        pass

    files = IO.list_file_to_strings(files)

    output = pd.DataFrame( index=np.arange(1,len(files),1), columns=['exp_name'])

    # Coordinates of KECK 1
    # https://www.ifa.hawaii.edu/mko/coordinates.shtml
    longitude = 155 + (28./60) + (28.98665/3600)
    latitude = 19+ (49./60)+  (33.40757/3600)
    altitude = 4159.581216 

    for i in xrange(len(files)):
        fname = files[i]
        thishdr, data, bs = IO.readmosfits(fname, options, extension=extension)

        try:
            Julian_date = thishdr['MJD-OBS'] + 2.4e6 #this is some magical number
        except KeyError:
            print "Header don't have information"
            continue

        RA  = thishdr['TARGRA']
        DEC = thishdr['TARGDEC']

        #corrections are calculated 
        corr, hjd = pyasl.helcorr(longitude, latitude, altitude, RA, DEC, Julian_date, debug=False)

        print fname, "--> Barycentric correction [km/s]: ", corr, " Heliocentric Julian day: ", hjd

        output.ix[i+1,'index']            = i+1
        output.ix[i+1,'exp_name']         = fname
        output.ix[i+1,'date']             = thishdr['DATE-OBS']
        output.ix[i+1,'UTC']              = thishdr['UTC']
        output.ix[i+1,'julian_date']      = Julian_date
        output.ix[i+1,'heliocentric_julian_date']      = hjd
        output.ix[i+1,'Barycentric_correction_kms-1']  = corr
        output.ix[i+1,'RA']   = RA
        output.ix[i+1,'DEC']  = DEC

    output.set_index('index', inplace=True)     
    output.to_csv('barycentric_corrections.csv')#fix path

    median_bary_corr = np.median(output['Barycentric_correction_kms-1'])

    #apply correction
    Lambda_corrected = np.asarray(lam[1] * (1.0+(float(median_bary_corr)/3e5)))
    wavename = wavename.rstrip(".fits")
    lam[0]['BARYCORR'] = (median_bary_corr,'Km/s,+->moving towards object')
    IO.writefits(Lambda_corrected, maskname, wavename, 
            options, header=lam[0], overwrite=True, rename=True)

    return

[lucarizzi]

I believe we could offer this as a standalone method, but probably it won't be implemented in the master release. Thank you for your contribution! I will try to document it on the instruction manual. Could you tell me at what point you think it should be called during the wavelength calculation?

[themiyan]

Sure, thanks. It should be added before background subtraction to avoid multiple interpolation of the data. eg: I've changed AutoDriver to generate something like this for my Driver file.

Wavelength.apply_lambda_simple(maskname, band, obsfiles, waveops)
Wavelength.apply_lambda_sky_and_arc(maskname, band, obsfiles,  'Ne.txt', LROIs, waveops)

Wavelength_file = 'merged_lambda_solution_wave_stack_K_m131224_0354-0372_and_wave_stack_K_m131224_0135-0137.fits'
Wavelength.bary_corr(obsfiles, Wavelength_file, maskname, band, waveops)

Background.handle_background(obsfiles,Wavelength_file,maskname,band,waveops)

[lucarizzi]

Thank you!

On May 3, 2016, at 9:13 PM, Themiya Nanayakkara notifications@github.com wrote:

Sure, thanks. It should be added before background subtraction to avoid multiple interpolation of the data. eg: I've changed AutoDriver to generate something like this for my Driver file.

Wavelength.apply_lambda_simple(maskname, band, obsfiles, waveops) Wavelength.apply_lambda_sky_and_arc(maskname, band, obsfiles, 'Ne.txt', LROIs, waveops)

Wavelength_file = 'merged_lambda_solution_wave_stack_K_m131224_0354-0372_and_wave_stack_K_m131224_0135-0137.fits' Wavelength.bary_corr(obsfiles, Wavelength_file, maskname, band, waveops)

Background.handle_background(obsfiles,Wavelength_file,maskname,band,waveops) — You are receiving this because you commented. Reply to this email directly or view it on GitHub https://github.com/Keck-DataReductionPipelines/MosfireDRP/issues/47#issuecomment-216763283