Overload Base.write for pyramids

src/PyramidScheme.jl

@@ -12,7 +12,7 @@ using DiskArrayEngine: create_outwindows, engine
 using DiskArrays: DiskArrays
 #using YAXArrays: savecube
 using YAXArrayBase: YAXArrayBase as YAB
-using YAXArrays: Cube, YAXArray, to_dataset, savedataset, setchunks, open_dataset
+using YAXArrays: Cube, YAXArray, to_dataset, savedataset, setchunks, open_dataset, savecube
 using Zarr: Zarr, zcreate, zopen, writeattrs
 using DimensionalData: DimensionalData as DD
 using DimensionalData.Dimensions: XDim, YDim
@@ -41,7 +41,7 @@ end
 
 function Pyramid(data::DD.AbstractDimArray; resampling_method= mean ∘ skipmissing, kwargs...)
     pyrdata, pyraxs = getpyramids(resampling_method, data; kwargs...)
-    levels = DD.DimArray.(pyrdata, pyraxs)
+    levels = DD.DimArray.(pyrdata, pyraxs, metadata=DD.Lookups.Metadata())
     meta = Dict(deepcopy(DD.metadata(data)))
     push!(meta, "resampling_method" => "mean_skipmissing")
     Pyramid(data, levels, meta)
@@ -320,7 +320,7 @@ function gen_output(t,s; path=tempname())
     if outsize > 100e6
         # This should be zgroup instead of zcreate, could use savedataset(skelton=true)
         # Dummy dataset with FillArrays with the shape of the pyramidlevel
-        zcreate(t,s...;path,chunks = (1024,1024),fill_value=zero(t))
+        zcreate(t,s...;path,chunks = (1024,1024),fill_value=typemax(t))
     else
         zeros(t,s...)
     end
@@ -463,7 +463,8 @@ function selectlevel(pyramid, ext;target_imsize=(1024, 512))
     minlevel = maximum(dimlevels)
     n_agg = min(max(ceil(Int,minlevel),0),nlevels(pyramid))
     @debug "Selected level $n_agg"
-    levels(pyramid)[n_agg][ext]
+    outlevel = levels(pyramid)[n_agg][ext]
+    outlevel
 end
 
 
@@ -502,19 +503,20 @@ function trans_bounds(
     return Extent(X = xlims, Y = ylims)
 end
 
-function write(path, pyramid::Pyramid; kwargs...)
-    savecube(parent(pyramid), path; kwargs...)
-
-    for (i,l) in enumerate(reverse(pyramid.levels))
-        outpath = joinpath(path, string(i-1))
-        savecube(l,outpath)
+function Base.write(path::AbstractString, pyramid::Pyramid; kwargs...)
+    println("saving base")
+    savecube(YAB.yaxconvert(YAXArray, parent(pyramid)), path; kwargs...)
+        for (i,l) in enumerate(pyramid.levels)
+        @show i
+        outpath = joinpath(path, string(i))
+        savecube(YAB.yaxconvert(YAXArray, l),outpath)
     end
 end
 
 """
     tms_json(dimarray)
 Construct a Tile Matrix Set json description from an AbstractDimArray.
-This assumes, that we use an ag3ycgregation of two by two pixels in the spatial domain to derive the underlying layers of the pyramids. 
+This assumes, that we use an aggregation of two by two pixels in the spatial domain to derive the underlying layers of the pyramids. 
 This returns a string representation of the json and is mainly used for writing the TMS definition to the metadata of the Zarr dataset.
 """
 function tms_json(pyramid)

test/runtests.jl

@@ -71,6 +71,7 @@ end
     using ArchGDAL: ArchGDAL as AG
     using PyramidScheme: PyramidScheme as PS
     using Rasters
+
     data = rand(2000,2000)
     r = Raster(data, (X(1:2000), Y(1:2000)))
     tname = tempname() * ".tif"
@@ -81,6 +82,19 @@ end
     @test PS.nlevels(ptif) == 2
     sub = ptif[1:10,1:10]
     @test sub isa PS.Pyramid
+
+    # Three dimensions
+    data = rand(2000,2000,3)
+    r = Raster(data, (X(1:2000), Y(1:2000), Band(1:3)))
+    tname = tempname() * ".tif"
+    write(tname, r, driver="cog", force=true)
+    ptif = Pyramid(tname)
+    #ras = Raster(path, lazy=true)
+    @test ptif isa PS.Pyramid
+    @test PS.nlevels(ptif) == 2
+    @test ndims(ptif) == 3
+    sub = ptif[1:10,1:10, 1]
+    @test sub isa PS.Pyramid
 end
 
 @testitem "Zarr build Pyramid inplace" begin
@@ -200,9 +214,11 @@ end
 
 @testitem "isequal of two pyramids" begin
     using DimensionalData
+    using PyramidScheme
     pyr1 = Pyramid(fill(1, X(1:128),Y(1:128)), tilesize=16)
     pyr2 = Pyramid(fill(1, X(1:128),Y(1:128)), tilesize=16, resampling_method=sum)
     @test !isequal(pyr1, pyr2)
+    @test isequal(pyr1, deepcopy(pyr1))
 
     a = fill(1, 100,100)
     a[1:2:end, 1:2:end] .= 2
@@ -217,7 +233,24 @@ end
     pyrb = Pyramid(ddb, tilesize=25)
     @test !isequal(pyra, pyrb)
     @test !isequal(pyra, pyr1)
+
+end
+
+@testitem "save pyramid" begin
+    using DimensionalData
+    using YAXArrays
+    using PyramidScheme
+    s = (2048, 1024)
+    a = rand(s...)
+    yax = YAXArray((X(1.:size(a,1)),Y(1.:size(a,2))), a)
+    pyr1 = Pyramid(yax)
+    path = tempname() * ".zarr"
+    write(path, pyr1)
+    pyrload = Pyramid(path)
+    @test isequal(pyr1, pyrload)
 end
+
+
 #=
 @testitem "Comparing zarr pyramid with tif pyramid" begin
     using PyramidScheme: PyramidScheme as PS