Determining 2-D or 3-D locations using the WHERE function

IDL's WHERE function returns the one-dimensional subscript locations of specific data values within an array. It does not return the column, row, and layer locations of a three-dimensional array, or the column and row locations of a two-dimensional array. However, these locations can be easily determined using IDL's MOD operator and integer division. 

Discussion
You can use IDL's SIZE function to compute the dimensions of a given three-dimensional array. With these dimensions and the MOD operator combined with integer division you can determine the X, Y, Z location.

The example code provides a function WHERE_XYZ, which takes the same arguments and keywords as the WHERE function, but in addition it has 3 keywords that return the X, Y, Z locations for 1-D, 2-D or 3-D arrays.

To use the function with a 2D array you could use:

IDL> data=findgen(3,3)
IDL> ind=WHERE_XYZ(data gt 5, XIND=xind ,YIND=yind)
IDL> print, data[xind, yind]
6.00000 7.00000 8.00000
IDL> print, data[ind]
6.00000 7.00000 8.00000


An example of how to use WHERE_XYZ with a 3D array:

IDL> data=findgen(3,3,3)
IDL> ind=WHERE_XYZ(data gt 5 and data lt 10, XIND=xind, YIND=yind, ZIND=zind)
IDL> print, data[ind]
6.00000 7.00000 8.00000 9.00000
IDL> print, data[xind, yind, zind]
6.00000 7.00000 8.00000 9.00000Solution:

function WHERE_XYZ, Array_expression, Count, XIND=xind, YIND=yind,ZIND=zind
  ; works for 1, 2 or 3 dimensional arrays
  ;
  ; Returns the 1D indices (same as WHERE)
  ;
  ; ARGUMENTS
  ;  - same as WHERE (see WHERE)
  ;
  ; KEYWORDS
  ; - Optionally returns X, Y, Z locations through:
  ;
  ; XIND: Output keyword, array of locations along the first dimension
  ; YIND: Output keyword, array of locations along the second dimension (if present)
  ; ZIND: Output keyword, array of locations along the third dimension (if present)
  ;
  ; If no matches where found, then XIND returns -1
  ;
  index_array=where(Array_expression, Count)
  dims=size(Array_expression,/dim)
  xind=index_array mod dims[0]
  case n_elements(dims) of
    2: yind=index_array / dims[0]
    3: begin
        yind=index_array / dims[0] mod dims[1]
        zind=index_array / dims[0] / dims[1]
      end
    else:
  endcase
  return, index_array
end