[Numpy-discussion] accuracy issues with numpy arrays?

Anne Archibald peridot.faceted at gmail.com
Tue Apr 29 18:10:26 EDT 2008


On 30/04/2008, eli bressert <bressert at gmail.com> wrote:
>  I'm writing a quick script to import a fits (astronomy) image that has
>  very low values for each pixel. Mostly on the order of 10^-9. I have
>  written a python script that attempts to take low values and put them
>  in integer format. I basically do this by taking the mean of the 1000
>  lowest pixel values, excluding zeros, and dividing the rest of the
>  image by that mean. Unfortunately, when I try this in practice, *all*
>  of the values in the image are being treated as zeros. But, if I use a
>  scipy.ndimage function, I get proper values. For example, I take the
>  pixel that I know has the highest value and do

I think the bug is something else:

>  import pyfits as p
>  import scipy as s
>  import scipy.ndimage as nd
>  import numpy as n
>
>  def flux2int(name):
>         d  = p.getdata(name)
>         x,y     = n.shape(d)
>         l = x*y
>         arr1    = n.array(d.reshape(x*y,1))
>         temp = n.unique(arr1[0]) # This is where the bug starts. All values
>  are treated as zeros. Hence only one value remains, zero.

Actually, since arr1 has shape (x*y, 1), arr1[0] has shape (1,), and
so it has only one entry:

In [82]: A = np.eye(3)

In [83]: A.reshape(9,1)
Out[83]:
array([[ 1.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 1.],
       [ 0.],
       [ 0.],
       [ 0.],
       [ 1.]])

In [84]: A.reshape(9,1)[0]
Out[84]: array([ 1.])

The python debugger is a good way to check this sort of thing out; if
you're using ipython, typing %pdb will start the debugger when an
exception is raised, at which point you can poke around in all your
local variables and evaluate expressions.

>         arr1    = arr1.sort()
>         arr1    = n.array(arr1)
>         arr1    = n.array(arr1[s.where(arr1 >= temp)])
>
>         val = n.mean(arr1[0:1000])
>
>         d               = d*(1.0/val)
>         d               = d.round()
>         p.writeto(name[0,]+'fixed.fits',d,h)

Good luck,
Anne



More information about the NumPy-Discussion mailing list