Sextractor【2】MAP_WEIGHT and HOW it works.

The parameters are only breifly mentioned in the manual of Sextractor. However, you should be careful if you really want to change the paramters.

(1) Choose a right type of WEIGHT_TYPE
MAP_RMS is the most safe way to use weight image. It is the map showing the rms (you can call it error or variance.). Sextractor will select pixels with value higher than threshold*rms as a source (the value I mentioned here is background subtracted and smoothed). The error are also calculated based on the counts of object and the rms value.
MAP_WEIGHT is not the same as define in Multidrizzle or somewhere else. It can be a flat map, a exposure map. The Sextractor will convert it to rms map whatever your input is. What is most important and often missed is that it will scale the weight image according to the rms measure from the background. I think it just multiplies the rms value which is show in the output, 78.5964 in the following output.
(M+D) Background: 0          RMS: 78.5964    / Threshold: 78.5964 
So the weight image is used to tell which part of the image is more important than other parts. The arbitrary value is meaningless. The program will scale it and only use the related values!
MAP_VAR In the manual, "The data are converted to variance units (by defi- nition variance ∝ 1/weight), and scaled as for MAP VAR" However,  I have not used it! 
BACKGROUND The program will created a rms map based on the background measure my the program. This method are really not good. If you want to do real science, you should not choose it.
None Then every part of the image will be treated as the same weight, which means the regions which are not observed and are set to value 0 will also be treated as other region. The errors of objects in the corner will also be under-estimated. Notice that even you have set the MAP_WEIGHT, the overall errors are also probably under-estimated, because the rms is calculate from the image and the information of read noises and dark noises are missed.

(2) Detect and measure with different weight map.
People usually know that in dual model, we can define the MAP_TYPE to the detection image and measurement image separately. We can also set this in one image run.
"

All the weighting options listed in §7.1 can be applied separately to detection and measurement images (§3), — even if some combinations may not always make sense. For instance, the following set of configuration lines:
WEIGHT_IMAGE rms.fits,weight.fits 
WEIGHT_TYPE MAP_RMS,MAP_WEIGHT
will load the FITS file rms.fits and use it as an RMS map for adjusting the detection threshold and CLEANing, while the weight.fits weight map will only be used for scaling the error estimates on measurements. This can be done in single- as well as in dual-image mode (§3). WEIGHT IMAGEs can be ignored for BACKGROUND WEIGHT TYPEs. 
"

(3) Created you own weight and rms image.
The details will come later, you can check my other blog and also give me an email.

   

Sextractor【1】The size measurement

Usually, we need to run sextractor for several bands with one detection image, which is called dual mode. But when you check the FWHM which is the half light radius for faint galaxies, you may find something different. How do these values come from? Does Sextractor use the same aperture for all bands?

The FWHM of one object in C1226 ( one cluster in CLASH)

The object is not detected in ACS/uvis bands and F775W F814W band, thus the size is only the largest radius of the source region which can be checked in segmentation map.
The object are clearly detected in F105W, so the size measured in F105W could be the most accuracy. The size in other bands are smaller due to the weak detection in these bands.

size measured by myself

Spectrum of this object

images of 16 bands

# photometry for this object

#a1423 nir  833    x    y         ra        dec
f225w           f275w           f336w           f390w        
f435w           f475w           f606w           f625w        
f775w           f814w          f850lp           f105w        
f110w           f125w           f140w           f160w           redshift

c1226   1 1388 2543 2764 186.754096  33.545302
99.0000 99.0000  99.0000 99.0000  29.2562 2.1399  99.0000 99.0000
27.1057 0.1604  27.2714 0.1534  27.2058 0.0910  26.8436 0.1355
27.3479 0.2488  27.2449 0.0892  27.6370 0.4182  26.9311 0.1483
26.5273 0.0706  26.1972 0.0958  26.4933 0.0874  26.9528 0.1281   1.7550

You can see that while the size is as large as the segmentation map, the object will be either not detected or deblend for bad pixels.

size 为何不同
有些波段的流量偏低，甚至无探测，导致sex测量的size不准。没有探测的时候，size将会设置成定值，即segmentation的大小。



Here I add another example:

# photometry for this object

#a1423 nir  833    x    y         ra        dec
f225w           f275w           f336w           f390w        
f435w           f475w           f606w           f625w        
f775w           f814w          f850lp           f105w        
f110w           f125w           f140w           f160w           redshift

c1226   1 2101 2303 2233 186.759314  33.535719  99.0000 99.0000  

30.7607 8.0779  27.8950 0.5042  28.3753 0.4528  27.1087 0.1303  

27.2161 0.1204  27.4034 0.0910  27.6560 0.2317  27.5740 0.2553  

27.9762 0.1456  27.2453 0.2345  27.3148 0.1553  27.2221 0.1081  

26.9002 0.1156  27.0465 0.1132  27.7699 0.2222   1.7550

Sextractor【2.2】RMS image and WEIGHT image problem in Hubble Image process

#!/usr/bin/env python
# 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