Splitchips

.gitignore

@@ -7,7 +7,6 @@ __pycache__
 a.out
 .cache
 .eggs
-*.fits
 *.cat
 
 # Ignore .c files by default to avoid including generated code. If you want to

lcogtgemini/__init__.py

@@ -1,9 +1,9 @@
 #!/usr/bin/env python
-'''
+"""
 Created on Nov 7, 2014
 
 @author: cmccully
-'''
+"""
 
 import os
 from pyraf import iraf
@@ -16,6 +16,7 @@
 iraf.onedspec()
 
 bluecut = 3450
+redcut = 10000
 
 iraf.gmos.logfile = "log.txt"
 iraf.gmos.mode = 'h'

lcogtgemini/bias.py

@@ -3,6 +3,7 @@
 from astropy.io import fits
 from pyraf import iraf
 import numpy as np
+import os
 
 
 def makebias(fs, obstypes, rawpath):
@@ -24,4 +25,4 @@ def makebias(fs, obstypes, rawpath):
             biastxtfile.close()
             iraf.gbias('@%s/bias{binning}.txt'.format(binning=binning) % os.getcwd(),
                        'bias{binning}'.format(binning=binning), rawpath=rawpath, fl_over=lcogtgemini.dooverscan,
-                       fl_vardq=lcogtgemini.dodq)
+                       fl_vardq=lcogtgemini.dodq)

lcogtgemini/bpm.py

@@ -2,14 +2,15 @@
 from pyraf import iraf
 from astropy.io import fits
 import numpy as np
+import pkg_resources
 
 
 def get_bad_pixel_mask(binnings, yroi):
     if lcogtgemini.detector == 'Hamamatsu':
         if lcogtgemini.is_GS:
-            bpm_file = 'bpm_gs.fits'
+            bpm_file = pkg_resources.resource_filename('lcogtgemini', 'bpm_gs.fits')
         else:
-            bpm_file = 'bpm_gn.fits'
+            bpm_file = pkg_resources.resource_filename('lcogtgemini', 'bpm_gn.fits')
 
         bpm_hdu = fits.open(bpm_file, uint16=True)
 
@@ -20,7 +21,7 @@ def get_bad_pixel_mask(binnings, yroi):
                 binning_list = binning.split('x')
                 binx, biny = int(binning_list[0]), int(binning_list[1])
                 iraf.unlearn('blkavg')
-                binned_bpm_filename = 'bpm.{i}.{x}x{y}.fits'.format(i=i,x=binx, y=biny)
+                binned_bpm_filename = 'bpm.{i}.{x}x{y}.fits'.format(i=i, x=binx, y=biny)
                 iraf.blkavg('bpm.{i}.unbinned.fits[1]'.format(i=i),
                             binned_bpm_filename, binx, biny)
                 averaged_bpm = fits.open(binned_bpm_filename)
@@ -39,4 +40,3 @@ def get_bad_pixel_mask(binnings, yroi):
                 bpm_data = np.zeros((2048 // int(biny), 1080 // int(binx)), dtype=np.uint8)
                 binned_bpm_filename = 'bpm.{i}.{x}x{y}.fits'.format(i=i, x=binx, y=biny)
                 fits.writeto(binned_bpm_filename, bpm_data)
-

lcogtgemini/combine.py

@@ -6,12 +6,11 @@
 from astropy.convolution import convolve, Gaussian1DKernel
 from lcogtgemini import utils
 
-from pyraf import iraf
 
 def find_bad_pixels(data, threshold=30.0):
     # Take the abs next pixel diff
     absdiff = np.abs(data[1:] - data[:-1])
-    scaled_diff =  absdiff / convolve(data, Gaussian1DKernel(stddev=20.0))[1:]
+    scaled_diff = absdiff / convolve(data, Gaussian1DKernel(stddev=20.0))[1:]
     # divide by square root of 2
     scatter = scaled_diff / np.sqrt(2.0)
     # Take the median
@@ -47,58 +46,51 @@ def speccombine(fs, outfile):
         wavelength_step = min([wavelength_step, wavelength[1] - wavelength[0]])
 
     min_w = np.max([min_w, lcogtgemini.bluecut])
+    max_w = np.min([max_w, lcogtgemini.redcut])
     first_hdu = fits.open(fs[0])
     first_wavelengths = fits_utils.fitshdr_to_wave(first_hdu['SCI'].header)
-    bad_pixels = find_bad_pixels(first_hdu['SCI'].data[0])
-    first_hdu['SCI'].data[0][bad_pixels] = 0.0
+    first_fluxes = first_hdu['SCI'].data[0, 0]
+    bad_pixels = find_bad_pixels(first_fluxes)
+    first_fluxes[bad_pixels] = 0.0
 
-    wavelength_grid = np.arange(min_w, max_w + wavelength_step, wavelength_step)
+    wavelength_grid = np.arange(min_w, max_w, wavelength_step)
     data_to_combine = np.zeros((len(fs), wavelength_grid.shape[0]))
     for i, f in enumerate(fs):
         hdu = fits.open(f)
         wavelengths = fits_utils.fitshdr_to_wave(hdu['SCI'].header)
+        fluxes = hdu['SCI'].data[0, 0]
         in_chips = np.zeros(wavelengths.shape, dtype=bool)
         basename = file_utils.get_base_name(f)
-        wavelengths_hdu = fits.open(basename +'.wavelengths.fits')
+        wavelengths_hdu = fits.open(basename + '.wavelengths.fits')
         chips = utils.get_wavelengths_of_chips(wavelengths_hdu)
         for chip in chips:
             in_chip = np.logical_and(wavelengths >= min(chip), wavelengths <= max(chip))
             in_chips[in_chip] = True
 
-        hdu['SCI'].data[0][~in_chips] = 0.0
-
+        fluxes[~in_chips] = 0.0
         overlap = np.logical_and(wavelengths >= overlap_min_w, wavelengths <= overlap_max_w)
 
         # Reject outliers
-        bad_pixels = find_bad_pixels(hdu['SCI'].data[0])
+        bad_pixels = find_bad_pixels(fluxes)
         in_telluric = np.logical_and(wavelengths >= 6640.0, wavelengths <= 7040.0)
         in_telluric = np.logical_or(in_telluric, np.logical_and(wavelengths >= 7550.0, wavelengths <= 7750.0))
         bad_pixels[in_telluric] = False
-        # Take the median of the ratio of each spectrum to the first to get the rescaling
 
-        first_fluxes = np.interp(wavelengths[overlap], first_wavelengths, first_hdu['SCI'].data[0], left=0.0, right=0.0)
-        good_pixels = np.ones(hdu['SCI'].data[0].shape[0], dtype=bool)
-        good_pixels[overlap] = np.logical_and(hdu['SCI'].data[0][overlap] != 0.0, first_fluxes != 0.0)
-        good_pixels = np.logical_and(good_pixels, ~(bad_pixels))
+        first_fluxes_interp = np.interp(wavelengths[overlap], first_wavelengths, first_fluxes, left=0.0, right=0.0)
+        good_pixels = (fluxes[overlap] * first_fluxes_interp != 0.) & ~bad_pixels[overlap]
 
-        scale = np.median(first_fluxes[good_pixels[overlap]] / hdu['SCI'].data[0][overlap][good_pixels[overlap]])
+        # Take the median of the ratio of each spectrum to the first to get the rescaling
+        scale = np.nanmedian(first_fluxes_interp[good_pixels] / fluxes[overlap][good_pixels])
         scales.append(scale)
-        hdu['SCI'].data[0][hdu['SCI'].data[0] == 0.0] = -1.e9
-        hdu['SCI'].data[0][bad_pixels] = -1.e9
-        data_to_combine[i] = np.interp(wavelength_grid, wavelengths, hdu['SCI'].data[0], left=0.0, right=0.0)
-        data_to_combine[data_to_combine < 0.0] = 0.0
+        fluxes[fluxes == 0.] = np.nan
+        fluxes[bad_pixels] = np.nan
+        data_to_combine[i] = np.interp(wavelength_grid, wavelengths, fluxes, left=0.0, right=0.0)
         data_to_combine[i] *= scale
 
-
     # write the scales into a file
     ascii.write({'scale': scales}, 'scales.dat', names=['scale'], format='fast_no_header')
 
-    combined_data = data_to_combine.sum(axis=0)
-    weights = (data_to_combine > 0).sum(axis=0)
-    weights[weights == 0.0] = 1.0
-    combined_data /= weights
-
-    first_hdu[0].data = combined_data
+    first_hdu[0].data = np.nanmedian(data_to_combine, axis=0)
     first_hdu[0].header['CRPIX1'] = 1
     first_hdu[0].header['CRVAL1'] = min_w
     first_hdu[0].header['CD1_1'] = wavelength_step
@@ -111,5 +103,3 @@ def speccombine(fs, outfile):
     first_hdu[0].header['APNUM1'] = first_hdu['SCI'].header['APNUM1']
 
     first_hdu[0].writeto(outfile)
-    iraf.unlearn(iraf.scombine)
-

lcogtgemini/cosmicrays.py

@@ -10,26 +10,20 @@
 def crreject(scifiles):
     for f in scifiles:
         # run lacosmicx
-        hdu = fits.open('st' + f.replace('.txt', '.fits'))
-
-        readnoise = float(hdu[2].header['RDNOISE'])
-        # figure out what pssl should be approximately
+        hdu = fits.open('t' + f.replace('.txt', '.fits'))
         d = hdu[2].data.copy()
-        dsort = np.sort(d.ravel())
-        nd = dsort.shape[0]
-        # Calculate the difference between the 16th and 84th percentiles to be
-        # robust against outliers
-        dsig = (dsort[int(round(0.84 * nd))] - dsort[int(round(0.16 * nd))]) / 2.0
-        pssl = (dsig * dsig - readnoise * readnoise)
+        readnoise = float(hdu[2].header['RDNOISE'])
+
+        hdu_skysub = fits.open('st' + f.replace('.txt', '.fits'))
+        bkg = d - hdu_skysub[2].data
 
         mask = d == 0.0
         mask = np.logical_or(mask, d > (50000.0 * float(hdu[0].header['GAINMULT'])))
         if lcogtgemini.dodq:
             mask = np.logical_or(mask,  hdu['DQ'].data)
 
-        crmask, _cleanarr = detect_cosmics(d, inmask=mask, sigclip=5.0,
-                                           objlim=10.0, sigfrac=0.2, gain=1.0,
-                                           readnoise=readnoise, pssl=pssl)
+        crmask, _cleanarr = detect_cosmics(d, inmask=mask, bkg=bkg, sigclip=5.0,
+                                           objlim=10.0, sigfrac=0.2, gain=1.0, readnoise=readnoise)
 
-        tofits(f[:-4] + '.lamask.fits', np.array(crmask, dtype=np.uint8), hdr=hdu['SCI'].header.copy())
+        tofits(f[:-4] + '.lamask.fits', np.array(crmask, dtype=np.uint8), hdr=hdu['SCI'].header.copy(), clobber=True)
         fixpix.run_fixpix('t' + f[:-4] + '.fits[2]', f[:-4] + '.lamask.fits')

lcogtgemini/file_utils.py

@@ -4,6 +4,7 @@
 import os
 from glob import glob
 from pyraf import iraf
+import pkg_resources
 
 
 def getobstypes(fs):
@@ -45,7 +46,7 @@ def getsetupname(f, calfile=False):
     return setupname
 
 
-def gettxtfiles(fs, objname):
+def gettxtfiles(objname):
 
     flatfiles = np.array(glob('*.flat.txt'))
 
@@ -58,7 +59,8 @@ def gettxtfiles(fs, objname):
     nonscifiles = []
     # remove the arcs and flats
     for f in scifiles:
-        if 'arc' in f or 'flat' in f: nonscifiles.append(f)
+        if 'arc' in f or 'flat' in f:
+            nonscifiles.append(f)
 
     for f in nonscifiles:
         scifiles.remove(f)
@@ -69,21 +71,19 @@ def gettxtfiles(fs, objname):
 
 def get_base_name(f):
     objname = getobjname(np.array([f]), np.array(['OBJECT']))
-    # Drop the raw/
-    fname = f.split('/')[-1]
     # red or blue setting
     redblue = fits.getval(f, 'GRATING')[0].lower()
     # central wavelength
     lamcentral = fits.getval(f, 'CENTWAVE')
 
-    return  '%s.%s%i' % (objname, redblue, lamcentral)
+    return '%s.%s%i' % (objname, redblue, lamcentral)
 
 
 def maketxtfiles(fs, obstypes, obsclasses, objname):
     # go through each of the files (Ignore bias and aquisition files)
     goodfiles = np.logical_and(obsclasses != 'acqCal', obsclasses != 'acq')
     goodfiles = np.logical_and(goodfiles, obstypes != 'BIAS')
-    goodfiles = np.logical_and(goodfiles, obstypes!='BPM')
+    goodfiles = np.logical_and(goodfiles, obstypes != 'BPM')
     goodfiles = np.logical_and(goodfiles, obsclasses != 'sensitivity')
     correct_names = np.logical_or([os.path.basename(f)[0] == 'S' for f in fs],
                                   [os.path.basename(f)[0] == 'N' for f in fs])
@@ -136,29 +136,38 @@ def get_images_from_txt_file(filename):
     return lines
 
 
-def get_standard_file(objname, base_stddir):
-    if os.path.exists(objname+'.std.dat'):
-        standard_file = objname+'.std.dat'
+def get_standard_file(objname):
+    places_to_look = [objname + '.std.dat',
+                      os.path.join(iraf.osfn('onedstds$spec50cal/'), objname + '.dat'),
+                      os.path.join(iraf.osfn('onedstds$ctionewcal/'), objname + '.dat'),
+                      os.path.join(iraf.osfn('onedstds$oke90/'), objname + '.dat'),
+                      os.path.join(iraf.osfn('onedstds$iidscal/'), objname + '.dat'),
+                      os.path.join(iraf.osfn('gmos$calib/'), objname + '.dat')]
+    for standard_file in places_to_look:
+        if os.path.exists(standard_file):
+            return standard_file
     else:
-        standard_file = os.path.join(iraf.osfn('gmisc$lib/onedstds/'), base_stddir, objname + '.dat')
-    return standard_file
+        raise IOError('could not find standard file for ' + objname)
 
 
 def read_standard_file(filename, maskname):
     standard = ascii.read(filename)
     standard['col2'] = smooth(maskname, standard['col2'])
     return standard
 
+
 def read_telluric_model(maskname):
     # Read in the telluric model
-    telluric = ascii.read('telluric_model.dat')
+    telluric_model_filename = pkg_resources.resource_filename('lcogtgemini', 'telluric_model.dat')
+    telluric = ascii.read(telluric_model_filename)
     telluric['col2'] = smooth(maskname, telluric['col2'])
     return telluric
 
+
 def smooth(maskname, data):
     # Smooth the spectrum to the size of the slit
     # I measured 5 angstroms FWHM for a 1 arcsecond slit
     # the 2.355 converts to sigma
     smoothing_scale = 5.0 * float(maskname.split('arc')[0]) / 2.355
 
-    return convolve(data, Gaussian1DKernel(stddev=smoothing_scale))
+    return convolve(data, Gaussian1DKernel(stddev=smoothing_scale))

lcogtgemini/fits_utils.py

@@ -49,7 +49,7 @@ def fitshdr_to_wave(hdr):
     else:
         cdelt = float(hdr['CD1_1'])
     npix = float(hdr['NAXIS1'])
-    lam = np.arange(crval - cdelt * crpix ,
+    lam = np.arange(crval - cdelt * crpix,
                     crval + cdelt * (npix - crpix) - 1e-4,
                     cdelt)
     return lam
@@ -77,23 +77,23 @@ def cut_gs_image(filename, output_filename, pixel_range, namps):
     :param output_filename:
     :param pixel_range: array-like, The range of pixels to keep, python indexed,
                         given in binned pixels
+    :param namps: int, Number of amplifiers (same as lcogtgemini.namps)
     :return:
     """
     hdu = fits.open(filename, unit16=True)
     for i in range(1, namps + 1):
         ccdsum = hdu[i].header['CCDSUM']
         ccdsum = np.array(ccdsum.split(), dtype=np.int)
 
