Python photUtils.Sed类代码示例

def mCalcs(airmass, bandName, ra, dec, expMJD,  FWHMeff, hwbpdict, photparams=None, sm=None):
    """
    sm : 
    """
    if photparams is None:
        photparams = PhotometricParameters()
    if sm is None:
        sm = SkyModel(observatory='LSST', mags=False, preciseAltAz=True)
    # Obtain full sky transmission at airmass
    # Note that this method is not interpolating but choosing the atmospheric transmission from
    # Modtran simulations of the closest airmass in a sequence of np.arange(1., 2.51, 0.1)
    fname = atmTransName(airmass)
    print(fname)
    atmTrans = np.loadtxt(fname)
    wave, trans = hwbpdict[bandName].multiplyThroughputs(atmTrans[:, 0], atmTrans[:, 1])
    bp = Bandpass(wavelen=wave, sb=trans)
    # Set the observing condition
    sm.setRaDecMjd(lon=[ra], lat=[dec], filterNames=[bandName],
                   mjd=expMJD, degrees=False, azAlt=False)
    # Get the sky sed
    wave, spec = sm.returnWaveSpec()
    sed = Sed(wavelen=wave, flambda=spec[0])
    sed.writeSED('skySED_laptop.csv')
    m5 = calcM5(sed, bp, hwbpdict[bandName], photparams, FWHMeff)
    # Get the sky magnitude only in the band concerned
    m = sm.returnMags(bandpasses=hwbpdict)[bandName][0]
    return m5, m

    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 returnMags(self, bandpass=None):
        """
        Convert the computed spectra to magnitudes using the supplied bandpasses,
        or, if self.mags=True, just return the mags in the LSST filters

        If mags=True when initialized, return mags returns an structured array with
        dtype names u,g,r,i,z,y.
        """
        if self.mags:
            if bandpass:
                warnings.warn('Ignoring set bandpasses and returning LSST ugrizy.')
            mags = -2.5*np.log10(self.spec)+np.log10(3631.)
            # Mask out high airmass
            mags[self.mask] *= np.nan
            # Convert to a structured array
            mags = np.core.records.fromarrays(mags.transpose(),
                                              names='u,g,r,i,z,y',
                                              formats='float,'*6)
        else:
            mags = np.zeros(self.npts, dtype=float)-666
            tempSed = Sed()
            isThrough = np.where(bandpass.sb > 0)
            minWave = bandpass.wavelen[isThrough].min()
            maxWave = bandpass.wavelen[isThrough].max()
            inBand = np.where((self.wave >= minWave) & (self.wave <= maxWave))
            for i, ra in enumerate(self.ra):
                # Check that there is flux in the band, otherwise calcMag fails
                if np.max(self.spec[i, inBand]) > 0:
                    tempSed.setSED(self.wave, flambda=self.spec[i, :])
                    mags[i] = tempSed.calcMag(bandpass)

            # Mask out high airmass
            mags[self.mask] *= np.nan
        return mags

    def calcBasicColors(self, sedList, bandpassDict, makeCopy = False):

        """
        This will calculate a set of colors from a list of SED objects when there is no need to redshift
        the SEDs.

        @param [in] sedList is the set of spectral objects from the models SEDs provided by loaders in
        rgStar or rgGalaxy. NOTE: Since this uses photometryBase.manyMagCalc_list the SED objects
        will be changed.

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

        @param [in] makeCopy indicates whether or not to operate on copies of the SED objects in sedList
        since this method will change the wavelength grid.

        @param [out] modelColors is the set of colors in the Bandpasses provided for the given sedList.
        """

        modelColors = []

        for specObj in sedList:
            if makeCopy==True:
                fileSED = Sed()
                fileSED.setSED(wavelen = specObj.wavelen, flambda = specObj.flambda)
                sEDMags = bandpassDict.magListForSed(fileSED)
            else:
                sEDMags = bandpassDict.magListForSed(specObj)
            colorInfo = []
            for filtNum in range(0, len(bandpassDict)-1):
                colorInfo.append(sEDMags[filtNum] - sEDMags[filtNum+1])
            modelColors.append(colorInfo)

        return modelColors

    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 _calcColors(self, sedname='C.dat'):
     """
     Calculate the colors for a moving object with sed 'sedname'.
     """
     # Do we need to read in the LSST bandpasses?
     try:
         self.lsst
     except AttributeError:
         filterdir = os.getenv('LSST_THROUGHPUTS_BASELINE')
         filterlist = ('u', 'g', 'r', 'i', 'z', 'y')
         self.lsst ={}
         for f in filterlist:
             self.lsst[f] = Bandpass()
             self.lsst[f].readThroughput(os.path.join(filterdir, 'total_'+f+'.dat'))
         self.vband = Bandpass()
         self.vband.readThroughput('harris_V.dat')
         self.colors = {}
     # See if the sed's colors are in memory already.
     if sedname not in self.colors:
         moSed = Sed()
         moSed.readSED_flambda(sedname)
         vmag = moSed.calcMag(self.vband)
         self.colors[sedname] = {}
         for f in filterlist:
             self.colors[sedname][f] = moSed.calcMag(self.lsst[f]) - vmag
     return self.colors[sedname]

    def loadwdSEDs(self, subset = None):

        """
        By default will load all seds in wd directory. The user can also define a subset of
        what's in the directory and load only those SEDs instead. Will skip over extraneous
        files in sed folder.

        @param [in] subset is the list of the subset of files wanted if one doesn't want all files
        in the kurucz directory.

        @param [out] sedListH is the set of model SED spectra objects for Hydrogen WDs to be passed onto
        the matching routines.

        @param [out] sedListHE is the set of model SED spectra objects for Helium WDs to be passed onto
        the matching routines.
        """

        if self.wdDir is None:
            try:
                self.wdDir = str(self.sEDDir + '/' +
                                 self.specMapDict['wd'] + '/')
            except:
                raise ValueError(str('self.wdDir is None. ' +
                                     'Add path to wddirectory.'))


        files = []

        if subset is None:
            for fileName in os.listdir(self.wdDir):
                files.append(fileName)
        else:
            for fileName in subset:
                files.append(fileName)

        numFiles = len(files)
        numOn = 0

        sedListH = []
        sedListHE = []

        for fileName in files:
            if numOn % 100 == 0:
                print('Loading %i of %i: WD SEDs' % (numOn, numFiles))

            try:
                spec = Sed()
                spec.readSED_flambda(str(self.wdDir + '/' + fileName))
                spec.name = fileName
                if fileName.split("_")[1] == 'He':
                    sedListHE.append(spec)
                else:
                    sedListH.append(spec)

            except:
                continue

            numOn += 1

        return sedListH, sedListHE

    def testNoSystematicUncertainty(self):
        """
        Test that systematic uncertainty is handled correctly when set to None.
        """
        m5 = [23.5, 24.3, 22.1, 20.0, 19.5, 21.7]
        photParams = PhotometricParameters(sigmaSys=0.0)

        obs_metadata = ObservationMetaData(
            unrefractedRA=23.0, unrefractedDec=45.0, m5=m5, bandpassName=self.filterNameList
        )

        magnitudes = []
        for bp in self.bpList:
            mag = self.starSED.calcMag(bp)
            magnitudes.append(mag)

        skySedList = []

        for bp, hardware, filterName in zip(self.bpList, self.hardwareList, self.filterNameList):
            skyDummy = Sed()
            skyDummy.readSED_flambda(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "darksky.dat"))
            normalizedSkyDummy = setM5(
                obs_metadata.m5[filterName],
                skyDummy,
                bp,
                hardware,
                seeing=LSSTdefaults().seeing(filterName),
                photParams=photParams,
            )

            skySedList.append(normalizedSkyDummy)

        sigmaList = snr.calcMagError_m5(numpy.array(magnitudes), numpy.array(self.bpList), numpy.array(m5), photParams)

        for i in range(len(self.bpList)):
            snrat = snr.calcSNR_sed(
                self.starSED,
                self.bpList[i],
                skySedList[i],
                self.hardwareList[i],
                seeing=LSSTdefaults().seeing(self.filterNameList[i]),
                photParams=PhotometricParameters(),
            )

            testSNR, gamma = snr.calcSNR_m5(
                numpy.array([magnitudes[i]]),
                [self.bpList[i]],
                numpy.array([m5[i]]),
                photParams=PhotometricParameters(sigmaSys=0.0),
            )

            self.assertAlmostEqual(
                snrat, testSNR[0], 10, msg="failed on calcSNR_m5 test %e != %e " % (snrat, testSNR[0])
            )

            control = snr.magErrorFromSNR(testSNR)

            msg = "%e is not %e; failed" % (sigmaList[i], control)

            self.assertAlmostEqual(sigmaList[i], control, 10, msg=msg)

    def __init__(self, m5Col='fiveSigmaDepth', units='mag', maps=['DustMap'],
                 lsstFilter='r', wavelen_min=None , wavelen_max=None , wavelen_step=1., **kwargs ):
        """
        Args:
            m5Col (str): Column name that ('fiveSigmaDepth')
            units (str): units of the metric ('mag')
            maps (list): List of maps to use with the metric (['DustMap'])
            lsstFilter (str): Which LSST filter to calculate m5 for
            wavelen_min (float): Minimum wavength of your filter (None)
            wavelen_max (float): (None)
            wavelen_step (float): (1.)
            **kwargs:
        """

        waveMins={'u':330.,'g':403.,'r':552.,'i':691.,'z':818.,'y':950.}
        waveMaxes={'u':403.,'g':552.,'r':691.,'i':818.,'z':922.,'y':1070.}

        if lsstFilter is not None:
            wavelen_min = waveMins[lsstFilter]
            wavelen_max = waveMaxes[lsstFilter]

        self.m5Col = m5Col
        super(ExgalM5, self).__init__(col=[self.m5Col],
                                      maps=maps, units=units, **kwargs)

        testsed = Sed()
        testsed.setFlatSED(wavelen_min = wavelen_min,
                           wavelen_max = wavelen_max, wavelen_step = 1)
        self.a,self.b = testsed.setupCCMab()
        self.R_v = 3.1
        self.Coaddm5Metric = Coaddm5Metric(m5Col=m5Col)

 def testRedshiftName(self):
     testsed = Sed(self.testsed.wavelen, self.testsed.flambda, name=self.testsed.name)
     redshift = .2
     testsed.redshiftSED(redshift=redshift)
     newname = testsed.name + '_Z' + '%.2f' % (redshift)
     testsed.name = newname
     self.assertEqual(testsed.name, newname)

    def magListForSed(self, sedobj, indices=None):
        """
        Return a list of magnitudes for a single Sed object.

        @param [in] sedobj is an Sed object.  Its wavelength grid can be arbitrary.  If necessary,
        a copy will be created and resampled onto the wavelength grid of the Bandpasses before
        magnitudes are calculated.  The original Sed will be unchanged.

        @param [in] indices is an optional list of indices indicating which bandpasses to actually
        calculate magnitudes for.  Other magnitudes will be listed as numpy.NaN (i.e. this method will
        return as many magnitudes as were loaded with the loadBandpassesFromFiles methods; it will
        just return numpy.NaN for magnitudes you did not actually ask for)

        @param [out] magList is a list of magnitudes in the bandpasses stored in this BandpassDict
        """

        if sedobj.wavelen is not None:

            # If the Sed's wavelength grid agrees with self._wavelen_match to one part in
            # 10^6, just use the Sed as-is.  Otherwise, copy it and resample it onto
            # self._wavelen_match
            if sedobj._needResample(wavelen_match=self._wavelen_match):
                dummySed = Sed(wavelen=sedobj.wavelen, flambda=sedobj.flambda)
                dummySed.resampleSED(force=True, wavelen_match=self._bandpassDict.values()[0].wavelen)
            else:
                dummySed = sedobj

            return numpy.array(self._magListForSed(dummySed, indices=indices))

        else:
            return numpy.array([numpy.NaN]*len(self._bandpassDict))

