Python MOSFIRE.IO類代碼示例

def imcombine(filelist, maskname, fname, options, sum_type):
    """ combine the images in file list into fname.

    Sum type:
        rate -- filelist is in cnt/s
        ivar-rate -- filelist is in s/cnt
        snr-rate -- filelist is in SNR
    """

    ARR = None
    hdr = None
    i = 1

    itime = 0

    for file in filelist:

        this_hdr, img = IO.readfits(file)
        cards = this_hdr.ascardlist()

        thisitime = this_hdr["truitime"]
        itime += thisitime

        if ARR is None:
            ARR = np.zeros(img.shape)

        if sum_type == "rate":
            ARR += img * thisitime
        if sum_type == "ivar-rate":
            ARR += thisitime / img
        if sum_type == "snr-rate":
            ARR += img * thisitime

        if hdr is None:
            hdr = this_hdr
        hdr.update("fno%2.2i" % i, file, "--")
        for card in cards:
            key, value, comment = (card.key, card.value, card.comment)

            if hdr.has_key(key) and hdr[key] != value:
                key = key + ("_img%2.2i" % i)

            if len(key) > 8:
                key = "HIERARCH " + key

            try:
                hdr.update(key, value, comment)
            except ValueError:
                pass

    hdr.update("itime", itime, "Itime for %i rectified images" % len(filelist))
    if sum_type == "rate":
        ARR /= itime
    if sum_type == "ivar-rate":
        ARR = itime / ARR
    if sum_type == "snr-rate":
        ARR /= itime

    IO.writefits(ARR, maskname, fname, options, header=hdr, overwrite=True, lossy_compress=True)

def combine(flatlist, maskname, band, options):
    '''
    combine list of flats into a flat file'''

    out = os.path.join("combflat_2d_%s.fits" 
                    % (band))
    if os.path.exists(out):
            os.remove(out)

    IO.imcombine(flatlist, out, options, reject="minmax", nlow=1, nhigh=1)

    def __init__(self, maskname, bandname, options, fig):

        self.fig = fig
        self.flat = IO.read_drpfits(maskname, "combflat_2d_%s.fits" % bandname,
                options)

        self.edges, meta = IO.load_edges(maskname, bandname, options)

        self.edgeno=2
        self.cid = self.fig.canvas.mpl_connect('key_press_event', self)

        self.draw()

def rectify(dname, lamdat, A, B, maskname, band, wavoptions, 
        longoptions):

    header, data = IO.readfits(dname)
    raw_img = data * Detector.gain / header['TRUITIME']

    dlam = Wavelength.grating_results(band)
    hpp = np.array(Filters.hpp[band]) 
    ll_fid = np.arange(hpp[0], hpp[1], dlam)

    rectified = np.zeros((2048, len(ll_fid)))

    from scipy.interpolate import interp1d

    for i in xrange(2048):
        ll = lamdat[i,:]
        ss = raw_img[i,:]
        ok = np.isfinite(ll) & np.isfinite(ss) & (ll < hpp[1]) & (ll >
                hpp[0])

        if len(np.where(ok)[0]) < 30:
            continue

        f = interp1d(ll[ok], ss[ok], bounds_error=False)
        rectified[i,:] = f(ll_fid)

    header.update("wat0_001", "system=world")
    header.update("wat1_001", "wtype=linear")
    header.update("wat2_001", "wtype=linear")
    header.update("dispaxis", 1)
    header.update("dclog1", "Transform")
    header.update("dc-flag", 0)
    header.update("ctype1", "AWAV")
    header.update("cunit1", "Angstrom")
    header.update("crval1", ll_fid[0])
    header.update("crval2", 0)
    header.update("crpix1", 1)
    header.update("crpix2", 1)
    header.update("cdelt1", 1)
    header.update("cdelt2", 1)
    header.update("cname1", "angstrom")
    header.update("cname2", "pixel")
    header.update("cd1_1", dlam)
    header.update("cd1_2", 0)
    header.update("cd2_1", 0)
    header.update("cd2_2", 1)


    header.update("object", "rectified [eps]")
    IO.writefits(rectified, maskname, "rectified_%s" % (dname), 
        wavoptions, header=header, overwrite=True, lossy_compress=True)

 def isEmpty(self,file):
     if not os.path.exists(file):
         return True
     fname = IO.list_file_to_strings(file)
     if len(fname):
         return False
     else:
         return True

