本文整理汇总了Python中lsst.sims.photUtils.Sed.Sed.multiplyFluxNorm方法的典型用法代码示例。如果您正苦于以下问题:Python Sed.multiplyFluxNorm方法的具体用法?Python Sed.multiplyFluxNorm怎么用?Python Sed.multiplyFluxNorm使用的例子?那么, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类lsst.sims.photUtils.Sed.Sed
的用法示例。
在下文中一共展示了Sed.multiplyFluxNorm方法的12个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: testAlternateBandpassesStars
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testAlternateBandpassesStars(self):
"""
This will test our ability to do photometry using non-LSST bandpasses.
It will first calculate the magnitudes using the getters in cartoonPhotometryStars.
It will then load the alternate bandpass files 'by hand' and re-calculate the magnitudes
and make sure that the magnitude values agree. This is guarding against the possibility
that some default value did not change and the code actually ended up loading the
LSST bandpasses.
"""
obs_metadata_pointed = ObservationMetaData(
mjd=2013.23, boundType="circle", unrefractedRA=200.0, unrefractedDec=-30.0, boundLength=1.0
)
bandpassDir = os.path.join(lsst.utils.getPackageDir("sims_photUtils"), "tests", "cartoonSedTestData")
cartoon_dict = BandpassDict.loadTotalBandpassesFromFiles(
["u", "g", "r", "i", "z"], bandpassDir=bandpassDir, bandpassRoot="test_bandpass_"
)
testBandPasses = {}
keys = ["u", "g", "r", "i", "z"]
bplist = []
for kk in keys:
testBandPasses[kk] = Bandpass()
testBandPasses[kk].readThroughput(os.path.join(bandpassDir, "test_bandpass_%s.dat" % kk))
bplist.append(testBandPasses[kk])
sedObj = Sed()
phiArray, waveLenStep = sedObj.setupPhiArray(bplist)
sedFileName = os.path.join(lsst.utils.getPackageDir("sims_sed_library"), "starSED", "kurucz")
sedFileName = os.path.join(sedFileName, "km20_5750.fits_g40_5790.gz")
ss = Sed()
ss.readSED_flambda(sedFileName)
controlBandpass = Bandpass()
controlBandpass.imsimBandpass()
ff = ss.calcFluxNorm(22.0, controlBandpass)
ss.multiplyFluxNorm(ff)
testMags = cartoon_dict.magListForSed(ss)
ss.resampleSED(wavelen_match=bplist[0].wavelen)
ss.flambdaTofnu()
mags = -2.5 * numpy.log10(numpy.sum(phiArray * ss.fnu, axis=1) * waveLenStep) - ss.zp
self.assertTrue(len(mags) == len(testMags))
self.assertTrue(len(mags) > 0)
for j in range(len(mags)):
self.assertAlmostEqual(mags[j], testMags[j], 10)
示例2: calcMagNorm
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
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.values()[0].wavelen)
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
示例3: testMatchToRestFrame
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testMatchToRestFrame(self):
"""Test that Galaxies with no effects added into catalog mags are matched correctly."""
np.random.seed(42)
galPhot = BandpassDict.loadTotalBandpassesFromFiles()
imSimBand = Bandpass()
imSimBand.imsimBandpass()
testMatching = selectGalaxySED(galDir = self.testSpecDir)
testSEDList = testMatching.loadBC03()
testSEDNames = []
testMags = []
testMagNormList = []
magNormStep = 1
for testSED in testSEDList:
getSEDMags = Sed()
testSEDNames.append(testSED.name)
getSEDMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
testMagNorm = np.round(np.random.uniform(20.0,22.0),magNormStep)
testMagNormList.append(testMagNorm)
fluxNorm = getSEDMags.calcFluxNorm(testMagNorm, imSimBand)
getSEDMags.multiplyFluxNorm(fluxNorm)
testMags.append(galPhot.magListForSed(getSEDMags))
#Also testing to make sure passing in non-default bandpasses works
#Substitute in nan values to simulate incomplete data.
testMags[0][1] = np.nan
testMags[0][2] = np.nan
testMags[0][4] = np.nan
testMags[1][1] = np.nan
testMatchingResults = testMatching.matchToRestFrame(testSEDList, testMags,
bandpassDict = galPhot)
self.assertEqual(None, testMatchingResults[0][0])
self.assertEqual(testSEDNames[1:], testMatchingResults[0][1:])
self.assertEqual(None, testMatchingResults[1][0])
np.testing.assert_almost_equal(testMagNormList[1:], testMatchingResults[1][1:], decimal = magNormStep)
#Test Match Errors
errMags = np.array((testMags[2], testMags[2], testMags[2], testMags[2]))
errMags[1,1] += 1. #Total MSE will be 2/(5 colors) = 0.4
errMags[2, 0:2] = np.nan
errMags[2, 3] += 1. #Total MSE will be 2/(3 colors) = 0.667
errMags[3, :] = None
errSED = testSEDList[2]
testMatchingResultsErrors = testMatching.matchToRestFrame([errSED], errMags,
bandpassDict = galPhot)
np.testing.assert_almost_equal(np.array((0.0, 0.4, 2./3.)), testMatchingResultsErrors[2][0:3],
decimal = 3)
self.assertEqual(None, testMatchingResultsErrors[2][3])
示例4: testSignalToNoise
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testSignalToNoise(self):
"""
Test that calcSNR_m5 and calcSNR_sed give similar results
"""
defaults = LSSTdefaults()
photParams = PhotometricParameters()
totalDict, hardwareDict = BandpassDict.loadBandpassesFromFiles()
skySED = Sed()
skySED.readSED_flambda(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "darksky.dat"))
m5 = []
for filt in totalDict:
m5.append(calcM5(skySED, totalDict[filt], hardwareDict[filt], photParams, seeing=defaults.seeing(filt)))
sedDir = lsst.utils.getPackageDir("sims_sed_library")
sedDir = os.path.join(sedDir, "starSED", "kurucz")
fileNameList = os.listdir(sedDir)
numpy.random.seed(42)
offset = numpy.random.random_sample(len(fileNameList)) * 2.0
for ix, name in enumerate(fileNameList):
if ix > 100:
break
spectrum = Sed()
spectrum.readSED_flambda(os.path.join(sedDir, name))
ff = spectrum.calcFluxNorm(m5[2] - offset[ix], totalDict.values()[2])
spectrum.multiplyFluxNorm(ff)
magList = []
controlList = []
magList = []
for filt in totalDict:
controlList.append(
calcSNR_sed(
spectrum, totalDict[filt], skySED, hardwareDict[filt], photParams, defaults.seeing(filt)
)
)
magList.append(spectrum.calcMag(totalDict[filt]))
testList, gammaList = calcSNR_m5(
numpy.array(magList), numpy.array(totalDict.values()), numpy.array(m5), photParams
)
for tt, cc in zip(controlList, testList):
msg = "%e != %e " % (tt, cc)
self.assertTrue(numpy.abs(tt / cc - 1.0) < 0.001, msg=msg)
示例5: uncertaintyUnitTest
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
class uncertaintyUnitTest(unittest.TestCase):
"""
Test the calculation of photometric uncertainties
"""
def setUp(self):
starName = os.path.join(lsst.utils.getPackageDir('sims_sed_library'),defaultSpecMap['km20_5750.fits_g40_5790'])
self.starSED = Sed()
self.starSED.readSED_flambda(starName)
imsimband = Bandpass()
imsimband.imsimBandpass()
fNorm = self.starSED.calcFluxNorm(22.0, imsimband)
self.starSED.multiplyFluxNorm(fNorm)
self.totalBandpasses = []
self.hardwareBandpasses = []
componentList = ['detector.dat', 'm1.dat', 'm2.dat', 'm3.dat',
'lens1.dat', 'lens2.dat', 'lens3.dat']
hardwareComponents = []
for c in componentList:
hardwareComponents.append(os.path.join(lsst.utils.getPackageDir('throughputs'),'baseline',c))
self.bandpasses = ['u', 'g', 'r', 'i', 'z', 'y']
for b in self.bandpasses:
filterName = os.path.join(lsst.utils.getPackageDir('throughputs'),'baseline','filter_%s.dat' % b)
components = hardwareComponents + [filterName]
bandpassDummy = Bandpass()
bandpassDummy.readThroughputList(components)
self.hardwareBandpasses.append(bandpassDummy)
components = components + [os.path.join(lsst.utils.getPackageDir('throughputs'),'baseline','atmos.dat')]
bandpassDummy = Bandpass()
bandpassDummy.readThroughputList(components)
self.totalBandpasses.append(bandpassDummy)
def tearDown(self):
del self.starSED
del self.bandpasses
del self.hardwareBandpasses
del self.totalBandpasses
示例6: testCalcMagNorm
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testCalcMagNorm(self):
"""Tests the calculation of magnitude normalization for an SED with the given magnitudes
in the given bandpasses."""
testUtils = matchBase()
bandpassDir = os.path.join(lsst.utils.getPackageDir('throughputs'), 'sdss')
testPhot = BandpassDict.loadTotalBandpassesFromFiles(self.filterList,
bandpassDir = bandpassDir,
bandpassRoot = 'sdss_')
unChangedSED = Sed()
unChangedSED.readSED_flambda(str(self.galDir + os.listdir(self.galDir)[0]))
imSimBand = Bandpass()
imSimBand.imsimBandpass()
testSED = Sed()
testSED.setSED(unChangedSED.wavelen, flambda = unChangedSED.flambda)
magNorm = 20.0
redVal = 0.1
testSED.redshiftSED(redVal)
fluxNorm = testSED.calcFluxNorm(magNorm, imSimBand)
testSED.multiplyFluxNorm(fluxNorm)
sedMags = testPhot.magListForSed(testSED)
stepSize = 0.001
testMagNorm = testUtils.calcMagNorm(sedMags, unChangedSED, testPhot, redshift = redVal)
# Test adding in mag_errors. If an array of np.ones is passed in we should get same result
testMagNormWithErr = testUtils.calcMagNorm(sedMags, unChangedSED, testPhot,
mag_error = np.ones(len(sedMags)), redshift = redVal)
# Also need to add in test for filtRange
sedMagsIncomp = sedMags
sedMagsIncomp[1] = None
filtRangeTest = [0, 2, 3, 4]
testMagNormFiltRange = testUtils.calcMagNorm(sedMagsIncomp, unChangedSED, testPhot,
redshift = redVal, filtRange = filtRangeTest)
self.assertAlmostEqual(magNorm, testMagNorm, delta = stepSize)
self.assertAlmostEqual(magNorm, testMagNormWithErr, delta = stepSize)
self.assertAlmostEqual(magNorm, testMagNormFiltRange, delta = stepSize)
示例7: loadGalfast
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
#.........这里部分代码省略.........
"""
Info about the following population cuts:
From Zeljko: "This color corresponds to the temperature (roughly spectral type M0) where
Kurucz models become increasingly bad, and thus we switch to empirical SEDs (the problem
is that for M and later stars, the effective surface temperature is low enough for
molecules to form, and their opacity is too complex to easily model, especially TiO)."
"""
mIn = np.where(((pop < 10) | (pop >= 20)) & (sDSSunred[:,2] - sDSSunred[:,3] > 0.59))
kIn = np.where(((pop < 10) | (pop >= 20)) & (sDSSunred[:,2] - sDSSunred[:,3] <= 0.59))
hIn = np.where((pop >= 10) & (pop < 15))
heIn = np.where((pop >= 15) & (pop < 20))
sEDNameK, magNormK, matchErrorK = selectStarSED0.findSED(listDict['kurucz'],
sDSSunred[kIn], ra[kIn], dec[kIn],
reddening = False,
colors = colorDict['kurucz'])
sEDNameM, magNormM, matchErrorM = selectStarSED0.findSED(listDict['mlt'],
sDSSunred[mIn], ra[mIn], dec[mIn],
reddening = False,
colors = colorDict['mlt'])
sEDNameH, magNormH, matchErrorH = selectStarSED0.findSED(listDict['H'],
sDSSunred[hIn], ra[hIn], dec[hIn],
reddening = False,
colors = colorDict['H'])
sEDNameHE, magNormHE, matchErrorHE = selectStarSED0.findSED(listDict['HE'],
sDSSunred[heIn],
ra[heIn], dec[heIn],
reddening = False,
colors = colorDict['HE'])
chunkNames = np.empty(readSize, dtype = 'S32')
chunkTypes = np.empty(readSize, dtype = 'S8')
chunkMagNorms = np.zeros(readSize)
chunkMatchErrors = np.zeros(readSize)
chunkNames[kIn] = sEDNameK
chunkTypes[kIn] = 'kurucz'
chunkMagNorms[kIn] = magNormK
chunkMatchErrors[kIn] = matchErrorK
chunkNames[mIn] = sEDNameM
chunkTypes[mIn] = 'mlt'
chunkMagNorms[mIn] = magNormM
chunkMatchErrors[mIn] = matchErrorM
chunkNames[hIn] = sEDNameH
chunkTypes[hIn] = 'H'
chunkMagNorms[hIn] = magNormH
chunkMatchErrors[hIn] = matchErrorH
chunkNames[heIn] = sEDNameHE
chunkTypes[heIn] = 'HE'
chunkMagNorms[heIn] = magNormHE
chunkMatchErrors[heIn] = matchErrorHE
lsstMagsUnred = []
for sedName, sedType, magNorm, matchError in zip(chunkNames.astype(str),
chunkTypes.astype(str),
chunkMagNorms,
chunkMatchErrors):
testSED = Sed()
testSED.setSED(listDict[sedType][positionDict[sedName]].wavelen,
flambda = listDict[sedType][positionDict[sedName]].flambda)
fluxNorm = testSED.calcFluxNorm(magNorm, imSimBand)
testSED.multiplyFluxNorm(fluxNorm)
lsstMagsUnred.append(lsstPhot.magListForSed(testSED))
#If the extinction value is negative then it will add the reddening back in
lsstMags = selectStarSED0.deReddenMags((-1.0*am), lsstMagsUnred,
lsstExtCoeffs)
distKpc = self.convDMtoKpc(DM)
ebv = am / 2.285 #From Schlafly and Finkbeiner 2011, (ApJ, 737, 103) for sdssr
ebvInf = amInf / 2.285
for line in range(0, readSize):
outFmt = '%i,%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,' +\
'%3.7f,%3.7f,%s,' +\
'%3.7f,%3.7f,' +\
'%3.7f,%3.7f,%3.7f,' +\
'%3.7f,%3.7f,%3.7f,' +\
'%3.7f,%3.7f,%3.7f,%3.7f,' +\
'%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,' +\
'%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,%3.7f,' +\
'%3.7f,%i,%3.7f,%3.7f,%3.7f\n'
if readSize == 1:
if inFits == True:
sDSS = sDSS[0]
outDat = (oID, ra[line], dec[line], gall, galb, coordX,
coordY, coordZ, chunkNames,
chunkMagNorms, chunkMatchErrors,
lsstMags[line][0], lsstMags[line][1], lsstMags[line][2],
lsstMags[line][3], lsstMags[line][4], lsstMags[line][5],
sDSS[0], sDSS[1], sDSS[2], sDSS[3],
sDSS[4], absSDSSr, pmRA, pmDec, vRad,
pml, pmb, vRadlb, vR, vPhi, vZ,
FeH, pop, distKpc, ebv, ebvInf)
else:
outDat = (oID[line], ra[line], dec[line], gall[line], galb[line], coordX[line],
coordY[line], coordZ[line], chunkNames[line],
chunkMagNorms[line], chunkMatchErrors[line],
lsstMags[line][0], lsstMags[line][1], lsstMags[line][2],
lsstMags[line][3], lsstMags[line][4], lsstMags[line][5],
sDSS[line][0], sDSS[line][1], sDSS[line][2], sDSS[line][3],
sDSS[line][4], absSDSSr[line], pmRA[line], pmDec[line], vRad[line],
pml[line], pmb[line], vRadlb[line], vR[line], vPhi[line], vZ[line],
FeH[line], pop[line], distKpc[line], ebv[line], ebvInf[line])
fOut.write(outFmt % outDat)
print('Chunk Num Done = %i out of %i' % (chunk+1, numChunks))
示例8: testAlternateBandpassesGalaxies
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testAlternateBandpassesGalaxies(self):
"""
the same as testAlternateBandpassesStars, but for galaxies
"""
obs_metadata_pointed = ObservationMetaData(mjd=50000.0,
boundType='circle',
pointingRA=0.0, pointingDec=0.0,
boundLength=10.0)
dtype = np.dtype([('galid', np.int),
('ra', np.float),
('dec', np.float),
('uTotal', np.float),
('gTotal', np.float),
('rTotal', np.float),
('iTotal', np.float),
('zTotal', np.float),
('uBulge', np.float),
('gBulge', np.float),
('rBulge', np.float),
('iBulge', np.float),
('zBulge', np.float),
('uDisk', np.float),
('gDisk', np.float),
('rDisk', np.float),
('iDisk', np.float),
('zDisk', np.float),
('uAgn', np.float),
('gAgn', np.float),
('rAgn', np.float),
('iAgn', np.float),
('zAgn', np.float),
('bulgeName', str, 200),
('bulgeNorm', np.float),
('bulgeAv', np.float),
('diskName', str, 200),
('diskNorm', np.float),
('diskAv', np.float),
('agnName', str, 200),
('agnNorm', np.float),
('redshift', np.float)])
test_cat = cartoonGalaxies(self.galaxy, obs_metadata=obs_metadata_pointed)
with lsst.utils.tests.getTempFilePath('.txt') as catName:
test_cat.write_catalog(catName)
catData = np.genfromtxt(catName, dtype=dtype, delimiter=', ')
self.assertGreater(len(catData), 0)
cartoonDir = getPackageDir('sims_photUtils')
cartoonDir = os.path.join(cartoonDir, 'tests', 'cartoonSedTestData')
sedDir = getPackageDir('sims_sed_library')
testBandpasses = {}
keys = ['u', 'g', 'r', 'i', 'z']
for kk in keys:
testBandpasses[kk] = Bandpass()
testBandpasses[kk].readThroughput(os.path.join(cartoonDir, "test_bandpass_%s.dat" % kk))
imsimBand = Bandpass()
imsimBand.imsimBandpass()
specMap = defaultSpecMap
ct = 0
for line in catData:
bulgeMagList = []
diskMagList = []
agnMagList = []
if line['bulgeName'] == 'None':
for bp in keys:
np.testing.assert_equal(line['%sBulge' % bp], np.NaN)
bulgeMagList.append(np.NaN)
else:
ct += 1
dummySed = Sed()
dummySed.readSED_flambda(os.path.join(sedDir, specMap[line['bulgeName']]))
fnorm = dummySed.calcFluxNorm(line['bulgeNorm'], imsimBand)
dummySed.multiplyFluxNorm(fnorm)
a_int, b_int = dummySed.setupCCM_ab()
dummySed.addDust(a_int, b_int, A_v=line['bulgeAv'])
dummySed.redshiftSED(line['redshift'], dimming=True)
dummySed.resampleSED(wavelen_match=testBandpasses['u'].wavelen)
for bpName in keys:
mag = dummySed.calcMag(testBandpasses[bpName])
self.assertAlmostEqual(mag, line['%sBulge' % bpName], 10)
bulgeMagList.append(mag)
if line['diskName'] == 'None':
for bp in keys:
np.assert_equal(line['%sDisk' % bp], np.NaN)
diskMagList.append(np.NaN)
else:
ct += 1
dummySed = Sed()
dummySed.readSED_flambda(os.path.join(sedDir, specMap[line['diskName']]))
fnorm = dummySed.calcFluxNorm(line['diskNorm'], imsimBand)
dummySed.multiplyFluxNorm(fnorm)
#.........这里部分代码省略.........
示例9: testFindSED
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testFindSED(self):
"""Pull SEDs from each type and make sure that each SED gets matched to itself.
Includes testing with extinction and passing in only colors."""
np.random.seed(42)
starPhot = BandpassDict.loadTotalBandpassesFromFiles(('u','g','r','i','z'),
bandpassDir = os.path.join(lsst.utils.getPackageDir('throughputs'),'sdss'),
bandpassRoot = 'sdss_')
imSimBand = Bandpass()
imSimBand.imsimBandpass()
testMatching = selectStarSED(sEDDir = self.testSpecDir, kuruczDir = self.testKDir,
mltDir = self.testMLTDir, wdDir = self.testWDDir)
testSEDList = []
testSEDList.append(testMatching.loadKuruczSEDs())
testSEDList.append(testMatching.loadmltSEDs())
testSEDListH, testSEDListHE = testMatching.loadwdSEDs()
testSEDList.append(testSEDListH)
testSEDList.append(testSEDListHE)
testSEDNames = []
testMags = []
testMagNormList = []
magNormStep = 1
for typeList in testSEDList:
if len(typeList) != 0:
typeSEDNames = []
typeMags = []
typeMagNorms = []
for testSED in typeList:
getSEDMags = Sed()
typeSEDNames.append(testSED.name)
getSEDMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
testMagNorm = np.round(np.random.uniform(20.0,22.0),magNormStep)
typeMagNorms.append(testMagNorm)
fluxNorm = getSEDMags.calcFluxNorm(testMagNorm, imSimBand)
getSEDMags.multiplyFluxNorm(fluxNorm)
typeMags.append(starPhot.magListForSed(getSEDMags))
testSEDNames.append(typeSEDNames)
testMags.append(typeMags)
testMagNormList.append(typeMagNorms)
fakeRA = np.ones(len(testSEDList[0]))
fakeDec = np.ones(len(testSEDList[0]))
#Since default bandpassDict should be SDSS ugrizy shouldn't need to specify it
#Substitute in nan values to simulate incomplete data.
for typeList, names, mags, magNorms in zip(testSEDList, testSEDNames, testMags, testMagNormList):
if len(typeList) > 2:
nanMags = np.array(mags)
nanMags[0][0] = np.nan
nanMags[0][2] = np.nan
nanMags[0][3] = np.nan
nanMags[1][1] = np.nan
testMatchingResults = testMatching.findSED(typeList, nanMags, reddening = False)
self.assertEqual(None, testMatchingResults[0][0])
self.assertEqual(names[1:], testMatchingResults[0][1:])
self.assertEqual(None, testMatchingResults[1][0])
np.testing.assert_almost_equal(magNorms[1:], testMatchingResults[1][1:],
decimal = magNormStep)
else:
testMatchingResults = testMatching.findSED(typeList, mags, reddening = False)
self.assertEqual(names, testMatchingResults[0])
np.testing.assert_almost_equal(magNorms, testMatchingResults[1], decimal = magNormStep)
#Test Null Values option
nullMags = np.array(testMags[0])
nullMags[0][0] = -99.
nullMags[0][4] = -99.
nullMags[1][0] = -99.
nullMags[1][1] = -99.
testMatchingResultsNull = testMatching.findSED(testSEDList[0], nullMags,
nullValues = -99., reddening = False)
self.assertEqual(testSEDNames[0], testMatchingResultsNull[0])
np.testing.assert_almost_equal(testMagNormList[0], testMatchingResultsNull[1],
decimal = magNormStep)
#Test Error Output
errMags = np.array((testMags[0][0], testMags[0][0], testMags[0][0], testMags[0][0]))
errMags[1,1] += 1. #Total MSE will be 2/(4 colors) = 0.5
errMags[2, 0:2] = np.nan
errMags[2, 3] += 1. #Total MSE will be 2/(2 colors) = 1.0
errMags[3, :] = None
errSED = testSEDList[0][0]
testMatchingResultsErrors = testMatching.findSED([errSED], errMags, reddening = False)
np.testing.assert_almost_equal(np.array((0.0, 0.5, 1.0)), testMatchingResultsErrors[2][0:3],
decimal = 3)
self.assertEqual(None, testMatchingResultsErrors[2][3])
#Now test what happens if we pass in a bandpassDict
testMatchingResultsNoDefault = testMatching.findSED(testSEDList[0], testMags[0],
bandpassDict = starPhot,
reddening = False)
self.assertEqual(testSEDNames[0], testMatchingResultsNoDefault[0])
np.testing.assert_almost_equal(testMagNormList[0], testMatchingResultsNoDefault[1],
decimal = magNormStep)
#Test Reddening
testRA = np.random.uniform(10,170,len(testSEDList[0]))
#.........这里部分代码省略.........
示例10: testMatchToObserved
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def testMatchToObserved(self):
"""Test that Galaxy SEDs with extinction or redshift are matched correctly"""
np.random.seed(42)
galPhot = BandpassDict.loadTotalBandpassesFromFiles()
imSimBand = Bandpass()
imSimBand.imsimBandpass()
testMatching = selectGalaxySED(galDir = self.testSpecDir)
testSEDList = testMatching.loadBC03()
testSEDNames = []
testRA = []
testDec = []
testRedshifts = []
testMagNormList = []
magNormStep = 1
extCoeffs = [1.8140, 1.4166, 0.9947, 0.7370, 0.5790, 0.4761]
testMags = []
testMagsRedshift = []
testMagsExt = []
for testSED in testSEDList:
#As a check make sure that it matches when no extinction and no redshift are present
getSEDMags = Sed()
testSEDNames.append(testSED.name)
getSEDMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
testMags.append(galPhot.magListForSed(getSEDMags))
#Check Extinction corrections
sedRA = np.random.uniform(10,170)
sedDec = np.random.uniform(10,80)
testRA.append(sedRA)
testDec.append(sedDec)
raDec = np.array((sedRA, sedDec)).reshape((2,1))
ebvVal = ebv().calculateEbv(equatorialCoordinates = raDec)
extVal = ebvVal*extCoeffs
testMagsExt.append(galPhot.magListForSed(getSEDMags) + extVal)
#Setup magnitudes for testing matching to redshifted values
getRedshiftMags = Sed()
testZ = np.round(np.random.uniform(1.1,1.3),3)
testRedshifts.append(testZ)
testMagNorm = np.round(np.random.uniform(20.0,22.0),magNormStep)
testMagNormList.append(testMagNorm)
getRedshiftMags.setSED(wavelen = testSED.wavelen, flambda = testSED.flambda)
getRedshiftMags.redshiftSED(testZ)
fluxNorm = getRedshiftMags.calcFluxNorm(testMagNorm, imSimBand)
getRedshiftMags.multiplyFluxNorm(fluxNorm)
testMagsRedshift.append(galPhot.magListForSed(getRedshiftMags))
#Will also test in passing of non-default bandpass
testNoExtNoRedshift = testMatching.matchToObserved(testSEDList, testMags, np.zeros(20),
reddening = False,
bandpassDict = galPhot)
testMatchingEbvVals = testMatching.matchToObserved(testSEDList, testMagsExt, np.zeros(20),
catRA = testRA, catDec = testDec,
reddening = True, extCoeffs = extCoeffs,
bandpassDict = galPhot)
#Substitute in nan values to simulate incomplete data and make sure magnorm works too.
testMagsRedshift[0][1] = np.nan
testMagsRedshift[0][3] = np.nan
testMagsRedshift[0][4] = np.nan
testMagsRedshift[1][1] = np.nan
testMatchingRedshift = testMatching.matchToObserved(testSEDList, testMagsRedshift, testRedshifts,
dzAcc = 3, reddening = False,
bandpassDict = galPhot)
self.assertEqual(testSEDNames, testNoExtNoRedshift[0])
self.assertEqual(testSEDNames, testMatchingEbvVals[0])
self.assertEqual(None, testMatchingRedshift[0][0])
self.assertEqual(testSEDNames[1:], testMatchingRedshift[0][1:])
self.assertEqual(None, testMatchingRedshift[1][0])
np.testing.assert_almost_equal(testMagNormList[1:], testMatchingRedshift[1][1:],
decimal = magNormStep)
#Test Match Errors
errMag = testMagsRedshift[2]
errRedshift = testRedshifts[2]
errMags = np.array((errMag, errMag, errMag, errMag))
errRedshifts = np.array((errRedshift, errRedshift, errRedshift, errRedshift))
errMags[1,1] += 1. #Total MSE will be 2/(5 colors) = 0.4
errMags[2, 0:2] = np.nan
errMags[2, 3] += 1. #Total MSE will be 2/(3 colors) = 0.667
errMags[3, :] = None
errSED = testSEDList[2]
testMatchingResultsErrors = testMatching.matchToObserved([errSED], errMags, errRedshifts,
reddening = False,
bandpassDict = galPhot,
dzAcc = 3)
np.testing.assert_almost_equal(np.array((0.0, 0.4, 2./3.)), testMatchingResultsErrors[2][0:3],
decimal = 2) #Give a little more leeway due to redshifting effects
self.assertEqual(None, testMatchingResultsErrors[2][3])
示例11: uncertaintyUnitTest
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
class uncertaintyUnitTest(unittest.TestCase):
"""
Test the calculation of photometric uncertainties
"""
def setUp(self):
starName = os.path.join(lsst.utils.getPackageDir("sims_sed_library"), defaultSpecMap["km20_5750.fits_g40_5790"])
self.starSED = Sed()
self.starSED.readSED_flambda(starName)
imsimband = Bandpass()
imsimband.imsimBandpass()
fNorm = self.starSED.calcFluxNorm(22.0, imsimband)
self.starSED.multiplyFluxNorm(fNorm)
self.totalBandpasses = []
self.hardwareBandpasses = []
componentList = ["detector.dat", "m1.dat", "m2.dat", "m3.dat", "lens1.dat", "lens2.dat", "lens3.dat"]
hardwareComponents = []
for c in componentList:
hardwareComponents.append(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", c))
self.bandpasses = ["u", "g", "r", "i", "z", "y"]
for b in self.bandpasses:
filterName = os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "filter_%s.dat" % b)
components = hardwareComponents + [filterName]
bandpassDummy = Bandpass()
bandpassDummy.readThroughputList(components)
self.hardwareBandpasses.append(bandpassDummy)
components = components + [os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "atmos.dat")]
bandpassDummy = Bandpass()
bandpassDummy.readThroughputList(components)
self.totalBandpasses.append(bandpassDummy)
def tearDown(self):
del self.starSED
del self.bandpasses
del self.hardwareBandpasses
del self.totalBandpasses
def testUncertaintyExceptions(self):
"""
Test that calcSNR_m5 raises exceptions when it needs to
"""
totalDict, hardwareDict = BandpassDict.loadBandpassesFromFiles()
magnitudes = numpy.array([22.0, 23.0, 24.0, 25.0, 26.0, 27.0])
shortMagnitudes = numpy.array([22.0])
photParams = PhotometricParameters()
shortGamma = numpy.array([1.0, 1.0])
self.assertRaises(RuntimeError, calcSNR_m5, magnitudes, totalDict.values(), shortMagnitudes, photParams)
self.assertRaises(RuntimeError, calcSNR_m5, shortMagnitudes, totalDict.values(), magnitudes, photParams)
self.assertRaises(
RuntimeError, calcSNR_m5, magnitudes, totalDict.values(), magnitudes, photParams, gamma=shortGamma
)
snr, gg = calcSNR_m5(magnitudes, totalDict.values(), magnitudes, photParams)
def testSignalToNoise(self):
"""
Test that calcSNR_m5 and calcSNR_sed give similar results
"""
defaults = LSSTdefaults()
photParams = PhotometricParameters()
totalDict, hardwareDict = BandpassDict.loadBandpassesFromFiles()
skySED = Sed()
skySED.readSED_flambda(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "darksky.dat"))
m5 = []
for filt in totalDict:
m5.append(calcM5(skySED, totalDict[filt], hardwareDict[filt], photParams, seeing=defaults.seeing(filt)))
sedDir = lsst.utils.getPackageDir("sims_sed_library")
sedDir = os.path.join(sedDir, "starSED", "kurucz")
fileNameList = os.listdir(sedDir)
numpy.random.seed(42)
offset = numpy.random.random_sample(len(fileNameList)) * 2.0
for ix, name in enumerate(fileNameList):
if ix > 100:
break
spectrum = Sed()
spectrum.readSED_flambda(os.path.join(sedDir, name))
ff = spectrum.calcFluxNorm(m5[2] - offset[ix], totalDict.values()[2])
spectrum.multiplyFluxNorm(ff)
magList = []
controlList = []
magList = []
for filt in totalDict:
controlList.append(
calcSNR_sed(
spectrum, totalDict[filt], skySED, hardwareDict[filt], photParams, defaults.seeing(filt)
)
)
magList.append(spectrum.calcMag(totalDict[filt]))
testList, gammaList = calcSNR_m5(
numpy.array(magList), numpy.array(totalDict.values()), numpy.array(m5), photParams
)
#.........这里部分代码省略.........
示例12: read_asteroids_reflectance
# 需要导入模块: from lsst.sims.photUtils.Sed import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed.Sed import multiplyFluxNorm [as 别名]
def read_asteroids_reflectance(dataDir='.'):
# Read the sun's spectrum.
sun = Sed()
sun.readSED_flambda('kurucz_sun')
# Read the asteroid reflectance spectra.
allfiles = os.listdir(dataDir)
asteroidDtype = numpy.dtype([
('wavelength', numpy.float),
('A', numpy.float),
('A_sig', numpy.float),
('B', numpy.float),
('B_sig', numpy.float),
('C', numpy.float),
('C_sig', numpy.float),
('Cb', numpy.float),
('Cb_sig', numpy.float),
('Cg', numpy.float),
('Cg_sig', numpy.float),
('Cgh', numpy.float),
('Cgh_sig', numpy.float),
('Ch', numpy.float),
('Ch_sig', numpy.float),
('D', numpy.float),
('D_sig', numpy.float),
('K', numpy.float),
('K_sig', numpy.float),
('L', numpy.float),
('L_sig', numpy.float),
('O', numpy.float),
('O_sig', numpy.float),
('Q', numpy.float),
('Q_sig', numpy.float),
('R', numpy.float),
('R_sig', numpy.float),
('S', numpy.float),
('S_sig', numpy.float),
('Sa', numpy.float),
('Sa_sig', numpy.float),
('Sq', numpy.float),
('Sq_sig', numpy.float),
('Sr', numpy.float),
('Sr_sig', numpy.float),
('Sv', numpy.float),
('Sv_sig', numpy.float),
('T', numpy.float),
('T_sig', numpy.float),
('V', numpy.float),
('V_sig', numpy.float),
('X', numpy.float),
('X_sig', numpy.float),
('Xc', numpy.float),
('Xc_sig', numpy.float),
('Xe', numpy.float),
('Xe_sig', numpy.float),
('Xk', numpy.float),
('Xk_sig', numpy.float),
])
data = numpy.loadtxt(os.path.join(dataDir,'meanspectra.tab'),
dtype = asteroidDtype)
data['wavelength'] *= 1000.0 #because spectra are in microns
wavelen_step = min(numpy.diff(data['wavelength']).min(), numpy.diff(sun.wavelen).min())
wavelen = numpy.arange(sun.wavelen[0], data['wavelength'][-1], wavelen_step)
ast_reflect = {}
for a in data.dtype.names:
if a == 'wavelength' or a[-3:] == 'sig':
continue
# Read the basic reflectance data
ast_reflect[a] = Sed(wavelen=data['wavelength'], flambda=data[a])
# And now add an extrapolation to the blue end.
# Binzel cuts off at 450nm.
condition = ((ast_reflect[a].wavelen >= 450) & (ast_reflect[a].wavelen < 700))
x = ast_reflect[a].wavelen[condition]
y = ast_reflect[a].flambda[condition]
p = numpy.polyfit(x, y, deg=2)
condition = (wavelen < 450)
flambda = numpy.zeros(len(wavelen), 'float')
interpolated = numpy.polyval(p, wavelen[condition])
flambda[condition] = [ii if ii > 0.0 else 0.0 for ii in interpolated]
condition = (wavelen >= 450)
flambda[condition] = numpy.interp(wavelen[condition], ast_reflect[a].wavelen, ast_reflect[a].flambda)
ast_reflect[a] = Sed(wavelen, flambda)
ast = {}
for a in ast_reflect:
ast[a] = sun.multiplySED(ast_reflect[a], wavelen_step=wavelen_step)
for a in ast:
name = a + '.dat'
normalizedSed = Sed(wavelen=ast[a].wavelen, flambda=ast[a].flambda)
norm = numpy.interp(500.0, normalizedSed.wavelen, normalizedSed.flambda)
normalizedSed.multiplyFluxNorm(1.0/norm)
normalizedSed.writeSED(name, print_header='21 April 2015; normalized to flambda=1 at 500nm')
#.........这里部分代码省略.........