comp.lang.idl-pvwave archive - RDF feed
type: channel
description: Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
items: Re: Using where() on slices of data cubes

Re: Using where() on slices of data cubes
type: item
description: Conor writes: > I feel like this should be an easy one, but I've never quite figured > it out. Let's say I got a data cube and I want to do something on > just a slice of it, say I want to turn certain values in a column into > something else: > > w = where( cube[1,*,*] lt 0 ) > > It seems like you should be able to do something like this: > > cube[1,w] = 1e24 > > But that doesn't work... Somehow I can't quite figure out the right > way to do this. It *seems* like there should be a simple way to do this, but if there is, I haven't found it. What would make sense to me is this: (cube[1,*,*])[w] = 1e24 But this gives the error message that you can't store into a temporary variable. (The ol' pass by reference/pass by value thing, I suppose.) I have always resorted to making a subset of the data, like this: sub = cube[1,*,*] w = where(sub lt 0) sub[w] = 1e24 cube[1,*,*] = sub Not at all elegant, I admit. :-( Cheers, David -- David Fanning, Ph.D. Fanning Software Consulting, Inc. Coyote's Guide to IDL Programming: http://www.dfanning.com/ Sepore ma de ni thui. ("Perhaps thou speakest truth.")
creator: David Fanning
date: 2009-10-20T13:03:22-00:00

Re: Using where() on slices of data cubes
type: item
description: Well it works, so I can't knock it that much :) I had resorted to using a for loop :( On Oct 20, 9:03 am, David Fanning <n...@dfanning.com> wrote: > Conor writes: >> I feel like this should be an easy one, but I've never quite figured >> it out. Let's say I got a data cube and I want to do something on >> just a slice of it, say I want to turn certain values in a column into >> something else: > >> w = where( cube[1,*,*] lt 0 ) > >> It seems like you should be able to do something like this: > >> cube[1,w] = 1e24 > >> But that doesn't work... Somehow I can't quite figure out the right >> way to do this. > > It *seems* like there should be a simple way to do this, > but if there is, I haven't found it. What would make > sense to me is this: > > (cube[1,*,*])[w] = 1e24 > > But this gives the error message that you can't store into > a temporary variable. (The ol' pass by reference/pass by value > thing, I suppose.) > > I have always resorted to making a subset of the data, like > this: > > sub = cube[1,*,*] > w = where(sub lt 0) > sub[w] = 1e24 > cube[1,*,*] = sub > > Not at all elegant, I admit. :-( > > Cheers, > > David > > -- > David Fanning, Ph.D. > Fanning Software Consulting, Inc. > Coyote's Guide to IDL Programming:http://www.dfanning.com/ > Sepore ma de ni thui. ("Perhaps thou speakest truth.")
creator: Conor
date: 2009-10-20T13:20:47-00:00

Re: Using where() on slices of data cubes
type: item
description: On 20 Okt., 14:22, Conor <cmanc...@gmail.com> wrote: > I feel like this should be an easy one, but I've never quite figured > it out. Let's say I got a data cube and I want to do something on > just a slice of it, say I want to turn certain values in a column into > something else: > > w = where( cube[1,*,*] lt 0 ) > > It seems like you should be able to do something like this: > > cube[1,w] = 1e24 > > But that doesn't work... Somehow I can't quite figure out the right > way to do this. This method avoids the need for any duplication: w=where(cube[1,*,*] lt 0) sz=size(cube,/dim) cube=reform(cube,[sz[0],sz[1]*sz[2]],/overwrite) cube[1,w]=1e24 cube=reform(cube,sz,/overwrite) regards, Greg
creator: greg.addr
date: 2009-10-20T14:20:46-00:00

Re: Using where() on slices of data cubes
type: item
description: On Tue, 20 Oct 2009, Conor wrote: > I feel like this should be an easy one, but I've never quite figured > it out. Let's say I got a data cube and I want to do something on > just a slice of it, say I want to turn certain values in a column into > something else: > > w = where( cube[1,*,*] lt 0 ) > > It seems like you should be able to do something like this: > > cube[1,w] = 1e24 > > But that doesn't work... Somehow I can't quite figure out the right > way to do this. > cube[where(cube[1,*,*] lt 0)*(size(cube, /dim))[0]+1] = 1e24 regards, lajos
creator: Foldy Lajos
date: 2009-10-20T14:26:52-00:00

Re: Using where() on slices of data cubes
type: item
description: > This method avoids the need for any duplication: > > w=where(cube[1,*,*] lt 0) > sz=size(cube,/dim) > cube=reform(cube,[sz[0],sz[1]*sz[2]],/overwrite) > cube[1,w]=1e24 > cube=reform(cube,sz,/overwrite) > > regards, > Greg Is'nt this more efficient ? cube[1,*,*] = 1e24*(cube[1,*,*] lt 0) + cube[1,*,*]*(cube[1,*,*] ge 0) alx.
creator: lecacheux.alain
date: 2009-10-20T14:34:34-00:00

Re: Using where() on slices of data cubes
type: item
description: On Oct 20, 10:22 am, Conor <cmanc...@gmail.com> wrote: > I feel like this should be an easy one, but I've never quite figured > it out. Let's say I got a data cube and I want to do something on > just a slice of it, say I want to turn certain values in a column into > something else: > > w = where( cube[1,*,*] lt 0 ) > > It seems like you should be able to do something like this: > > cube[1,w] = 1e24 > > But that doesn't work... Somehow I can't quite figure out the right > way to do this. It does not apply to this case, but for future reference, it may be useful to know that it could be done directly if the slice was the other way around (a subset of indices on the left): w=where(cube[*,*,i] lt 0) nel=n_elements(cube[*,*,i]) cube[i*nel+w]=1e24
creator: penteado
date: 2009-10-20T14:39:18-00:00

Re: Using where() on slices of data cubes
type: item
description: David Fanning wrote: > Conor writes: > >> I feel like this should be an easy one, but I've never quite figured >> it out. Let's say I got a data cube and I want to do something on >> just a slice of it, say I want to turn certain values in a column into >> something else: >> >> w = where( cube[1,*,*] lt 0 ) >> >> It seems like you should be able to do something like this: >> >> cube[1,w] = 1e24 >> >> But that doesn't work... Somehow I can't quite figure out the right >> way to do this. > > It *seems* like there should be a simple way to do this, > but if there is, I haven't found it. What would make > sense to me is this: > > (cube[1,*,*])[w] = 1e24 > > But this gives the error message that you can't store into > a temporary variable. (The ol' pass by reference/pass by value > thing, I suppose.) Is it possible at all to use IDL pointer for aliasing? E.g. in Fortran I would do something like the following real, target :: cube(:,:,:) real, pointer :: alias(:,:) => null() alias => cube(1,:,:) where ( alias < 0.0 ) alias = 1e24 end where nullify(alias) Is there an IDL equivalent to the "alias => cube(1,:,:)" statement in Fortran? If there is, then that might get around the tired old "can't store into a temporary variable" issue that plagues IDL due to it sticking to pass by reference only. (Man, I wish they would fix that.) The only way I can figure out to do the aliasing ends up copying the data. Anyway... cheers, paulv
creator: Paul Van Delst[1]
date: 2009-10-20T15:00:32-00:00

Re: Using where() on slices of data cubes
type: item
description: > Is'nt this more efficient ? > cube[1,*,*] = 1e24*(cube[1,*,*] lt 0) + cube[1,*,*]*(cube[1,*,*] ge > 0) I reckon you're copying out the slice twice, processing through it 5 times, and copying it back. That's compared to one scan through the slice and one indexed substitution (possibly many fewer items than one slice). The 'reform' with /overwrite costs nothing. regards, Greg
creator: greg.addr
date: 2009-10-20T16:18:43-00:00

Re: Using where() on slices of data cubes
type: item
description: Lajos writes: > > On Tue, 20 Oct 2009, Conor wrote: > >> I feel like this should be an easy one, but I've never quite figured >> it out. Let's say I got a data cube and I want to do something on >> just a slice of it, say I want to turn certain values in a column into >> something else: >> >> w = where( cube[1,*,*] lt 0 ) >> >> It seems like you should be able to do something like this: >> >> cube[1,w] = 1e24 >> >> But that doesn't work... Somehow I can't quite figure out the right >> way to do this. >> > > > cube[where(cube[1,*,*] lt 0)*(size(cube, /dim))[0]+1] = 1e24 Humm. Well, I suppose, for this *specific* problem. I can't wait to see the general solution, though, for a slice in any direction. Do you suppose it will fit on a page of paper? :-) Cheers, David -- David Fanning, Ph.D. Coyote's Guide to IDL Programming (www.dfanning.com) Sepore ma de ni thui. ("Perhaps thou speakest truth.")
creator: David Fanning
date: 2009-10-20T20:03:29-00:00

Re: Using where() on slices of data cubes
type: item
description: On Oct 20, 8:22 am, Conor <cmanc...@gmail.com> wrote: > I feel like this should be an easy one, but I've never quite figured > it out. Let's say I got a data cube and I want to do something on > just a slice of it, say I want to turn certain values in a column into > something else: > > w = where( cube[1,*,*] lt 0 ) > > It seems like you should be able to do something like this: > > cube[1,w] = 1e24 > > But that doesn't work... Somehow I can't quite figure out the right > way to do this. It's impossible for IDL to know the dimensionality of the original array which from which you extracted a given set of elements. WHERE simply gives you the zero-based running index into a given pile of data, independent of its form and (much less) the form of the parent array from which it was derived. So it's not at all surprising that the returned set of indices doesn't "plug right in" in some convenient way. The (likely) quickest solution in this case is the REFORM,/ OVERWRITE method given by Greg, since it doesn't actually copy any data. But unfortunately, it's not at all general. For example, suppose you had been interested in longitudinal slices instead, ala: w=where(cube[*,1,*] lt 0) A solution similar to Greg's would require first TRANSPOSE'ing the cube such that the elements of the slice of interest could be put in order along some dimension, e.g.: w=where(cube[*,1,*] lt 0) sz=size(cube,/dim) cube=reform(transpose(cube,[0,2,1]),[sz[0]*sz[2],sz[1]]) cube[w,1]=1e24 cube=transpose(reform(cube,sz[0],sz[2],sz[1],/overwrite),[0, 2,1]) This obviously creates two transposed copies of the cube, which is memory and CPU inefficient, not to mention difficult to remember. For operating on small cubes, or a large fraction of the cube elements, it might be reasonable. If, however, you have only a few elements to access in a very large array (say just a sprinkling of negatives, in your example), this would be very wasteful indeed. In my opinion, a far better general solution to these sorts…
creator: JDS
date: 2009-10-20T20:03:30-00:00

Re: Using where() on slices of data cubes
type: item
description: JD Smith writes: > you should easily be able to generalize the above arguments to access > these elements I think in this case the word "easily" might be too subtly sarcastic to be easily appreciated by the vast majority of this newsgroup. :-) Cheers, David -- David Fanning, Ph.D. Coyote's Guide to IDL Programming (www.dfanning.com) Sepore ma de ni thui. ("Perhaps thou speakest truth.")
creator: David Fanning
date: 2009-10-20T20:32:35-00:00

Re: Using where() on slices of data cubes
type: item
description: On Oct 20, 4:32 pm, David Fanning <n...@dfanning.com> wrote: > JD Smith writes: >> you should easily be able to generalize the above arguments to access >> these elements > > I think in this case the word "easily" might be > too subtly sarcastic to be easily appreciated by > the vast majority of this newsgroup. :-) (Almost) no sarcasm was intended. Suppose you have this: w=where(cube[1,5:*,10:1024] lt 0) The "slice" is no longer as large as the cube in the yz dimensions, and is offset by [5,10] too. So y_full_cube = slice_column + 5 z_full_cube = slice_row + 10 and since the slice is smaller than the cube by 5 columns, to convert our WHERE index vector w into col,row in the slice, we use slice_column = w mod (sz[1]-5) slice_row = w/(sz[1]-5) Putting it all together we have: ind = 1 + sz[0] * (5 + w mod (sz[1]-5) + (10 + w/(sz[1]-5)) * sz[1]) JD
creator: JDS
date: 2009-10-21T21:52:45-00:00