def combine(flatlist, maskname, band, options, lampsOff=False):
    '''
    combine list of flats into a flat file'''

    if lampsOff:
        out = os.path.join("combflat_lamps_off_2d_%s.fits" 
                    % (band))
    else:
        out = os.path.join("combflat_2d_%s.fits" 
                    % (band))
    if os.path.exists(out):
            os.remove(out)

    if len(flatlist)>1:
        IO.imcombine(flatlist, out, options, reject="minmax", nlow=1, nhigh=1)
    else:
        IO.imcombine(flatlist, out, options, reject="none", nlow=1, nhigh=1)

    def set_header(self, header, ssl=None, msl=None, asl=None, targs=None):
        '''Passed "header" a FITS header dictionary and converts to a Barset'''
        self.pos = np.array(IO.parse_header_for_bars(header))
        self.set_pos_pix()

        self.ssl = ssl
        self.msl = msl
        self.asl = asl
        self.targs = targs

        def is_alignment_slit(slit):
            return (np.float(slit["Target_Priority"]) < 0)

        # If len(ssl) == 0 then the header is for a long slit
        if (header['MASKNAME'] == 'long2pos'):
            info("long2pos mode in CSU slit determination")
            self.long2pos_slit = True

        if (len(ssl) == 0):
        
            self.long_slit = True

            start = np.int(msl[0]["Slit_Number"])
            stop = np.int(msl[-1]["Slit_Number"])


            for mech_slit in msl:
                mech_slit["Target_in_Slit"] = "long"

            self.ssl = np.array([("1", "??", "??", "??", "??", "??", "??", msl[0]['Slit_width'],
                (stop-start+1)*7.6, "0", "long", "0")],
                dtype= [ ('Slit_Number', '|S2'), 
                ('Slit_RA_Hours', '|S2'), ('Slit_RA_Minutes', '|S2'), ('Slit_RA_Seconds', '|S5'),
                ('Slit_Dec_Degrees', '|S3'), ('Slit_Dec_Minutes', '|S2'), ('Slit_Dec_Seconds', '|S5'), 
                ('Slit_width', '|S5'), ('Slit_length', '|S5'), ('Target_to_center_of_slit_distance', '|S5'), 
                ('Target_Name', '|S80'), ('Target_Priority', '|S1')])
            self.scislit_to_slit = [ np.arange(start,stop) ]
            ssl = None

        # Create a map between scislit number and mechanical slit
        # recall that slits count from 1
        if ssl is not None:
            prev = self.msl[0]["Target_in_Slit"]

            v = []

            for science_slit in ssl:
                targ = science_slit["Target_Name"]
                v.append([int(x) for x in self.msl.field("Slit_Number")[np.where(self.msl.field("Target_in_Slit").rstrip() == targ)[0]]])
            self.scislit_to_slit = v

            if (len(self.scislit_to_slit) != len(ssl)) and not (self.long_slit
                    and len(self.scislit_to_slit) == 1):
                error("SSL should match targets in slit")
                raise Exception("SSL should match targets in slit")

    def test_trace_edge(self):
            (header, data1, targs, ssl, msl, asl) = \
                            IO.readfits_all("/users/npk/desktop/c9/m110326_3242.fits")
            data = data1

            ssl = ssl[ssl["Slit_Number"] != ' ']
            numslits = np.round(np.array(ssl["Slit_length"], 
                    dtype=np.float) / 7.02)

            for i in range(len(ssl)):
                    print ssl[i]["Target_Name"], numslits[i]

def audit(filename):
    
    header, data = IO.readfits(filename)


    ll0 = header['crval1']
    dlam = header['cd1_1']

    ls = ll0 + dlam * np.arange(data.shape[1])

    linelist = Wavelength.pick_linelist(header)


    deltas = []
    sigs = []
    xpos = []
    ys = []
    for y in np.linspace(750, 1100, 30):
    #for y in np.linspace(5, 640, 50):
        sp = np.ma.array(data[y,:])

        xs, sxs, sigmas = Wavelength.find_known_lines(linelist, ls, sp,
                Options.wavelength)

        xpos.append(xs)
        ys.append([y] * len(xs))

        deltas.append(xs - (linelist - ll0)/dlam)
        sigs.append(sxs)


    xpos, ys, deltas, sigs = map(np.array, [xpos, ys, deltas, sigs])

    deltas[np.abs(deltas) > .75] = np.nan
    sigs[np.abs(sigs) > .001] = np.nan

    pl.clf()
    size = 0.003/sigs 
    size[size > 30] = 30
    size[size < 1] = 1
    pl.scatter( xpos, ys, c=deltas, s=size)
    pl.xlim([0, data.shape[1]])
    pl.ylim([0, data.shape[0]])
    pl.xlabel("Spectral pixel")
    pl.ylabel("Spatial pixel")
    pl.title("Night sky line deviation from solution [pixel]")
    pl.colorbar()

    pl.savefig("audit.pdf")
    
    pdb.set_trace()

