# anton 与 aplus weight rms image的区别
Used data:
F140W and F160W plus the detection image of a2744 from aplus pipeline and the anton
F160W exptime = 68202.
F140W exptime = 28140.213257
notice that SExtractor use "weight image" in a different way, it scaled the weight map based on the sextracted image.
(*) create the weight image based on the rms image and weight image by my self.
Guess:
(*) relation of weight images between APLUS and Anton’s
anton: weight = 1/rms*2*(exposure map)**2
APLUS: weight = exposure map
_rnVarSuppFac_ = 1.38 # boost N*rn^2 by this factor for RMS images
_exVarPerSec_ = 0.02222 # add this much extra var per sec for RMS ims
dark ＝ 0.048 # WFC3 http://www.stsci.edu/hst/wfc3/documents/handbooks/currentIHB/c05_detector8.html
########################################################################################################
Var = [f*(dark+back)/g+sig**2*f**2] / t**2
weight = 1 / (Var * scale*4)
########################################################################################################
# 1. construct variance from sky, sci[], wght[]
# 2. clip zeros/negatives and infinities to values that will work w/sqrt(),
# and so that sqrt(val) = 2e30 > 1e30, => zero weight in SExtractor;
# have to work in Float64 to avoid Inf's.
# 3. take sqrt() and cast as float32
# 4. tidy header
# 5. write it out
################################################# APSIS ##################################################
skyval *= area_ratio
readVariance = Ncombed*(rn/gain)*(rn/gain)
totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime) #### here
totInstVar *= area_ratio
# maybe doing arithmetic in two steps will help conserve memory...
newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio)
# newDat[] is now variance *in counts* times maximum expTime; divide by expTime map...
newDat /= wgtfits[0].data
rmsfits[0].data = numpy.sqrt(newDat).astype(numpy.float32)
################################################# APLUS ##################################################
# about rms image
dark = 0.00008 #WZ, UVIS dark rate
# skyval *= area_ratio #WZ redundant
readVariance = Ncombed*(rn/gain)*(rn/gain)
totInstVar = readVariance + dark * exptime #WZ #### here
newDat = ((skyval+ nullDat)/gain + totInstVar) * (exptime * area_ratio) # This is an array
newDat /= wgtfits[0].data
## now fix up problem values...
newDat = numpy.where(numpy.logical_or(numpy.greater_equal(newDat,1e38),\
numpy.less_equal(newDat,0.)),4e38,newDat) #WZ Mar 2013
# about the weight image for detection
def mkIrafIm(self):
""" make the detection image. This method creates the detection
Image using the CHI^2 algorithm.
"""
def mkChi2Im(self):
"""
Make the detection image using Numeric package, not IRAF.
This method creates the detection Image using the CHI^2 algorithm.
"""
def mkInvarIm(self):
"""
Make the detection image using Numeric package, not IRAF.
This method creates a detection image using the inverse variance
algorithm.
"""
def mkWghtIm(self, doMedFilt=0, medfiltsize=5):
""" Makes the weight image to be used by Sextractor when finding
sources in the (inverse variance-weighted) detection image.
The method for constructing the weight image follows analogously
to that used in making the (inverse-variance) detection image.
"""
for im in self.statsList:
noise = float(im[2]) # [1] is background, [2] is noise based on the sextractor.
x = (ffi[0].data)/noise/noise
wgtFF[0].data += x.astype(numpy.float32)
weightNorm += 1.0/noise/noise
wgtFF[0].data = wgtFF[0].data/weightNorm
#########################################################################################################
####################################################
IN APSIS
####################################################
def makeRmsImage(self):
"""turns the drizzled weight images into RMS images that
SExtractor will use.
"""
self.logfile.write("starting makeRmsImage . . .")
if not self.weightImageList:
errtxt="No Weight Images present."
self.errorList.append((self.modName,errtxt))
raise Exception, errtxt
# reset rms image list
while self.rmsImageList:
del self.rmsImageList[0]
curdir = os.getcwd()
os.chdir(self.obsFits)
for im_wgt in self.weightImageList:
im_sci = im_wgt[:-12]+'.fits'
if im_sci not in self.sciImageList:
errtxt = 'makeRmsImage: '+im_sci+' not in sciImageList[]!'
self.errorList.append((self.modName,errtxt))
self.logfile.write(errtxt)
try:
wgtfits = pyfits.open(im_wgt)
scifits = pyfits.open(im_sci)
except:
self.errorList.append((self.modName,"Cannot make a FITS object out of file "+im_wgt))
raise Exception,"Cannot make a FITS object out of file "+im_wgt
if len(wgtfits) > 1 or len(scifits) > 1:
self.errorList.append((self.modName,"image file is not simple fits."))
raise Exception,"image file is not simple fits."
# build rms image name and open as a new file.
rmsfile = im_wgt.split("_drz")[0]+'_RMS.fits'
self.rmsImageList.append(rmsfile)
self.outputList[rmsfile] = [im_wgt]
# make new fits obj and copy WGT/SCI hdr/data to RMS image initially
rmsfitsobj = pyfits.HDUList()
rmsfitsobj.append(pyfits.PrimaryHDU())
try:
del rmsfitsobj[0].header.ascard["EXTEND"]
except KeyError:
pass
rmsfitsobj[0].header = wgtfits[0].header
rmsfitsobj[0].data = scifits[0].data
if os.path.isfile(rmsfile):
os.remove(rmsfile)
rmsfitsobj.writeto(rmsfile)
del rmsfitsobj
# reopen the rms image for editing.
rmsfits = pyfits.open(rmsfile,'update')
# ratio of default to specified output scales
area_ratio = (self.asecpix / self.origscale)**2
if abs(1-area_ratio) < 1e-4: area_ratio = 1
self.logfile.write('Using area_ratio = %.6f in makeRmsImage' %(area_ratio))
skyval = scifits[0].header.get('ALIGNSKY') # this rescaled below
exptime = scifits[0].header.get('EXPTIME')
Ncombed = scifits[0].header.get('NCOMBINE')
if Ncombed == None:
errtxt='Error: NCOMBINE not in '+im_sci+' header. _sumSubSkyNcombine() not run?'
self.logfile.write(errtxt)
self.errorList.append((self.modName,errtxt))
raise Exception,errtxt
gain,rn = pyblot._gain_rn(scifits, self.logfile, ext=0)
# if not told to use header gain, then use 1.0 (data in electrons)
if not self.hdrGain:
gain = 1.0
self.logfile.write(im_sci+": gain,rn = "+str(gain)+","+str(rn)+\
" NCOMBINE = "+str(Ncombed)+" EXPTIME = "+str(exptime))
if not exptime:
raise Exception,"No EXPTIME in "+im_sci
if (skyval == None or skyval < 1):
warntxt = 'WARNING: found ALIGNSKY of '+str(skyval)+' in '+im_sci+\
' : RMS image may be in WRONG!'
self.logfile.write(warntxt)
self.errorList.append((self.modName,warntxt))
del warntxt
if skyval == None: skyval=0
skyval *= area_ratio
# 1. construct variance from sky, sci[], wght[]
# 2. clip zeros/negatives and infinities to values that will work w/sqrt(),
# and so that sqrt(val) = 2e30 > 1e30, => zero weight in SExtractor;
# have to work in Float64 to avoid Inf's.
# 3. take sqrt() and cast as float32
# 4. tidy header
# 5. write it out
readVariance = Ncombed*(rn/gain)*(rn/gain)
self.logfile.write("total read variance = "+str(readVariance)+" for "+im_sci)
if self.suppInstVar:
# supplement factor for reference bais subtraction, etc
# extra var per sec for dark subtraction, repaired cosmic rays, etc.
totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
self.logfile.write("adjusted instrumental variance = "+str(totInstVar)+" for "+im_sci)
else:
totInstVar = readVariance
totInstVar *= area_ratio
# maybe doing arithmetic in two steps will help conserve memory...
newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio)
# newDat[] is now variance *in counts* times maximum expTime; divide by expTime map...
newDat /= wgtfits[0].data
scifits.close()
wgtfits.close()
del scifits, wgtfits, im_wgt, im_sci, readVariance, totInstVar, area_ratio
## now fix up problem values...
newDat = numpy.where(numpy.logical_or(numpy.greater_equal(newDat,1e60),\
numpy.less_equal(newDat,0.)),4e60,newDat)
rmsfits[0].data = numpy.sqrt(newDat).astype(numpy.float32)
# a few token updates to the header, then write it out
rmsfits[0].header.update('FILENAME',rmsfile)
rmsfits[0].header.update('FILETYPE','RMS')
rmsfits.close()
self.logfile.write('Made rms image '+rmsfile)
del newDat, rmsfile, rmsfits
os.chdir(curdir)
return
####################################################
IN APLUS
####################################################
def makeRmsImage(self):
"""turns the drizzled weight images into RMS images that
SExtractor will use.
"""
sn0 = 20. # WZ Threshold for bright pixels
self.logfile.write("starting makeRmsImage . . .")
if not self.weightImageList:
errtxt="No Weight Images present."
self.errorList.append((self.modName,errtxt))
raise Exception, errtxt
# reset rms image list
while self.rmsImageList:
del self.rmsImageList[0]
curdir = os.getcwd()
os.chdir(self.obsFits)
for im_wgt in self.weightImageList:
#im_sci = im_wgt[:-9]+'.fits' #WZ
#pdb.set_trace()
im_sci = im_wgt[:-12]+'.fits'
if im_sci not in self.sciImageList:
errtxt = 'makeRmsImage: '+im_sci+' not in sciImageList[]!'
self.errorList.append((self.modName,errtxt))
self.logfile.write(errtxt)
try:
wgtfits = pyfits.open(im_wgt)
scifits = pyfits.open(im_sci)
except:
self.errorList.append((self.modName,"Cannot make a FITS object out of file "+im_wgt))
raise Exception,"Cannot make a FITS object out of file "+im_wgt
if len(wgtfits) > 1 or len(scifits) > 1:
self.errorList.append((self.modName,"image file is not simple fits."))
raise Exception,"image file is not simple fits."
# build rms image name and open as a new file.
rmsfile = im_wgt.split("_drz")[0]+'_RMS.fits'
self.rmsImageList.append(rmsfile)
self.outputList[rmsfile] = [im_wgt]
# make new fits obj and copy WGT/SCI hdr/data to RMS image initially
rmsfitsobj = pyfits.HDUList()
rmsfitsobj.append(pyfits.PrimaryHDU())
try:
del rmsfitsobj[0].header.ascard["EXTEND"]
except KeyError:
pass
rmsfitsobj[0].header = wgtfits[0].header
rmsfitsobj[0].data = scifits[0].data
if os.path.isfile(rmsfile):
os.remove(rmsfile)
rmsfitsobj.writeto(rmsfile)
del rmsfitsobj
# reopen the rms image for editing.
rmsfits = pyfits.open(rmsfile,'update')
# ratio of default to specified output scales
area_ratio = (self.asecpix / self.origscale)**2
if abs(1-area_ratio) < 1e-4: area_ratio = 1
self.logfile.write('Using area_ratio = %.6f in makeRmsImage' %(area_ratio))
skyval = scifits[0].header.get('ALIGNSKY') # this rescaled below
exptime = scifits[0].header.get('EXPTIME')
Ncombed = scifits[0].header.get('NCOMBINE')
if Ncombed == None:
errtxt='Error: NCOMBINE not in '+im_sci+' header. _sumSubSkyNcombine() not run?'
self.logfile.write(errtxt)
self.errorList.append((self.modName,errtxt))
raise Exception,errtxt
gain,rn = pyblot._gain_rn(scifits, self.logfile, ext=0)
# if not told to use header gain, then use 1.0 (data in electrons)
if not self.hdrGain:
gain = 1.54 # WZ was 1.0 # XX really?
self.logfile.write(im_sci+": gain,rn = "+str(gain)+","+str(rn)+\
" NCOMBINE = "+str(Ncombed)+" EXPTIME = "+str(exptime))
if not exptime:
raise Exception,"No EXPTIME in "+im_sci
if (skyval == None or skyval < 1):
warntxt = 'WARNING: found ALIGNSKY of '+str(skyval)+' in '+im_sci+\
' : RMS image may be in WRONG!'
self.logfile.write(warntxt)
self.errorList.append((self.modName,warntxt))
del warntxt
if skyval == None: skyval=0
# skyval *= area_ratio #WZ redundant
# 1. construct variance from sky, sci[], wght[]
# 2. clip zeros/negatives and infinities to values that will work w/sqrt(),
# and so that sqrt(val) = 2e30 > 1e30, => zero weight in SExtractor;
# have to work in Float64 to avoid Inf's.
# 3. take sqrt() and cast as float32
# 4. tidy header
# 5. write it out
readVariance = Ncombed*(rn/gain)*(rn/gain)
self.logfile.write("total read variance = "+str(readVariance)+" for "+im_sci)
if self.suppInstVar:
# supplement factor for reference bais subtraction, etc
# extra var per sec for dark subtraction, repaired cosmic rays, etc.
dark = 0.00008 #WZ, UVIS dark rate
totInstVar = readVariance + dark * exptime #WZ
#totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
self.logfile.write("adjusted instrumental variance = "+str(totInstVar)+" for "+im_sci)
else:
totInstVar = readVariance
# totInstVar *= area_ratio # WZ redundant
# maybe doing arithmetic in two steps will help conserve memory...
# newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio
factor = 0.
nullDat = factor * scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
nullDat = scifits[0].data.astype(numpy.float64) - scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
# Set up an array that does not include source counts WZ
# newDat = ((skyval+ nullDat)/gain + totInstVar) * area_ratio
newDat = ((skyval+ nullDat)/gain + totInstVar) * (exptime * area_ratio) # This is an array
# newDat[] is now variance *in counts* times maximum expTime; To be divided by expTime map...
#sn = newDat/rmsfits[0].data #WZ Mark bright pixels
#idx = numpy.where(numpy.logical_and(numpy.greater(sn,sn0),numpy.less(wgtfits[0].data,0.)))
#wgtfits[0].data[idx] = numpy.abs(wgtfits[0].data[idx])
#for i in range(len(wgtfits[0].data)):
# for j in range(len(wgtfits[0].data[i])):
# if (wgtfits[0].data[i,j]< 1e-10):
# wgtfits[0].data[i,j]= 1e-10
# idx = numpy.where(numpy.less(wgtfits[0].data,1e-10))
# wgtfits[0].data[idx] = 1e-10
# wgtfits[0].data = wgtfits[0].data + 1e-10
newDat /= wgtfits[0].data
## now fix up problem values...
newDat = numpy.where(numpy.logical_or(numpy.greater_equal(newDat,1e38),\
numpy.less_equal(newDat,0.)),4e38,newDat) #WZ Mar 2013
scifits.close()
wgtfits.close()
rmsfits[0].data = numpy.sqrt(newDat).astype(numpy.float32)
del scifits, wgtfits, im_wgt, im_sci, readVariance, totInstVar, area_ratio
# a few token updates to the header, then write it out
rmsfits[0].header.update('FILENAME',rmsfile)
rmsfits[0].header.update('FILETYPE','RMS')
rmsfits.close()
self.logfile.write('Made rms image '+rmsfile)
del newDat, rmsfile, rmsfits
os.chdir(curdir)
return
#########################################################################################################
different:
#########################################################################################################
if not self.hdrGain:
gain = 1.54 # WZ was 1.0 # XX really?
#skyval *= area_ratio #WZ #<-------------
readVariance = Ncombed*(rn/gain)*(rn/gain) # same
if self.suppInstVar:
# supplement factor for reference bais subtraction, etc
# extra var per sec for dark subtraction, repaired cosmic rays, etc.
dark = 0.00008 #WZ, UVIS dark rate #<-------------
totInstVar = readVariance + dark * exptime #WZ #<-------------
#totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
factor = 0.
nullDat = factor * scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
nullDat = scifits[0].data.astype(numpy.float64) - scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
# Set up an array that does not include source counts WZ
# newDat = ((skyval+ nullDat)/gain + totInstVar) * area_ratio
newDat = ((skyval+ nullDat)/gain + totInstVar) * (exptime * area_ratio) # This is an array #<-------------
#########################################################################################################
In apsis:
#########################################################################################################
_rnVarSuppFac_ = 1.38 # boost N*rn^2 by this factor for RMS images
_exVarPerSec_ = 0.02222 # add this much extra var per sec for RMS ims
totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio)
# newDat[] is now variance *in counts* times maximum expTime; divide by expTime map...
newDat /= wgtfits[0].data
different:
if not self.hdrGain:
gain = 1.54 # WZ was 1.0 # XX really?
#skyval *= area_ratio #WZ #<-------------
readVariance = Ncombed*(rn/gain)*(rn/gain) # same
if self.suppInstVar:
# supplement factor for reference bais subtraction, etc
# extra var per sec for dark subtraction, repaired cosmic rays, etc.
dark = 0.00008 #WZ, UVIS dark rate #<-------------
totInstVar = readVariance + dark * exptime #WZ #<-------------
#totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
factor = 0.
nullDat = factor * scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
nullDat = scifits[0].data.astype(numpy.float64) - scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
# Set up an array that does not include source counts WZ
# newDat = ((skyval+ nullDat)/gain + totInstVar) * area_ratio
newDat = ((skyval+ nullDat)/gain + totInstVar) * (exptime * area_ratio) # This is an array #<-------------
In apsis:
_rnVarSuppFac_ = 1.38 # boost N*rn^2 by this factor for RMS images
_exVarPerSec_ = 0.02222 # add this much extra var per sec for RMS ims
totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio)
# newDat[] is now variance *in counts* times maximum expTime; divide by expTime map...
newDat /= wgtfits[0].data
APLUS
APSIS
aplus (CombDither.py 1641):
def makeRmsImage(self):
"""turns the drizzled weight images into RMS images that
SExtractor will use.
"""
sn0 = 20. # WZ Threshold for bright pixels
self.logfile.write("starting makeRmsImage . . .")
if not self.weightImageList:
errtxt="No Weight Images present."
self.errorList.append((self.modName,errtxt))
raise Exception, errtxt
# reset rms image list
while self.rmsImageList:
del self.rmsImageList[0]
curdir = os.getcwd()
os.chdir(self.obsFits)
for im_wgt in self.weightImageList:
#im_sci = im_wgt[:-9]+'.fits' #WZ
#pdb.set_trace()
im_sci = im_wgt[:-12]+'.fits'
if im_sci not in self.sciImageList:
errtxt = 'makeRmsImage: '+im_sci+' not in sciImageList[]!'
self.errorList.append((self.modName,errtxt))
self.logfile.write(errtxt)
try:
wgtfits = pyfits.open(im_wgt)
scifits = pyfits.open(im_sci)
except:
self.errorList.append((self.modName,"Cannot make a FITS object out of file "+im_wgt))
raise Exception,"Cannot make a FITS object out of file "+im_wgt
if len(wgtfits) > 1 or len(scifits) > 1:
self.errorList.append((self.modName,"image file is not simple fits."))
raise Exception,"image file is not simple fits."
# build rms image name and open as a new file.
rmsfile = im_wgt.split("_drz")[0]+'_RMS.fits'
self.rmsImageList.append(rmsfile)
self.outputList[rmsfile] = [im_wgt]
# make new fits obj and copy WGT/SCI hdr/data to RMS image initially
rmsfitsobj = pyfits.HDUList()
rmsfitsobj.append(pyfits.PrimaryHDU())
try:
del rmsfitsobj[0].header.ascard["EXTEND"]
except KeyError:
pass
rmsfitsobj[0].header = wgtfits[0].header
rmsfitsobj[0].data = scifits[0].data
if os.path.isfile(rmsfile):
os.remove(rmsfile)
rmsfitsobj.writeto(rmsfile)
del rmsfitsobj
# reopen the rms image for editing.
rmsfits = pyfits.open(rmsfile,'update')
# ratio of default to specified output scales
area_ratio = (self.asecpix / self.origscale)**2
if abs(1-area_ratio) < 1e-4: area_ratio = 1
self.logfile.write('Using area_ratio = %.6f in makeRmsImage' %(area_ratio))
skyval = scifits[0].header.get('ALIGNSKY') # this rescaled below
exptime = scifits[0].header.get('EXPTIME')
Ncombed = scifits[0].header.get('NCOMBINE')
if Ncombed == None:
errtxt='Error: NCOMBINE not in '+im_sci+' header. _sumSubSkyNcombine() not run?'
self.logfile.write(errtxt)
self.errorList.append((self.modName,errtxt))
raise Exception,errtxt
gain,rn = pyblot._gain_rn(scifits, self.logfile, ext=0)
# if not told to use header gain, then use 1.0 (data in electrons)
if not self.hdrGain:
gain = 1.54 # WZ was 1.0
self.logfile.write(im_sci+": gain,rn = "+str(gain)+","+str(rn)+\
" NCOMBINE = "+str(Ncombed)+" EXPTIME = "+str(exptime))
if not exptime:
raise Exception,"No EXPTIME in "+im_sci
if (skyval == None or skyval < 1):
warntxt = 'WARNING: found ALIGNSKY of '+str(skyval)+' in '+im_sci+\
' : RMS image may be in WRONG!'
self.logfile.write(warntxt)
self.errorList.append((self.modName,warntxt))
del warntxt
if skyval == None: skyval=0
# skyval *= area_ratio #WZ redundant
# 1. construct variance from sky, sci[], wght[]
# 2. clip zeros/negatives and infinities to values that will work w/sqrt(),
# and so that sqrt(val) = 2e30 > 1e30, => zero weight in SExtractor;
# have to work in Float64 to avoid Inf's.
# 3. take sqrt() and cast as float32
# 4. tidy header
# 5. write it out
readVariance = Ncombed*(rn/gain)*(rn/gain)
self.logfile.write("total read variance = "+str(readVariance)+" for "+im_sci)
if self.suppInstVar:
# supplement factor for reference bais subtraction, etc
# extra var per sec for dark subtraction, repaired cosmic rays, etc.
dark = 0.00008 #WZ, UVIS dark rate
totInstVar = readVariance + dark * exptime #WZ
#totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
self.logfile.write("adjusted instrumental variance = "+str(totInstVar)+" for "+im_sci)
else:
totInstVar = readVariance
# totInstVar *= area_ratio # WZ redundant
# maybe doing arithmetic in two steps will help conserve memory...
# newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio
factor = 0.
nullDat = factor * scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
nullDat = scifits[0].data.astype(numpy.float64) - scifits[0].data.astype(numpy.float64)
mean=numpy.mean(nullDat)
std=numpy.std(nullDat)
min = numpy.min(nullDat)
max = numpy.max(nullDat)
# Set up an array that does not include source counts WZ
# newDat = ((skyval+ nullDat)/gain + totInstVar) * area_ratio
newDat = ((skyval+ nullDat)/gain + totInstVar) * (exptime * area_ratio) # This is an array
# newDat[] is now variance *in counts* times maximum expTime; To be divided by expTime map...
#sn = newDat/rmsfits[0].data #WZ Mark bright pixels
#idx = numpy.where(numpy.logical_and(numpy.greater(sn,sn0),numpy.less(wgtfits[0].data,0.)))
#wgtfits[0].data[idx] = numpy.abs(wgtfits[0].data[idx])
#for i in range(len(wgtfits[0].data)):
# for j in range(len(wgtfits[0].data[i])):
# if (wgtfits[0].data[i,j]< 1e-10):
# wgtfits[0].data[i,j]= 1e-10
# idx = numpy.where(numpy.less(wgtfits[0].data,1e-10))
# wgtfits[0].data[idx] = 1e-10
# wgtfits[0].data = wgtfits[0].data + 1e-10
newDat /= wgtfits[0].data
## now fix up problem values...
newDat = numpy.where(numpy.logical_or(numpy.greater_equal(newDat,1e38),\
numpy.less_equal(newDat,0.)),4e38,newDat) #WZ Mar 2013
scifits.close()
wgtfits.close()
rmsfits[0].data = numpy.sqrt(newDat).astype(numpy.float32)
del scifits, wgtfits, im_wgt, im_sci, readVariance, totInstVar, area_ratio
# a few token updates to the header, then write it out
rmsfits[0].header.update('FILENAME',rmsfile)
rmsfits[0].header.update('FILETYPE','RMS')
rmsfits.close()
self.logfile.write('Made rms image '+rmsfile)
del newDat, rmsfile, rmsfits
os.chdir(curdir)
return
apsis (CombDither.py 1641):
def makeRmsImage(self):
"""turns the drizzled weight images into RMS images that
SExtractor will use.
"""
self.logfile.write("starting makeRmsImage . . .")
if not self.weightImageList:
errtxt="No Weight Images present."
self.errorList.append((self.modName,errtxt))
raise Exception, errtxt
# reset rms image list
while self.rmsImageList:
del self.rmsImageList[0]
curdir = os.getcwd()
os.chdir(self.obsFits)
for im_wgt in self.weightImageList:
im_sci = im_wgt[:-12]+'.fits'
if im_sci not in self.sciImageList:
errtxt = 'makeRmsImage: '+im_sci+' not in sciImageList[]!'
self.errorList.append((self.modName,errtxt))
self.logfile.write(errtxt)
try:
wgtfits = pyfits.open(im_wgt)
scifits = pyfits.open(im_sci)
except:
self.errorList.append((self.modName,"Cannot make a FITS object out of file "+im_wgt))
raise Exception,"Cannot make a FITS object out of file "+im_wgt
if len(wgtfits) > 1 or len(scifits) > 1:
self.errorList.append((self.modName,"image file is not simple fits."))
raise Exception,"image file is not simple fits."
# build rms image name and open as a new file.
rmsfile = im_wgt.split("_drz")[0]+'_RMS.fits'
self.rmsImageList.append(rmsfile)
self.outputList[rmsfile] = [im_wgt]
# make new fits obj and copy WGT/SCI hdr/data to RMS image initially
rmsfitsobj = pyfits.HDUList()
rmsfitsobj.append(pyfits.PrimaryHDU())
try:
del rmsfitsobj[0].header.ascard["EXTEND"]
except KeyError:
pass
rmsfitsobj[0].header = wgtfits[0].header
rmsfitsobj[0].data = scifits[0].data
if os.path.isfile(rmsfile):
os.remove(rmsfile)
rmsfitsobj.writeto(rmsfile)
del rmsfitsobj
# reopen the rms image for editing.
rmsfits = pyfits.open(rmsfile,'update')
# ratio of default to specified output scales
area_ratio = (self.asecpix / self.origscale)**2
if abs(1-area_ratio) < 1e-4: area_ratio = 1
self.logfile.write('Using area_ratio = %.6f in makeRmsImage' %(area_ratio))
skyval = scifits[0].header.get('ALIGNSKY') # this rescaled below
exptime = scifits[0].header.get('EXPTIME')
Ncombed = scifits[0].header.get('NCOMBINE')
if Ncombed == None:
errtxt='Error: NCOMBINE not in '+im_sci+' header. _sumSubSkyNcombine() not run?'
self.logfile.write(errtxt)
self.errorList.append((self.modName,errtxt))
raise Exception,errtxt
gain,rn = pyblot._gain_rn(scifits, self.logfile, ext=0)
# if not told to use header gain, then use 1.0 (data in electrons)
if not self.hdrGain:
gain = 1.0
self.logfile.write(im_sci+": gain,rn = "+str(gain)+","+str(rn)+\
" NCOMBINE = "+str(Ncombed)+" EXPTIME = "+str(exptime))
if not exptime:
raise Exception,"No EXPTIME in "+im_sci
if (skyval == None or skyval < 1):
warntxt = 'WARNING: found ALIGNSKY of '+str(skyval)+' in '+im_sci+\
' : RMS image may be in WRONG!'
self.logfile.write(warntxt)
self.errorList.append((self.modName,warntxt))
del warntxt
if skyval == None: skyval=0
skyval *= area_ratio
# 1. construct variance from sky, sci[], wght[]
# 2. clip zeros/negatives and infinities to values that will work w/sqrt(),
# and so that sqrt(val) = 2e30 > 1e30, => zero weight in SExtractor;
# have to work in Float64 to avoid Inf's.
# 3. take sqrt() and cast as float32
# 4. tidy header
# 5. write it out
readVariance = Ncombed*(rn/gain)*(rn/gain)
self.logfile.write("total read variance = "+str(readVariance)+" for "+im_sci)
if self.suppInstVar:
# supplement factor for reference bais subtraction, etc
# extra var per sec for dark subtraction, repaired cosmic rays, etc.
totInstVar = (_rnVarSuppFac_ * readVariance) + (_exVarPerSec_ * exptime)
self.logfile.write("adjusted instrumental variance = "+str(totInstVar)+" for "+im_sci)
else:
totInstVar = readVariance
totInstVar *= area_ratio
# maybe doing arithmetic in two steps will help conserve memory...
newDat = ((skyval + scifits[0].data.astype(numpy.float64))/gain + totInstVar) * (exptime * area_ratio)
# newDat[] is now variance *in counts* times maximum expTime; divide by expTime map...
newDat /= wgtfits[0].data
scifits.close()
wgtfits.close()
del scifits, wgtfits, im_wgt, im_sci, readVariance, totInstVar, area_ratio
## now fix up problem values...
newDat = numpy.where(numpy.logical_or(numpy.greater_equal(newDat,1e60),\
numpy.less_equal(newDat,0.)),4e60,newDat)
rmsfits[0].data = numpy.sqrt(newDat).astype(numpy.float32)
# a few token updates to the header, then write it out
rmsfits[0].header.update('FILENAME',rmsfile)
rmsfits[0].header.update('FILETYPE','RMS')
rmsfits.close()
self.logfile.write('Made rms image '+rmsfile)
del newDat, rmsfile, rmsfits
os.chdir(curdir)
return
没有评论:
发表评论