From b56101bb1c167db21770fe361ea103f6fe7b6712 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Tue, 2 Jul 2024 11:18:18 +0200 Subject: [PATCH 01/13] Add test for constructing a pyramid with additional dimension --- test/runtests.jl | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/runtests.jl b/test/runtests.jl index 114521b..3015cd8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -123,3 +123,19 @@ end end =# +@testitem "Building pyramid with additional dimension in memory" begin + # The aim of this test is to check whether we can build a pyramid from a data cube with an extra dimension. + # We will only build the pyramids on the spatial dimensions and keep the other dimensions as is. + using YAXArrays + using Zarr + using PyramidScheme + orgpath = joinpath(@__DIR__, "data", "world.zarr") + path = tempname() * ".zarr" + cp(orgpath, path) + c = Cube(path) + @time "Building world pyramid" pyr = Pyramid(c) + @test pyr isa Pyramid + @test length(dims(pyr)) == 3 + @test size(pyr.levels)[end] = (256,128,3) + +end \ No newline at end of file From 2831eb89ab675116ce0332e5dea0c303c4cad869 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Tue, 2 Jul 2024 11:19:27 +0200 Subject: [PATCH 02/13] Try to make fill_pyramids handle additional dimension This is currently work in progress. The output arrays need to be constructed accordingly and the additional dimensions need to be added to the windows. --- src/PyramidScheme.jl | 49 ++++++++++++++++++++++++++++++++------------ 1 file changed, 36 insertions(+), 13 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 10528aa..655e937 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -207,13 +207,40 @@ Fill the pyramids generated from the `data` with the aggregation function `func` `recursive` indicates whether higher tiles are computed from lower tiles or directly from the original data. This is an optimization which for functions like median might lead to misleading results. """ -function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, kwargs...) +function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose=false, outtype=:mem, kwargs...) + t = typeof(func(zeros(eltype(data), 2,2))) + n_level = compute_nlevels(data) + input_axes = pyramidedaxes(data) + nonpyramiddims = DD.otherdims(data, input_axes) + @show nonpyramiddims + if length(input_axes) != 2 + throw(ArgumentError("Expected two spatial dimensions got $input_axes")) + end + verbose && println("Constructing output arrays") + spatialsize = size(data)[collect(DD.dimnum(data, input_axes))] + pyramid_sizes = [ceil.(Int, spatialsize ./ 2^i) for i in 1:n_level] + allsizes = [spatialsize..., [1 for o in nonpyramiddims]...] + sizeperm = [DD.dimnum(data, input_axes)..., DD.dimnum(data, nonpyramiddims)...] + permute!(allsizes, sizeperm) + @show allsizes + outputs = if outtype == :zarr + [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] + elseif outtype == :mem + outmin = output_arrays(pyramid_sizes, t) + else + throw(ArgumentError("Output type not valied got $outtype expected :mem or :zarr")) + end + + verbose && println("Start computation") n_level = length(outputs) + @show typeof(data) + @show n_level + @show size.(outputs) pixel_base_size = 2^n_level pyramid_sizes = size.(outputs) tmp_sizes = [ceil(Int,pixel_base_size / 2^i) for i in 1:n_level] - - ia = InputArray(data, windows = arraywindows(size(data),pixel_base_size)) + windows = arraywindows(allsizes,pixel_base_size) + ia = InputArray(data;windows) oa = ntuple(i->create_outwindows(pyramid_sizes[i],windows = arraywindows(pyramid_sizes[i],tmp_sizes[i])),n_level) @@ -255,6 +282,8 @@ Construct a list of `RegularWindows` for the size list in `s` for windows `w`. ?? """ function arraywindows(s,w) + @show s + @show w map(s) do l RegularWindows(1,l,window=w) end @@ -308,6 +337,9 @@ Union of Dimensions which are assumed to be in space and are therefore used in t """ SpatialDim = Union{DD.Dimensions.XDim, DD.Dimensions.YDim} +pyramidedaxes(input) = filter(x-> x isa SpatialDim, DD.dims(input)) + + """ buildpyramids(path; resampling_method=mean) Build the pyramids for the zarr dataset at `path` and write the pyramid layers into the zarr folder. @@ -324,15 +356,6 @@ function buildpyramids(path; resampling_method=mean, recursive=true, runner=Loca # 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))) - n_level = compute_nlevels(org) - input_axes = filter(x-> x isa SpatialDim, DD.dims(org)) - if length(input_axes) != 2 - throw(ArgumentError("Expected two spatial dimensions got $input_axes")) - end - verbose && println("Constructing output arrays") - outarrs = [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] - verbose && println("Start computation") fill_pyramids(org, outarrs, resampling_method, recursive;runner) pyraxs = [agg_axis.(input_axes, 2^n) for n in 1:n_level] pyrlevels = DD.DimArray.(outarrs, pyraxs) @@ -376,7 +399,7 @@ 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) - input_axes = DD.dims(ras) + input_axes = pyramidedaxes(ras) n_level = compute_nlevels(ras) if iszero(n_level) @info "Array is smaller than the tilesize no pyramids are computed" From 6fb9b1c50551bf94192eff5cf6a30479e5e7563a Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 5 Aug 2024 16:42:34 +0200 Subject: [PATCH 03/13] Make compute levels take a size tuple instead of the full dataset This makes it easier to call it only for the pyramided axes. --- src/PyramidScheme.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 655e937..ed7a101 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -209,7 +209,7 @@ This is an optimization which for functions like median might lead to misleading """ function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose=false, outtype=:mem, kwargs...) t = typeof(func(zeros(eltype(data), 2,2))) - n_level = compute_nlevels(data) + n_level = compute_nlevels(size(data)) input_axes = pyramidedaxes(data) nonpyramiddims = DD.otherdims(data, input_axes) @show nonpyramiddims @@ -293,9 +293,9 @@ end """ compute_nlevels(data, tilesize=1024) -Compute the number of levels for the aggregation based on the size of `data`. +Compute the number of levels for the aggregation based on the tuple `datasize`. """ -compute_nlevels(data, tilesize=256) = max(0,ceil(Int,log2(maximum(size(data))/tilesize))) +compute_nlevels(datasize, tilesize=256) = max(0,ceil(Int,log2(maximum(datasize)/tilesize))) function agg_axis(d,n) # TODO this might be problematic for explicitly set axes From daae85838aea9ac0272e79936a101a62c82267de Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 5 Aug 2024 16:45:50 +0200 Subject: [PATCH 04/13] Make comparison in test on the size of the ouput pyramid --- test/runtests.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 3015cd8..827d423 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -136,6 +136,6 @@ end @time "Building world pyramid" pyr = Pyramid(c) @test pyr isa Pyramid @test length(dims(pyr)) == 3 - @test size(pyr.levels)[end] = (256,128,3) + @test size(pyr.levels)[end] == (256,128,3) end \ No newline at end of file From 3e8da5c90d67808ffac0e0b0702b8f84dfc08edd Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 5 Aug 2024 16:47:39 +0200 Subject: [PATCH 05/13] Try to get windows use unpyramided dims in GMWOP --- src/PyramidScheme.jl | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index ed7a101..6cc3b82 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -223,13 +223,13 @@ function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose= sizeperm = [DD.dimnum(data, input_axes)..., DD.dimnum(data, nonpyramiddims)...] permute!(allsizes, sizeperm) @show allsizes - outputs = if outtype == :zarr - [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] - elseif outtype == :mem - outmin = output_arrays(pyramid_sizes, t) - else - throw(ArgumentError("Output type not valied got $outtype expected :mem or :zarr")) - end + #outputs = if outtype == :zarr + # [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] + #elseif outtype == :mem + # output_arrays(pyramid_sizes, t) + #else + # throw(ArgumentError("Output type not valid got $outtype, but expected :mem or :zarr")) + #end verbose && println("Start computation") n_level = length(outputs) @@ -238,9 +238,11 @@ function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose= @show size.(outputs) pixel_base_size = 2^n_level pyramid_sizes = size.(outputs) + # What is tmp_sizes supposed to be? + # Does this have to be tmp_sizes = [ceil(Int,pixel_base_size / 2^i) for i in 1:n_level] windows = arraywindows(allsizes,pixel_base_size) - ia = InputArray(data;windows) + ia = InputArray(data;windows = Tuple(windows)) oa = ntuple(i->create_outwindows(pyramid_sizes[i],windows = arraywindows(pyramid_sizes[i],tmp_sizes[i])),n_level) @@ -284,9 +286,10 @@ Construct a list of `RegularWindows` for the size list in `s` for windows `w`. function arraywindows(s,w) @show s @show w - map(s) do l + windows = map(s) do l RegularWindows(1,l,window=w) end + windows end @@ -400,15 +403,16 @@ This returns the data of the pyramids and the dimension values of the aggregated """ function getpyramids(reducefunc, ras;recursive=true) input_axes = pyramidedaxes(ras) - n_level = compute_nlevels(ras) + n_level = compute_nlevels(length.(input_axes)) if iszero(n_level) @info "Array is smaller than the tilesize no pyramids are computed" [ras], [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] - + pyramid_axes = [map(x-> in(x, input_axes) ? agg_axis(x, 2^l) : x , DD.dims(ras)) for l in 1:n_level] + pyramid_sizes = size.(pyramid_axes) + @show pyramid_sizes outmin = output_arrays(pyramid_sizes, Float32) + @show size.(outmin) fill_pyramids(ras,outmin,reducefunc,recursive; threaded=true) outmin, pyramid_axes From b6c6bbbfd84735217e21eeef641725e3189869ac Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Wed, 7 Aug 2024 17:51:13 +0200 Subject: [PATCH 06/13] Make DimArray before filling it This makes it easier to handle the axes position in the fill_pyramids function --- src/PyramidScheme.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 6cc3b82..c1e6781 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -45,8 +45,7 @@ struct Pyramid{T,N,D,A,B<:DD.AbstractDimArray{T,N,D,A},L, Me} <: DD.AbstractDimA end function Pyramid(data::DD.AbstractDimArray) - pyrdata, pyraxs = getpyramids(mean ∘ skipmissing, data, recursive=false) - levels = DD.DimArray.(pyrdata, pyraxs) + levels = getpyramids(mean ∘ skipmissing, data, recursive=false) meta = Dict(deepcopy(DD.metadata(data))) push!(meta, "resampling_method" => "mean_skipmissing") Pyramid(data, levels, meta) @@ -413,9 +412,10 @@ function getpyramids(reducefunc, ras;recursive=true) @show pyramid_sizes outmin = output_arrays(pyramid_sizes, Float32) @show size.(outmin) - fill_pyramids(ras,outmin,reducefunc,recursive; threaded=true) + levels = DD.DimArray.(outmin, pyramid_axes) + fill_pyramids(ras,levels,reducefunc,recursive; threaded=true) - outmin, pyramid_axes + levels end """ From d80b20cc0bad8b392182067153d1389360778439 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Wed, 7 Aug 2024 17:52:30 +0200 Subject: [PATCH 07/13] Try to make windows correct size for non pyramided dims --- src/PyramidScheme.jl | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index c1e6781..cef3122 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -218,9 +218,9 @@ function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose= verbose && println("Constructing output arrays") spatialsize = size(data)[collect(DD.dimnum(data, input_axes))] pyramid_sizes = [ceil.(Int, spatialsize ./ 2^i) for i in 1:n_level] - allsizes = [spatialsize..., [1 for o in nonpyramiddims]...] + allsizes = [spatialsize...,] sizeperm = [DD.dimnum(data, input_axes)..., DD.dimnum(data, nonpyramiddims)...] - permute!(allsizes, sizeperm) + #permute!(allsizes, sizeperm) @show allsizes #outputs = if outtype == :zarr # [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] @@ -240,17 +240,32 @@ function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose= # What is tmp_sizes supposed to be? # Does this have to be tmp_sizes = [ceil(Int,pixel_base_size / 2^i) for i in 1:n_level] - windows = arraywindows(allsizes,pixel_base_size) + windows = Any[Base.OneTo.(size(data))...] + #pseudocode + windows[DD.dimnum(data, XDim)] = RegularWindows(1,size(data, XDim),window=pixel_base_size) + windows[DD.dimnum(data, YDim)] = RegularWindows(1,size(data, YDim),window=pixel_base_size) +@show size.(windows) + ia = InputArray(data;windows = Tuple(windows)) + # Replace + function outwindows(i, outputs, tmpsize) + outwins = Any[Base.OneTo.(size(data))...] + outwins[DD.dimnum(outputs[i], XDim)] = RegularWindows(1,size(outputs[i], XDim),window=pixel_base_size) + outwins[DD.dimnum(outputs[i], YDim)] = RegularWindows(1,size(outputs[i], YDim),window=pixel_base_size) + Tuple(outwins) + end - oa = ntuple(i->create_outwindows(pyramid_sizes[i],windows = arraywindows(pyramid_sizes[i],tmp_sizes[i])),n_level) + for i in 1:1 + @show outwindows(1, outputs, tmp_sizes[i]) + end + oa = ntuple(i->create_outwindows(pyramid_sizes[i],windows = outwindows(i,outputs, tmp_sizes[i])),n_level) func = DiskArrayEngine.create_userfunction(gen_pyr,ntuple(_->eltype(first(outputs)),length(outputs));is_mutating=true,kwargs = (;recursive),args = (func,)) op = GMDWop((ia,), oa, func) lr = DiskArrayEngine.optimize_loopranges(op,5e8,tol_low=0.2,tol_high=0.05,max_order=2) - r = runner(op,lr,outputs;kwargs...) + r = runner(op,lr,getproperty.(outputs, :data);kwargs...) run(r) end From 1731f0f47a671846d01e0ce535f52ff8a8c215d6 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Thu, 8 Aug 2024 11:50:07 +0200 Subject: [PATCH 08/13] Get pyramid building with additional dims working This sets the chunking to 256 in all dimensions. Just because it was easier than having to find out what are the pyramided dimensions. --- src/PyramidScheme.jl | 9 +++------ test/runtests.jl | 2 +- 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index cef3122..67886c5 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -250,14 +250,11 @@ function fill_pyramids(data, outputs,func,recursive;runner=LocalRunner, verbose= # Replace function outwindows(i, outputs, tmpsize) outwins = Any[Base.OneTo.(size(data))...] - outwins[DD.dimnum(outputs[i], XDim)] = RegularWindows(1,size(outputs[i], XDim),window=pixel_base_size) - outwins[DD.dimnum(outputs[i], YDim)] = RegularWindows(1,size(outputs[i], YDim),window=pixel_base_size) + outwins[DD.dimnum(outputs[i], XDim)] = RegularWindows(1,size(outputs[i], XDim),window=tmpsize) + outwins[DD.dimnum(outputs[i], YDim)] = RegularWindows(1,size(outputs[i], YDim),window=tmpsize) Tuple(outwins) end - for i in 1:1 - @show outwindows(1, outputs, tmp_sizes[i]) - end oa = ntuple(i->create_outwindows(pyramid_sizes[i],windows = outwindows(i,outputs, tmp_sizes[i])),n_level) func = DiskArrayEngine.create_userfunction(gen_pyr,ntuple(_->eltype(first(outputs)),length(outputs));is_mutating=true,kwargs = (;recursive),args = (func,)) @@ -331,7 +328,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=ntuple(i->min(256, s[i]), length(s)),fill_value=zero(t)) else zeros(t,s...) end diff --git a/test/runtests.jl b/test/runtests.jl index 827d423..c1fe1bc 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -136,6 +136,6 @@ end @time "Building world pyramid" pyr = Pyramid(c) @test pyr isa Pyramid @test length(dims(pyr)) == 3 - @test size(pyr.levels)[end] == (256,128,3) + @test size(pyr.levels[end]) == (256,128,3) end \ No newline at end of file From 0e0f5cf44a9b5ad5be83dbac3cd86fa028235bee Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Thu, 8 Aug 2024 14:40:40 +0200 Subject: [PATCH 09/13] Fix compute_nlevels tests after change of input --- test/runtests.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index c1fe1bc..790d79a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -45,9 +45,9 @@ end end @testitem "Helper functions" begin using PyramidScheme: PyramidScheme as PS - @test PS.compute_nlevels(zeros(1000)) == 2 - @test PS.compute_nlevels(zeros(1000,1025)) == 3 - @test PS.compute_nlevels(zeros(10000, 8000)) == 6 + @test PS.compute_nlevels((1000,)) == 2 + @test PS.compute_nlevels((1000,1025)) == 3 + @test PS.compute_nlevels((10000, 8000)) == 6 end @testitem "ArchGDAL Loading of Geotiff Overviews" begin @@ -129,6 +129,7 @@ end using YAXArrays using Zarr using PyramidScheme + using DimensionalData orgpath = joinpath(@__DIR__, "data", "world.zarr") path = tempname() * ".zarr" cp(orgpath, path) From 3ea60076eca3fc4dd8f96d43730f43ecb9fd384e Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Thu, 8 Aug 2024 14:41:22 +0200 Subject: [PATCH 10/13] Fix inplace building for zarr --- src/PyramidScheme.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 67886c5..2895e50 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -305,7 +305,7 @@ end """ - compute_nlevels(data, tilesize=1024) + compute_nlevels(data, tilesize=256) Compute the number of levels for the aggregation based on the tuple `datasize`. """ @@ -370,6 +370,15 @@ function buildpyramids(path; resampling_method=mean, recursive=true, runner=Loca # 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))) + n_level = compute_nlevels(org) + input_axes = filter(x-> x isa SpatialDim, DD.dims(org)) + if length(input_axes) != 2 + throw(ArgumentError("Expected two spatial dimensions got $input_axes")) + end + verbose && println("Constructing output arrays") + outarrs = [output_zarr(n, input_axes, t, joinpath(path, string(n))) for n in 1:n_level] + verbose && println("Start computation") fill_pyramids(org, outarrs, resampling_method, recursive;runner) pyraxs = [agg_axis.(input_axes, 2^n) for n in 1:n_level] pyrlevels = DD.DimArray.(outarrs, pyraxs) From 237db3d5bf4759a7b1b3ad882ee60a1e4c7c0cd8 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Wed, 14 Aug 2024 16:15:23 +0200 Subject: [PATCH 11/13] Rename presentation quarto in docs --- .gitignore | 4 +++- docs/Manifest.toml | 2 +- docs/{juliacon.qmd => presentation.qmd} | 27 +++++++++++++++++++------ 3 files changed, 25 insertions(+), 8 deletions(-) rename docs/{juliacon.qmd => presentation.qmd} (89%) diff --git a/.gitignore b/.gitignore index b202fac..723dd62 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,4 @@ /Manifest.toml -test/data \ No newline at end of file +test/data +docs/.quarto +docs/presentation_files \ No newline at end of file diff --git a/docs/Manifest.toml b/docs/Manifest.toml index 670ab53..0370f03 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -1785,7 +1785,7 @@ uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d" version = "1.2.0" [[deps.PyramidScheme]] -deps = ["Colors", "DimensionalData", "DiskArrayEngine", "DiskArrayTools", "DiskArrays", "Extents", "FillArrays", "GeoInterface", "Makie", "MakieCore", "Observables", "OffsetArrays", "Proj", "Statistics", "WGLMakie", "YAXArrayBase", "YAXArrays", "Zarr"] +deps = ["Colors", "DimensionalData", "DiskArrayEngine", "DiskArrayTools", "DiskArrays", "Extents", "FillArrays", "GeoInterface", "Makie", "MakieCore", "Observables", "OffsetArrays", "Proj", "Statistics", "YAXArrayBase", "YAXArrays", "Zarr"] path = ".." uuid = "ec211b67-1c2c-4319-878f-eaee078ee145" version = "0.1.0" diff --git a/docs/juliacon.qmd b/docs/presentation.qmd similarity index 89% rename from docs/juliacon.qmd rename to docs/presentation.qmd index 99acaff..9051d04 100644 --- a/docs/juliacon.qmd +++ b/docs/presentation.qmd @@ -1,12 +1,23 @@ --- -engine: julia title: "PyramidScheme.jl \n interactively plotting large raster data" author: "Felix Cremer, Fabian Gans" format: revealjs +engine: julia execute: echo: true --- +## Why use pyramids? + +julia> agbdiffbase = agb2020 .- agb2017 + +╭───────────────────────────────────╮ + +│ 405000×157500 YAXArray{Float32,2} │ + +- julia> savecube(agbdiffbase, tempname()*".zarr") +Progress: 1%|▍ | ETA: 3:15:43 + ## What are pyramids? ```{julia} @@ -65,11 +76,12 @@ plotpyramid(europyr) ``` + ## Compute the pyramids for your data - In Memory data ```{julia} -#| eval: false +#| eval: false using Rasters using RasterDataSources ras = Raster(WorldClim{Elevation},:elev, res="30s", lazy=true) @@ -101,13 +113,13 @@ agb2020 = replacenan.(Pyramid(joinpath(@__DIR__,"..", "test/data/ESACCI-BIOMASS- # Interactive visualisation ```{julia} -#| eval: false +#| eval: false plot(agb2020, colormap=:speed, colorscale=sqrt) ``` # Computations on Pyramids ```{julia} -#| eval: false +#| eval: false replacenan(data) = data <= 0 ? NaN32 : Float32(data) agb2020 = replacenan.(Pyramid(joinpath(@__DIR__,"..", "test/data/ESACCI-BIOMASS-L4-AGB-MERGED-100m-2020-fv4.0.zarr/"))) agb2017 = replacenan.(Pyramid(joinpath(@__DIR__,"..", "test/data/ESACCI-BIOMASS-L4-AGB-MERGED-100m-2017-fv4.0.zarr/"))) @@ -118,16 +130,19 @@ plot(agbdiff, colormap=:bukavu, colorrange=(-200, 200)) # Outlook -- pyramids in multiple dimensions +- pyramids in multiple dimensions (WIP) - Understanding more formats +- Expose data as WMTS/WMS for online usage + - Integrate DiskArrayEngine computations - Open it up to more use cases +## Example of cloud access ```{julia} -#| eval:false +#| eval: false using YAXArrays, PyramidScheme using GLMakie import PyramidScheme as PS From 594dd9aac5572b9f35f51550b4569270437814a0 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 23 Sep 2024 16:16:59 +0200 Subject: [PATCH 12/13] Fix Pyramid Constructor from DimArray --- src/PyramidScheme.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 93446b2..55fbe72 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -45,8 +45,7 @@ struct Pyramid{T,N,D,A,B<:DD.AbstractDimArray{T,N,D,A},L, Me} <: DD.AbstractDimA end function Pyramid(data::DD.AbstractDimArray; resampling_method=mean ∘ skipmissing, kwargs...) - pyrdata, pyraxs = getpyramids(resampling_method, data; kwargs...) - levels = DD.DimArray.(pyrdata, pyraxs) + levels = getpyramids(resampling_method, data; kwargs...) meta = Dict(deepcopy(DD.metadata(data))) push!(meta, "resampling_method" => "mean_skipmissing") Pyramid(data, levels, meta) From 93f2b9ca66bb3e105689a249597334a477093b3d Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Tue, 24 Sep 2024 15:03:55 +0200 Subject: [PATCH 13/13] Determine eltype of output from eltype of the input --- src/PyramidScheme.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 55fbe72..f45f79d 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -338,9 +338,9 @@ function DiskArrayEngine.engine(dimarr::DD.AbstractDimArray) end DiskArrayEngine.engine(pyr::Pyramid) = Pyramid(engine(parent(pyr)), engine.(pyr.levels), DD.metadata(pyr)) """ - output_arrays(pyramid_sizes) + output_arrays(pyramid_sizes, T) -Create the output arrays for the given `pyramid_sizes` +Create the output arrays for the given `pyramid_sizes` with the element type T. """ output_arrays(pyramid_sizes, T) = [gen_output(T,p) for p in pyramid_sizes] @@ -430,7 +430,9 @@ function getpyramids(reducefunc, ras;recursive=true) pyramid_axes = [map(x-> in(x, input_axes) ? agg_axis(x, 2^l) : x , DD.dims(ras)) for l in 1:n_level] pyramid_sizes = size.(pyramid_axes) @show pyramid_sizes - outmin = output_arrays(pyramid_sizes, Float32) + t = typeof(reducefunc(zeros(eltype(ras), 2,2))) + + outmin = output_arrays(pyramid_sizes, t) @show size.(outmin) levels = DD.DimArray.(outmin, pyramid_axes) fill_pyramids(ras,levels,reducefunc,recursive; threaded=true)