cannot read data after `permutedims`

[tiemvanderdeure]

MWE:

using NCDatasets
ds = NCDataset("file.nc","c")
defDim(ds,"lon",100)
defDim(ds,"lat",110)
v = defVar(ds,"temperature",Float32,("lon","lat"), chunksizes = (10,10))
data = [Float32(i+j) for i = 1:100, j = 1:110];
v[:,:] = data;
permutedims(v.var) |> Array

What happens here is that DiskArray has its own lazy PermuteDiskArray type. When you call readblock! it wraps the destination array in aPermutedDimsArray, which C can't handle.

I'm not entire sure if this is best addressed here or in DiskArrays.

Stacktrace:

ERROR: conversion to pointer not defined for PermutedDimsArray{Float32, 2, (2, 1), (2, 1), Matrix{Float32}}
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] unsafe_convert(::Type{Ptr{Float32}}, a::PermutedDimsArray{Float32, 2, (2, 1), (2, 1), Matrix{Float32}})
    @ Base .\pointer.jl:68
  [3] unsafe_convert(::Type{Ptr{Nothing}}, a::PermutedDimsArray{Float32, 2, (2, 1), (2, 1), Matrix{Float32}})
    @ Base .\pointer.jl:66
  [4] nc_get_vars!(ncid::Int32, varid::Int32, startp::Vector{…}, countp::Vector{…}, stridep::Vector{…}, ip::PermutedDimsArray{…})
    @ NCDatasets C:\Users\tsh371\.julia\packages\NCDatasets\LjXBn\src\netcdf_c.jl:1048
  [5] _read_data_from_nc!(::NCDatasets.Variable{…}, ::PermutedDimsArray{…}, ::UnitRange{…}, ::Vararg{…})
    @ NCDatasets C:\Users\tsh371\.julia\packages\NCDatasets\LjXBn\src\variable.jl:407
  [6] readblock!(::NCDatasets.Variable{…}, ::PermutedDimsArray{…}, ::UnitRange{…}, ::Vararg{…})
    @ NCDatasets C:\Users\tsh371\.julia\packages\NCDatasets\LjXBn\src\variable.jl:395
  [7] readblock!(::DiskArrays.PermutedDiskArray{…}, ::Matrix{…}, ::UnitRange{…}, ::Vararg{…})
    @ DiskArrays C:\Users\tsh371\.julia\dev\DiskArrays\src\permute.jl:29
  [8] readblock_checked!
    @ C:\Users\tsh371\.julia\dev\DiskArrays\src\indexing.jl:290 [inlined]
  [9] getindex_disk_nobatch!(out::Nothing, a::DiskArrays.PermutedDiskArray{…}, i::Tuple{…})
    @ DiskArrays C:\Users\tsh371\.julia\dev\DiskArrays\src\indexing.jl:102
 [10] getindex_disk!
    @ C:\Users\tsh371\.julia\dev\DiskArrays\src\indexing.jl:65 [inlined]
 [11] getindex_disk
    @ C:\Users\tsh371\.julia\dev\DiskArrays\src\indexing.jl:58 [inlined]
 [12] getindex
    @ C:\Users\tsh371\.julia\dev\DiskArrays\src\indexing.jl:308 [inlined]
 [13] _disk_collect
    @ C:\Users\tsh371\.julia\dev\DiskArrays\src\array.jl:2 [inlined]
 [14] Array
    @ C:\Users\tsh371\.julia\dev\DiskArrays\src\array.jl:54 [inlined]
 [15] |>(x::DiskArrays.PermutedDiskArray{Float32, 2, PermutedDimsArray{…}}, f::Type{Array})
    @ Base .\operators.jl:926
 [16] top-level scope

[tiemvanderdeure]

This is a little hacky, but it does work:

function DiskArrays.getindex_disk(
    a::DiskArrays.PermutedDiskArray{<:Any,<:Any,<:Base.PermutedDimsArray{<:Any,<:Any,<:Any,<:Any,<:Variable}}, i...
)
    iperm = DiskArrays._getiperm(a)
    # Permute the indices
    inew = Base.PermutedDimsArrays.genperm(i, iperm)
    DiskArrays.getindex_disk!(nothing, a.a.parent, inew...)
end

[Alexander-Barth]

I am wondering if you would have the same issue if you use Zarr.jl (or NetCDF.jl). Can you check?

[tiemvanderdeure]

Just tried and with Zarr.jl it just works (I think it's pure julia so you don't run into the same issue?)

With NetCDF it errors because Base.iscontiguous is not defined on PermutedDimsArray, but after manually defining Base.iscontiguous(::PermutedDimsArray) = true it works (but it seems to just loop through every index?)