Greetings,
I have an I,J,K 3D volume of amplitude values at regularly sampled time intervals. I have an I,J 2D slice which contains a time (K) value at each I, J location. What I would like to do is extract a subvolume at a constant +/ K window around the slice. Is there an easy way to do this using advanced indexing or some other method? Thanks in advanced for your help.
 Jack
Kind Regards,
Jack Cook Jack Cook
On Wed, Jul 16, 2008 at 16:45, Jack.Cook@shell.com wrote:
Greetings,
I have an I,J,K 3D volume of amplitude values at regularly sampled time intervals. I have an I,J 2D slice which contains a time (K) value at each I, J location. What I would like to do is extract a subvolume at a constant +/ K window around the slice. Is there an easy way to do this using advanced indexing or some other method? Thanks in advanced for your help.
cube[:,:,Khalf_width:K+half_width]
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 twodimensional 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,khalf_width:k+half_width) # shove trace in sub volume
What am I missing?
 Jack
Original Message From: numpydiscussionbounces@scipy.org [mailto:numpydiscussionbounces@scipy.org]On Behalf Of Robert Kern Sent: Wednesday, July 16, 2008 4:56 PM To: Discussion of Numerical Python Subject: Re: [Numpydiscussion] Numpy Advanced Indexing Question
On Wed, Jul 16, 2008 at 16:45, Jack.Cook@shell.com wrote:
Greetings,
I have an I,J,K 3D volume of amplitude values at regularly sampled time intervals. I have an I,J 2D slice which contains a time (K) value at each I, J location. What I would like to do is extract a subvolume at a constant +/ K window around the slice. Is there an easy way to do this using advanced indexing or some other method? Thanks in advanced for your help.
cube[:,:,Khalf_width:K+half_width]
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 twodimensional 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,khalf_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 fancyindex 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.
On Wed, Jul 16, 2008 at 4:08 PM, Robert Kern robert.kern@gmail.com wrote:
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:
+1 for TOW (Tip of the Week).
I'm sure some helpful soul will be kind enough to put this in one of the example/tutorial/cookbook pages.
Cheers,
f
Hi Robert
2008/7/17 Robert Kern robert.kern@gmail.com:
In [42]: smallcube = cube[idx_i,idx_j,idx_k]
Fantastic  a good way to warm up the braincircuit in the morning! Is there an easytoremember rule that predicts the output shape of the operation above? I'm trying to imaging how the output would change if I altered the dimensions of idx_i or idx_j, but it's hard.
It looks like you can do all sorts of interesting things by manipulation the indices. For example, if I take
In [137]: x = np.arange(12).reshape((3,4))
I can produce either
In [138]: x[np.array([[0,1]]), np.array([[1, 2]])] Out[138]: array([[1, 6]])
or
In [140]: x[np.array([[0],[1]]), np.array([[1], [2]])] Out[140]: array([[1], [6]])
and even
In [141]: x[np.array([[0],[1]]), np.array([[1, 2]])] Out[141]: array([[1, 2], [5, 6]])
or its transpose
In [143]: x[np.array([[0,1]]), np.array([[1], [2]])] Out[143]: array([[1, 5], [2, 6]])
Is it possible to separate the indexing in order to understand it better? My thinking was
cube_i = cube[idx_i,:,:].squeeze() cube_j = cube_i[:,idx_j,:].squeeze() cube_k = cube_j[:,:,idx_k].squeeze()
Not sure what would happen if the original array had single dimensions, though.
Back to the original problem:
In [127]: idx_i.shape Out[127]: (10, 1, 1)
In [128]: idx_j.shape Out[128]: (1, 15, 1)
In [129]: idx_k.shape Out[129]: (10, 15, 7)
For the constant slice case, I guess idx_k also have been (1, 1, 7)?
The construction of the cube could probably be done using only
cube.flat = np.arange(nk)
Fernando is right: this is good food for thought and excellent cookbook material!
Regards Stéfan
On Thu, Jul 17, 2008 at 03:16, Stéfan van der Walt stefan@sun.ac.za wrote:
Hi Robert
2008/7/17 Robert Kern robert.kern@gmail.com:
In [42]: smallcube = cube[idx_i,idx_j,idx_k]
Fantastic  a good way to warm up the braincircuit in the morning! Is there an easytoremember rule that predicts the output shape of the operation above? I'm trying to imaging how the output would change if I altered the dimensions of idx_i or idx_j, but it's hard.
Like I said, they all get broadcasted against each other. The final output is the shape of the broadcasted index arrays and takes values found by iterating in parallel over those broadcasted index arrays.
It looks like you can do all sorts of interesting things by manipulation the indices. For example, if I take
In [137]: x = np.arange(12).reshape((3,4))
I can produce either
In [138]: x[np.array([[0,1]]), np.array([[1, 2]])] Out[138]: array([[1, 6]])
or
In [140]: x[np.array([[0],[1]]), np.array([[1], [2]])] Out[140]: array([[1], [6]])
and even
In [141]: x[np.array([[0],[1]]), np.array([[1, 2]])] Out[141]: array([[1, 2], [5, 6]])
or its transpose
In [143]: x[np.array([[0,1]]), np.array([[1], [2]])] Out[143]: array([[1, 5], [2, 6]])
Is it possible to separate the indexing in order to understand it better? My thinking was
cube_i = cube[idx_i,:,:].squeeze() cube_j = cube_i[:,idx_j,:].squeeze() cube_k = cube_j[:,:,idx_k].squeeze()
Not sure what would happen if the original array had single dimensions, though.
You'd have a problem.
So the way fancy indexing interacts with slices is a bit tricky, and this is why we couldn't use the nicer syntax of cube[:,:,idx_k]. All axes with fancy indices are collected together. Their index arrays are broadcasted and iterated over. *For each iterate*, all of the slices are collected, and those sliced axes are *added* to the output array. If you had used fancy indexing on all of the axes, then the iterate would be a scalar value pulled from the original array. If you mix fancy indexing and slices, the iterate is the *array* formed by the remaining slices.
So if idx_k is shaped (ni,nj,3), for example, cube[:,:,idx_k] will have the shape (ni,nj,ni,nj,3). So smallcube[:,:,i,j,k]==cube[:,:,idx_k[i,j,k]].
Is that clear, or am I obfuscating the subject more?
Back to the original problem:
In [127]: idx_i.shape Out[127]: (10, 1, 1)
In [128]: idx_j.shape Out[128]: (1, 15, 1)
In [129]: idx_k.shape Out[129]: (10, 15, 7)
For the constant slice case, I guess idx_k also have been (1, 1, 7)?
The construction of the cube could probably be done using only
cube.flat = np.arange(nk)
Yes, but only due to a weird feature of assigning to .flat. If the RHS is too short, it gets repeated. Since the last axis is contiguous, repeating arange(nk) happily coincides with the desired result of cube[i,j] == arange(nk) for all i,j. This won't check the size, though. If I give it cube.flat=np.arange(nk+1), it will repeat that array just fine, although it doesn't line up.
cube[:,:,:]=np.arange(nk), on the other hand broadcasts the RHS to the shape of cube, then does the assignment. If the RHS cannot be broadcasted to the right shape (in this case because it is not the same length as the final axis of the LHS), an error is raised. I find the reuse of the broadcasting concept to be more memorable, and robust over the (mostly) ad hoc use of plain repetition with .flat.
2008/7/17 Robert Kern robert.kern@gmail.com:
So the way fancy indexing interacts with slices is a bit tricky, and this is why we couldn't use the nicer syntax of cube[:,:,idx_k]. All axes with fancy indices are collected together. Their index arrays are broadcasted and iterated over. *For each iterate*, all of the slices are collected, and those sliced axes are *added* to the output array. If you had used fancy indexing on all of the axes, then the iterate would be a scalar value pulled from the original array. If you mix fancy indexing and slices, the iterate is the *array* formed by the remaining slices.
So if idx_k is shaped (ni,nj,3), for example, cube[:,:,idx_k] will have the shape (ni,nj,ni,nj,3). So smallcube[:,:,i,j,k]==cube[:,:,idx_k[i,j,k]].
Is that clear, or am I obfuscating the subject more?
Crystal, thank you for taking the time to explain! This is such valuable information; we should consider adding a section to numpy.doc.indexing or wherever is more appropriate.
For the constant slice case, I guess idx_k also have been (1, 1, 7)?
The construction of the cube could probably be done using only
cube.flat = np.arange(nk)
[...]
If the RHS cannot be broadcasted to the right shape (in this case because it is not the same length as the final axis of the LHS), an error is raised. I find the reuse of the broadcasting concept to be more memorable, and robust over the (mostly) ad hoc use of plain repetition with .flat.
I've become used to exploiting the repeating property of .flat, and forgot its dangers. Thanks for the reminder!
Cheers Stéfan
participants (4)

Fernando Perez

Jack.Cook＠shell.com

Robert Kern

Stéfan van der Walt