Merge branch 'main' into fg/multidim

Author: Felix Cremer

CHANGELOG.md

@@ -0,0 +1,29 @@
+
+File filter
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+Links are updated with Changelog.jl using the command:
+
+```julia
+Changelog.generate(
+    Changelog.CommonMark(),             # output type
+    "CHANGELOG.md";                     # input and output file
+    repo = "JuliaDataCubes/PyramidScheme.jl", # default repository for links
+)
+```
+
+## [0.1.2]
+
+### Added 
+
+- Implemented proper MakieRecipies this should make Observable(::Pyramid) work
+- Added `cat` of Pyramids to concatenate on all levels
+
+### Fixed
+
+- Improved computation of the dimensions of the aggregated levels

Project.toml

@@ -1,7 +1,7 @@
 name = "PyramidScheme"
 uuid = "ec211b67-1c2c-4319-878f-eaee078ee145"
 authors = ["Felix Cremer <[email protected]> and contributors"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
@@ -42,7 +42,7 @@ Extents = "0.1.3"
 FillArrays = "1.11"
 Flatten = "0.4.3"
 GeoInterface = "1.3.4"
-Makie = "0.21.2, 0.22, 0.23, 0.24"
+Makie = "0.24"
 OffsetArrays = "1.14"
 Proj = "1.7.1"
 Rasters = "0.11.8, 0.12, 0.13, 0.14"

README.md

@@ -1,4 +1,6 @@
 # PyramidScheme
+[![codecov](https://codecov.io/gh/JuliaDataCubes/PyramidScheme.jl/graph/badge.svg?token=EBhwEBXlnx)](https://codecov.io/gh/JuliaDataCubes/PyramidScheme.jl)
+
 
 PyramidScheme.jl is a package to easily and efficiently compute pyramids of a given datacube which might be larger than RAM.
 It uses DiskArrayEngine.jl as the computational backend.

ext/PyramidSchemeMakieExt.jl

@@ -1,62 +1,150 @@
 module PyramidSchemeMakieExt
 using Makie: Axis, Colorbar, DataAspect, Figure, FigureAxisPlot, Observable, Relative
-using Makie: on, heatmap!, image!
-import Makie: plot, plot!
+using Makie: on, heatmap!, map!
+using Makie.ComputePipeline: add_input!
+using Makie
+
 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, size=(1155, 1043), kwargs...)
-    #This should be converted into a proper recipe for Makie but this would depend on a pyramid type.
-    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...)
-    if colorbar
-        Colorbar(fig[1,2], hmap, height = Relative(3.5 / 4))
+
+# hacks to get around DD hacks that get around Makie issues
+for p in (Heatmap, Image, Contour, Contourf, Contour3d, Spy, Surface) 
+    f = Makie.plotkey(p)
+    @eval begin
+        function Makie.$f(A::Pyramid; kwargs...)
+            invoke(Makie.$f, Tuple{AbstractMatrix{<: Any}}, A; kwargs...)
+        end
+        function Makie.$f(A::Observable{<: Pyramid}; kwargs...)
+            invoke(Makie.$f, Tuple{<:Union{<:AbstractMatrix{<: Any}, <:Observable{<: AbstractMatrix{<:Any}}}}, A; kwargs...)
+        end
+        Makie.expand_dimensions(::Type{$p}, p::Pyramid) = (p,)
+        Makie.convert_arguments(::Type{$p}, p::Pyramid) = (p,)
     end
-    ax.autolimitaspect = 1
-    FigureAxisPlot(fig, ax, hmap)
 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)
