python类PrimaryHDU()的实例源码

def regrid_h2(nh3_image,h2_image):
    # Edit to write out regridded image - glue won't work if files not on same grid
    h2fits  = fits.open(h2_image)
    nh3_hdr = fits.getheader(nh3_image)
    new_h2  = FITS_tools.hcongrid.hcongrid(h2fits[0].data,h2fits[0].header,nh3_hdr)
    new_h2_hdu = fits.PrimaryHDU(new_h2,nh3_hdr)
    return new_h2_hdu

def log10_h2(h2_image):
    log_h2_data = np.log10(h2_image.data)
    log_h2_hdu = fits.PrimaryHDU(log_h2_data,h2_image.header)
    return log_h2_hdu

def make_obs_mask(region_list,file_ext):
    for region in region_list:
        rms_hdu = fits.open('{0}/{0}_NH3_11_{1}_mom0_QA_trim.fits'.format(region,file_ext))
        # Set nan pixels to zero to create mask
        obs_mask = mask_obs_area(rms_hdu[0].data)
        new_hdu = fits.PrimaryHDU(obs_mask,rms_hdu[0].header)
        new_hdu.writeto('{0}/{0}_NH3_11_{1}_obsMask.fits'.format(region,file_ext),clobber=True)

def to_fits(self):
        primary_hdu = fits.PrimaryHDU()
        image_hdu = fits.ImageHDU(self.data, self.header)
        hdulist = fits.HDUList([primary_hdu, image_hdu])
        return hdulist

def create_mock_fits():
    x = np.ones((5, 5))
    prihdu = fits.PrimaryHDU(x)

    # Single extension FITS
    img = fits.ImageHDU(data=x)
    singlehdu = fits.HDUList([prihdu, img])
    singlehdu.writeto('image.fits', clobber=True)

def wraptable2fits(cat_columns, extname):
    tbhdu = fits_from_columns(pyfits.ColDefs(cat_columns))
    hdu = pyfits.PrimaryHDU()
    import datetime, time
    now = datetime.datetime.fromtimestamp(time.time())
    nowstr = now.isoformat()
    nowstr = nowstr[:nowstr.rfind('.')]
    hdu.header['DATE'] = nowstr
    hdu.header['ANALYSIS'] = 'NWAY matching'
    tbhdu.header['EXTNAME'] = extname
    hdulist = pyfits.HDUList([hdu, tbhdu])
    return hdulist

def load_data(self):
        hdulist = fits.open(self.filename)
        print("MuseCube: Loading the cube fluxes and variances...")

        # import pdb; pdb.set_trace()
        self.cube = hdulist[1].data
        self.stat = hdulist[2].data


        # for ivar weighting ; consider creating it in init ; takes long
        # self.flux_over_ivar = self.cube / self.stat

        self.header_1 = hdulist[1].header  # Necesito el header para crear una buena copia del white.
        self.header_0 = hdulist[0].header

        if self.filename_white is None:
            print("MuseCube: No white image given, creating one.")

            w_data = self.create_white(save=False)

            w_header_0 = copy.deepcopy(self.header_0)
            w_header_1 = copy.deepcopy(self.header_1)

            # These loops remove the third dimension from the header's keywords. This is neccesary in order to
            # create the white image and preserve the cube astrometry
            for i in w_header_0.keys():
                if '3' in i:
                    del w_header_0[i]
            for i in w_header_1.keys():
                if '3' in i:
                    del w_header_1[i]

            # prepare the header
            hdu = fits.HDUList()
            hdu_0 = fits.PrimaryHDU(header=w_header_0)
            hdu_1 = fits.ImageHDU(data=w_data, header=w_header_1)
            hdu.append(hdu_0)
            hdu.append(hdu_1)
            hdu.writeto('new_white.fits', clobber=True)
            self.filename_white = 'new_white.fits'
            print("MuseCube: `new_white.fits` image saved to disk.")

