Python IO.writefits方法代码示例

# 需要导入模块: from MOSFIRE import IO [as 别名]
# 或者: from MOSFIRE.IO import writefits [as 别名]
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)

# 需要导入模块: from MOSFIRE import IO [as 别名]
# 或者: from MOSFIRE.IO import writefits [as 别名]
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)

# 需要导入模块: from MOSFIRE import IO [as 别名]
# 或者: from MOSFIRE.IO import writefits [as 别名]
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)

# 需要导入模块: from MOSFIRE import IO [as 别名]
# 或者: from MOSFIRE.IO import writefits [as 别名]

#.........这里部分代码省略.........
    el_per_sec = electrons / itimes

    output = np.zeros((2048, 2048))
    exptime = np.zeros((2048, 2048))

    numreads = header["READS0"]
    RN_adu = Detector.RN / np.sqrt(numreads) / Detector.gain
    RN = Detector.RN / np.sqrt(numreads)

    # Cosmic ray rejection code begins here. This code construction the
    # electrons and itimes arrays.
    if len(files) >= 9:
        print "Sigclip CRR"
        srt = np.argsort(electrons, axis=0, kind='quicksort')
        shp = el_per_sec.shape
        sti = np.ogrid[0:shp[0], 0:shp[1], 0:shp[2]]

        electrons = electrons[srt, sti[1], sti[2]]
        el_per_sec = el_per_sec[srt, sti[1], sti[2]]
        itimes = itimes[srt, sti[1], sti[2]]

        # Construct the mean and standard deviation by dropping the top and bottom two 
        # electron fluxes. This is temporary.
        mean = np.mean(el_per_sec[2:-2,:,:], axis = 0)
        std = np.std(el_per_sec[2:-2,:,:], axis = 0)

        drop = np.where( (el_per_sec > (mean+std*4)) | (el_per_sec < (mean-std*4)) )
        print "dropping: ", len(drop[0])
        electrons[drop] = 0.0
        itimes[drop] = 0.0

        electrons = np.sum(electrons, axis=0)
        itimes = np.sum(itimes, axis=0)
        Nframe = len(files) 

    elif len(files) > 5:
        print "WARNING: Drop min/max CRR"
        srt = np.argsort(el_per_sec,axis=0)
        shp = el_per_sec.shape
        sti = np.ogrid[0:shp[0], 0:shp[1], 0:shp[2]]

        electrons = electrons[srt, sti[1], sti[2]]
        itimes = itimes[srt, sti[1], sti[2]]

        electrons = np.sum(electrons[1:-1,:,:], axis=0)
        itimes = np.sum(itimes[1:-1,:,:], axis=0)

        Nframe = len(files) - 2

    else:
        print "WARNING: CRR through median filtering"
        for i in xrange(len(files)):
            el = electrons[i,:,:]
            it = itimes[i,:,:]
            el_mf = scipy.signal.medfilt(el, 5)

            bad = np.abs(el - el_mf) / np.abs(el) > 10.0
            el[bad] = 0.0
            it[bad] = 0.0

            electrons[i,:,:] = el
            itimes[i,:,:] = it

        electrons = np.sum(electrons, axis=0)
        itimes = np.sum(itimes, axis=0)
        Nframe = len(files) 


    ''' Now handle variance '''
    numreads = header["READS0"]
    RN_adu = Detector.RN / np.sqrt(numreads) / Detector.gain
    RN = Detector.RN / np.sqrt(numreads)

    var = (electrons + RN**2) 

    ''' Now mask out bad pixels '''
    electrons[data.mask] = np.nan
    var[data.mask] = np.inf

    if "RN" in header: raise Exception("RN Already populated in header")
    header['RN'] = ("%1.3f" , "Read noise in e-")
    header['NUMFRM'] = (Nframe, 'Typical number of frames in stack')

    if outname is not None:
        header['BUNIT'] = 'ELECTRONS/SECOND'
        IO.writefits(np.float32(electrons/itimes), maskname, "eps_%s" % (outname),
                options, header=header, overwrite=True)

        # Update itimes after division in order to not introduce nans
        itimes[data.mask] = 0.0

        header['BUNIT'] = 'ELECTRONS^2'
        IO.writefits(var, maskname, "var_%s" % (outname),
                options, header=header, overwrite=True, lossy_compress=True)

        header['BUNIT'] = 'SECOND'
        IO.writefits(np.float32(itimes), maskname, "itimes_%s" % (outname),
                options, header=header, overwrite=True, lossy_compress=True)

    return header, electrons, var, bs, itimes, Nframe