-    on(ax.scene.viewport) do viewport
-        limext = extent(ax.finallimits[])
-
-        datalimit = switchkeys(limext, rasext)
-        data.val = miss2nan.(selectlevel(pyramid, datalimit, target_imsize=viewport.widths))
-        notify(data)
-    end
-    on(ax.finallimits) do limits
-        limext = extent(limits)
-        # Compute limit in raster projection
-        #trans = Proj.Transformation(plotcrs, rastercrs, always_xy=true)
-        #datalimit = trans_bounds(trans, limext)
-        datalimit = switchkeys(limext, rasext)
-        if intersects(rasdataext, limext)
-            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 = miss2nan.(rasdata)
+struct PyramidConversion <: Makie.ConversionTrait end
+
+function Makie.conversion_trait(::Type{<:Heatmap}, ::Pyramid)
+    return PyramidConversion()
+end
+
+function Makie.types_for_plot_arguments(::Type{<:Heatmap}, ::PyramidConversion)
+    return Tuple{<:Pyramid}
+end
+
+Makie.expand_dimensions(::PyramidConversion, p::Pyramid) = (p,)
+Makie.convert_arguments(::PyramidConversion, p::Pyramid) = (p,)
+
+
+
+Makie.plottype(::Pyramid) = Makie.Heatmap
+
+function Makie.plot!(plot::Heatmap{<: Tuple{<: Pyramid}})
+    #=
+    Go from a relative space rectangle 
+    to the rectangle in data space and pixel space,
+    thus getting `ax.finallimits` and `ax.viewport`
+    respectively.
+
+    This is a bit arcane but probably the simplest thing in the long run.
+    =#
+    inputpositions = [Point2f(0, 0), Point2f(1, 1)]
+    add_input!(plot.attributes, :__pyramid_input_positions, inputpositions)
+    Makie.register_positions_projected!(
+        plot; input_space = :relative, output_space = :space,
+        input_name = :__pyramid_input_positions,
+        output_name = :__pyramid_dataspace_positions,
+    )
+    Makie.register_positions_projected!(
+        plot; input_space = :relative, output_space = :pixel,
+        input_name = :__pyramid_input_positions,
+        output_name = :__pyramid_pixelspace_positions,
+    )
+
+    data = Observable{DD.AbstractDimMatrix}(levels(plot.arg1[])[end-2][:, :])
+    onany(plot, plot.arg1, plot.__pyramid_dataspace_positions, plot.__pyramid_pixelspace_positions) do pyramid, datapos, pixelpos
+        xyext = values.(extent(pyramid, (XDim, YDim)))
+        xval, yval = first(xyext), last(xyext)
+        pyramid_data_ext = Extent(X=xval, Y=yval)
+        pyramid_ext = extent(pyramid)
+
+        data_limits_ext = Extent(X = extrema(first, datapos), Y = extrema(x -> x[2], datapos))
+        pixel_widths = Point2f(abs.(pixelpos[2] .- pixelpos[1]))
+
+        datalimit = switchkeys(data_limits_ext, pyramid_ext)
+
+        if intersects(pyramid_data_ext, data_limits_ext)
+            data[] = miss2nan.(
+                selectlevel(pyramid, datalimit, target_imsize = pixel_widths)
+            )
         end
-        notify(data)
+        nothing
     end
-    #@show typeof(data)
-    hmap = heatmap!(ax, data; interpolate=interp, kwargs...)
+
+    heatmap!(plot, plot.attributes, data)
 end