def spec_to_redmonster_format(spec, fitsname, n_id=None, mag=None):
    """
    Function used to create a spectrum in the REDMONSTER software format
    :param spec: XSpectrum1D object
    :param mag: List containing 2 elements, the first is the keyword in the header that will contain the magnitud saved in the second element
    :param fitsname:  Name of the fitsfile that will be created
    :return:
    """
    from scipy import interpolate
    wave = spec.wavelength.value
    wave_log = np.log10(wave)
    n = len(wave)
    spec.wavelength = wave_log * u.angstrom
    new_wave_log = np.arange(wave_log[1], wave_log[n - 2], 0.0001)
    spec_rebined = spec.rebin(new_wv=new_wave_log * u.angstrom)
    flux = spec_rebined.flux.value
    f = interpolate.interp1d(wave_log, spec.sig.value)
    sig = f(new_wave_log)
    inv_sig = 1. / np.array(sig) ** 2
    inv_sig = np.where(np.isinf(inv_sig), 0, inv_sig)
    inv_sig = np.where(np.isnan(inv_sig), 0, inv_sig)
    hdu1 = fits.PrimaryHDU([flux])
    hdu2 = fits.ImageHDU([inv_sig])
    hdu1.header['COEFF0'] = new_wave_log[0]
    hdu1.header['COEFF1'] = new_wave_log[1] - new_wave_log[0]
    if n_id != None:
        hdu1.header['ID'] = n_id
    if mag != None:
        hdu1.header[mag[0]] = mag[1]
    hdulist_new = fits.HDUList([hdu1, hdu2])
    hdulist_new.writeto(fitsname, clobber=True)

def save(self, filename, clobber=False):
        hdulist = pf.HDUList()
        ext = '.oif.fits'
        if filename.find(ext)==-1: filename += ext
        hdu = pf.PrimaryHDU()
        hdu.header.set('EXT', 'MODEL', comment='Type of information in the HDU')
        hdu.header.set('DATE', strftime('%Y%m%dT%H%M%S'), comment='Creation Date')
        hdu.header.set('NOBJ', self.nobj, comment='Number of objects in the model')
        hdu.header.set('NPARAMS', self.nparams, comment='Total number of free parameters')
        for i, item in enumerate(self._objs):
            hdu.header.set('OBJ'+str(i), item.name, comment='Name of object '+str(i))
        allmodes = ['vis2', 't3phi', 't3amp', 'visphi', 'visamp']
        for mode in allmodes:
            number = []
            if getattr(self.oidata, mode):
                number = [len(i) for i in getattr(self.oidata, "_input_"+mode) if i is not None]
            hdu.header.set('N'+mode, int(np.sum(number)), comment='Total number of '+mode+' measurements')
            hdu.header.set('F_'+mode, len(number), comment='Number of '+mode+' files')

        hdu.header.add_comment('Written by Guillaume SCHWORER')
        hdulist.append(hdu)

        hdulist.writeto(filename, clobber=clobber)

        for item in self._objs:
            item.save(filename, append=True, clobber=clobber)

        if self._hasdata: self.oidata.save(filename, append=True, clobber=clobber)

        return filename

def test_succeeds_func_fits_hdu():
    from astropy.io import fits
    return fits.PrimaryHDU(np.arange(3 * 5).reshape((3, 5)))

def write(filename, data, **kwargs):
        from astropy.io import fits
        if isinstance(data, np.ndarray):
            data = fits.PrimaryHDU(data)
        return data.writeto(filename, **kwargs)

def tofits(outfilename, pixelarray, hdr = None, verbose = True):
    """
    Takes a 2D numpy array and write it into a FITS file.
    If you specify a header (pyfits format, as returned by fromfits()) it will be used for the image.
    You can give me boolean numpy arrays, I will convert them into 8 bit integers.
    """
    pixelarrayshape = pixelarray.shape
    if verbose :
        print "FITS export shape : (%i, %i)" % (pixelarrayshape[0], pixelarrayshape[1])

    if pixelarray.dtype.name == "bool":
        pixelarray = np.cast["uint8"](pixelarray)

    if os.path.isfile(outfilename):
        os.remove(outfilename)

    if hdr == None: # then a minimal header will be created 
        hdu = fits.PrimaryHDU(pixelarray.transpose())
    else: # this if else is probably not needed but anyway ...
        hdu = fits.PrimaryHDU(pixelarray.transpose(), hdr)

    hdu.writeto(outfilename)

    if verbose :
        print "Wrote %s" % outfilename


# Array manipulation

