Handle `PermutedDimsArray`

[tiemvanderdeure]

The chunks of a PermutedDimsArray of a chunked array is not chunked - this is really annoying for Rasters in particular, because the @d macro internally uses PermutedDimsArray

using Rasters, NCDatasets
import Rasters.DiskArrays: eachchunk
A = Rasters.create(
    "myraster.nc", Float64, (X(1:1000), Y(1:1000)); 
    force = true, chunks = (128, 128)
)

Aperm = @d A (Y, X)
eachchunk(Aperm) # 1×1 DiskArrays.GridChunks
eachchunk(A) # 8x8 DiskArrays.GridChunks

Or a more minimal example

using DiskArrays
chunkedA = DiskArrays.TestTypes.ChunkedDiskArray(rand(100,100), (10, 10))
PermutedDimsArray(chunkedA, (2,1)) |> eachchunk

1×1 DiskArrays.GridChunks{2, Tuple{DiskArrays.RegularChunks, DiskArrays.RegularChunks}}:
 (1:100, 1:100)

I'm not entirely sure how this should work - probably we need to get the chunks of the parent and then apply the same permutation, right?

[tiemvanderdeure]

Okay, silly of me, I just found out about PermutedDiskArray and this old issue #42 - but we probably need to handle this in Rasters.

[rafaqz]

No its a real issue - Julia base needs a function rather than a constructor for PermutedDimsArray so we can return what we like, or DiskArrays.jl needs to lower its standards and do what DD does and just return something else from PermutedDimsArray constructor.

Rasters already has a bunch of special casing for isdisk and PermutedDimsArray because the fast path for memory is PermutedDimsArray but for disk arrays its permutedims.

I've pushed for both fixes but got nowhere concrete so far. Maybe chime in on my Base.lazypermutedims issue?

[meggart]

It looks like this is indeed an issue, so I would be ok with being a bit less strict from the DiskArrays side to allow PermutedDimsArray constructors to return something else.