+
+function Makie.data_limits(p::Heatmap{<: Tuple{<: Pyramid}})
+    extX, extY = extent(p.arg1[], (XDim, YDim))
+    rect = Rect3f((extX[1], extY[1],0), (extX[2] - extX[1], extY[2] - extY[1], 0))
+    return rect
+end
+Makie.boundingbox(p::Heatmap{<: Tuple{<: Pyramid}}, space::Symbol = :data) = Makie.apply_transform_and_model(p, Makie.data_limits(p))
+
+
+# """
+#     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, size=(1155, 1043), kwargs...)
+#     #This should be converted into a proper recipe for Makie but this would depend on a pyramid type.
+#     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...)
+#     if colorbar
+#         Colorbar(fig[1,2], hmap, height = Relative(3.5 / 4))
+#     end
+#     ax.autolimitaspect = 1
+#     FigureAxisPlot(fig, ax, hmap)
+# 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)
+#     xyext = values.(extent(pyramid, (XDim, YDim)))
+#     xval = first(xyext)
+#     yval = last(xyext)
+#     rasdataext = Extent(X=xval, Y=yval)
+#     rasext = extent(pyramid)
+#     on(ax.scene.viewport) do viewport
+#         limext = extent(ax.finallimits[])
+
+#         datalimit = switchkeys(limext, rasext)
+#         data[] = miss2nan.(selectlevel(pyramid, datalimit, target_imsize=viewport.widths))
+#     end
+#     on(ax.finallimits) do limits
+#         limext = extent(limits)
+#         # Compute limit in raster projection
+#         #trans = Proj.Transformation(plotcrs, rastercrs, always_xy=true)
+#         #datalimit = trans_bounds(trans, limext)
+#         datalimit = switchkeys(limext, rasext)
+#         if intersects(rasdataext, limext)
+#             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 = miss2nan.(rasdata)
+#         end
+#         notify(data)
+#     end
+#     #@show typeof(data)
+#     hmap = heatmap!(ax, data; interpolate=interp, kwargs...)
+# end
 end

src/PyramidScheme.jl

@@ -124,10 +124,10 @@ end
     Pyramid(newbase, newlevels, DD.metadata(A))
 end
 
-function Base.map(f, A::Pyramid)
-    newbase = map(f, parent(A))
-    newlevels = [map(f, levels(A,i)) for i in 1:nlevels(A)]
-    Pyramid(newbase, newlevels, DD.metadata(A)) # This should handle metadata better.
+function Base.map(f, A1::Pyramid, As::Pyramid...)
+    newbase = map(f, parent(A1), parent.(As)...)
+    newlevels = [map(f, levels(A1,i), levels.(As, i)...) for i in 1:nlevels(A1)]
+    Pyramid(newbase, newlevels, DD.metadata(A1)) # This should handle metadata better.
 end
 
 function DD.show_after(io::IO, mime, A::Pyramid)
@@ -458,8 +458,8 @@ end
 
 
 
-xkey(keyext) = DD.dim2key(DD.dims(keyext, XDim))
-ykey(keyext) = DD.dim2key(DD.dims(keyext, YDim))
+xkey(keyext) = DD.name(DD.dims(keyext, XDim))
+ykey(keyext) = DD.name(DD.dims(keyext, YDim))
 #TODO write test, move to utils.jl
 
 """

test/runtests.jl

@@ -42,6 +42,9 @@ end
     @test size(subpyramid) == (10,10)
     @test parent(subpyramid) == data[1:10,1:10]
     fig, axis, h = plot(pyramid)
+    @test length(axis.scene.plots) == 1
+    plot!(axis, pyramid)
+    @test length(axis.scene.plots) == 2
 end
 
 #= building an RGB pyramid doesn't work, need to think more about it.
@@ -161,6 +164,23 @@ end
     end
 
 end
+@testitem "Map of pyramids" begin
+    using DimensionalData
+    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)
+    pyr1_neg = map(x-> x-1, pyr1)
+    @test all(all.(iszero, pyr1_neg.levels))
+    @test iszero(pyr1_neg.base)
+    pyr2_neg = map(x-> x-1, pyr2)
+    @test pyr2_neg.levels[1][1,1] == 3
+
+    pyrsum = map((x,y) -> x + y, pyr1, pyr2)
+    @test pyrsum[100,30] == 2
+    @test pyrsum.levels[1][10,10] == 5
+    @test pyrsum.levels[2][10,10] == 17
+end
+
+
 #=
 @testitem "Comparing zarr pyramid with tif pyramid" begin
     using PyramidScheme: PyramidScheme as PS