def make_fiber_image(Fe, header, outname, args, amp):
    print("Making Fiberextract image for %s" %op.basename(outname))
    a,b = Fe.shape
    hdu = fits.PrimaryHDU(np.array(Fe, dtype='float32'), header=header)
    hdu.header.remove('BIASSEC')
    hdu.header.remove('TRIMSEC')
    hdu.header['DATASEC'] = '[%i:%i,%i:%i]' %(1,b,1,a)
    hdu.header['CRVAL1'] = args.wvl_dict[amp][0]
    hdu.header['CDELT1'] = args.disp[amp]
    hdu.header['CD1_1'] = args.disp[amp]
    hdu.header['CRPIX1'] = 1
    write_to_fits(hdu, outname)

def write_spectrograph_image(self, side, ext, prefix):
        if not self.check_side(side):
            return None
        outname = self.build_outname(side, prefix)
        self.log.info('Making spectrograph image for %s' %op.basename(outname))
        image = []
        for i, amp in enumerate(self.side_dict[side]):
            F = self.load_file(amp)
            if F is not None:
                image.append(F[ext].data)
            else:
                image.append(np.zeros((self.N, self.D)))
        new = np.zeros((self.N*2, self.D))
        new[:self.N,:] = image[0]
        new[self.N:,:] = image[1]
        hdu = fits.PrimaryHDU(np.array(new, dtype='float32'), 
                              header=self.header)
        hdu.header.remove('BIASSEC')
        hdu.header.remove('TRIMSEC')
        hdu.header.remove('AMPSEC')
        hdu.header.remove('DETSEC')
        hdu.header.remove('CCDSEC')
        hdu.header['DATASEC']='[%i:%i,%i:%i]'%(1,self.D,1,2*self.N)
        hdu.header['CCDPOS']=side
        hdu.header.remove('CCDHALF')
        hdu.header.remove('AMPLIFIE')
        self.write_to_fits(hdu, outname)

def save(self, image_list=[], spec_list=[]):
        '''
        Save the entire amplifier include the list of fibers.  
        This property is not used often as "amp*.pkl" is large and typically
        the fibers can be loaded and the other amplifier properties quickly
        recalculated.
        '''
        self.log.info('Saving images/properties to %s' %self.path)
        fn = op.join(self.path, 'multi_%s_%s_%s_%s.fits' %(self.specid, 
                                                           self.ifuslot,
                                                           self.ifuid,
                                                           self.amp))
        fits_list = []
        for i,image in enumerate(image_list):
            if i==0:
                fits_list.append(fits.PrimaryHDU(np.array(getattr(self, image), dtype='float32')))
            else:
                fits_list.append(fits.ImageHDU(np.array(getattr(self, image), dtype='float32')))

            fits_list[-1].header['EXTNAME'] = image

        for i, spec in enumerate(spec_list):
            try:
                s = np.array([getattr(fiber, spec) for fiber in self.fibers], dtype='float32')
                fits_list.append(fits.ImageHDU(s))
                fits_list[-1].header['EXTNAME'] = spec
            except AttributeError:
                self.log.warning('Attribute %s does not exist to save' %spec)
        if fits_list:
            fits_list[0] = self.write_header(fits_list[0])
            hdu = fits.HDUList(fits_list)
            self.write_to_fits(hdu, fn)

def make_fiber_image(Fe, header, outname, args, amp):
    print("Making Fiberextract image for %s" %op.basename(outname))
    a,b = Fe.shape
    hdu = fits.PrimaryHDU(np.array(Fe, dtype='float32'), header=header)
    hdu.header.remove('BIASSEC')
    hdu.header.remove('TRIMSEC')
    hdu.header['DATASEC'] = '[%i:%i,%i:%i]' %(1,b,1,a)
    hdu.header['CRVAL1'] = args.wvl_dict[amp][0]
    hdu.header['CDELT1'] = args.disp[amp]
    hdu.header['CD1_1'] = args.disp[amp]
    hdu.header['CRPIX1'] = 1
    hdu.writeto(outname, overwrite=True)

