Python photUtils.Bandpass类代码示例

    def testStellarPhotometricUncertainties(self):
        """
        Test in the case of a catalog of stars
        """
        lsstDefaults = LSSTdefaults()
        starDB = testStarsDBObj(driver=self.driver, host=self.host, database=self.dbName)
        starCat = testStarCatalog(starDB, obs_metadata=self.obs_metadata)
        phot = PhotometryStars()

        ct = 0
        for line in starCat.iter_catalog():
            starSed = Sed()
            starSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                                 defaultSpecMap[line[14]]))
            imsimband = Bandpass()
            imsimband.imsimBandpass()
            fNorm = starSed.calcFluxNorm(line[15], imsimband)
            starSed.multiplyFluxNorm(fNorm)

            aV = numpy.float(line[16])
            a_int, b_int = starSed.setupCCMab()
            starSed.addCCMDust(a_int, b_int, A_v=aV)

            for i in range(len(self.bandpasses)):
                controlSigma = calcMagError_sed(starSed, self.totalBandpasses[i],
                                             self.skySeds[i],
                                             self.hardwareBandpasses[i],
                                             FWHMeff=lsstDefaults.FWHMeff(self.bandpasses[i]),
                                             photParams=PhotometricParameters())

                testSigma = line[8+i]
                self.assertAlmostEqual(controlSigma, testSigma, 4)
                ct += 1
        self.assertGreater(ct, 0)

    def testLoadTotalBandpassesFromFiles(self):
        """
        Test that the class method loadTotalBandpassesFromFiles produces the
        expected result
        """

        bandpassDir = os.path.join(getPackageDir('sims_photUtils'), 'tests', 'cartoonSedTestData')
        bandpassNames = ['g', 'r', 'u']
        bandpassRoot = 'test_bandpass_'

        bandpassDict = BandpassDict.loadTotalBandpassesFromFiles(bandpassNames=bandpassNames,
                                                                 bandpassDir=bandpassDir,
                                                                 bandpassRoot = bandpassRoot)

        controlBandpassList = []
        for bpn in bandpassNames:
            dummyBp = Bandpass()
            dummyBp.readThroughput(os.path.join(bandpassDir,bandpassRoot+bpn+'.dat'))
            controlBandpassList.append(dummyBp)

        wMin = controlBandpassList[0].wavelen[0]
        wMax = controlBandpassList[0].wavelen[-1]
        wStep = controlBandpassList[0].wavelen[1]-controlBandpassList[0].wavelen[0]

        for bp in controlBandpassList:
            bp.resampleBandpass(wavelen_min=wMin, wavelen_max=wMax, wavelen_step=wStep)

        for test, control in zip(bandpassDict.values(), controlBandpassList):
            numpy.testing.assert_array_almost_equal(test.wavelen, control.wavelen, 19)
            numpy.testing.assert_array_almost_equal(test.sb, control.sb, 19)

    def loadTotalBandpassesFromFiles(cls,
                                    bandpassNames=['u', 'g', 'r', 'i', 'z', 'y'],
                                    bandpassDir = os.path.join(getPackageDir('throughputs'),'baseline'),
                                    bandpassRoot = 'total_'):
        """
        This will take the list of band passes named by bandpassNames and load them into
        a BandpassDict

        The bandpasses loaded this way are total bandpasses: they account for instrumental
        and atmospheric transmission.

        @param [in] bandpassNames is a list of names identifying each filter.
        Defaults to ['u', 'g', 'r', 'i', 'z', 'y']

        @param [in] bandpassDir is the name of the directory where the bandpass files are stored

        @param [in] bandpassRoot contains the first part of the bandpass file name, i.e., it is assumed
        that the bandpasses are stored in files of the type

        bandpassDir/bandpassRoot_bandpassNames[i].dat

        if we want to load bandpasses for a telescope other than LSST, we would do so
        by altering bandpassDir and bandpassRoot

        @param [out] bandpassDict is a BandpassDict containing the loaded throughputs
        """

        bandpassList = []

        for w in bandpassNames:
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughput(os.path.join(bandpassDir,"%s.dat" % (bandpassRoot + w)))
            bandpassList.append(bandpassDummy)

        return cls(bandpassList, bandpassNames)

    def test_norm(self):
        """
        Test that the special test case getImsimFluxNorm
        returns the same value as calling calcFluxNorm actually
        passing in the imsim Bandpass
        """

        bp = Bandpass()
        bp.imsimBandpass()

        rng = np.random.RandomState(1123)
        wavelen = np.arange(300.0, 2000.0, 0.17)

        for ix in range(10):
            flux = rng.random_sample(len(wavelen))*100.0
            sed = Sed()
            sed.setSED(wavelen=wavelen, flambda=flux)
            magmatch = rng.random_sample()*5.0 + 10.0

            control = sed.calcFluxNorm(magmatch, bp)
            test = getImsimFluxNorm(sed, magmatch)

            # something about how interpolation is done in Sed means
            # that the values don't come out exactly equal.  They come
            # out equal to 8 seignificant digits, though.
            self.assertEqual(control, test)

    def testMags(self):
        """
        Test that the interpolated mags are similar to mags computed from interpolated spectra
        """

        throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
        filters = ['u', 'g', 'r', 'i', 'z', 'y']

        bps = []
        for filterName in filters:
            bp = np.loadtxt(os.path.join(throughPath, 'filter_%s.dat' % filterName),
                            dtype=zip(['wave', 'trans'], [float]*2))
            lsst_bp = Bandpass()
            lsst_bp.setBandpass(bp['wave'], bp['trans'])
            bps.append(lsst_bp)

        sm1 = self.sm_spec
        sm1.setRaDecMjd([36.], [-68.], 49353.18, degrees=True)
        mags1 = []
        for bp in bps:
            mags1.append(sm1.returnMags(bandpass=bp))
        mags1 = np.array(mags1)

        sm2 = self.sm_mags
        sm2.setRaDecMjd([36.], [-68.], 49353.18, degrees=True)
        mag2 = sm2.returnMags()
        for i, filtername in enumerate(filters):
            np.testing.assert_allclose(mags1[i, :], mag2[filtername], rtol=1e-4)

