Fix building pyramids and add cat method (#75)

* Fix building pyramids

* Convert missing to nan before plotting

Co-authored-by: Anshul Singhvi <[email protected]>

* Update dependencies

* Add cat for Pyramids

This also improves the computation of the aggregated dimensions for points
Co-authored-by: Anshul Singhvi <[email protected]>

---------

Co-authored-by: Felix Cremer <[email protected]>
Co-authored-by: Anshul Singhvi <[email protected]>

Author: 3 people

Project.toml

@@ -35,22 +35,22 @@ ArchGDAL = "0.10.4"
 CairoMakie = "0.12"
 Colors = "0.12"
 DimensionalData = "0.28, 0.29"
-DiskArrayEngine = "0.1.2"
+DiskArrayEngine = "0.1.2, 0.2"
 DiskArrayTools = "0.1.10"
 DiskArrays = "0.4.2"
 Extents = "0.1.3"
 FillArrays = "1.11"
 Flatten = "0.4.3"
 GeoInterface = "1.3.4"
-Makie = "0.21.2"
+Makie = "0.21.2, 0.22, 0.23, 0.24"
 OffsetArrays = "1.14"
 Proj = "1.7.1"
 Rasters = "0.11.8, 0.12, 0.13"
 Statistics = "1.10"
 Test = "1.10"
 TestItemRunner = "0.2.3"
 YAXArrayBase = "0.6.1, 0.7"
-YAXArrays = "0.5.6"
+YAXArrays = "0.5.6, 0.6"
 Zarr = "0.9.3"
 julia = "1.9"
 

ext/PyramidSchemeMakieExt.jl

@@ -6,15 +6,16 @@ using PyramidScheme: Pyramid, switchkeys, levels, selectlevel, xkey, ykey
 using DimensionalData: DimensionalData as DD
 using DimensionalData.Dimensions: XDim, YDim
 using Extents: Extent, extent, intersects
+miss2nan(x) = ismissing(x) ? NaN : x
 """
     plot(pyramids)
 Plot a Pyramid. 
 This will plot the coarsest resolution level at the beginning and will plot higher resolutions after zooming into the plot.
 This is expected to be used with interactive Makie backends.
 """
-function plot(pyramid::Pyramid;colorbar=true, kwargs...)
+function plot(pyramid::Pyramid;colorbar=true, size=(1155, 1043), kwargs...)
     #This should be converted into a proper recipe for Makie but this would depend on a pyramid type.
-    fig = Figure()
+    fig = Figure(;size)
     lon, lat = DD.dims(DD.parent(pyramid))
     ax = Axis(fig[1,1], limits=(extrema(lon), extrema(lat)), aspect=DataAspect())
     hmap = plot!(ax, pyramid;kwargs...)
@@ -28,19 +29,17 @@ end
 function plot!(ax, pyramid::Pyramid;interp=false, kwargs...)#; rastercrs=crs(parent(pyramid)),plotcrs=EPSG(3857), kwargs...)
     tip = levels(pyramid)[end-2][:,:]
     #@show typeof(tip)
+    subtypes = Union{typeof.(pyramid.levels)..., typeof(parent(pyramid))}
     data = Observable{DD.AbstractDimMatrix}(tip)
     xval = only(values(extent(pyramid, XDim)))
     yval = only(values(extent(pyramid, YDim)))
     rasdataext = Extent(X=xval, Y=yval)
     rasext = extent(pyramid)
-    xk = xkey(rasext)
-    yk = ykey(rasext)
     on(ax.scene.viewport) do viewport
         limext = extent(ax.finallimits[])
 
         datalimit = switchkeys(limext, rasext)
-
-        data.val = selectlevel(pyramid, datalimit, target_imsize=viewport.widths)
+        data.val = miss2nan.(selectlevel(pyramid, datalimit, target_imsize=viewport.widths))
         notify(data)
     end
     on(ax.finallimits) do limits
@@ -53,7 +52,7 @@ function plot!(ax, pyramid::Pyramid;interp=false, kwargs...)#; rastercrs=crs(par
             rasdata = selectlevel(pyramid, datalimit, target_imsize=ax.scene.viewport[].widths)
             # Project selected data to plotcrs
             #data.val = Rasters.resample(rasdata, crs=plotcrs, method=:bilinear )
-            data.val = rasdata
+            data.val = miss2nan.(rasdata)
         end
         notify(data)
     end

src/PyramidScheme.jl

@@ -16,7 +16,7 @@ using YAXArrays: Cube, YAXArray, to_dataset, savedataset, setchunks, open_datase
 using Zarr: Zarr, zcreate, zopen, writeattrs
 using DimensionalData: DimensionalData as DD
 using DimensionalData.Dimensions: XDim, YDim
-using Extents: Extent, extent
+using Extents: Extent, extent, intersects
 using FillArrays: Zeros
 using Proj
 using OffsetArrays
@@ -39,7 +39,7 @@ struct Pyramid{T,N,D,A,B<:DD.AbstractDimArray{T,N,D,A},L, Me} <: DD.AbstractDimA
     metadata::Me
 end
 
-function Pyramid(data::DD.AbstractDimArray; resampling_method=mean ∘ skipmissing, kwargs...)
+function Pyramid(data::DD.AbstractDimArray; resampling_method= mean ∘ skipmissing, kwargs...)
     pyrdata, pyraxs = getpyramids(resampling_method, data; kwargs...)
     levels = DD.DimArray.(pyrdata, pyraxs)
     meta = Dict(deepcopy(DD.metadata(data)))
@@ -103,7 +103,7 @@ function DD.modify(f, pyr::Pyramid)
     Pyramid(pbase, plevels, pyr.metadata)
 end
 Base.read(pyr::Pyramid) = DD.modify(Array, pyr)
-@inline function DD.rebuild(A::Pyramid, data, dims::Tuple=dims(A), refdims=refdims(A), name=name(A))
+@inline function DD.rebuild(A::Pyramid, data, dims::Tuple=DD.dims(A), refdims=DD.refdims(A), name=DD.name(A))
     Pyramid(DD.rebuild(parent(A), data, dims, refdims, name, nothing), A.levels, A.metadata)
 end
 
@@ -271,9 +271,20 @@ Compute the number of levels for the aggregation based on the size of `data`.
 compute_nlevels(data, tilesize=256) = max(0,ceil(Int,log2(maximum(size(data))/tilesize)))
 
 function agg_axis(d,n)
-    # TODO this might be problematic for explicitly set axes
-    DD.set(d, LinRange(first(d), last(d), cld(length(d), n)))
+    # TODO this might be problematic for Reversed axes
+    # TODO this is only correct for points not intervals
+    npoints = cld(length(d), n)
+    half_stepsize = step(d) * (n-1) / 2 
+    @show half_stepsize
+    sgn = DD.isreverse(d) ? -1 : 1
+    DD.set(d, LinRange(first(d) + sgn * half_stepsize, last(d) - sgn * half_stepsize, npoints))
 end
+#=
+.   .
+  .    <--- new point
+.   .
+=#
+
 """
     gen_output(t,s)
 
@@ -326,7 +337,11 @@ function buildpyramids(path::AbstractString; resampling_method=mean, recursive=t
     # Build a loop for all variables in a dataset?
     org = Cube(path)
     # We run the method once to derive the output type
-    t = typeof(resampling_method(zeros(eltype(org), 2,2)))
+    #tfunc = typeof(resampling_method(zeros(eltype(org), 2,2)))
+    #t = Missing <: eltype(org) ? Union{Missing, tfunc} : tfunc
+
+    t = Base.infer_return_type(resampling_method, (Matrix{nonmissingtype(eltype(org))},))
+
     n_level = compute_nlevels(org)            
     input_axes = filter(x-> x isa SpatialDim,  DD.dims(org))
     if length(input_axes) != 2
@@ -377,16 +392,17 @@ end
 Compute the data of the pyramids of a given data cube `ras`.
 This returns the data of the pyramids and the dimension values of the aggregated axes.
 """
-function getpyramids(reducefunc, ras;recursive=true)
+function getpyramids(reducefunc, ras;recursive=true, tilesize=256)
     input_axes = DD.dims(ras)
-    n_level = compute_nlevels(ras)            
+    n_level = compute_nlevels(ras, tilesize)
     if iszero(n_level)
-        @info "Array is smaller than the tilesize no pyramids are computed"
-        [ras], [dims(ras)]
+        @info "Array is smaller than the tilesize no pyramidlevels are computed"
+        [ras], [DD.dims(ras)]
     end 
     pyramid_sizes =  [ceil.(Int, size(ras) ./ 2^i) for i in 1:n_level]
     pyramid_axes = [agg_axis.(input_axes,2^i) for i in 1:n_level]
-    outtype = typeof(reducefunc(rand(eltype(ras),2,2)))
+    outtype = Base.infer_return_type(reducefunc, (Matrix{eltype(ras)},))
+    #outtype = Missing <: eltype(ras) ? Union{Missing, outtypefunc} : outtypefunc
     outmin = output_arrays(pyramid_sizes, outtype)
     fill_pyramids(ras,outmin,reducefunc,recursive; threaded=true)
 
@@ -402,6 +418,7 @@ Internal function to select the raster data that should be plotted on screen.
 """
 function selectlevel(pyramid, ext;target_imsize=(1024, 512))
     pyrext = extent(pyramid)
+    intersects(pyrext, ext) || return zero(pyramid.levels[end])
     basepixels = map(keys(pyrext)) do bb
         pyrspan = pyrext[bb][2] - pyrext[bb][1]
         imsize = ext[bb][2] - ext[bb][1]
@@ -473,5 +490,15 @@ function tms_json(pyramid)
     push!(tms, "orderedAxes" => pyramidaxes())
     return tms
 end
+function Base.cat(A1::Pyramid, As::Pyramid...;dims)
+    println("Inside pyr cat")
+    @show typeof(levels.(As, 1))
+    catlevels = [cat(A1.levels[i], levels.(As, i)...; dims) for i in eachindex(A1.levels)]
+    catbase = cat(parent(A1), parent.(As)...; dims)
+    Pyramid(catbase, catlevels, merge(DD.metadata(A1), DD.metadata.(As)...))
+end
+
+
+
 include("broadcast.jl")
 end

test/runtests.jl

@@ -122,6 +122,21 @@ end
 end
 
 
+@testitem "cat Pyramids" begin
+    using PyramidScheme: PyramidScheme as PS
+    using DimensionalData
+    pyr1 = Pyramid(rand(X(1:128),Y(1:128)), tilesize=16)
+    pyrcat = cat(pyr1, pyr1, dims=Dim{:new}([1,2]))
+    pyrcat3 = cat(pyr1, pyr1, pyr1, dims=Dim{:new}([1,2,3]))
+    pyr2 = Pyramid(rand(X(129:256), Y(1:128)), tilesize=16)
+    pyrcat2 = cat(pyr1, pyr2, dims=X)
+    for l in 1:3
+        @test PS.levels(pyrcat,l) == cat(PS.levels(pyr1, l), PS.levels(pyr1, l), dims=Dim{:new}([1,2]))
+        @test PS.levels(pyrcat3, l) == cat(PS.levels(pyr1, l), PS.levels(pyr1, l), PS.levels(pyr1, l), dims=Dim{:new}([1,2,3]))
+        @test PS.levels(pyrcat2, l) == cat(PS.levels(pyr1, l), PS.levels(pyr2, l), dims=X)
+    end
+end
+
 #=
 @testitem "Comparing zarr pyramid with tif pyramid" begin
     using PyramidScheme: PyramidScheme as PS