From f877d688da6fbcf15cac791db97ccd9928ccfdd9 Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 29 Jul 2024 11:57:42 +0200 Subject: [PATCH 1/2] Add Juliacon talk --- docs/Manifest.toml | 82 +++++++++++++++++++++++- docs/Project.toml | 4 ++ docs/juliacon.qmd | 153 +++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 236 insertions(+), 3 deletions(-) create mode 100644 docs/juliacon.qmd diff --git a/docs/Manifest.toml b/docs/Manifest.toml index f9d4001..670ab53 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "ac78e054db06bc8c9510d748ed232ad5aec1caea" +project_hash = "5be710695837568d8d3ab6b3b0a0760cf10b84ed" [[deps.ADTypes]] git-tree-sha1 = "ae44a0c3d68a420d4ac0fa1f7e0c034ccecb6dc5" @@ -211,6 +211,12 @@ git-tree-sha1 = "19b98ee7e3db3b4eff74c5c9c72bf32144e24f10" uuid = "0b7ba130-8d10-5ba8-a3d6-c5182647fed9" version = "1.21.5+0" +[[deps.Bonito]] +deps = ["Base64", "CodecZlib", "Colors", "Dates", "Deno_jll", "HTTP", "Hyperscript", "LinearAlgebra", "Markdown", "MsgPack", "Observables", "RelocatableFolders", "SHA", "Sockets", "Tables", "ThreadPools", "URIs", "UUIDs", "WidgetsBase"] +git-tree-sha1 = "afaccf851288791c6a05102d3d677632a9690421" +uuid = "824d6782-a2ef-11e9-3a09-e5662e0c26f8" +version = "3.1.1" + [[deps.Bzip2_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "9e2a6b69137e6969bab0152632dcb3bc108c8bdd" @@ -434,6 +440,12 @@ git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae" uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab" version = "1.9.1" +[[deps.Deno_jll]] +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "cd6756e833c377e0ce9cd63fb97689a255f12323" +uuid = "04572ae6-984a-583e-9378-9577a1c2574d" +version = "1.33.4+0" + [[deps.DiffResults]] deps = ["StaticArraysCore"] git-tree-sha1 = "782dd5f4561f5d267313f23853baaaa4c52ea621" @@ -448,7 +460,9 @@ version = "1.15.1" [[deps.DimensionalData]] deps = ["Adapt", "ArrayInterface", "ConstructionBase", "DataAPI", "Dates", "Extents", "Interfaces", "IntervalSets", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "PrecompileTools", "Random", "RecipesBase", "SparseArrays", "Statistics", "TableTraits", "Tables"] -path = "../dev/DimensionalData" +git-tree-sha1 = "8c5b7b1446faa10b8fe152a380276bc70c728bb3" +repo-rev = "main" +repo-url = "https://github.com/rafaqz/DimensionalData.jl.git" uuid = "0703355e-b756-11e9-17c0-8b28908087d0" version = "0.27.3" @@ -774,6 +788,16 @@ git-tree-sha1 = "fb1156076f24f1dfee45b3feadb31d05730a49ac" uuid = "0329782f-3d07-4b52-b9f6-d3137cf03c7a" version = "1.0.2" +[[deps.GeoJSON]] +deps = ["Extents", "GeoFormatTypes", "GeoInterface", "GeoInterfaceMakie", "GeoInterfaceRecipes", "JSON3", "StructTypes", "Tables"] +git-tree-sha1 = "e2ae0c6d4f6b8c49eccc261fef29c290998e44a5" +uuid = "61d90e0f-e114-555e-ac52-39dfb47a3ef9" +version = "0.8.1" +weakdeps = ["Makie"] + + [deps.GeoJSON.extensions] + GeoJSONMakieExt = "Makie" + [[deps.GeometryBasics]] deps = ["EarCut_jll", "Extents", "GeoInterface", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"] git-tree-sha1 = "b62f2b2d76cee0d61a2ef2b3118cd2a3215d3134" @@ -867,6 +891,12 @@ git-tree-sha1 = "f218fe3736ddf977e0e772bc9a586b2383da2685" uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a" version = "0.3.23" +[[deps.Hyperscript]] +deps = ["Test"] +git-tree-sha1 = "179267cfa5e712760cd43dcae385d7ea90cc25a4" +uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91" +version = "0.0.5" + [[deps.ICU_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] git-tree-sha1 = "20b6765a3016e1fca0c9c93c80d50061b94218b7" @@ -1042,6 +1072,16 @@ git-tree-sha1 = "31e996f0a15c7b280ba9f76636b3ff9e2ae58c9a" uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" version = "0.21.4" +[[deps.JSON3]] +deps = ["Dates", "Mmap", "Parsers", "PrecompileTools", "StructTypes", "UUIDs"] +git-tree-sha1 = "eb3edce0ed4fa32f75a0a11217433c31d56bd48b" +uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" +version = "1.14.0" +weakdeps = ["ArrowTypes"] + + [deps.JSON3.extensions] + JSON3ArrowExt = ["ArrowTypes"] + [[deps.JpegTurbo]] deps = ["CEnum", "FileIO", "ImageCore", "JpegTurbo_jll", "TOML"] git-tree-sha1 = "fa6d0bcff8583bac20f1ffa708c3913ca605c611" @@ -1394,6 +1434,12 @@ version = "0.3.4" uuid = "14a3606d-f60d-562e-9121-12d972cd8159" version = "2023.1.10" +[[deps.MsgPack]] +deps = ["Serialization"] +git-tree-sha1 = "f5db02ae992c260e4826fe78c942954b48e1d9c2" +uuid = "99f44e22-a591-53d1-9472-aa23ef4bd671" +version = "1.2.1" + [[deps.MutableArithmetics]] deps = ["LinearAlgebra", "SparseArrays", "Test"] git-tree-sha1 = "898c56fbf8bf71afb0c02146ef26f3a454e88873" @@ -1412,6 +1458,12 @@ git-tree-sha1 = "0877504529a3e5c3343c6f8b4c0381e57e4387e4" uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3" version = "1.0.2" +[[deps.NaturalEarth]] +deps = ["Downloads", "GeoJSON", "Pkg", "Scratch"] +git-tree-sha1 = "3f75210ac08fe4496a55f9694b95859c40b8eaea" +uuid = "436b0209-26ab-4e65-94a9-6526d86fea76" +version = "0.1.0" + [[deps.NetCDF_jll]] deps = ["Artifacts", "Blosc_jll", "Bzip2_jll", "HDF5_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "OpenMPI_jll", "XML2_jll", "Zlib_jll", "Zstd_jll", "libzip_jll"] git-tree-sha1 = "a8af1798e4eb9ff768ce7fdefc0e957097793f15" @@ -1733,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", "YAXArrayBase", "YAXArrays", "Zarr"] +deps = ["Colors", "DimensionalData", "DiskArrayEngine", "DiskArrayTools", "DiskArrays", "Extents", "FillArrays", "GeoInterface", "Makie", "MakieCore", "Observables", "OffsetArrays", "Proj", "Statistics", "WGLMakie", "YAXArrayBase", "YAXArrays", "Zarr"] path = ".." uuid = "ec211b67-1c2c-4319-878f-eaee078ee145" version = "0.1.0" @@ -2099,6 +2151,12 @@ weakdeps = ["Adapt", "GPUArraysCore", "SparseArrays", "StaticArrays"] StructArraysSparseArraysExt = "SparseArrays" StructArraysStaticArraysExt = "StaticArrays" +[[deps.StructTypes]] +deps = ["Dates", "UUIDs"] +git-tree-sha1 = "ca4bccb03acf9faaf4137a9abc1881ed1841aa70" +uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4" +version = "1.10.0" + [[deps.SuiteSparse]] deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"] uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" @@ -2163,6 +2221,12 @@ version = "0.1.7" deps = ["InteractiveUtils", "Logging", "Random", "Serialization"] uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +[[deps.ThreadPools]] +deps = ["Printf", "RecipesBase", "Statistics"] +git-tree-sha1 = "50cb5f85d5646bc1422aa0238aa5bfca99ca9ae7" +uuid = "b189fb0b-2eb5-4ed4-bc0c-d34c51242431" +version = "2.1.1" + [[deps.Thrift_jll]] deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "boost_jll"] git-tree-sha1 = "fd7da49fae680c18aa59f421f0ba468e658a2d7a" @@ -2234,6 +2298,12 @@ git-tree-sha1 = "6331ac3440856ea1988316b46045303bef658278" uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f" version = "0.2.1" +[[deps.WGLMakie]] +deps = ["Bonito", "Colors", "FileIO", "FreeTypeAbstraction", "GeometryBasics", "Hyperscript", "LinearAlgebra", "Makie", "Observables", "PNGFiles", "PrecompileTools", "RelocatableFolders", "ShaderAbstractions", "StaticArrays"] +git-tree-sha1 = "ba21aa663fb95ebf522377b62a9b4d85213f01e6" +uuid = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" +version = "0.10.3" + [[deps.WeightedOnlineStats]] deps = ["LinearAlgebra", "OnlineStats", "OnlineStatsBase", "Statistics", "StatsBase"] git-tree-sha1 = "97af6ba86935292d5ed4a76cfee6e47d7ce02366" @@ -2246,6 +2316,12 @@ version = "0.6.3" [deps.WeightedOnlineStats.weakdeps] MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411" +[[deps.WidgetsBase]] +deps = ["Observables"] +git-tree-sha1 = "30a1d631eb06e8c868c559599f915a62d55c2601" +uuid = "eead4739-05f7-45a1-878c-cee36b57321c" +version = "0.1.4" + [[deps.WoodburyMatrices]] deps = ["LinearAlgebra", "SparseArrays"] git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511" diff --git a/docs/Project.toml b/docs/Project.toml index e6a3957..1ae17e3 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,8 +1,12 @@ [deps] ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" +DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" +NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76" PyramidScheme = "ec211b67-1c2c-4319-878f-eaee078ee145" RasterDataSources = "3cb90ccd-e1b6-4867-9617-4276c8b2ca36" Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689" +WGLMakie = "276b4fcb-3e11-5398-bf8b-a0c2d153d008" YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c" diff --git a/docs/juliacon.qmd b/docs/juliacon.qmd new file mode 100644 index 0000000..99acaff --- /dev/null +++ b/docs/juliacon.qmd @@ -0,0 +1,153 @@ +--- +engine: julia +title: "PyramidScheme.jl \n interactively plotting large raster data" +author: "Felix Cremer, Fabian Gans" +format: revealjs +execute: + echo: true +--- + +## What are pyramids? + +```{julia} +#| echo: false +using WGLMakie +using PyramidScheme +using DimensionalData.Dimensions +using DimensionalData +using Extents +@dim lat YDim "Latitude" +@dim lon XDim "Longitude" +replacenan(data) = data <= 0 ? NaN32 : Float32(data) +fullpyr = replacenan.(Pyramid(joinpath(@__DIR__,"..", "test/data/ESACCI-BIOMASS-L4-AGB-MERGED-100m-2020-fv4.0.zarr/"))) +europe= Extent(lon=(-10.0,35.0),lat=(33.0,70.0)) +europyr = Pyramid(fullpyr.base[lon = Between(europe.lon...), lat=Between(europe.lat...)], [l[lon = Between(europe.lon...), lat=Between(europe.lat...)] for l in fullpyr.levels], metadata(fullpyr)) + +lim = Extent(X = (-63.36854472609895, -57.18529373390659), Y = (-2.666626089016638, -1.9161481184310643)) +subpyr = Pyramid(fullpyr.levels[3][lon = Between(lim.X...), lat=Between(lim.Y...)], [l[lon = Between(lim.X...), lat=Between(lim.Y...)] for l in fullpyr.levels[4:7]], metadata(fullpyr)) +function plotpyramid(fullpyr) +ps = Pyramid(fullpyr.levels[end-2], fullpyr.levels[end-1:end], metadata(fullpyr)) +n=4 +levelsizes = collect(size.(PyramidScheme.levels(ps))) +numlevels = PyramidScheme.nlevels(ps) +offsets = [(0.,0.)] +for i in 1:numlevels + s = size(PyramidScheme.levels(ps, i)) ./2 + push!(offsets, offsets[i] .+ s) +end + + + + fig = Figure() + #globax = Axis(fig[1,1]) + #plot!(globax, fullpyr, colormap=:speed, colorscale=sqrt) + ax = Axis3(fig[1:2,1]) + + for (i, offset) in enumerate(offsets) + axn = Axis(fig[3-i,2]) + l = levelsizes[i] + x = first(offset):(first(offset)+first(l)-1) + y = last(offset):last(offset) + last(l)-1 + r = Rect3(Point3f( offset..., 2*i), Vec3f(l..., 1)) + mesh!(ax, r, color=:grey80) + heatmap!(ax, reverse(x), y,PyramidScheme.levels(ps,i-1), transformation=(:xy, 2*i+1.2), colormap=:speed, colorscale=sqrt) + heatmap!(axn, PyramidScheme.levels(ps, i-1), colormap=:speed, colorscale=sqrt) + hidexdecorations!(axn) + axn.aspect = DataAspect() +end +hidedecorations!(ax) +hidespines!(ax) +colsize!(fig.layout, 2, Relative(0.3)) +fig +end + +plotpyramid(europyr) +``` + + +## Compute the pyramids for your data +- In Memory data + +```{julia} +#| eval: false +using Rasters +using RasterDataSources +ras = Raster(WorldClim{Elevation},:elev, res="30s", lazy=true) +pyr = Pyramid(ras, resampling_method=mean) +Progress: 100%|████████████████████████████████████████████| Time: 0:00:47 +``` + +- Larger data on disk + +```{julia} +#| eval: false +using YAXArrays +elevpath = getraster(WorldClim{Elevation},:elev, res="30s") +elev = Cube(elevpath) +zarrpath = tempname() * ".zarr" +savecube(elev, zarrpath) +PyramidScheme.buildpyramids(zarrpath) +``` + +## Access Pyramiddata + +- Zarr data +- Geotiff +- Locally or in the cloud +```{julia} +agb2020 = replacenan.(Pyramid(joinpath(@__DIR__,"..", "test/data/ESACCI-BIOMASS-L4-AGB-MERGED-100m-2020-fv4.0.zarr/"))) +``` + +# Interactive visualisation + +```{julia} +#| eval: false +plot(agb2020, colormap=:speed, colorscale=sqrt) +``` + +# Computations on Pyramids +```{julia} +#| 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/"))) +agbdiff = agb2020 .- agb2017 +plot(agbdiff, colormap=:bukavu, colorrange=(-200, 200)) +``` + + +# Outlook + +- pyramids in multiple dimensions + +- Understanding more formats + +- Integrate DiskArrayEngine computations + +- Open it up to more use cases + +```{julia} +#| eval:false +using YAXArrays, PyramidScheme +using GLMakie +import PyramidScheme as PS +using DiskArrays +import DimensionalData as DD +using DimensionalData.Dimensions +@dim lat YDim "Latitude" +@dim lon XDim "Longitude" + +p2020 = PS.Pyramid("https://s3.bgc-jena.mpg.de:9000/pyramids/ESACCI-BIOMASS-L4-AGB-MERGED-100m-2020-fv4.0.zarr") +function DiskArrays.cache(a::YAXArray;maxsize=1000) + DD.rebuild(a,DiskArrays.cache(a.data;maxsize)) +end +function DiskArrays.cache(p::Pyramid;maxsize = 1000) + maxsize = maxsize ÷ (length(p.levels)+1) + Pyramid(DiskArrays.cache(p.base;maxsize),DiskArrays.cache.(p.levels;maxsize),p.metadata) +end +p2020 = DiskArrays.cache(p2020) + +replacenan(nanval) = data -> <=(nanval)(data) ? NaN32 : Float32(data) +p2020nan = replacenan(0).(p2020) +plot(p2020nan, colormap=:speed, colorscale=sqrt) +``` \ No newline at end of file From 19660b36512dd43ba9b9dfb8b5c253b718dbc1ad Mon Sep 17 00:00:00 2001 From: Felix Cremer Date: Mon, 29 Jul 2024 12:05:13 +0200 Subject: [PATCH 2/2] Forward kwargs from construtor to get pyramids Forward metadata from base to the pyramid struct change plotting to heatmap instead of image! --- src/PyramidScheme.jl | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/PyramidScheme.jl b/src/PyramidScheme.jl index 10528aa..5969f82 100644 --- a/src/PyramidScheme.jl +++ b/src/PyramidScheme.jl @@ -44,8 +44,8 @@ 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) - pyrdata, pyraxs = getpyramids(mean ∘ skipmissing, data, recursive=false) +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))) push!(meta, "resampling_method" => "mean_skipmissing") @@ -72,10 +72,10 @@ function _pyramid_gdal end function _pyramid_zarr(path) g = zopen(path) allkeys = collect(keys(g.groups)) - base = Cube(path)[Ti=1] # This getindex should be unnecessary and I should rather fix my data on disk + base = Cube(path) # This getindex should be unnecessary and I should rather fix my data on disk levavail = extrema(parse.(Int,allkeys[contains.(allkeys, r"\d")])) clevels = [Cube(open_dataset(g[string(l)])) for l in 1:last(levavail)] - Pyramid(base, clevels, Dict()) + Pyramid(base[Ti=1], clevels, DD.metadata(base)) end # refdims # name @@ -124,7 +124,6 @@ end function Base.map(f, A::Pyramid) newbase = map(f, parent(A)) newlevels = [map(f, levels(A,i)) for i in 1:nlevels(A)] - @show typeof(newlevels) Pyramid(newbase, newlevels, DD.metadata(A)) # This should handle metadata better. end @@ -315,7 +314,7 @@ The different scales are written according to the GeoZarr spec and a multiscales The data is aggregated with the specified `resampling_method`. Keyword arguments are forwarded to the `fill_pyramids` function. """ -function buildpyramids(path; resampling_method=mean, recursive=true, runner=LocalRunner, verbose=false) +function buildpyramids(path::AbstractString; resampling_method=mean, recursive=true, runner=LocalRunner, verbose=false) if YAB.backendfrompath(path) != YAB.ZarrDataset throw(ArgumentError("$path is not a Zarr dataset therefore we can't build the Pyramids inplace")) end @@ -441,7 +440,7 @@ ykey(keyext) = DD.dim2key(DD.dims(keyext, YDim)) Internal function ### Extended help - Return an Extent with the limits from `dataext` on the keys of `keyext`. + Return an Extent with the limits from `dataext` and the keys of `keyext`. We assume that dataext has keys X, and Y and the keys of keyext are XDim and YDim from DimensionalData """ function switchkeys(dataext, keyext) @@ -484,7 +483,7 @@ function plot!(ax, pyramid::Pyramid;interp=false, kwargs...)#; rastercrs=crs(par notify(data) end #@show typeof(data) - hmap = image!(ax, data; interpolate=interp, kwargs...) + hmap = heatmap!(ax, data; interpolate=interp, kwargs...) end """