test_fits.py

from astropy.io import fits
import numpy as np

hdul = fits.HDUList(fits.PrimaryHDU(np.random.random((10, 10))))
hdul.writeto('tmp1.fits', overwrite=True)

from __future__ import division
import numpy as np
obj1 = fits.open('tmp1.fits') 
data_obj1 = obj1[0].data
meanpix = np.mean(data_obj1)
noise = np.linspace(-meanpix,meanpix,100000)
shape = data_obj1.shape
result = np.zeros(shape)
for y in range(0,shape[0]):
    for x in range(0,shape[1]):       
        if 0.5 > data_obj1[x,y] > 0.4:
            random_noise = np.random.choice(noise,1) 
            result[x,y] = 2
        else:
            result[x,y] = data_obj1[x,y]
out = obj1
out[0].data = result
out.writeto('tmp2.fits', overwrite=True)

print(fits.FITSDiff('tmp1.fits', 'tmp2.fits').report())