def setM5(m5target, skysed, totalBandpass, hardware,
          photParams,
          FWHMeff = None):
    """
    Take an SED representing the sky and normalize it so that
    m5 (the magnitude at which an object is detected in this
    bandpass at 5-sigma) is set to some specified value.

    The 5-sigma limiting magnitude (m5) for an observation is
    determined by a combination of the telescope and camera parameters
    (such as diameter of the mirrors and the readnoise) together with the
    sky background. This method (setM5) scales a provided sky background
    Sed so that an observation would have a target m5 value, for the
    provided hardware parameters. Using the resulting Sed in the
    'calcM5' method will return this target value for m5.

    @param [in] the desired value of m5

    @param [in] skysed is an instantiation of the Sed class representing
    sky emission

    @param [in] totalBandpass is an instantiation of the Bandpass class
    representing the total throughput of the telescope (instrumentation
    plus atmosphere)

    @param [in] hardware is an instantiation of the Bandpass class representing
    the throughput due solely to instrumentation.

    @param [in] photParams is an instantiation of the
    PhotometricParameters class that carries details about the
    photometric response of the telescope.

    @param [in] FWHMeff in arcseconds

    @param [out] returns an instantiation of the Sed class that is the skysed renormalized
    so that m5 has the desired value.

    Note that the returned SED will be renormalized such that calling the method
    self.calcADU(hardwareBandpass) on it will yield the number of counts per square
    arcsecond in a given bandpass.
    """

    #This is based on the LSST SNR document (v1.2, May 2010)
    #www.astro.washington.edu/users/ivezic/Astr511/LSST_SNRdoc.pdf

    if FWHMeff is None:
        FWHMeff = LSSTdefaults().FWHMeff('r')

    skyCountsTarget = calcSkyCountsPerPixelForM5(m5target, totalBandpass, FWHMeff=FWHMeff,
                                             photParams=photParams)

    skySedOut = Sed(wavelen=numpy.copy(skysed.wavelen),
                    flambda=numpy.copy(skysed.flambda))

    skyCounts = skySedOut.calcADU(hardware, photParams=photParams) \
                    * photParams.platescale * photParams.platescale
    skySedOut.multiplyFluxNorm(skyCountsTarget/skyCounts)

    return skySedOut

    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 testFluxListForSedList(self):
        """
        Test that fluxListForSedList calculates the correct fluxes
        """

        nBandpasses = 7
        bpNameList, bpList = self.getListOfBandpasses(nBandpasses)
        testBpDict = BandpassDict(bpList, bpNameList)

        nSed = 20
        sedNameList = self.getListOfSedNames(nSed)
        magNormList = self.rng.random_sample(nSed)*5.0 + 15.0
        internalAvList = self.rng.random_sample(nSed)*0.3 + 0.1
        redshiftList = self.rng.random_sample(nSed)*5.0
        galacticAvList = self.rng.random_sample(nSed)*0.3 + 0.1

        # first, test on an SedList without a wavelenMatch
        testSedList = SedList(sedNameList, magNormList,
                              fileDir=self.sedDir,
                              internalAvList=internalAvList,
                              redshiftList=redshiftList,
                              galacticAvList=galacticAvList)

        fluxList = testBpDict.fluxListForSedList(testSedList)
        self.assertEqual(fluxList.shape[0], nSed)
        self.assertEqual(fluxList.shape[1], nBandpasses)

        for ix, sedObj in enumerate(testSedList):
            dummySed = Sed(wavelen=copy.deepcopy(sedObj.wavelen),
                           flambda=copy.deepcopy(sedObj.flambda))

            for iy, bp in enumerate(testBpDict):
                flux = dummySed.calcFlux(bpList[iy])
                self.assertAlmostEqual(flux/fluxList[ix][iy], 1.0, 2)

        # now use wavelenMatch
        testSedList = SedList(sedNameList, magNormList,
                              fileDir=self.sedDir,
                              internalAvList=internalAvList,
                              redshiftList=redshiftList,
                              galacticAvList=galacticAvList,
                              wavelenMatch=testBpDict.wavelenMatch)

        fluxList = testBpDict.fluxListForSedList(testSedList)
        self.assertEqual(fluxList.shape[0], nSed)
        self.assertEqual(fluxList.shape[1], nBandpasses)

        for ix, sedObj in enumerate(testSedList):
            dummySed = Sed(wavelen=copy.deepcopy(sedObj.wavelen),
                           flambda=copy.deepcopy(sedObj.flambda))

            for iy, bp in enumerate(testBpDict):
                flux = dummySed.calcFlux(bpList[iy])
                self.assertAlmostEqual(flux/fluxList[ix][iy], 1.0, 2)

    def setUpClass(cls):
        lsstDefaults=LSSTdefaults()
        cls.dbName = 'uncertaintyTestDB.db'
        if os.path.exists(cls.dbName):
            os.unlink(cls.dbName)

        default_obs_metadata = makePhoSimTestDB(filename=cls.dbName, size=10, radius = 5.0)
        bandpass = ['u', 'g', 'r', 'i', 'z', 'y']
        m5 = lsstDefaults._m5.values()

        cls.obs_metadata = ObservationMetaData(
                                              pointingRA = default_obs_metadata.pointingRA,
                                              pointingDec = default_obs_metadata.pointingDec,
                                              rotSkyPos = default_obs_metadata.rotSkyPos,
                                              bandpassName = bandpass,
                                              m5 = m5
                                              )

        cls.obs_metadata.setBandpassM5andSeeing(bandpassName=bandpass, m5=m5)
        cls.driver = 'sqlite'
        cls.host = ''

        cls.skySeds = []
        cls.hardwareBandpasses = []
        cls.totalBandpasses = []
        cls.bandpasses = ['u', 'g', 'r', 'i', 'z', 'y']

        components = ['detector.dat', 'm1.dat', 'm2.dat', 'm3.dat',
                      'lens1.dat', 'lens2.dat', 'lens3.dat']

        for b in cls.bandpasses:
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughput(os.path.join(lsst.utils.getPackageDir('throughputs'),
                                                      'baseline', 'total_%s.dat' % b))
            cls.totalBandpasses.append(bandpassDummy)

        for b in cls.bandpasses:
            finalComponents = []
            for c in components:
                finalComponents.append(os.path.join(lsst.utils.getPackageDir('throughputs'), 'baseline', c))
            finalComponents.append(os.path.join(lsst.utils.getPackageDir('throughputs'), 'baseline', 'filter_%s.dat' %b))
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughputList(finalComponents)
            cls.hardwareBandpasses.append(bandpassDummy)

        for i in range(len(cls.bandpasses)):
            sedDummy = Sed()
            sedDummy.readSED_flambda(os.path.join(lsst.utils.getPackageDir('throughputs'), 'baseline', 'darksky.dat'))
            normalizedSedDummy = setM5(cls.obs_metadata.m5[cls.bandpasses[i]], sedDummy,
                                       cls.totalBandpasses[i], cls.hardwareBandpasses[i],
                                       FWHMeff=lsstDefaults.FWHMeff(cls.bandpasses[i]),
                                       photParams=PhotometricParameters())

            cls.skySeds.append(normalizedSedDummy)