Skip to content

Commit fce9f83

Browse files
author
Felix Cremer
authored
Merge branch 'main' into fc/upyax
2 parents 91e3baf + 226a525 commit fce9f83

File tree

5 files changed

+228
-66
lines changed

5 files changed

+228
-66
lines changed

CHANGELOG.md

Lines changed: 29 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,29 @@
1+
2+
File filter
3+
# Changelog
4+
5+
All notable changes to this project will be documented in this file.
6+
7+
The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
8+
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
9+
10+
Links are updated with Changelog.jl using the command:
11+
12+
```julia
13+
Changelog.generate(
14+
Changelog.CommonMark(), # output type
15+
"CHANGELOG.md"; # input and output file
16+
repo = "JuliaDataCubes/PyramidScheme.jl", # default repository for links
17+
)
18+
```
19+
20+
## [0.1.2]
21+
22+
### Added
23+
24+
- Implemented proper MakieRecipies this should make Observable(::Pyramid) work
25+
- Added `cat` of Pyramids to concatenate on all levels
26+
27+
### Fixed
28+
29+
- Improved computation of the dimensions of the aggregated levels

Project.toml

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "PyramidScheme"
22
uuid = "ec211b67-1c2c-4319-878f-eaee078ee145"
33
authors = ["Felix Cremer <fcremer@bgc-jena.mpg.de> and contributors"]
4-
version = "0.1.1"
4+
version = "0.1.2"
55

66
[deps]
77
ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
@@ -32,7 +32,7 @@ PyramidSchemeMakieExt = "Makie"
3232
[compat]
3333
Aqua = "0.8"
3434
ArchGDAL = "0.10.4"
35-
CairoMakie = "0.12"
35+
CairoMakie = "0.12,0.13,0.14,0.15"
3636
Colors = "0.12"
3737
DimensionalData = "0.28, 0.29"
3838
DiskArrayEngine = "0.1.2, 0.2"
@@ -42,13 +42,13 @@ Extents = "0.1.3"
4242
FillArrays = "1.11"
4343
Flatten = "0.4.3"
4444
GeoInterface = "1.3.4"
45-
Makie = "0.21.2"
45+
Makie = "0.24"
4646
OffsetArrays = "1.14"
4747
Proj = "1.7.1"
48-
Rasters = "0.11.8, 0.12, 0.13"
48+
Rasters = "0.11.8, 0.12, 0.13, 0.14"
4949
Statistics = "1.10"
5050
Test = "1.10"
51-
TestItemRunner = "0.2.3"
51+
TestItemRunner = "1"
5252
YAXArrayBase = "0.6.1, 0.7"
5353
YAXArrays = "0.5.6, 0.6, 0.7"
5454
Zarr = "0.9.3"

ext/PyramidSchemeMakieExt.jl

