Add GDAL multidimensional interface
I can't use NCDatasets.jl because of DiskArrays/Zarr version mismatches.
using Rasters, ArchGDAL
ras = Raster(rand(LinRange(0, 10, 100), X(1:100), Y(5:150), Ti(1:12)))
write("test.nc", ras; source = Rasters.GDALsource(), force = true)
julia> write("test.nc", ras; source = Rasters.GDALsource(), force = true)
ERROR: ArgumentError: no valid permutation of dimensions
Stacktrace:
[1] permutedims(B::Array{Float64, 3}, perm::Tuple{Int64, Int64})
@ Base ./multidimensional.jl:1668
[2] permutedims
@ ~/.julia/dev/DimensionalData/src/array/methods.jl:230 [inlined]
[3] _maybe_permute_to_gdal(A::Raster{Float64, 3, Tuple{X{…}, Y{…}, Ti{…}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing}, dims::Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{…}}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}})
@ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:572
[4] _maybe_permute_to_gdal(A::Raster{Float64, 3, Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{…}}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}, Ti{DimensionalData.Dimensions.Lookups.Sampled{…}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing})
@ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:569
[5] |>(x::Raster{Float64, 3, Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{…}}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}, Ti{DimensionalData.Dimensions.Lookups.Sampled{…}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing}, f::typeof(RastersArchGDALExt._maybe_permute_to_gdal))
@ Base ./operators.jl:972
[6] _maybe_correct_to_write(lookup::DimensionalData.Dimensions.Lookups.Sampled{…}, A::Raster{…}, args::Rasters.NoKW)
@ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:329
[7] _maybe_correct_to_write(A::Raster{Float64, 3, Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{…}}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}, Ti{DimensionalData.Dimensions.Lookups.Sampled{…}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing}, args::Rasters.NoKW)
@ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:323
[8] write(filename::String, ::Rasters.GDALsource, A::Raster{Float64, 3, Tuple{X{…}, Y{…}, Ti{…}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing}; force::Bool, verbose::Bool, missingval::Rasters.NoKW, kw::@Kwargs{})
@ RastersArchGDALExt ~/.julia/dev/Rasters/ext/RastersArchGDALExt/gdal_source.jl:56
[9] write(filename::String, A::Raster{Float64, 3, Tuple{X{DimensionalData.Dimensions.Lookups.Sampled{…}}, Y{DimensionalData.Dimensions.Lookups.Sampled{…}}, Ti{DimensionalData.Dimensions.Lookups.Sampled{…}}}, Tuple{}, Array{Float64, 3}, DimensionalData.NoName, DimensionalData.Dimensions.Lookups.NoMetadata, Nothing}; source::Rasters.GDALsource, kw::@Kwargs{force::Bool})
@ Rasters ~/.julia/dev/Rasters/src/write.jl:78
[10] top-level scope
@ REPL[50]:1
Some type information was truncated. Use `show(err)` to see complete types.
If I slice this by Ti(At(1)), it just works. Similarly, if I replace Ti with Band, it works.
Yeah, you cant save a time dimension with GDAL in Rasters. It's only using the band interface.
There's not much point implementing that whole multidimensional gdal interface when we have DiskArrays/CommonDataModel etc. But I do think someone recently added it to ArchGDAL at least.
So if you can be bothered, but I won't.
(And yes @meggart @Alexander-Barth this version split with DiskArrays is getting pretty dire)
Yeah I just found that PR, looks like it's still open and under review but close to the finish line. Not at all required since I found the Band hack. But it would be nice to have.
The error could be better though, might improve on that.
Also we should make a better error here, I though we checked for X/Y/Band...