Can't plot Raster based on ArchGDAL overview

[felixcremer]

This is most likely an ArchGDAL limitation but I wanted to also report it here, because it comes from somewhere in the Raster code base using PermuteDims.

julia> using Rasters, CairoMakie
julia> import ArchGDAL as AG, PyramidScheme as PS
julia> rastop = AG.readraster("pyramidmiddle.tif")
GDAL Dataset (Driver: GTiff/GeoTIFF)
File(s): 
  pyramidmiddle.tif

Dataset (width x height): 938 x 938 (pixels)
Number of raster bands: 1
  [GA_ReadOnly] Band 1 (Gray): 938 x 938 (Int16)


julia> b =AG.getband(rastop, 1)
[GA_ReadOnly] Band 1 (Gray): 938 x 938 (Int16)
    blocksize: 938×128, nodata: -32768.0, units: 1.0px + 0.0
    overviews: (0) 469x469 (1) 235x235 

julia> o = AG.getoverview(b,1)
[GA_ReadOnly] Band 1 (Gray): 235 x 235 (Int16)
    blocksize: 128×128, nodata: -32768.0, units: 1.0px + 0.0
    overviews: 

julia> r = Raster("pyramidmiddle.tif")
938×938 Raster{Int16,2} with dimensions: 
  X Projected{Float64} LinRange{Float64}(4.8e6, 5.09998e6, 938) ForwardOrdered Regular Points crs: WellKnownText,
  Y Projected{Float64} LinRange{Float64}(2.09998e6, 1.8e6, 938) ReverseOrdered Regular Points crs: WellKnownText
and reference dimensions: 
  Band Categorical{Int64} 1:1 ForwardOrdered