Lines changed: 137 additions & 50 deletions
Original file line numberDiff line numberDiff line change
@@ -1,63 +1,150 @@
11
module PyramidSchemeMakieExt
22
using Makie: Axis, Colorbar, DataAspect, Figure, FigureAxisPlot, Observable, Relative
3-
using Makie: on, heatmap!, image!
4-
import Makie: plot, plot!
3+
using Makie: on, heatmap!, map!
4+
using Makie.ComputePipeline: add_input!
5+
using Makie
6+
57
using PyramidScheme: Pyramid, switchkeys, levels, selectlevel, xkey, ykey
68
using DimensionalData: DimensionalData as DD
79
using DimensionalData.Dimensions: XDim, YDim
810
using Extents: Extent, extent, intersects
9-
"""
10-
plot(pyramids)
11-
Plot a Pyramid.
12-
This will plot the coarsest resolution level at the beginning and will plot higher resolutions after zooming into the plot.
13-
This is expected to be used with interactive Makie backends.
14-
"""
15-
function plot(pyramid::Pyramid;colorbar=true, kwargs...)
16-
#This should be converted into a proper recipe for Makie but this would depend on a pyramid type.
17-
fig = Figure()
18-
lon, lat = DD.dims(DD.parent(pyramid))
19-
ax = Axis(fig[1,1], limits=(extrema(lon), extrema(lat)), aspect=DataAspect())
20-
hmap = plot!(ax, pyramid;kwargs...)
21-
if colorbar
22-
Colorbar(fig[1,2], hmap, height = Relative(3.5 / 4))
11+
miss2nan(x) = ismissing(x) ? NaN : x
12+
13+
# hacks to get around DD hacks that get around Makie issues
14+
for p in (Heatmap, Image, Contour, Contourf, Contour3d, Spy, Surface)
15+
f = Makie.plotkey(p)
16+
@eval begin
17+
function Makie.$f(A::Pyramid; kwargs...)
18+
invoke(Makie.$f, Tuple{AbstractMatrix{<: Any}}, A; kwargs...)
19+
end
20+
function Makie.$f(A::Observable{<: Pyramid}; kwargs...)
21+
invoke(Makie.$f, Tuple{<:Union{<:AbstractMatrix{<: Any}, <:Observable{<: AbstractMatrix{<:Any}}}}, A; kwargs...)
22+
end
23+
Makie.expand_dimensions(::Type{$p}, p::Pyramid) = (p,)
24+
Makie.convert_arguments(::Type{$p}, p::Pyramid) = (p,)
2325
end
24-
ax.autolimitaspect = 1
25-
FigureAxisPlot(fig, ax, hmap)
2626
end
2727