-        y_ccdsec = [(pixel_range[0]  * ccdsum[1]) + 1,
-                    (pixel_range[1]) * ccdsum[1]]
+        y_ccdsec = [(pixel_range[0] * ccdsum[1]) + 1,
+                    pixel_range[1] * ccdsum[1]]
 
         x_detector_section = get_x_pixel_range(hdu[i].header['DETSEC'])
         hdu[i].header['DETSEC'] = hdr_pixel_range(int(x_detector_section[0]), int(x_detector_section[1]), y_ccdsec[0], y_ccdsec[1])
 
         x_ccd_section = get_x_pixel_range(hdu[i].header['CCDSEC'])
         hdu[i].header['CCDSEC'] = hdr_pixel_range(int(x_ccd_section[0]), int(x_ccd_section[1]), y_ccdsec[0], y_ccdsec[1])
 
-
         numpix = pixel_range[1] - pixel_range[0]
 
         x_bias_section = get_x_pixel_range(hdu[i].header['BIASSEC'])
@@ -115,13 +115,14 @@ def updatecomheader(extractedfiles, filename):
         exptimes.append(float(fits.getval(f, 'EXPTIME')))
 
     fits.setval(filename, 'AIRMASS', value=np.mean(airmasses))
+    fits.setval(filename, 'EXPTIME', value=np.sum(exptimes))
     fits.setval(filename, 'SLIT', value=fits.getval(extractedfiles[0], 'MASKNAME').replace('arcsec', ''))
 
     comhdu = fits.open(filename, mode='update')
 
     extractedhdu = fits.open(extractedfiles[0])
     for k in extractedhdu[0].header.keys():
-        if not k in comhdu[0].header.keys():
+        if k not in comhdu[0].header.keys():
             extractedhdu[0].header.cards[k].verify('fix')
             comhdu[0].header.append(extractedhdu[0].header.cards[k])
 
@@ -131,6 +132,7 @@ def updatecomheader(extractedfiles, filename):
     dateobs = fits.getval(filename, 'DATE-OBS')
     dateobs += 'T' + fits.getval(filename, 'TIME-OBS')
     fits.setval(filename, 'DATE-OBS', value=dateobs)
+    fits.setval(filename, 'BUNIT', value='1E-15 erg / (s cm2 angstrom)')
 
 
 def spectoascii(infilename, outfilename):
@@ -145,4 +147,4 @@ def spectoascii(infilename, outfilename):
     d = np.zeros((2, len(lam)))
     d[0] = lam
     d[1] = flux
-    np.savetxt(outfilename, d.transpose())
+    np.savetxt(outfilename, d.transpose())