def calcEffWavelen(hardware, title, throughputDir=None):
    photParams = PhotometricParameters()
    lsstDefaults = LSSTdefaults()
    atmos = {}
    stdAtmoFile = os.path.join(os.getenv('SYSENG_THROUGHPUTS_DIR'), 'siteProperties/pachonModtranAtm_12.dat')
    atmos['std'] = Bandpass()
    atmos['std'].readThroughput(stdAtmoFile)
    multiAtmos = False
    if throughputDir is None:
        throughputDir = os.getenv('THROUGHPUTS_DIR')
    if throughputDir is not None:
        multiAtmos = True
        Xarr = np.arange(1.0, 2.55, 0.1)
        for X in Xarr:
            atmos['%.1f' %X] = Bandpass()
            atmos['%.1f' %X].readThroughput(os.path.join(throughputDir,
                                                         'atmos/atmos_%d.dat' %(int(X*10))))

    atmoskeys = sorted(atmos.keys())
    print title
    print '    %s' %('    '.join(atmoskeys))
    system = Bandpass()
    effsb = {}
    for f in ['u', 'g', 'r', 'i', 'z', 'y']:
        writestring = '%s ' %f
        for k in atmoskeys:
            system.wavelen, system.sb = hardware[f].multiplyThroughputs(atmos[k].wavelen, atmos[k].sb)
            effphi, effsb[k] = system.calcEffWavelen()
            writestring += '%.2f ' %(effsb[k])
        print writestring

    def testMags(self):
        """
        Test that the interpolated mags are similar to mags computed from interpolated spectra
        """

        throughPath = os.path.join(getPackageDir("throughputs"), "baseline")
        filters = ["u", "g", "r", "i", "z", "y"]

        bps = []
        for filterName in filters:
            bp = np.loadtxt(
                os.path.join(throughPath, "filter_%s.dat" % filterName), dtype=zip(["wave", "trans"], [float] * 2)
            )
            lsst_bp = Bandpass()
            lsst_bp.setBandpass(bp["wave"], bp["trans"])
            bps.append(lsst_bp)

        sm1 = sb.SkyModel()
        sm1.setRaDecMjd([36.0], [-68.0], 49353.18, degrees=True)
        sm1.computeSpec()
        mags1 = []
        for bp in bps:
            mags1.append(sm1.computeMags(bandpass=bp))
        mags1 = np.array(mags1)

        sm2 = sb.SkyModel(mags=True)
        sm2.setRaDecMjd([36.0], [-68.0], 49353.18, degrees=True)
        sm2.computeSpec()
        mag2 = sm2.computeMags()
        np.testing.assert_allclose(mags1, mag2.T, rtol=1e-4)

 def setUp(self):
     self.normband = Bandpass()
     self.normband.imsimBandpass()
     self.uband = Bandpass()
     self.uband.readThroughput(os.path.join(getPackageDir('sims_photUtils'), 'tests',
                               'cartoonSedTestData', 'test_bandpass_u.dat'))
     self.gband = Bandpass()
     self.gband.readThroughput(os.path.join(getPackageDir('sims_photUtils'), 'tests',
                               'cartoonSedTestData', 'test_bandpass_g.dat'))

