rafaqz / Rasters.jl

Raster manipulation for the Julia language
MIT License
215 stars 36 forks source link

Add GDAL multidimensional interface #784

Open asinghvi17 opened 1 month ago

asinghvi17 commented 1 month ago

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.

rafaqz commented 1 month ago

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)

asinghvi17 commented 1 month ago

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.

rafaqz commented 1 month ago

Also we should make a better error here, I though we checked for X/Y/Band...