# 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