def fit_invsnr(mags, floor, m5, bandpass_name='r'):
    filterdir = os.getenv('LSST_THROUGHPUTS_BASELINE')
    bandpass = Bandpass()
    bandpass.readThroughput(os.path.join(filterdir, 'total_%s.dat'%bandpass_name))
    invsnr_arr = []
    for mag in mags:
        snr = flat_sed_snr(mag, bandpass, m5)
        invsnr_arr.append(math.hypot(1./snr, floor))
    return invsnr_arr

    def get_phosimVars(self):
        """
        Obtain variables sedFilepath to be used to obtain unique filenames
        for each SED for phoSim and MagNorm which is also used. Note that aside
        from acting as a getter, this also writes spectra to 
        `self.sn_sedfile_prefix`snid_mjd_band.dat for each observation of
        interest
        """
        # construct the unique filename
        # method: snid_mjd(to 4 places of decimal)_bandpassname
        mjd = "_{:0.4f}_".format(self.mjdobs)
        mjd += self.obs_metadata.bandpass + '.dat'
        fnames = np.array([self.sn_sedfile_prefix + str(int(elem)) + mjd
                           if isinstance(elem, numbers.Number) else
                           self.sn_sedfile_prefix + str(elem) + mjd
                           for elem in self.column_by_name('snid')], dtype='str')

        c, x1, x0, t0, z = self.column_by_name('c'),\
            self.column_by_name('x1'),\
            self.column_by_name('x0'),\
            self.column_by_name('t0'),\
            self.column_by_name('redshift')

        bp = Bandpass()
        bp.imsimBandpass()

        magNorms = np.zeros(len(fnames))

        snobject = SNObject()
        snobject.rectifySED = True
        for i in range(len(self.column_by_name('snid'))):
            # if t0 is nan, this was set by the catalog for dim SN, or SN
            #   outside redshift range, We will not provide a SED file for these
            if np.isnan(t0[i]):
                magNorms[i] = np.nan
                fnames[i] = None

            else:
                snobject.set(c=c[i], x1=x1[i], x0=x0[i], t0=t0[i],
                             z=z[i])
                if snobject.modelOutSideTemporalRange == 'zero':
                    if self.mjdobs > snobject.maxtime() or self.mjdobs < snobject.mintime():
                        magNorms[i] = np.nan
                        fnames[i] = None

                # SED in rest frame
                sed = snobject.SNObjectSourceSED(time=self.mjdobs)
                try:
                    magNorms[i] = sed.calcMag(bandpass=bp)
                except:
                    # sed.flambda = 1.0e-20
                    magNorms[i] = 1000.  # sed.calcMag(bandpass=bp)

                if self.writeSedFile:
                    sed.writeSED(fnames[i])

        return (fnames, magNorms)

