-2

I have a small problem. I know there are some similar questions but I do not know that I am doing wrong because they are not working for me, I would appreciate any help. I want to change some pixel values in a fits file. They are basically empty spots and I want to fill them with ~the mean pixels value of the image. I do like this:

from __future__ import division
import pyfits as fits
import numpy as np
obj1 = fits.open(raw_input('Name of the image to be improved? ')) 
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 x in range(0,shape[0]):
    for y in range(0,shape[1]):       
        if data_obj1[x,y] > -5.48e-14 and data_obj1[x,y] < -5.46e-14:
            random_noise = np.random.choice(noise,1) 
            result[x,y] = random_noise
        else:
            result[x,y] = data_obj1[x,y]
out = obj1
out[0].data = result
out.writeto(raw_input('Name of the output file? '), clobber=True)

I know it is doing the operation I want to do, because if I print result[x,y] it is how it is supposed to be. Nevertheless, when I open the generated fits file, it is exactly the same as it was at the beginning. So probably I do not understand i) how to properly save the fits file or ii) how to build my new image correctly. Can someone help me?

Enrique
  • 55
  • 7
  • Have you considered switching to `astropy`? `pyfits` has been unsupported for several years now. You just need to install `astropy` and replace the `import pyfits as fits` with `from astropy.io import fits`. I'm using astropy 2.0 and your script works for me. How do you compare the files? You can easily compare files using `print(fits.FITSDiff(filename1, filename2).report())` – MSeifert Aug 17 '17 at 22:15
  • 1
    You have typos/syntax errors in your code sample. You're also confused about array indices/axes. – Iguananaut Aug 17 '17 at 22:30
  • @Enrique I don't know why it doesn't work for you. But I used this [script](https://gist.github.com/MSeifert04/2e8096b67242b6b44d9fa00350e9d05b) and that reports several differences for me, so I assume it's working correctly. – MSeifert Aug 17 '17 at 22:34
  • Hello @MSerifert, thanks! I changed to astropy, but something weird is happening. I look at both files with ds9 (via IRAF), and they look exactly the same, with the same pixel values at the same positions. Nevertheless, when I use your command to see the differences, it says that there are differences in one value of the header (in BITPIX, for one 32 and 64 for the other, for some unknown reason) and in some pixel values (in those where I want!). So I do not understand if something is wrong or if it is a visualization problem. Do you have any idea? Also no idea about the change in header – Enrique Aug 17 '17 at 22:42
  • @Iguananaut, I only found one typo in my code, and it was absolutely irrelevant, I think. Anyway, now that you are pointing out I am confused with something, it would be nice if you could clarify this with one line of text, instead of just mention it. :) – Enrique Aug 17 '17 at 22:50
  • The BITPIX difference is because you you create a `float64` array (`result = np.zeros(shape)`) if you want to keep the data type "identical" then you could use `result = np.zeros_like(data_obj1)`. No idea why DS9 doesn't show differences. Are you looking at the correct pixels? It could be that DS9 switches the x and y axis and adds 1 to the indexes (that's fits convention but astropy + pyfits use the numpy convention for indexing). – MSeifert Aug 17 '17 at 22:51
  • 1
    @Enrique Your `else` has no indentation and you should iterate over `for y in range(0,shape[1]):` (note that it's the second item in `shape`) and `for x in range(0,shape[0]):` (first element) that's probably what he meant with "confusing the axis". Or keep your loops but then you have to index with `[y, x]` instead of `[x, y]`. – MSeifert Aug 17 '17 at 22:52
  • @MSeifert, oh, I see, many thanks again! And yes, the indentation was also a typo, as well as the [x,y], since the image is a square I did not realize, but thanks for explaining me! Now I know that the code works and I just have to find out what is going on with the visualization. Thanks! – Enrique Aug 17 '17 at 23:04

1 Answers1

-2

Apart from typos explained by @MSeifert, it is just a visualization problem. See the comments for clarifications!

Enrique
  • 55
  • 7
  • 1
    This is not an answer. You can self-answer of course and that's highly encouraged, but please post your actual answer here. – Iguananaut Aug 18 '17 at 10:22