def make_spectrograph_image(image1, image2, header, outname):
    print("Making spectrograph image for %s" %op.basename(outname))
    a,b = image1.shape
    new = np.zeros((a*2,b))
    new[:a,:] = image1
    new[a:,:] = image2
    hdu = fits.PrimaryHDU(np.array(new, dtype='float32'), header=header)
    hdu.header.remove('BIASSEC')
    hdu.header.remove('TRIMSEC')
    hdu.header['DATASEC'] = '[%i:%i,%i:%i]' %(1,b,1,2*a)
    hdu.writeto(outname, overwrite=True)

def make_amplifier_image(image, header, outname):
    print("Making amplifier image for %s" %op.basename(outname))
    a,b = image.shape
    hdu = fits.PrimaryHDU(np.array(image, dtype='float32'), header=header)
    hdu.header.remove('BIASSEC')
    hdu.header.remove('TRIMSEC')
    hdu.header['DATASEC'] = '[%i:%i,%i:%i]' %(1,b,1,a)
    hdu.writeto(outname, overwrite=True)

def make_library_image(amp_image, header, outname, fits_list, for_bias=True):
    a,b = amp_image.shape
    overscan = []
    trimsec = []
    bias = re.split('[\[ \] \: \,]', header['BIASSEC'])[1:-1]
    biassec = [int(t)-((i+1)%2) for i,t in enumerate(bias)]  
    trim = re.split('[\[ \] \: \,]', header['TRIMSEC'])[1:-1]
    trimsec = [int(t)-((i+1)%2) for i,t in enumerate(trim)]  
    for F in fits_list:
        overscan.append(biweight_location(F[0].data[biassec[2]:biassec[3],
                                                    biassec[0]:biassec[1]]))
        del F[0].data
    A = []                
    for j,hdu in enumerate(fits_list):
        if for_bias:
            A.append(biweight_filter2d(hdu[0].data[:,i], (25,5),(5,1)) 
            - overscan[j])
        else:
            A.append(hdu[0].data - overscan[j])
    amp_image[:,i] = biweight_location(A, axis=(0,))
    good = np.isfinite(amp_image[:,i])
    amp_image[:,i] = np.interp(np.arange(a), np.arange(a)[good], 
                                   amp_image[good,i])
    for hdu in fits_list:
        del hdu[0].data
        hdu.close()

    hdu = fits.PrimaryHDU(np.array(amp_image[trimsec[2]:trimsec[3],
                                             trimsec[0]:trimsec[1]], 
                                   dtype='float32'), 
                          header=header)
    hdu.header.remove('BIASSEC')
    hdu.header.remove('TRIMSEC')
    hdu.header['DATASEC'] = '[%i:%i,%i:%i]' %(1,trimsec[1]-trimsec[0],1,a)
    hdu.writeto(outname, overwrite=True)

def check_bias(args, amp, folder, edge=3, width=10):
    # Create empty lists for the left edge jump, right edge jump, and structure
    left_edge, right_edge, structure, overscan = [], [], [], []

    # Select only the bias frames that match the input amp, e.g., "RU"
    sel = [i for i, v in enumerate(args.bia_list) if v.amp == amp]
    log = args.bia_list[sel[0]].log

    overscan_list = [[v.overscan_value for i, v in enumerate(args.bia_list)
                      if v.amp == amp]]
    overscan = biweight_location(overscan_list)
    log.info('Overscan value for %s: %0.3f' % (amp, overscan))
    # Loop through the bias list and measure the jump/structure
    big_array = np.array([v.image for v in itemgetter(*sel)(args.bia_list)])
    if args.quick:
        func = np.median
    else:
        func = biweight_location
    masterbias = func(big_array, axis=(0,))

    a, b = masterbias.shape
    hdu = fits.PrimaryHDU(np.array(masterbias, dtype='float32'))

    log.info('Writing masterbias_%s.fits' % (amp))
    write_fits(hdu,
               op.join(folder, 'masterbias_%s_%s.fits' % (args.specid, amp)))

    left_edge = func(masterbias[:, edge:edge+width])
    right_edge = func(masterbias[:, (b-width-edge):(b-edge)])
    structure = func(masterbias[:, edge:(b-edge)], axis=(0,))

    log.info('Left edge - Overscan, Right edge - Overscan: %0.3f, %0.3f'
             % (left_edge, right_edge))
    return left_edge, right_edge, structure, overscan, masterbias