Adding poisson noise to a fits file

Question

I'm trying to add incremental amounts of poisson noise to a .fits file. I know how to do it for a regular filetype, but I can't seem to read in the fits and then add in the poisson noise. Does anybody know how to go about doing this?

Here's the code. Most of it isn't particularly relevant.

s=str(raw_input("filter name: "))
t=str(raw_input("sci or wht: "))
poisson = str(raw_input("Poisson noise amount: "))
for i in range(0,len(ra_new)):
    ra_new2=cat['ra'][z2&lmass2&ra2&dec2][i]
    dec_new2=cat['dec'][z2&lmass2&ra2&dec2][i]
    id_new=cat['id'][z2&lmass2&ra2&dec2][i]
    target_pixel_x = ((ra_new2-ra_ref)/(pixel_size_x))+reference_pixel_x       
    target_pixel_y = ((dec_new2-dec_ref)/(pixel_size_y))+reference_pixel_y   
    fig = plt.figure(figsize=(5.,5.))
    timage=img[target_pixel_y-65:target_pixel_y+65,target_pixel_x-65:target_pixel_x+65]
    plt.imshow(img[target_pixel_y-65:target_pixel_y+65,target_pixel_x-65:target_pixel_x+65], vmin=-0.01, vmax=0.1, cmap='Greys')

galimage =  pf.writeto(t+'PHOTO'+s+str(i)+'.fits',timage,clobber=True,header=hdr)
        imagea = (scipy.misc.imread(galimage)).astype(float)
        poissonNoise = numpy.random.poisson(poisson,imagea.shape).astype(float)
        noisyImage = imagea + poissonNoise
        pf.writeto(t+'POISSONPHOTO'+s+str(i)+poisson+'.fits',timage,clobber=True,header=hdr)
        lmass3=cat['lmass'][z2&lmass2&ra2&dec2][i]
        print id_new, ra_new2,dec_new2

Answer 1

Am I interpreting your question (in)correctly? Using pyfits.

import numpy
import scipy
import pyfits

# Use Pyfits to read in a file.
im = pyfits.open("example.fits")

# pyfits opens a "list" of all header extensions in the file.
isinstance(im,list)

# my example is a simple 2d image. so I want the first header unit
im0 = im[0]

# you access the data shape this way
print im0.data.shape

# simple check of the image variance
print numpy.var(im0.data)

# now I'm just repeating the same as your example
poisson = str(raw_input("Poisson noise amount: "))
poissonNoise = numpy.random.poisson(poisson, im0.data.shape).astype(float)
test = im0.data + poissonNoise

# check things out and write to new file
print numpy.var(test)

# you could do this.
im0.data = test

# write that new image back out with the old header. 
pyfits.writeto("/tmp/test.fits", data=test, header=im0.header)

# prove it worked. 
check = pyfits.open("/tmp/test.fits")

numpy.var(check[0].data)