extent: Extent(X = (4.8e6, 5.09998e6), Y = (1.8e6, 2.09998e6))
missingval: -32768
crs: LOCAL_CS["unnamed",UNIT["unknown",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
parent:
 ⋮      ⋱  

julia> ras = Raster(AG.getoverview(band, 1), PS.agg_axis.(dims(base), 2^(2)))
3750×3750 Raster{Int16,2} with dimensions: 
  X Sampled{Float64} LinRange{Float64}(4.8e6, 5.09998e6, 3750) ForwardOrdered Regular Points,
  Y Sampled{Float64} LinRange{Float64}(2.09998e6, 1.8e6, 3750) ReverseOrdered Regular Points
extent: Extent(X = (4.8e6, 5.09998e6), Y = (1.8e6, 2.09998e6))
missingval: missing


julia> ras = Raster(AG.getoverview(b, 1), PS.agg_axis.(dims(r), 2^(2)))
235×235 Raster{Int16,2} with dimensions: 
  X Sampled{Float64} LinRange{Float64}(4.8e6, 5.09998e6, 235) ForwardOrdered Regular Points,
  Y Sampled{Float64} LinRange{Float64}(2.09998e6, 1.8e6, 235) ReverseOrdered Regular Points
extent: Extent(X = (4.8e6, 5.09998e6), Y = (1.8e6, 2.09998e6))
missingval: missing


julia> plot(ras)
ERROR: MethodError: no method matching read!(::ArchGDAL.IRasterBand{…}, ::PermutedDimsArray{…}, ::Int64, ::Int64, ::Int64, ::Int64)

Closest candidates are:
  read!(::ArchGDAL.AbstractRasterBand, ::T, ::Integer, ::Integer, ::Integer, ::Integer) where T<:(Matrix)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/raster/rasterio.jl:263
  read!(::ArchGDAL.AbstractDataset, ::T, ::Any, ::Integer, ::Integer, ::Integer, ::Integer) where T<:(Array{<:Any, 3})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/raster/rasterio.jl:371
  read!(::ArchGDAL.AbstractDataset, ::T, ::Integer, ::Integer, ::Integer, ::Integer, ::Integer) where T<:(Matrix)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/raster/rasterio.jl:358
  ...

Stacktrace:
  [1] readblock!(band::ArchGDAL.IRasterBand{…}, buffer::PermutedDimsArray{…}, x::UnitRange{…}, y::UnitRange{…})
    @ ArchGDAL ~/.julia/packages/ArchGDAL/lgE4A/src/raster/array.jl:170
  [2] readblock!(::DiskArrays.PermutedDiskArray{…}, ::Matrix{…}, ::UnitRange{…}, ::UnitRange{…})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/permute.jl:24
  [3] getindex_disk(::DiskArrays.PermutedDiskArray{…}, ::UnitRange{…}, ::Vararg{…})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/diskarray.jl:40
  [4] getindex(::DiskArrays.PermutedDiskArray{Int16, 2, PermutedDimsArray{…}}, ::UnitRange{Int64}, ::UnitRange{Int64})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/diskarray.jl:211
  [5] subsetarg(x::DiskArrays.PermutedDiskArray{Int16, 2, PermutedDimsArray{…}}, a::Tuple{UnitRange{…}, UnitRange{…}})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/broadcast.jl:145
  [6] (::DiskArrays.var"#62#64"{Tuple{…}})(i::DiskArrays.PermutedDiskArray{Int16, 2, PermutedDimsArray{…}})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/broadcast.jl:41
  [7] map
    @ ./tuple.jl:291 [inlined]
  [8] (::DiskArrays.var"#61#63"{Matrix{…}, Base.Broadcast.Broadcasted{…}})(cnow::Tuple{UnitRange{…}, UnitRange{…}})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/broadcast.jl:41
  [9] foreach(f::DiskArrays.var"#61#63"{Matrix{…}, Base.Broadcast.Broadcasted{…}}, itr::DiskArrays.GridChunks{2, Tuple{…}})
    @ Base ./abstractarray.jl:3094
 [10] copyto!(dest::Matrix{…}, bc::Base.Broadcast.Broadcasted{…})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/broadcast.jl:38
 [11] materialize!
    @ ./broadcast.jl:914 [inlined]
 [12] materialize!
    @ ./broadcast.jl:911 [inlined]
 [13] _copyto!(dest::Matrix{Int16}, source::DiskArrays.PermutedDiskArray{Int16, 2, PermutedDimsArray{…}})
    @ DiskArrays ~/.julia/packages/DiskArrays/QlfRF/src/array.jl:57
 [14] copyto!
    @ ~/.julia/packages/DiskArrays/QlfRF/src/array.jl:7 [inlined]
 [15] copymutable
    @ ./abstractarray.jl:1199 [inlined]
 [16] _reverse
    @ ./arraymath.jl:60 [inlined]
 [17] reverse
    @ ./arraymath.jl:59 [inlined]
 [18] reverse(A::Raster{…}; dims::Y{…})
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/array/methods.jl:536
 [19] reverse
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/array/methods.jl:534 [inlined]
 [20] _reorder(::Type{…}, x::Raster{…}, dim::Y{…})
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/utils.jl:54
 [21] reorder(x::Raster{…}, orderdim::Y{…})
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/utils.jl:49
 [22] _reorder
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/utils.jl:45 [inlined]
 [23] _reorder(x::Raster{…}, orderdims::Tuple{…})
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/utils.jl:45
 [24] reorder
    @ DimensionalData ~/Documents/PyramidScheme/dev/DimensionalData/src/utils.jl:38 [inlined]
 [25] _reorder
    @ DimensionalDataMakie ~/Documents/PyramidScheme/dev/DimensionalData/ext/DimensionalDataMakie.jl:487 [inlined]
 [26] |>(x::Raster{…}, f::typeof(DimensionalDataMakie._reorder))
    @ Base ./operators.jl:917
 [27] _prepare_for_makie(A::Raster{…}, replacements::Tuple{})
    @ DimensionalDataMakie ~/Documents/PyramidScheme/dev/DimensionalData/ext/DimensionalDataMakie.jl:441
 [28] _surface2(A::Raster{…}, attributes::@Kwargs{}, replacements::Tuple{})
    @ DimensionalDataMakie ~/Documents/PyramidScheme/dev/DimensionalData/ext/DimensionalDataMakie.jl:169
 [29] #plot#20
    @ DimensionalDataMakie ~/Documents/PyramidScheme/dev/DimensionalData/ext/DimensionalDataMakie.jl:138 [inlined]
 [30] plot(A::Raster{…})
    @ DimensionalDataMakie ~/Documents/PyramidScheme/dev/DimensionalData/ext/DimensionalDataMakie.jl:134
 [31] top-level scope
    @ REPL[259]:1
Some type information was truncated. Use `show(err)` to see complete types.

pyramidmiddle.zip

[rafaqz]

Does this work now? No file so I cant check ;)

[felixcremer]

The file is attached to the first comment I will test this again next week.

[rafaqz]

Ugh right thanks

[asinghvi17]

Huh!

ras = Raster(AG.getoverview(b, 1), PS.agg_axis.(dims(r), 2^(2)))
heatmap(read(ras))

heatmap(ras)

This looks like good old DiskArrays collect again.

[rafaqz]

Oh no