[Numpy-discussion] Numpy Advanced Indexing Question

Robert Kern robert.kern@gmail....
Wed Jul 16 18:08:12 CDT 2008


On Wed, Jul 16, 2008 at 17:12,  <Jack.Cook@shell.com> wrote:
> Robert,
>
> I can understand how this works if K is a constant time value but in my case K varies at each location in the two-dimensional slice. In other words, if I was doing this in a for loop I would do something like this
>
> for i in range(numI):
>   for j in range(numJ):
>      k = slice(i,j)
>      trace = cube(i,j,k-half_width:k+half_width)
>      # shove trace in sub volume
>
> What am I missing?

Ah, okay. It's a bit tricky, though. Yes, you need to use fancy
indexing. Since axis you want to be index fancifully is not the first
one, you have to be more explicit than you might otherwise want. For
example, it would be great if you could just use slices for the first
two axes:

  cube[:,:,slice + numpy.arange(-half_width,half_width)]

but the semantics of that are a bit different for reasons I can
explain later, if you want. Instead, you have to have explicit
fancy-index arrays for the first two axes. Further, the arrays for
each axis need to be broadcastable to each other. Fancy indexing will
iterate over these broadcasted arrays in parallel with each other to
form the new array. The liberal application of numpy.newaxis will help
us achieve that. So this is the complete recipe:


In [29]: import numpy
In [30]: ni, nj, nk = (10, 15, 20)

# Make a fake data cube such that cube[i,j,k] == k for all i,j,k.
In [31]: cube = numpy.empty((ni,nj,nk), dtype=int)
In [32]: cube[:,:,:] = numpy.arange(nk)[numpy.newaxis,numpy.newaxis,:]

# Pick out a random fake horizon in k.
In [34]: kslice = numpy.random.randint(5, 15, size=(ni, nj))
In [35]: kslice
Out[35]:
array([[13, 14,  9, 12, 12, 11,  8, 14, 11, 13, 13, 13,  8, 11,  8],
       [ 7, 12, 12,  6, 10, 12,  9, 11, 13,  9, 14, 11,  5, 12, 12],
       [ 7,  5, 10,  9,  6,  5,  5, 14,  5,  6,  7, 10,  6, 10, 11],
       [ 6,  9, 11, 14,  7, 11, 10,  6,  6,  9,  9, 11,  5,  5, 14],
       [12,  8, 11,  6, 10,  8,  5,  9,  8, 10,  7,  5,  9,  9, 14],
       [ 9,  8, 10,  9, 10, 12, 10, 10,  6, 10, 11,  6,  8,  7,  7],
       [11, 12,  7, 13,  5,  5,  8, 14,  5, 14,  9, 10, 12,  7, 14],
       [ 7,  7,  7, 12, 10,  6, 13, 13, 11, 13,  8, 11, 13, 14, 14],
       [ 6, 13, 13, 10, 10, 14, 10,  8,  9, 14, 13, 12,  9,  9,  5],
       [13, 14, 10,  8, 11, 11, 10,  6, 12, 11, 12, 12, 13, 11,  7]])

In [36]: half_width = 3

# These two replace the empty slices for the first two axes.
In [37]: idx_i = numpy.arange(ni)[:,numpy.newaxis,numpy.newaxis]
In [38]: idx_j = numpy.arange(nj)[numpy.newaxis,:,numpy.newaxis]

# This is the substantive part that actually picks out our window.
In [41]: idx_k = kslice[:,:,numpy.newaxis] +
numpy.arange(-half_width,half_width+1)
In [42]: smallcube = cube[idx_i,idx_j,idx_k]
In [43]: smallcube.shape
Out[43]: (10, 15, 7)

# Now verify that our window is centered on kslice everywhere:
In [47]: smallcube[:,:,3]
Out[47]:
array([[13, 14,  9, 12, 12, 11,  8, 14, 11, 13, 13, 13,  8, 11,  8],
       [ 7, 12, 12,  6, 10, 12,  9, 11, 13,  9, 14, 11,  5, 12, 12],
       [ 7,  5, 10,  9,  6,  5,  5, 14,  5,  6,  7, 10,  6, 10, 11],
       [ 6,  9, 11, 14,  7, 11, 10,  6,  6,  9,  9, 11,  5,  5, 14],
       [12,  8, 11,  6, 10,  8,  5,  9,  8, 10,  7,  5,  9,  9, 14],
       [ 9,  8, 10,  9, 10, 12, 10, 10,  6, 10, 11,  6,  8,  7,  7],
       [11, 12,  7, 13,  5,  5,  8, 14,  5, 14,  9, 10, 12,  7, 14],
       [ 7,  7,  7, 12, 10,  6, 13, 13, 11, 13,  8, 11, 13, 14, 14],
       [ 6, 13, 13, 10, 10, 14, 10,  8,  9, 14, 13, 12,  9,  9,  5],
       [13, 14, 10,  8, 11, 11, 10,  6, 12, 11, 12, 12, 13, 11,  7]])

In [50]: (smallcube[:,:,3] == kslice).all()
Out[50]: True


Clear as mud? I can go into more detail if you like particularly about
how newaxis works.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
 -- Umberto Eco


More information about the Numpy-discussion mailing list