diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts index 622c99e0..6ecb5280 100644 --- a/docs/src/.vitepress/config.mts +++ b/docs/src/.vitepress/config.mts @@ -26,6 +26,7 @@ const navTemp = { { text: 'Create', link: '/UserGuide/create' }, { text: 'Select', link: '/UserGuide/select' }, { text: 'Compute', link: '/UserGuide/compute' }, + { text: 'Plot', link: '/UserGuide/plot' }, { text: 'Chunk', link: '/UserGuide/chunk' }, { text: 'Cache', link: '/UserGuide/cache' }, { text: 'Group', link: '/UserGuide/group' }, diff --git a/docs/src/UserGuide/plot.md b/docs/src/UserGuide/plot.md new file mode 100644 index 00000000..0806305a --- /dev/null +++ b/docs/src/UserGuide/plot.md @@ -0,0 +1,72 @@ +# Plot YAXArrays + +This section describes how to visualize YAXArrays. +See also the [Plotting maps tutorial](/tutorials/plottingmaps.html) to plot geospatial data. +All [plotting capabilities](https://rafaqz.github.io/DimisensionalData.jl/dev/plots) of `AbstractDimArray` apply to a `YAXArrays` as well, because every `YAXArray` is also an `AbstractDimArray`. + +## Plot a YAXArrray + +Create a simple YAXArray: + +```@example plot +using CairoMakie +using YAXArrays +using DimensionalData + +data = collect(reshape(1:20, 4, 5)) +axlist = (X(1:4), Y(1:5)) +a = YAXArray(axlist, data) +``` + +Plot the entire array: + +```@example plot +plot(a) +``` + +This will plot a heatmap, because the array is a matrix. + +Plot the first column: + +```@example plot +plot(a[Y=1]) +``` + +This results in a scatter plot, because the subarray is a vector. + +## Plot a YAXArrray with CF conventions + +[Climate and Forecast Metadata Conventions](https://cfconventions.org/) are used to generate appropriate labels for the plot whenever possible. +This requires the YAXArray to have metadata properties like `standard_name` and `units`. + +Get a `Dataset` with CF meta data: + +```@example plot +using NetCDF +using Downloads: download + +path = download("https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc", "example.nc") +ds = open_dataset(path) +``` + +Plot the first time step of the sea surface temperature with CF metadata: + +```@example plot +plot(ds.tos[time=1]) +``` + +Time in [Climate and Forecasting (CF) conventions](https://cfconventions.org/Data/cf-conventions/cf-conventions-1.12/cf-conventions.html#time-coordinate-units) requires conversion before plotting, e.g., to plot the sea surface temperature over time at a given location (e.g. the null island): + +```@example plot +a = ds.tos[lon = Near(0), lat = Near(0)] +times = Ti(map(x -> DateTime(string(x)), a.time.val)) +a = YAXArray((times,), collect(a.data), a.properties) + +plot(a) +``` + +Or as as an explicit line plot: + +```@example plot +lines(a) +``` diff --git a/src/Cubes/Cubes.jl b/src/Cubes/Cubes.jl index fd0b6b55..ae046b32 100644 --- a/src/Cubes/Cubes.jl +++ b/src/Cubes/Cubes.jl @@ -15,7 +15,7 @@ using DiskArrayTools: CFDiskArray using DocStringExtensions using Tables: istable, schema, columns using DimensionalData: DimensionalData as DD, AbstractDimArray, NoName -import DimensionalData: name +import DimensionalData: name, label export concatenatecubes, caxes, subsetcube, readcubedata, renameaxis!, YAXArray, setchunks, cache @@ -91,7 +91,7 @@ It can wrap normal arrays or, more typically DiskArrays. $(FIELDS) """ -struct YAXArray{T,N,A<:AbstractArray{T,N}, D, Me} <: AbstractDimArray{T,N,D,A} +struct YAXArray{T,N,A<:AbstractArray{T,N},D,Me} <: AbstractDimArray{T,N,D,A} "`Tuple` of Dimensions containing the Axes of the Cube" axes::D "length(axes)-dimensional array which holds the data, this can be a lazy DiskArray" @@ -120,7 +120,7 @@ struct YAXArray{T,N,A<:AbstractArray{T,N}, D, Me} <: AbstractDimArray{T,N,D,A} throw(ArgumentError("Can not construct YAXArray, supplied chunk dimension is $(ndims(chunks)) while the number of dims is $(length(axes))")) else axes = DD.format(axes, data) - return new{eltype(data),ndims(data),typeof(data),typeof(axes), typeof(properties)}( + return new{eltype(data),ndims(data),typeof(data),typeof(axes),typeof(properties)}( axes, data, properties, @@ -131,16 +131,34 @@ struct YAXArray{T,N,A<:AbstractArray{T,N}, D, Me} <: AbstractDimArray{T,N,D,A} end end -name(::YAXArray) = NoName() -YAXArray(axes, data, properties = Dict{String,Any}(); cleaner = CleanMe[], chunks = eachchunk(data)) = +"Name of an YAXArray using CF conventions. Searches first for metadata `long_name`, followed by `standard_name` and finally `name`" +function name(a::YAXArray) + # as implemented in python xarray + isempty(a.properties) && return NoName() + haskey(a.properties, "long_name") && return a.properties["long_name"] |> Symbol + haskey(a.properties, "standard_name") && return a.properties["standard_name"] |> Symbol + haskey(a.properties, "name") && return a.properties["name"] |> Symbol + return NoName() +end + +"Label of an YAXArrray using CF conventions. Includes the name and the unit, if present. Searches first for metadata `units` followed by `unit``" +function label(a::YAXArray) + # as implemented in python xarray + isempty(a.properties) && return "" + haskey(a.properties, "units") && return "$(name(a)) [$(a.properties["units"])]" + haskey(a.properties, "unit") && return "$(name(a)) [$(a.properties["unit"])]" + return string(name(a)) +end + +YAXArray(axes, data, properties=Dict{String,Any}(); cleaner=CleanMe[], chunks=eachchunk(data)) = YAXArray(axes, data, properties, chunks, cleaner) -YAXArray(axes,data,properties,cleaner) = YAXArray(axes,data,properties,eachchunk(data),cleaner) +YAXArray(axes, data, properties, cleaner) = YAXArray(axes, data, properties, eachchunk(data), cleaner) function YAXArray(x::AbstractArray) ax = caxes(x) props = getattributes(x) chunks = eachchunk(x) - YAXArray(ax, x, props,chunks=chunks) + YAXArray(ax, x, props, chunks=chunks) end @@ -172,20 +190,20 @@ end Base.Generator(f, A::YAXArray) = Base.Generator(f, parent(A)) Base.ndims(a::YAXArray{<:Any,N}) where {N} = N Base.eltype(a::YAXArray{T}) where {T} = T -function Base.permutedims(c::YAXArray, p) +function Base.permutedims(c::YAXArray, p) newdims = DD.sortdims(DD.dims(c), Tuple(p)) dimnums = map(d -> DD.dimnum(c, d), p) newdata = permutedims(getdata(c), dimnums) newchunks = DiskArrays.GridChunks(eachchunk(c).chunks[collect(dimnums)]) YAXArray(newdims, newdata, c.properties, newchunks, c.cleaner) end -DiskArrays.cache(a::YAXArray;maxsize=1000) = DD.rebuild(a,cache(a.data;maxsize)) +DiskArrays.cache(a::YAXArray; maxsize=1000) = DD.rebuild(a, cache(a.data; maxsize)) # DimensionalData overloads -DD.dims(x::YAXArray) = getfield(x,:axes) +DD.dims(x::YAXArray) = getfield(x, :axes) DD.refdims(::YAXArray) = () -DD.metadata(x::YAXArray) = getfield(x,:properties) +DD.metadata(x::YAXArray) = getfield(x, :properties) function DD.rebuild(A::YAXArray, data::AbstractArray, dims::Tuple, refdims::Tuple, name, metadata) #chunks = map(dims, eachchunk(data).chunks) do d, chunk @@ -205,14 +223,14 @@ function DD.rebuild(A::YAXArray; data=parent(A), dims=DD.dims(A), metadata=DD.me end function caxes(x) - #@show x - #@show typeof(x) - dims = map(enumerate(dimnames(x))) do a - index, symbol = a - values = YAXArrayBase.dimvals(x, index) - DD.Dim{symbol}(values) - end - (dims... ,) + #@show x + #@show typeof(x) + dims = map(enumerate(dimnames(x))) do a + index, symbol = a + values = YAXArrayBase.dimvals(x, index) + DD.Dim{symbol}(values) + end + (dims...,) end caxes(x::DD.AbstractDimArray) = collect(DD.dims(x)) @@ -238,18 +256,18 @@ function readcubedata(x) YAXArray(caxes(x), getindex_all(x), getattributes(x)) end -interpret_cubechunks(cs::NTuple{N,Int},cube) where N = DiskArrays.GridChunks(getdata(cube),cs) -interpret_cubechunks(cs::DiskArrays.GridChunks,_) = cs -interpret_dimchunk(cs::Integer,s) = DiskArrays.RegularChunks(cs,0,s) +interpret_cubechunks(cs::NTuple{N,Int}, cube) where {N} = DiskArrays.GridChunks(getdata(cube), cs) +interpret_cubechunks(cs::DiskArrays.GridChunks, _) = cs +interpret_dimchunk(cs::Integer, s) = DiskArrays.RegularChunks(cs, 0, s) interpret_dimchunk(cs::DiskArrays.ChunkVector, _) = cs -function interpret_cubechunks(cs,cube) +function interpret_cubechunks(cs, cube) oldchunks = DiskArrays.eachchunk(cube).chunks for k in keys(cs) - i = findAxis(k,cube) + i = findAxis(k, cube) if i !== nothing - dimchunk = interpret_dimchunk(cs[k],size(cube.data,i)) - oldchunks = Base.setindex(oldchunks,dimchunk,i) + dimchunk = interpret_dimchunk(cs[k], size(cube.data, i)) + oldchunks = Base.setindex(oldchunks, dimchunk, i) end end GridChunks(oldchunks) @@ -267,7 +285,7 @@ of this chunking, use `savecube` on the resulting array. The `chunks` argument c - an AbstractDict or NamedTuple mapping one or more axis names to chunk sizes """ -setchunks(c::YAXArray,chunks) = YAXArray(c.axes,c.data,c.properties,interpret_cubechunks(chunks,c),c.cleaner) +setchunks(c::YAXArray, chunks) = YAXArray(c.axes, c.data, c.properties, interpret_cubechunks(chunks, c), c.cleaner) cubechunks(c) = approx_chunksize(eachchunk(c)) DiskArrays.eachchunk(c::YAXArray) = c.chunks getindex_all(a) = getindex(a, ntuple(_ -> Colon(), ndims(a))...).data @@ -283,45 +301,45 @@ end =# -function batchextract(x,i) +function batchextract(x, i) # This function should be documented and moved to DimensionalData sch = schema(i) axinds = map(sch.names) do n - findAxis(n,x) + findAxis(n, x) end tcols = columns(i) #Try to find a column denoting new axis name and values newaxcol = nothing - - if any(isnothing,axinds) - allnothings = findall(isnothing,axinds) + + if any(isnothing, axinds) + allnothings = findall(isnothing, axinds) if length(allnothings) == 1 newaxcol = allnothings[1] end - tcols = (;[p[1:2] for p in zip(keys(tcols), values(tcols), axinds) if !isnothing(last(p))]...) - axinds = filter(!isnothing,axinds) + tcols = (; [p[1:2] for p in zip(keys(tcols), values(tcols), axinds) if !isnothing(last(p))]...) + axinds = filter(!isnothing, axinds) end - + allax = 1:ndims(x) - axrem = setdiff(allax,axinds) + axrem = setdiff(allax, axinds) ai1, ai2 = extrema(axinds) - - if !all(diff(sort(collect(axinds))).==1) + + if !all(diff(sort(collect(axinds))) .== 1) #Axes to be extracted from are not consecutive in cube -> permute - p = [1:(ai1-1);collect(axinds);filter(!in(axinds),ai1:ai2);(ai2+1:ndims(x))] - x_perm = permutedims(x,p) - return batchextract(x_perm,i) + p = [1:(ai1-1); collect(axinds); filter(!in(axinds), ai1:ai2); (ai2+1:ndims(x))] + x_perm = permutedims(x, p) + return batchextract(x_perm, i) end - cartinds = map(axinds,tcols) do iax,col + cartinds = map(axinds, tcols) do iax, col axcur = caxes(x)[iax] map(col) do val - axVal2Index(axcur,val) + axVal2Index(axcur, val) end end - - before = ntuple(_->Colon(),ai1-1) - after = ntuple(_->Colon(),ndims(x)-ai2) + + before = ntuple(_ -> Colon(), ai1 - 1) + after = ntuple(_ -> Colon(), ndims(x) - ai2) sp = issorted(axinds) ? nothing : sortperm(collect(axinds)) function makeindex(sp, inds...) if sp === nothing @@ -330,37 +348,37 @@ function batchextract(x,i) CartesianIndex(inds[sp]...) end end - indlist = makeindex.(Ref(sp),cartinds...) - d = getdata(x)[before...,indlist,after...] + indlist = makeindex.(Ref(sp), cartinds...) + d = getdata(x)[before..., indlist, after...] cax = caxes(x) newax = if newaxcol == nothing - outaxis_from_data(cax,axinds,indlist) + outaxis_from_data(cax, axinds, indlist) else - outaxis_from_column(i,newaxcol) + outaxis_from_column(i, newaxcol) end outax = Tuple([axcopy(a) for a in cax][axrem]...) - insert!(outax,minimum(axinds),newax) - YAXArray(outax,d,x.properties) + insert!(outax, minimum(axinds), newax) + YAXArray(outax, d, x.properties) end -function outaxis_from_column(tab,icol) +function outaxis_from_column(tab, icol) axdata = columns(tab)[icol] axname = schema(tab).names[icol] if eltype(axdata) <: AbstractString || - (!issorted(axdata) && !issorted(axdata, rev = true)) + (!issorted(axdata) && !issorted(axdata, rev=true)) DD.rebuild(DD.name2dim(Symbol(axname)), axdata) else DD.rebuild(DD.name2dim(Symbol(axname)), axdata) end end -function outaxis_from_data(cax,axinds,indlist) - mergeaxes = getindex.(Ref(cax),axinds) +function outaxis_from_data(cax, axinds, indlist) + mergeaxes = getindex.(Ref(cax), axinds) mergenames = axname.(mergeaxes) - newname = join(mergenames,'_') + newname = join(mergenames, '_') minai = minimum(axinds) mergevals = map(indlist) do i - broadcast(mergeaxes,axinds) do ax,ai + broadcast(mergeaxes, axinds) do ax, ai ax.values[i[ai-minai+1]] end end @@ -471,14 +489,14 @@ function _subsetcube(z, subs; kwargs...) end -function Base.getindex(a::YAXArray, args::DD.Dimension...; kwargs...) +function Base.getindex(a::YAXArray, args::DD.Dimension...; kwargs...) kwargsdict = Dict{Any,Any}(kwargs...) for ext in YAXDefaults.subsetextensions ext(kwargsdict) end d2 = Dict() - for (k,v) in kwargsdict - d = getAxis(k,a) + for (k, v) in kwargsdict + d = getAxis(k, a) if d !== nothing d2[DD.name(d)] = v else @@ -491,7 +509,7 @@ end Base.read(d::YAXArray) = getindex_all(d) function formatbytes(x) - exts = ["bytes", "KB", "MB", "GB", "TB","PB"] + exts = ["bytes", "KB", "MB", "GB", "TB", "PB"] i = 1 while x >= 1024 i = i + 1 @@ -515,7 +533,7 @@ function DD.show_after(io::IO, mime, c::YAXArray) # ? sizeof : Check if the element type is a bitstype or a union of bitstypes if (isconcretetype(eltype(c)) && isbitstype(eltype(c))) || - (eltype(c) isa Union && all(isbitstype, Base.uniontypes(eltype(c)))) + (eltype(c) isa Union && all(isbitstype, Base.uniontypes(eltype(c)))) println(io, "\n data size: ", formatbytes(cubesize(c))) else # fallback diff --git a/test/Cubes/cubes.jl b/test/Cubes/cubes.jl index 7b436cdc..24779aac 100644 --- a/test/Cubes/cubes.jl +++ b/test/Cubes/cubes.jl @@ -18,7 +18,7 @@ using DimensionalData @testset "Basic array Functions" begin @test size(a) == (4, 5) - @test size(a, 1) == 4 + @test size(a, 1) == 4 @test size(a, :X) == 4 @test size(a, 2) == 5 @test size(a, :YVals) == 5 @@ -51,7 +51,7 @@ using DimensionalData cs = CartesianIndices.(collect(Iterators.product([1:2, 3:5, 6:6], [1:2, 3:4]))) - a2 = renameaxis!(a2, :X => :Ax1) + a2 = renameaxis!(a2, :X => :Ax1) @test YAXArrayBase.dimname(a2, 1) == :Ax1 @@ -104,28 +104,41 @@ using DimensionalData @test endswith(Cubes.formatbytes(1205), "KB") @test endswith(Cubes.formatbytes(1200000), "MB") end -#= - @testset "Subsets" begin - s = YAXArrays.Cubes.subsetcube(a, X = 1.5..3.5) - @test s.data == [2 6 10 14 18; 3 7 11 15 19] - @test s.axes[1] == X(2.0:3.0) - @test s.axes[2] == Y([1, 2, 3, 4, 5]) - ax = a.axes[1] - @test YAXArrays.Cubes.interpretsubset(CartesianIndices((1:2,)), ax) == 1:2 - @test YAXArrays.Cubes.interpretsubset(CartesianIndex((2,)), ax) == 2 - @test YAXArrays.Cubes.interpretsubset(2.1, ax) == 2 - @test YAXArrays.Cubes.interpretsubset((3.5, 1.5), ax) == 2:3 - @test YAXArrays.Cubes.interpretsubset(0.8..2.2, ax) == 1:2 - tax = RangeAxis("ADate", Date(2001):Day(1):Date(2003, 2, 28)) - @test YAXArrays.Cubes.interpretsubset((Date(2001, 1, 2), Date(2001, 1, 5)), tax) == - 2:4 - @test YAXArrays.Cubes.interpretsubset(2001:2002, tax) == 1:730 - @test YAXArrays.Cubes.interpretsubset([1, 3, 5], a.axes[2]) == [1, 3, 5] - - s2 = a[X = 0.5..3.5, Y = [1, 5, 4]] - @test s2.data == [1 17 13; 2 18 14; 3 19 15] - @test s2.axes[1] == X(1.0:3.0) - @test s2.axes[2] == Y([1, 5, 4]) + + @testset "YAXArrays with CF conventions" begin + attrs = Dict( + "standard_name" => "air_temperature", + "long_name" => "global mean air temperature", + "units" => "K", + ) + a_cf = YAXArray(a.axes, a.data, attrs) + @test DimensionalData.name(a) == DimensionalData.NoName() + @test DimensionalData.label(a) == "" + @test DimensionalData.name(a_cf) == Symbol("global mean air temperature") + @test DimensionalData.label(a_cf) == "global mean air temperature [K]" end -=# + #= + @testset "Subsets" begin + s = YAXArrays.Cubes.subsetcube(a, X = 1.5..3.5) + @test s.data == [2 6 10 14 18; 3 7 11 15 19] + @test s.axes[1] == X(2.0:3.0) + @test s.axes[2] == Y([1, 2, 3, 4, 5]) + ax = a.axes[1] + @test YAXArrays.Cubes.interpretsubset(CartesianIndices((1:2,)), ax) == 1:2 + @test YAXArrays.Cubes.interpretsubset(CartesianIndex((2,)), ax) == 2 + @test YAXArrays.Cubes.interpretsubset(2.1, ax) == 2 + @test YAXArrays.Cubes.interpretsubset((3.5, 1.5), ax) == 2:3 + @test YAXArrays.Cubes.interpretsubset(0.8..2.2, ax) == 1:2 + tax = RangeAxis("ADate", Date(2001):Day(1):Date(2003, 2, 28)) + @test YAXArrays.Cubes.interpretsubset((Date(2001, 1, 2), Date(2001, 1, 5)), tax) == + 2:4 + @test YAXArrays.Cubes.interpretsubset(2001:2002, tax) == 1:730 + @test YAXArrays.Cubes.interpretsubset([1, 3, 5], a.axes[2]) == [1, 3, 5] + + s2 = a[X = 0.5..3.5, Y = [1, 5, 4]] + @test s2.data == [1 17 13; 2 18 14; 3 19 15] + @test s2.axes[1] == X(1.0:3.0) + @test s2.axes[2] == Y([1, 5, 4]) + end + =# end