def apply_flat(scifilename, maskname, band):
    ''' Divides the contents of scifilename by the flat field and
        overwrites scifilename with the same file divided by the flat

        Args:
            scifilename: Path to science file name.
            maskname: The mask name
            band: The filter bands

        Results:
            Overwrites scifilename where the data contents of the file
                are divided by the pixel flat
    '''

    
    flat = IO.load_flat(maskname, band, {})
    flat_data = flat[1].filled(1.0)

    header, data = IO.readfits(scifilename)
    
    print("Applying flat to file {0}".format(scifilename))
    IO.writefits(data/flat_data, maskname, scifilename, {}, header=header,
        overwrite=True)

    def printMaskAndBand(self):
        offsetfile = self.offsetFiles[0]
        fname = IO.list_file_to_strings(offsetfile)

        if os.path.isabs(fname[0]): path = fname[0]
        else: path = os.path.join(fname_to_path(fname[0]), fname[0])
        hdulist = pf.open(path)
        header = hdulist[0].header

        self.maskName = header['maskname']
        self.band = header['filter']
        self.addLine("maskname = '"+str(self.maskName)+"'")
        self.addLine("band = '"+str(self.band)+"'")
        self.addLine("")

def rename_files(wavenames, maskname, band, wavops):

    lamname = Wavelength.filelist_to_wavename(wavenames[0], band, maskname, wavops).rstrip(".fits")

    suffix = lamname.lstrip("wave_stack_%s_" % band)

    path = os.path.join(wavops["outdir"], maskname)

    fnames = ["rectified_%s%s.fits", "rectified_ivar_%s%s.fits", "rectified_sn_%s%s.fits"]

    for fname in fnames:
        try:
            a = get_path(os.path.join(path, fname % (band, "_" + suffix)))
            b = os.path.join(path, fname % (band, "")) + gz(a)
            os.rename(a, b)
        except:
            print "Ignoring renaming of: ", fname
            pass

    edges = IO.load_edges(maskname, band, wavops)
    n_slits = len(edges[0])

    for i in xrange(1, n_slits + 1):
        S = "S%2.2i" % (i)
        a = get_path(os.path.join(path, "eps_%s_%s_%s.fits" % (band, suffix, S)))

        a_h = pf.open(a)[0].header
        obj = a_h["object"]

        b = os.path.join(path, "%s_%s_%s_eps.fits" % (maskname, band, obj)) + gz(a)
        os.rename(a, b)

        a = get_path(os.path.join(path, "ivar_%s_%s_%s.fits" % (band, suffix, S)))
        a_h = pf.open(a)[0].header
        obj = a_h["object"]

        b = os.path.join(path, "%s_%s_%s_ivar.fits" % (maskname, band, obj)) + gz(a)
        os.rename(a, b)

    a = get_path(os.path.join(path, "eps_%s_%s_%s.fits" % (maskname, suffix, band)))
    b = os.path.join(path, "%s_%s_eps.fits" % (maskname, band)) + gz(a)
    os.rename(a, b)

    a = get_path(os.path.join(path, "snrs_%s_%s_%s.fits" % (maskname, suffix, band)))
    b = os.path.join(path, "%s_%s_snrs.fits" % (maskname, band)) + gz(a)
    os.rename(a, b)

    a = get_path(os.path.join(path, "ivars_%s_%s_%s.fits" % (maskname, suffix, band)))
    b = os.path.join(path, "%s_%s_ivars.fits" % (maskname, band)) + gz(a)
    os.rename(a, b)