# 需要导入模块: from MOSFIRE import IO [as 别名]
# 或者: from MOSFIRE.IO import writefits [as 别名]
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,
#.........这里部分代码省略.........

# 需要导入模块: from MOSFIRE import IO [as 别名]
# 或者: from MOSFIRE.IO import writefits [as 别名]

#.........这里部分代码省略.........
        header["cunit1"] = "Angstrom"
        header["crval1"] = ll[0]
        header["crval2"] = -solution["eps_img"].shape[0]/2 - pixel_dist
        header["crpix1"] = 1
        header["crpix2"] = 1
        header["cdelt1"] = 1
        header["cdelt2"] = 1
        header["cname1"] = "angstrom"
        header["cname2"] = "pixel"
        header["cd1_1"] = ll[1]-ll[0]
        header["cd1_2"] = 0
        header["cd2_1"] = 0
        header["cd2_2"] = 1


        S = output.shape

        img = solution["eps_img"]
        std = solution["sd_img"]
        tms = solution["itime_img"]


        for i_solution in xrange(1,len(all_solutions)):
            print "Combining solution %i" %i_solution
            solution = all_solutions[i_solution][i_slit]
            img += solution["eps_img"]
            std += solution["sd_img"]
            tms += solution["itime_img"]

        output = np.append(output, img, 0)
        output = np.append(output, np.nan*np.zeros((3,S[1])), 0)
        snrs = np.append(snrs, img*tms/std, 0)
        snrs = np.append(snrs, np.nan*np.zeros((3,S[1])), 0)
        sdout = np.append(sdout, std, 0)
        sdout = np.append(sdout, np.nan*np.zeros((3,S[1])), 0)
        itout = np.append(itout, tms, 0)
        itout = np.append(itout, np.nan*np.zeros((3,S[1])), 0)

        header['bunit'] = ('electron/second', 'electron power')
        IO.writefits(img, maskname,
            "{0}_{1}_{2}_eps.fits".format(maskname, band, target_name), options,
            overwrite=True, header=header, lossy_compress=False)

        header['bunit'] = ('electron/second', 'sigma/itime')
        IO.writefits(std/tms, maskname,
            "{0}_{1}_{2}_sig.fits".format(maskname, band, target_name), options,
            overwrite=True, header=header, lossy_compress=False)

        header['bunit'] = ('second', 'exposure time')
        IO.writefits(tms, maskname,
            "{0}_{1}_{2}_itime.fits".format(maskname, band, target_name), options,
            overwrite=True, header=header, lossy_compress=False)

        header['bunit'] = ('', 'SNR')
        IO.writefits(img*tms/std, maskname,
            "{0}_{1}_{2}_snrs.fits".format(maskname, band, target_name), options,
            overwrite=True, header=header, lossy_compress=False)

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


    header["bunit"] = "ELECTRONS/SECOND"
    IO.writefits(output, maskname, "{0}_{1}_eps.fits".format(maskname,
        band), options, overwrite=True, header=header,
        lossy_compress=False)

    header["bunit"] = ""
    IO.writefits(snrs, maskname, "{0}_{1}_snrs.fits".format(maskname,
        band), options, overwrite=True, header=header,
        lossy_compress=False)

    header["bunit"] = "ELECTRONS/SECOND"
    IO.writefits(sdout/itout, maskname, "{0}_{1}_sig.fits".format(maskname,
        band), options, overwrite=True, header=header,
        lossy_compress=False)

    header["bunit"] = "SECOND"
    IO.writefits(itout, maskname, "{0}_{1}_itime.fits".format(maskname,
        band), options, overwrite=True, header=header,
        lossy_compress=False)