def read_stdatmo(override_filename = None):
    # read the standard atmosphere bandpass file, precomputed by MODTRAN & DaveBurke.
    # this is closely equivalent to atmos_12.dat
    atmosdir = os.getenv("LSST_THROUGHPUTS_ATMOS")
    atmos_bp = Bandpass()
    if override_filename == None:
        atmos_bp.readThroughput(os.path.join(atmosdir, "atmos_std.dat"))
    else:
        atmos_bp.readThroughput(os.path.join(atmosdir, override_filename))
    return atmos_bp

    def setUp(self):

        starFileName = os.path.join(lsst.utils.getPackageDir("sims_sed_library"), "starSED")
        starFileName = os.path.join(starFileName, "kurucz", "km20_5750.fits_g40_5790.gz")
        starName = os.path.join(lsst.utils.getPackageDir("sims_sed_library"), starFileName)
        self.starSED = Sed()
        self.starSED.readSED_flambda(starName)
        imsimband = Bandpass()
        imsimband.imsimBandpass()
        fNorm = self.starSED.calcFluxNorm(22.0, imsimband)
        self.starSED.multiplyFluxNorm(fNorm)

        hardwareDir = os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline")
        componentList = ["detector.dat", "m1.dat", "m2.dat", "m3.dat", "lens1.dat", "lens2.dat", "lens3.dat"]
        self.skySed = Sed()
        self.skySed.readSED_flambda(os.path.join(hardwareDir, "darksky.dat"))

        totalNameList = ["total_u.dat", "total_g.dat", "total_r.dat", "total_i.dat", "total_z.dat", "total_y.dat"]

        self.bpList = []
        self.hardwareList = []
        for name in totalNameList:
            dummy = Bandpass()
            dummy.readThroughput(os.path.join(hardwareDir, name))
            self.bpList.append(dummy)

            dummy = Bandpass()
            hardwareNameList = [os.path.join(hardwareDir, name)]
            for component in componentList:
                hardwareNameList.append(os.path.join(hardwareDir, component))
            dummy.readThroughputList(hardwareNameList)
            self.hardwareList.append(dummy)

        self.filterNameList = ["u", "g", "r", "i", "z", "y"]

    def getFilesAndBandpasses(
        self,
        catalog,
        nameRoot=None,
        bandpassDir=os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline"),
        bandpassRoot="total_",
    ):

        """
        Take a GalSimCatalog.  Return a list of fits files and and OrderedDict of bandpasses associated
        with that catalog

        @param [in] catalog is a GalSimCatalog instantiation

        @param [in] nameRoot is the nameRoot prepended to the fits files output by that catalog

        @param[in] bandpassDir is the directory where bandpass files can be found

        @param [in] bandpassRoot is the root of the name of the bandpass files

        @param [out] listOfFiles is a list of the names of the fits files written by this catalog

        @param [out] bandpassDict is an OrderedDict of Bandpass instantiations corresponding to the
        filters in this catalog.
        """

        # write the fits files
        catalog.write_images(nameRoot=nameRoot)

        # a list of bandpasses over which we are integraging
        listOfFilters = []
        listOfFiles = []

        # read in the names of all of the written fits files directly from the
        # InstanceCatalog's GalSimInterpreter
        # Use AFW to read in the FITS files and calculate the ADU
        for name in catalog.galSimInterpreter.detectorImages:
            if nameRoot is not None:
                name = nameRoot + "_" + name

            listOfFiles.append(name)

            if name[-6] not in listOfFilters:
                listOfFilters.append(name[-6])

        bandpassDict = OrderedDict()
        for filterName in listOfFilters:
            bandpassName = os.path.join(bandpassDir, bandpassRoot + filterName + ".dat")
            bandpass = Bandpass()
            bandpass.readThroughput(bandpassName)
            bandpassDict[filterName] = bandpass

        return listOfFiles, bandpassDict

    def setUpClass(cls):
        cls.scratchDir = os.path.join(getPackageDir('sims_GalSimInterface'), 'tests', 'scratchSpace')
        cls.obs = ObservationMetaData(pointingRA=122.0, pointingDec=-29.1,
                                      mjd=57381.2, rotSkyPos=43.2)


        cls.camera = camTestUtils.CameraWrapper().camera

        cls.dbFileName = os.path.join(cls.scratchDir, 'allowed_chips_test_db.txt')
        if os.path.exists(cls.dbFileName):
            os.unlink(cls.dbFileName)

        cls.controlSed = Sed()
        cls.controlSed.readSED_flambda(os.path.join(getPackageDir('sims_sed_library'),
                                               'flatSED','sed_flat.txt.gz'))
        cls.magNorm = 18.1
        imsim = Bandpass()
        imsim.imsimBandpass()
        ff = cls.controlSed.calcFluxNorm(cls.magNorm, imsim)
        cls.controlSed.multiplyFluxNorm(ff)
        a_x, b_x = cls.controlSed.setupCCMab()
        cls.controlSed.addCCMDust(a_x, b_x, A_v=0.1, R_v=3.1)
        bpd = BandpassDict.loadTotalBandpassesFromFiles()
        pp = PhotometricParameters()
        cls.controlADU = cls.controlSed.calcADU(bpd['u'], pp)
        cls.countSigma = np.sqrt(cls.controlADU/pp.gain)

        cls.x_pix = 50
        cls.y_pix = 50

        x_list = []
        y_list = []
        name_list = []
        for dd in cls.camera:
            x_list.append(cls.x_pix)
            y_list.append(cls.y_pix)
            name_list.append(dd.getName())

        x_list = np.array(x_list)
        y_list = np.array(y_list)

        ra_list, dec_list = raDecFromPixelCoords(x_list, y_list, name_list,
                                                 camera=cls.camera, obs_metadata=cls.obs,
                                                 epoch=2000.0)

        dra_list = 3600.0*(ra_list-cls.obs.pointingRA)
        ddec_list = 3600.0*(dec_list-cls.obs.pointingDec)

        create_text_catalog(cls.obs, cls.dbFileName, dra_list, ddec_list,
                            mag_norm=[cls.magNorm]*len(dra_list))

        cls.db = allowedChipsFileDBObj(cls.dbFileName, runtable='test')

    def calcMagNorm(self, objectMags, sedObj, bandpassDict, mag_error = None,
                    redshift = None, filtRange = None):

        """
        This will find the magNorm value that gives the closest match to the magnitudes of the object
        using the matched SED. Uses scipy.optimize.leastsq to find the values of fluxNorm that minimizes
        the function: ((flux_obs - (fluxNorm*flux_model))/flux_error)**2.

        @param [in] objectMags are the magnitude values for the object with extinction matching that of
        the SED object. In the normal case using the selectSED routines above it will be dereddened mags.

        @param [in] sedObj is an Sed class instance that is set with the wavelength and flux of the
        matched SED

        @param [in] bandpassDict is a BandpassDict class instance with the Bandpasses set to those
        for the magnitudes given for the catalog object

        @param [in] mag_error are provided error values for magnitudes in objectMags. If none provided
        then this defaults to 1.0. This should be an array of the same length as objectMags.

        @param [in] redshift is the redshift of the object if the magnitude is observed

        @param [in] filtRange is a selected range of filters specified by their indices in the bandpassList
        to match up against. Used when missing data in some magnitude bands.

        @param [out] bestMagNorm is the magnitude normalization for the given magnitudes and SED
        """

        import scipy.optimize as opt

        sedTest = Sed()
        sedTest.setSED(sedObj.wavelen, flambda = sedObj.flambda)
        if redshift is not None:
            sedTest.redshiftSED(redshift)
        imSimBand = Bandpass()
        imSimBand.imsimBandpass()
        zp = -2.5*np.log10(3631)  #Note using default AB zeropoint
        flux_obs = np.power(10,(objectMags + zp)/(-2.5))
        sedTest.resampleSED(wavelen_match=bandpassDict.wavelenMatch)
        sedTest.flambdaTofnu()
        flux_model = sedTest.manyFluxCalc(bandpassDict.phiArray, bandpassDict.wavelenStep)
        if filtRange is not None:
            flux_obs = flux_obs[filtRange]
            flux_model = flux_model[filtRange]
        if mag_error is None:
            flux_error = np.ones(len(flux_obs))
        else:
            flux_error = np.abs(flux_obs*(np.log(10)/(-2.5))*mag_error)
        bestFluxNorm = opt.leastsq(lambda x: ((flux_obs - (x*flux_model))/flux_error), 1.0)[0][0]
        sedTest.multiplyFluxNorm(bestFluxNorm)
        bestMagNorm = sedTest.calcMag(imSimBand)
        return bestMagNorm