28-
function plot!(ax, pyramid::Pyramid;interp=false, kwargs...)#; rastercrs=crs(parent(pyramid)),plotcrs=EPSG(3857), kwargs...)
29-
tip = levels(pyramid)[end-2][:,:]
30-
#@show typeof(tip)
31-
data = Observable{DD.AbstractDimMatrix}(tip)
32-
xval = only(values(extent(pyramid, XDim)))
33-
yval = only(values(extent(pyramid, YDim)))
34-
rasdataext = Extent(X=xval, Y=yval)
35-
rasext = extent(pyramid)
36-
xk = xkey(rasext)
37-
yk = ykey(rasext)
38-
on(ax.scene.viewport) do viewport
39-
limext = extent(ax.finallimits[])
40-
41-
datalimit = switchkeys(limext, rasext)
42-
43-
data.val = selectlevel(pyramid, datalimit, target_imsize=viewport.widths)
44-
notify(data)
45-
end
46-
on(ax.finallimits) do limits
47-
limext = extent(limits)
48-
# Compute limit in raster projection
49-
#trans = Proj.Transformation(plotcrs, rastercrs, always_xy=true)
50-
#datalimit = trans_bounds(trans, limext)
51-
datalimit = switchkeys(limext, rasext)
52-
if intersects(rasdataext, limext)
53-
rasdata = selectlevel(pyramid, datalimit, target_imsize=ax.scene.viewport[].widths)
54-
# Project selected data to plotcrs
55-
#data.val = Rasters.resample(rasdata, crs=plotcrs, method=:bilinear )
56-
data.val = rasdata
28+
struct PyramidConversion <: Makie.ConversionTrait end
29+
30+
function Makie.conversion_trait(::Type{<:Heatmap}, ::Pyramid)
31+
return PyramidConversion()
32+
end
33+
34+
function Makie.types_for_plot_arguments(::Type{<:Heatmap}, ::PyramidConversion)
35+
return Tuple{<:Pyramid}
36+
end
37+
38+
Makie.expand_dimensions(::PyramidConversion, p::Pyramid) = (p,)
39+
Makie.convert_arguments(::PyramidConversion, p::Pyramid) = (p,)
40+
41+
42+
43+
Makie.plottype(::Pyramid) = Makie.Heatmap
44+
45+
function Makie.plot!(plot::Heatmap{<: Tuple{<: Pyramid}})
46+
#=
47+
Go from a relative space rectangle
48+
to the rectangle in data space and pixel space,
49+
thus getting `ax.finallimits` and `ax.viewport`
50+
respectively.
51+
52+
This is a bit arcane but probably the simplest thing in the long run.
53+
=#
54+
inputpositions = [Point2f(0, 0), Point2f(1, 1)]
55+
add_input!(plot.attributes, :__pyramid_input_positions, inputpositions)
56+
Makie.register_positions_projected!(
57+
plot; input_space = :relative, output_space = :space,
58+
input_name = :__pyramid_input_positions,
59+
output_name = :__pyramid_dataspace_positions,
60+
)
61+
Makie.register_positions_projected!(
62+
plot; input_space = :relative, output_space = :pixel,
63+
input_name = :__pyramid_input_positions,
64+
output_name = :__pyramid_pixelspace_positions,
65+
)
66+
67+
data = Observable{DD.AbstractDimMatrix}(levels(plot.arg1[])[end-2][:, :])
68+
onany(plot, plot.arg1, plot.__pyramid_dataspace_positions, plot.__pyramid_pixelspace_positions) do pyramid, datapos, pixelpos
69+
xyext = values.(extent(pyramid, (XDim, YDim)))
70+
xval, yval = first(xyext), last(xyext)
71+
pyramid_data_ext = Extent(X=xval, Y=yval)
72+
pyramid_ext = extent(pyramid)
73+
74+
data_limits_ext = Extent(X = extrema(first, datapos), Y = extrema(x -> x[2], datapos))
75+
pixel_widths = Point2f(abs.(pixelpos[2] .- pixelpos[1]))
76+
77+
datalimit = switchkeys(data_limits_ext, pyramid_ext)
78+
79+
if intersects(pyramid_data_ext, data_limits_ext)
80+
data[] = miss2nan.(
81+
selectlevel(pyramid, datalimit, target_imsize = pixel_widths)
82+
)
5783
end
58-
notify(data)
84+
nothing
5985
end
60-
#@show typeof(data)
61-
hmap = heatmap!(ax, data; interpolate=interp, kwargs...)
86+
87+
heatmap!(plot, plot.attributes, data)
6288
end
89+
90+
function Makie.data_limits(p::Heatmap{<: Tuple{<: Pyramid}})
91+
extX, extY = extent(p.arg1[], (XDim, YDim))
92+
rect = Rect3f((extX[1], extY[1],0), (extX[2] - extX[1], extY[2] - extY[1], 0))
93+
return rect
94+
end
95+
Makie.boundingbox(p::Heatmap{<: Tuple{<: Pyramid}}, space::Symbol = :data) = Makie.apply_transform_and_model(p, Makie.data_limits(p))
96+
97+
98+
# """
99+
# plot(pyramids)
100+
# Plot a Pyramid.
101+
# This will plot the coarsest resolution level at the beginning and will plot higher resolutions after zooming into the plot.
102+
# This is expected to be used with interactive Makie backends.
103+
# """
104+
# function plot(pyramid::Pyramid;colorbar=true, size=(1155, 1043), kwargs...)
105+
# #This should be converted into a proper recipe for Makie but this would depend on a pyramid type.
106+
# fig = Figure(;size)
107+
# lon, lat = DD.dims(DD.parent(pyramid))
108+
# ax = Axis(fig[1,1], limits=(extrema(lon), extrema(lat)), aspect=DataAspect())
109+
# hmap = plot!(ax, pyramid;kwargs...)
110+
# if colorbar
111+
# Colorbar(fig[1,2], hmap, height = Relative(3.5 / 4))
112+
# end
113+
# ax.autolimitaspect = 1
114+
# FigureAxisPlot(fig, ax, hmap)
115+
# end
116+
117+
# function plot!(ax, pyramid::Pyramid;interp=false, kwargs...)#; rastercrs=crs(parent(pyramid)),plotcrs=EPSG(3857), kwargs...)
118+
# tip = levels(pyramid)[end-2][:,:]
119+
# #@show typeof(tip)
120+
# subtypes = Union{typeof.(pyramid.levels)..., typeof(parent(pyramid))}
121+
# data = Observable{DD.AbstractDimMatrix}(tip)
122+
# xyext = values.(extent(pyramid, (XDim, YDim)))
123+
# xval = first(xyext)
124+
# yval = last(xyext)
125+
# rasdataext = Extent(X=xval, Y=yval)
126+
# rasext = extent(pyramid)
127+
# on(ax.scene.viewport) do viewport
128+
# limext = extent(ax.finallimits[])
129+
130+
# datalimit = switchkeys(limext, rasext)
131+
# data[] = miss2nan.(selectlevel(pyramid, datalimit, target_imsize=viewport.widths))
132+
# end
133+
# on(ax.finallimits) do limits
134+
# limext = extent(limits)
135+
# # Compute limit in raster projection
136+
# #trans = Proj.Transformation(plotcrs, rastercrs, always_xy=true)
137+
# #datalimit = trans_bounds(trans, limext)
138+
# datalimit = switchkeys(limext, rasext)
139+
# if intersects(rasdataext, limext)
140+
# rasdata = selectlevel(pyramid, datalimit, target_imsize=ax.scene.viewport[].widths)
141+
# # Project selected data to plotcrs
142+
# #data.val = Rasters.resample(rasdata, crs=plotcrs, method=:bilinear )
143+
# data.val = miss2nan.(rasdata)
144+
# end
145+
# notify(data)
146+
# end
147+
# #@show typeof(data)
148+
# hmap = heatmap!(ax, data; interpolate=interp, kwargs...)
149+
# end
63150
end