def imdiff(A, B, maskname, band, header, options):
    s = "[0]"

    targname = A[1]["targname"].rstrip(" ")
    if targname == "":
        objname = A[1]["object"].replace(" ", "_")
    else:
        objname = targname.replace(" ", "_")

    operand1 = A[0] + '[0]'
    operand2 = B[0] + '[0]'

    imnumA = A[0].split('_')[-1].rstrip(".fits")
    imnumB = B[0].split('_')[-1].rstrip(".fits")

    dname = "{0}_{1}_{2}_{3}-{4}_{5}-{6}.fits".format(maskname, objname, band,
        A[1]["frameid"], B[1]["frameid"], imnumA, imnumB)

    try: os.remove(dname)
    except:pass
    print "Data Diff {0}-{1}".format(operand1,operand2)
    IO.imarith(operand1, '-', operand2, dname)

    ''' Now handle variance '''
    numreads = header["READS0"]
    RN_adu = Detector.RN / np.sqrt(numreads) / Detector.gain
    varname = "var_{0}_{1}_{2}_{3}+{4}_{5}+{6}.fits".format(maskname, objname, band,
        A[1]["frameid"], B[1]["frameid"], imnumA, imnumB)

    
    print "Var Sum {0}+{1}".format(operand1,operand2)
    IO.imarith(operand1, '+', operand2, "tmp_" + varname)
    try: os.remove(varname)
    except: pass
    print "Var add RN {0}+{1}".format(operand1,RN_adu**2)
    IO.imarith("tmp_" + varname, '+', RN_adu**2, varname)

    try: os.remove("tmp_" + varname)
    except: pass


    return dname, varname

def imcombine(files, maskname, options, flat, outname=None, shifts=None,
    extension=None):
    '''
    From a list of files it imcombine returns the imcombine of several values.
    The imcombine code also estimates the readnoise ad RN/sqrt(numreads) so
    that the variance per frame is equal to (ADU + RN^2) where RN is computed
    in ADUs.

    Arguments:
        files[]: list of full path to files to combine
        maskname: Name of mask
        options: Options dictionary
        flat[2048x2048]: Flat field (values should all be ~ 1.0)
        outname: If set, will write (see notes below for details)
            eps_[outname].fits: electron/sec file
            itimes_[outname].fits: integration time
            var_[outname].fits: Variance files
        shifts[len(files)]: If set, will "roll" each file by the 
            amount in the shifts vector in pixels. This argument
            is used when telescope tracking is poor. If you need
            to use this, please notify Keck staff about poor 
            telescope tracking.

    Returns 6-element tuple:
        header: The combined header
        electrons [2048x2048]:  e- (in e- units)
        var [2048x2048]: electrons + RN**2 (in e-^2 units)
        bs: The MOSFIRE.Barset instance
        itimes [2048x2048]: itimes (in s units)
        Nframe: The number of frames that contribute to the summed
            arrays above. If Nframe > 5 I use the sigma-clipping
            Cosmic Ray Rejection tool. If Nframe < 5 then I drop
            the max/min elements.

    Notes:

        header -- fits header
        ADUs -- The mean # of ADUs per frame
        var -- the Variance [in adu] per frame. 
        bs -- Barset
        itimes -- The _total_ integration time in second
        Nframe -- The number of frames in a stack.

        
        Thus the number of electron per second is derived as: 
            e-/sec = (ADUs * Gain / Flat) * (Nframe/itimes)

        The total number of electrons is:
            el = ADUs * Gain * Nframe


    '''

    ADUs = np.zeros((len(files), 2048, 2048))
    itimes = np.zeros((len(files), 2048, 2048))
    prevssl = None
    prevmn = None
    patternid = None
    maskname = None

    header = None

    if shifts is None:
        shifts = np.zeros(len(files))

    warnings.filterwarnings('ignore')
    for i in xrange(len(files)):
        fname = files[i]
        thishdr, data, bs = IO.readmosfits(fname, options, extension=extension)
        itimes[i,:,:] = thishdr["truitime"]

        base = os.path.basename(fname).rstrip(".fits")
        fnum = int(base.split("_")[1])
        
        if shifts[i] == 0:
            ADUs[i,:,:] = data.filled(0.0) / flat
        else:
            ADUs[i,:,:] = np.roll(data.filled(0.0) / flat, np.int(shifts[i]), axis=0)

        ''' Construct Header'''
        if header is None:
            header = thishdr

        header["imfno%3.3i" % (fnum)] =  (fname, "img%3.3i file name" % fnum)

        map(lambda x: rem_header_key(header, x), ["CTYPE1", "CTYPE2", "WCSDIM",
            "CD1_1", "CD1_2", "CD2_1", "CD2_2", "LTM1_1", "LTM2_2", "WAT0_001",
            "WAT1_001", "WAT2_001", "CRVAL1", "CRVAL2", "CRPIX1", "CRPIX2",
            "RADECSYS"])

        for card in header.cards:
            if card == '': continue
            key,val,comment = card
            
            if key in thishdr:
                if val != thishdr[key]:
                    newkey = key + ("_img%2.2i" % fnum)
                    try: header[newkey.rstrip()] = (thishdr[key], comment)
                    except: pass

