Feature request: way to track offsets from additional dimensions

[sefffal]

Over at https://github.com/JuliaAstro/AstroImages.jl/issues/29#issuecomment-1024632322 we have a use case that OffsetArrays is almost perfect for. In astronomy, it's common to have physical coordinates associated with the axes of a data cube. For instance, angular X/Y coordinates along the first and second axes and wavelength along the third axes (or polarization, or time, etc.). This is represented by something called a WCS transform for which we have the library WCSTransforms.jl. The WCS transform has functions that allow us to map from image indices to physical coordinates and back.

What's unusual is that a WCSTransform can also contain an affine transformation matrix that describes for example, rotation between these axes. That means that in general, moving along the a third (say, spectral) axis might change the physical coordinates associated with indices along the first and second axes.

OffsetArrays could almost solve this for us, since we can track the image coordinates through slicing etc. But if a user selects say a 2D slice of a 3+D cube, there's no way at the moment for an OffsetArray to track the offset along that third "droped" dimension.

What I propose is that the definition of OffsetArrays be expanded slightly so that

o = OffsetArray(rand(100,100), (1:100, 1:100, 5))
# where 5 is the position along a third axis

becomes valid.

Calling axes(o) would then return something like

(OffsetArrays.IdOffsetRange(values=1:100, indices=1:100), OffsetArrays.IdOffsetRange(values=1:100, indices=1:100), 5)

But ndims(o) would remain 2. In general we would want this to work for more than 3 dimensions with integer offsets occurring at any position, not just at the end.

Apologies if this request is completely out there or out of scope for OffsetArrays, just curious to hear what people think.

[timholy]

Does this already do what you want?

julia> A = OffsetArray([1 3; 2 4], 0:1, -11:-10)
2×2 OffsetArray(::Matrix{Int64}, 0:1, -11:-10) with eltype Int64 with indices 0:1×-11:-10:
 1  3
 2  4

julia> v = view(A, :, -10)
2-element view(OffsetArray(::Matrix{Int64}, 0:1, -11:-10), :, -10) with eltype Int64 with indices 0:1:
 3
 4

julia> dump(v)
SubArray{Int64, 1, OffsetMatrix{Int64, Matrix{Int64}}, Tuple{Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}, Int64}, true}
  parent: OffsetMatrix{Int64, Matrix{Int64}}
    parent: Array{Int64}((2, 2)) [1 3; 2 4]
    offsets: Tuple{Int64, Int64}
      1: Int64 -1
      2: Int64 -12
  indices: Tuple{Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}, Int64}
    1: Base.Slice{OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}}
      indices: OffsetArrays.IdOffsetRange{Int64, Base.OneTo{Int64}}
        parent: Base.OneTo{Int64}
          stop: Int64 2
        offset: Int64 -1
    2: Int64 -10
  offset1: Int64 3
  stride1: Int64 1

[jishnub]

In particular, in the above example, you may obtain the index along the dropped dimension using

julia> parentindices(v)
(Base.Slice(OffsetArrays.IdOffsetRange(values=0:1, indices=0:1)), -10)

Do views work for you? Or do you require slices?

[sefffal]

Thank you these are great points -- yes this is almost perfect except for the fact that it creates a view instead of a copy. Sometimes we would want to create a view, and other times a slice.