{"action":"opened","issue":{"url":"https://api.github.com/repos/Keck-DataReductionPipelines/MosfireDRP/issues/47","repository_url":"https://api.github.com/repos/Keck-DataReductionPipelines/MosfireDRP","labels_url":"https://api.github.com/repos/Keck-DataReductionPipelines/MosfireDRP/issues/47/labels{/name}","comments_url":"https://api.github.com/repos/Keck-DataReductionPipelines/MosfireDRP/issues/47/comments","events_url":"https://api.github.com/repos/Keck-DataReductionPipelines/MosfireDRP/issues/47/events","html_url":"https://github.com/Keck-DataReductionPipelines/MosfireDRP/issues/47","id":141748282,"number":47,"title":"Applying barycentric corrections as a part of the DRP?","user":{"login":"themiyan","id":5274688,"avatar_url":"https://avatars.githubusercontent.com/u/5274688?v=3","gravatar_id":"","url":"https://api.github.com/users/themiyan","html_url":"https://github.com/themiyan","followers_url":"https://api.github.com/users/themiyan/followers","following_url":"https://api.github.com/users/themiyan/following{/other_user}","gists_url":"https://api.github.com/users/themiyan/gists{/gist_id}","starred_url":"https://api.github.com/users/themiyan/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/themiyan/subscriptions","organizations_url":"https://api.github.com/users/themiyan/orgs","repos_url":"https://api.github.com/users/themiyan/repos","events_url":"https://api.github.com/users/themiyan/events{/privacy}","received_events_url":"https://api.github.com/users/themiyan/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":0,"created_at":"2016-03-18T01:25:23Z","updated_at":"2016-03-18T01:25:23Z","closed_at":null,"body":"Hi,\r\n\r\nFollowing 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. \r\nThis meant that for some of our galaxies, the emission line centroid could be offset by 1-2 pixels. \r\nIt took us a while to figure out why this happened, but once we did it was a quick fix. \r\n\r\nOne 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. \r\n\r\nSo 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. \r\n \r\nI've written a function within Wavelength.py to perform these corrections shown at the end. \r\nThis requires some minor changes to other functions as as well: https://github.com/themiyan/MosfireDRP_Themiyan/tree/master/MOSFIRE\r\n\r\nAnyways, just wanted to check if anyone would notice any obvious flaws in doing this using this method?\r\n\r\n```python\r\ndef bary_corr(files,wavename,options,maskname,extension=None):\r\n    \"\"\"This function will perfom barycentric corrections to the data. \r\n    Uses PyAstronomy module which may not be a part of the Ureka installation. Pandas can be replaced with\r\n    Astropy ascii or numpy if required.\r\n    \r\n    Will calculate the median barycentric velocity for all the exposures in files and apply the \r\n    velocity correction to the wavename file.\r\n    \r\n    Args: \r\n        files: list of strings with file names\r\n        wavename: path (relative or full) to the wavelength stack file, string\r\n        options: passed across from Options file\r\n        maskname: string with name of the mask\r\n        extension: path to file that contains a well formated fits header\r\n            this should be used only when the detector server fails\r\n            to write the full FITS header\r\n    Returns:\r\n        1. modifies the existing wavename file to wavename_before_correction.fits\r\n        2. saves the barycentric corrected wavelength solution as wavename\r\n        3. writes a csv file with barycentric velocities for the exposures in files. \r\n         \r\n    \"\"\"\r\n    try:\r\n        from PyAstronomy import pyasl\r\n    except ImportError:\r\n        import pip\r\n        pip.main(['install', 'PyAstronomy'])\r\n        from PyAstronomy import pyasl\r\n    import pandas as pd\r\n\r\n    lam = IO.readfits(wavename, options)\r\n\r\n    try:\r\n        assert lam[0]['BARYCORR']==None, 'The lambda solutions file has already been corrected for barycentric velocity. Comment this step from the Driver file.'\r\n    except KeyError:\r\n        pass\r\n    \r\n    files = IO.list_file_to_strings(files)\r\n        \r\n    output = pd.DataFrame( index=np.arange(1,len(files),1), columns=['exp_name'])\r\n    \r\n    # Coordinates of KECK 1\r\n    # https://www.ifa.hawaii.edu/mko/coordinates.shtml\r\n    longitude = 155 + (28./60) + (28.98665/3600)\r\n    latitude = 19+ (49./60)+  (33.40757/3600)\r\n    altitude = 4159.581216 \r\n    \r\n    \r\n    for i in xrange(len(files)):\r\n        fname = files[i]\r\n        thishdr, data, bs = IO.readmosfits(fname, options, extension=extension)\r\n        \r\n        try:\r\n            Julian_date = thishdr['MJD-OBS'] + 2.4e6 #this is some magical number\r\n        except KeyError:\r\n            print \"Header don't have information\"\r\n            continue\r\n            \r\n        RA  = thishdr['TARGRA']\r\n        DEC = thishdr['TARGDEC']\r\n        \r\n        #corrections are calculated \r\n        corr, hjd = pyasl.helcorr(longitude, latitude, altitude, RA, DEC, Julian_date, debug=False)\r\n\r\n        print fname, \"--> Barycentric correction [km/s]: \", corr, \" Heliocentric Julian day: \", hjd\r\n        \r\n        output.ix[i+1,'index']            = i+1\r\n        output.ix[i+1,'exp_name']         = fname\r\n        output.ix[i+1,'date']             = thishdr['DATE-OBS']\r\n        output.ix[i+1,'UTC']              = thishdr['UTC']\r\n        output.ix[i+1,'julian_date']      = Julian_date\r\n        output.ix[i+1,'heliocentric_julian_date']      = hjd\r\n        output.ix[i+1,'Barycentric_correction_kms-1']  = corr\r\n        output.ix[i+1,'RA']   = RA\r\n        output.ix[i+1,'DEC']  = DEC\r\n        \r\n        \r\n    output.set_index('index', inplace=True)     \r\n    output.to_csv('barycentric_corrections.csv')#fix path\r\n    \r\n    median_bary_corr = np.median(output['Barycentric_correction_kms-1'])\r\n    \r\n    #apply correction\r\n    Lambda_corrected = np.asarray(lam[1] * (1.0+(float(median_bary_corr)/3e5)))\r\n    wavename = wavename.rstrip(\".fits\")\r\n    lam[0]['BARYCORR'] = (median_bary_corr,'Km/s,+->moving towards object')\r\n    IO.writefits(Lambda_corrected, maskname, wavename, \r\n            options, header=lam[0], overwrite=True, rename=True)\r\n                    \r\n    return\r\n```"}}