#.........這裏部分代碼省略.........

def write_outputs(solutions, itime, header, maskname, band_name, plan, options):
    sky_sub_out = np.zeros((2048, 2048), dtype=np.float)
    sky_model_out = np.zeros((2048, 2048), dtype=np.float)

    p0 = plan[0].replace("'", "p")
    p1 = plan[1].replace("'", "p")
    suffix = "%s-%s" % (p0,p1)
    xroi = slice(0,2048)

    for sol in solutions:
        if not sol["ok"]: 
            continue

        yroi = slice(sol["bottom"], sol["top"])
        sky_sub_out[yroi, xroi] = sol["output"]
        sky_model_out[yroi, xroi] = sol["model"]
    
    header['BUNIT'] = 'SECOND'
    IO.writefits(itime, maskname, "itime_%s_%s_%s.fits" % (maskname, band,
        suffix), options, header=header, overwrite=True, lossy_compress=True)


    header['BUNIT'] = 'ELECTRONS/SECOND'
    IO.writefits(data, maskname, "sub_%s_%s_%s.fits" % (maskname, band,
        suffix), options, header=header, overwrite=True, lossy_compress=True)

    header['BUNIT'] = 'ELECTRONS/SECOND'
    IO.writefits(sky_sub_out, maskname, "bsub_%s_%s_%s.fits" % (maskname, band,
        suffix), options, header=header, overwrite=True)

    header['BUNIT'] = 'ELECTRONS'
    IO.writefits(Var, maskname, "var_%s_%s_%s.fits" % (maskname, band,
        suffix), options, header=header, overwrite=True, lossy_compress=True)

    header['BUNIT'] = 'ELECTRONS/SECOND'
    IO.writefits(sky_model_out, maskname, "bmod_%s_%s_%s.fits" % (maskname,
        band, suffix), options, header=header, overwrite=True,
        lossy_compress=True)

    '''Now create rectified solutions'''
    dlam = Wavelength.grating_results(band)
    hpp = np.array(Filters.hpp[band]) 
    ll_fid = np.arange(hpp[0], hpp[1], dlam)
    nspec = len(ll_fid)


    rectified = np.zeros((2048, nspec), dtype=np.float32)
    rectified_var = np.zeros((2048, nspec), dtype=np.float32)
    rectified_itime = np.zeros((2048, nspec), dtype=np.float32)

    from scipy.interpolate import interp1d
    for i in xrange(2048):
        ll = lam[1][i,:]
        ss = sky_sub_out[i,:]

        ok = np.isfinite(ll) & np.isfinite(ss) & (ll < hpp[1]) & (ll >
                hpp[0])

        if len(np.where(ok)[0]) < 100:
            continue
        f = interp1d(ll[ok], ss[ok], bounds_error=False)
        rectified[i,:] = f(ll_fid)

        f = interp1d(ll, Var[i,:], bounds_error=False)
        rectified_var[i,:] = f(ll_fid)

        f = interp1d(ll, itime[i,:], bounds_error=False)
        rectified_itime[i,:] = f(ll_fid)

    header["wat0_001"] = "system=world"
    header["wat1_001"] = "type=linear"
    header["wat2_001"] = "type=linear"
    header["dispaxis"] = 1
    header["dclog1"] = "Transform"
    header["dc-flag"] = 0
    header["type1"] = "AWAV"
    header["cunit1"] = "Angstrom"
    header["crval1"] = (ll_fid[0], "Starting wavelength Angstrom")
    header["crval2"] = 0
    header["crpix1"] = 1
    header["crpix2"] = 1
    header["cdelt1"] = 1
    header["cdelt2"] = 1
    header["cname1"] = "angstrom"
    header["cname2"] = "pixel"
    header["cd1_1"] = (dlam, "Angstrom/pixel")
    header["cd1_2"] = 0
    header["cd2_1"] = 0
    header["cd2_2"] = (1, "pixel/pixel")

    IO.writefits(rectified_itime, maskname,
        "%s_rectified_itime_%s_%s.fits" % (maskname, band_name,
        suffix), options, header=header, overwrite=True, lossy_compress=True)

    IO.writefits(rectified, maskname, "%s_rectified_%s_%s.fits" % (maskname,
        band_name, suffix), options, header=header, overwrite=True,
        lossy_compress=True)

    IO.writefits(rectified_var, maskname, "%s_rectified_var_%s_%s.fits" %
        (maskname, band_name, suffix), options, header=header, overwrite=True,
#.........這裏部分代碼省略.........