src/PyramidScheme.jl

Lines changed: 37 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@ using YAXArrays: Cube, YAXArray, to_dataset, savedataset, setchunks, open_datase
1616
using Zarr: Zarr, zcreate, zopen, writeattrs
1717
using DimensionalData: DimensionalData as DD
1818
using DimensionalData.Dimensions: XDim, YDim
19-
using Extents: Extent, extent
19+
using Extents: Extent, extent, intersects
2020
using FillArrays: Zeros
2121
using Proj
2222
using OffsetArrays
@@ -39,7 +39,7 @@ struct Pyramid{T,N,D,A,B<:DD.AbstractDimArray{T,N,D,A},L, Me} <: DD.AbstractDimA
3939
metadata::Me
4040
end
4141

42-
function Pyramid(data::DD.AbstractDimArray; resampling_method=mean skipmissing, kwargs...)
42+
function Pyramid(data::DD.AbstractDimArray; resampling_method= mean skipmissing, kwargs...)
4343
pyrdata, pyraxs = getpyramids(resampling_method, data; kwargs...)
4444
levels = DD.DimArray.(pyrdata, pyraxs)
4545
meta = Dict(deepcopy(DD.metadata(data)))
@@ -103,7 +103,7 @@ function DD.modify(f, pyr::Pyramid)
103103
Pyramid(pbase, plevels, pyr.metadata)
104104
end
105105
Base.read(pyr::Pyramid) = DD.modify(Array, pyr)
106-
@inline function DD.rebuild(A::Pyramid, data, dims::Tuple=dims(A), refdims=refdims(A), name=name(A))
106+
@inline function DD.rebuild(A::Pyramid, data, dims::Tuple=DD.dims(A), refdims=DD.refdims(A), name=DD.name(A))
107107
Pyramid(DD.rebuild(parent(A), data, dims, refdims, name, nothing), A.levels, A.metadata)
108108
end
109109

@@ -271,9 +271,19 @@ Compute the number of levels for the aggregation based on the size of `data`.
271271
compute_nlevels(data, tilesize=256) = max(0,ceil(Int,log2(maximum(size(data))/tilesize)))
272272

