# Measures the rms of an image. Designed for single-channel data cubes, i.e. for
# measuring rms across the spatial bandpass. Data MUST have only two axes !
# MinMed style approach. At each y pixel, splits the "spectrum" into five boxes, measures the
# rms in each, then finds the minimum and median of the boxes.
# Also reports the rms of a specified chunk of the image.
import math
import numpy
import pyfits
import os
# Changes to the directory where the script is located
abspath = os.path.abspath(__file__)
dname = os.path.dirname(abspath)
os.chdir(dname)
infile = 'VirgoCalTrim.fits'
outfile = str(infile.split('.fits')[0]+str('_SimMinMed.fits'))
# Specify number of boxes
n = 5
FitsFile = pyfits.open(infile)
image= FitsFile[0].data
header=FitsFile[0].header
#sizez = image.shape[0]
sizey = image.shape[0]
sizex = image.shape[1]
print sizex,sizey#,sizez
for yp in range(0,sizey):
#for zp in range(0,sizez):
#print yp,zp
vlist = []
medlist = []
bpass = numpy.array(image[yp,:])
# No. boxes examined so far
nb = 0
# No. values in each box
lb = int(float(len(bpass))/float(n))
for i in range(0,n):
for j in range(nb*lb,(nb*lb)+lb):
vlist.append(bpass[j])
# Remove nans from vlist. Do this now, rather than removing from the bandpass, because
# this makes it easier to iterate over the (fixed) length of the bandpass (rather than
# what would be the variable length of vlist).
# Have to convert vlist into numpy array first
vlist = numpy.array(vlist)
vlist = vlist[numpy.logical_not(numpy.isnan(vlist))]
#med = numpy.std(vlist)
med = numpy.sqrt(numpy.mean(numpy.square(vlist)))
medlist.append(med)
vlist = []
nb = nb + 1
minmed = numpy.min(medlist)
medmed = numpy.median(medlist)
#bpass = bpass - minmed
#image[zp,yp,:] = bpass
#FitsFile = pyfits.writeto(outfile,image,header=header)
print 'Spatial bandpass rms measurements'
print 'Minimum of rms :',minmed
print 'Median of rms :,',medmed
# Now get the rms of a specific chunk of the image
chunk = numpy.array(image[50:100,30:130])
crms = numpy.std(chunk)
print 'Image chunk rms :',crms