273273
function agg_axis(d,n)
274-
# TODO this might be problematic for explicitly set axes
275-
DD.set(d, LinRange(first(d), last(d), cld(length(d), n)))
274+
# TODO this might be problematic for Reversed axes
275+
# TODO this is only correct for points not intervals
276+
npoints = cld(length(d), n)
277+
half_stepsize = step(d) * (n-1) / 2
278+
sgn = DD.isreverse(d) ? -1 : 1
279+
DD.set(d, LinRange(first(d) + sgn * half_stepsize, last(d) - sgn * half_stepsize, npoints))
276280
end
281+
#=
282+
. .
283+
. <--- new point
284+
. .
285+
=#
286+
277287
"""
278288
gen_output(t,s)
279289
@@ -326,7 +336,11 @@ function buildpyramids(path::AbstractString; resampling_method=mean, recursive=t
326336
# Build a loop for all variables in a dataset?
327337
org = Cube(path)
328338
# We run the method once to derive the output type
329-
t = typeof(resampling_method(zeros(eltype(org), 2,2)))
339+
#tfunc = typeof(resampling_method(zeros(eltype(org), 2,2)))
340+
#t = Missing <: eltype(org) ? Union{Missing, tfunc} : tfunc
341+
342+
t = Base.infer_return_type(resampling_method, (Matrix{nonmissingtype(eltype(org))},))
343+
330344
n_level = compute_nlevels(org)
331345
input_axes = filter(x-> x isa SpatialDim, DD.dims(org))
332346
if length(input_axes) != 2
@@ -377,16 +391,17 @@ end
377391
Compute the data of the pyramids of a given data cube `ras`.
378392
This returns the data of the pyramids and the dimension values of the aggregated axes.
379393
"""
380-
function getpyramids(reducefunc, ras;recursive=true)
394+
function getpyramids(reducefunc, ras;recursive=true, tilesize=256)
381395
input_axes = DD.dims(ras)
382-
n_level = compute_nlevels(ras)
396+
n_level = compute_nlevels(ras, tilesize)
383397
if iszero(n_level)
384-
@info "Array is smaller than the tilesize no pyramids are computed"
385-
[ras], [dims(ras)]
398+
@info "Array is smaller than the tilesize no pyramidlevels are computed"
399+
[ras], [DD.dims(ras)]
386400
end
387401
pyramid_sizes = [ceil.(Int, size(ras) ./ 2^i) for i in 1:n_level]
388402
pyramid_axes = [agg_axis.(input_axes,2^i) for i in 1:n_level]
389-
outtype = typeof(reducefunc(rand(eltype(ras),2,2)))
403+
outtype = Base.infer_return_type(reducefunc, (Matrix{eltype(ras)},))
404+
#outtype = Missing <: eltype(ras) ? Union{Missing, outtypefunc} : outtypefunc
390405
outmin = output_arrays(pyramid_sizes, outtype)
391406
fill_pyramids(ras,outmin,reducefunc,recursive; threaded=true)
392407

@@ -402,6 +417,7 @@ Internal function to select the raster data that should be plotted on screen.
402417
"""
403418
function selectlevel(pyramid, ext;target_imsize=(1024, 512))
404419
pyrext = extent(pyramid)
420+
intersects(pyrext, ext) || return zero(pyramid.levels[end])
405421
basepixels = map(keys(pyrext)) do bb
406422
pyrspan = pyrext[bb][2] - pyrext[bb][1]
407423
imsize = ext[bb][2] - ext[bb][1]
@@ -473,5 +489,15 @@ function tms_json(pyramid)
473489
push!(tms, "orderedAxes" => pyramidaxes())
474490
return tms
475491
end
492+
function Base.cat(A1::Pyramid, As::Pyramid...;dims)
493+
println("Inside pyr cat")
494+
@show typeof(levels.(As, 1))
495+
catlevels = [cat(A1.levels[i], levels.(As, i)...; dims) for i in eachindex(A1.levels)]
496+
catbase = cat(parent(A1), parent.(As)...; dims)
497+
Pyramid(catbase, catlevels, merge(DD.metadata(A1), DD.metadata.(As)...))
498+
end
499+
500+
501+
476502
include("broadcast.jl")
477503
end

0 commit comments

Comments
 (0)