diff --git a/Project.toml b/Project.toml index cfc1eef..6e01d03 100644 --- a/Project.toml +++ b/Project.toml @@ -39,6 +39,7 @@ Infiltrator = "1.8.8" OrderedCollections = "1.8.1" StaticArrays = "1.9.13" YAXArrays = "0.7.3" +Zarr = "0.10.1" julia = "1.6.7" [extras] diff --git a/docs/src/public b/docs/src/public deleted file mode 120000 index f40fe05..0000000 --- a/docs/src/public +++ /dev/null @@ -1 +0,0 @@ -assets \ No newline at end of file diff --git a/docs/src/public/cube-shape.drawio.svg b/docs/src/public/cube-shape.drawio.svg new file mode 100644 index 0000000..d33b185 --- /dev/null +++ b/docs/src/public/cube-shape.drawio.svg @@ -0,0 +1,152 @@ + + + + + + + + + + + + +
+
+
+ + Longitude + +
+
+
+ + + Long... + + + + + + + +
+
+
+ + Time + +
+
+
+ + + Time + + + + + + + +
+
+
+ + Latitude +
+ +
+
+
+ + + Lati... + + + + + + + +
+
+
+ + Geo Cube + +
+
+
+ + + Geo Cube + + + + + + + +
+
+
+ + Cell Cube + +
+
+
+ + + Cell Cube + + + + + + + +
+
+
+ + Cell ID +
+ +
+
+
+ + + Cell... + + + + + + +
+
+
+ + Time + +
+
+
+ + + Time + + + + + + + + + Text is not SVG - cannot display + + + + \ No newline at end of file diff --git a/docs/src/public/dest.png b/docs/src/public/dest.png new file mode 100644 index 0000000..f50f2ec Binary files /dev/null and b/docs/src/public/dest.png differ diff --git a/docs/src/public/dggrid-grids-multi-levels.png b/docs/src/public/dggrid-grids-multi-levels.png new file mode 100644 index 0000000..3e53d2f Binary files /dev/null and b/docs/src/public/dggrid-grids-multi-levels.png differ diff --git a/docs/src/public/dggs-distortion.png b/docs/src/public/dggs-distortion.png new file mode 100644 index 0000000..9f8970a Binary files /dev/null and b/docs/src/public/dggs-distortion.png differ diff --git a/docs/src/public/grid-levels.png b/docs/src/public/grid-levels.png new file mode 100644 index 0000000..0294fd4 Binary files /dev/null and b/docs/src/public/grid-levels.png differ diff --git a/docs/src/public/hex-pattern.png b/docs/src/public/hex-pattern.png new file mode 100644 index 0000000..e09a7df Binary files /dev/null and b/docs/src/public/hex-pattern.png differ diff --git a/docs/src/public/hexagon-children-aperture.png b/docs/src/public/hexagon-children-aperture.png new file mode 100644 index 0000000..0f762bd Binary files /dev/null and b/docs/src/public/hexagon-children-aperture.png differ diff --git a/docs/src/public/horizontal-neighbors.drawio.svg b/docs/src/public/horizontal-neighbors.drawio.svg new file mode 100644 index 0000000..b12c497 --- /dev/null +++ b/docs/src/public/horizontal-neighbors.drawio.svg @@ -0,0 +1,114 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/docs/src/public/icon.drawio.png b/docs/src/public/icon.drawio.png new file mode 100644 index 0000000..8ca2b1a Binary files /dev/null and b/docs/src/public/icon.drawio.png differ diff --git a/docs/src/public/icon.drawio.svg b/docs/src/public/icon.drawio.svg new file mode 100644 index 0000000..e134b01 --- /dev/null +++ b/docs/src/public/icon.drawio.svg @@ -0,0 +1,10 @@ + + + + + + + + + + \ No newline at end of file diff --git a/docs/src/public/index.drawio.svg b/docs/src/public/index.drawio.svg new file mode 100644 index 0000000..279ba73 --- /dev/null +++ b/docs/src/public/index.drawio.svg @@ -0,0 +1,14 @@ + + + + + + + + + + + + + + \ No newline at end of file diff --git a/docs/src/public/logo-eu.png b/docs/src/public/logo-eu.png new file mode 100644 index 0000000..e10a6e6 Binary files /dev/null and b/docs/src/public/logo-eu.png differ diff --git a/docs/src/public/logo-mpi-bgc.svg b/docs/src/public/logo-mpi-bgc.svg new file mode 100644 index 0000000..8e35546 --- /dev/null +++ b/docs/src/public/logo-mpi-bgc.svg @@ -0,0 +1,224 @@ + + + + + + image/svg+xml + + + + + + + + + + + + minerva + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/public/logo-open-earth-monitor.png b/docs/src/public/logo-open-earth-monitor.png new file mode 100644 index 0000000..d615644 Binary files /dev/null and b/docs/src/public/logo-open-earth-monitor.png differ diff --git a/docs/src/public/logo.drawio.svg b/docs/src/public/logo.drawio.svg new file mode 100644 index 0000000..31bd9b3 --- /dev/null +++ b/docs/src/public/logo.drawio.svg @@ -0,0 +1,54 @@ + + + + + + + +
+
+
+ + DGGS.jl + +
+
+
+ + + DGGS.jl + + + + + + +
+
+
+ + Discrete Global Grid System + +
+
+
+ + + Discrete Global Grid System + + + + + + + + + + + + + Text is not SVG - cannot display + + + + \ No newline at end of file diff --git a/docs/src/public/modis-ndvi-level6.png b/docs/src/public/modis-ndvi-level6.png new file mode 100644 index 0000000..b790896 Binary files /dev/null and b/docs/src/public/modis-ndvi-level6.png differ diff --git a/docs/src/public/pentacube-overview.png b/docs/src/public/pentacube-overview.png new file mode 100644 index 0000000..2cb983e Binary files /dev/null and b/docs/src/public/pentacube-overview.png differ diff --git a/docs/src/public/pentagon-children-q2di.png b/docs/src/public/pentagon-children-q2di.png new file mode 100644 index 0000000..5f63e30 Binary files /dev/null and b/docs/src/public/pentagon-children-q2di.png differ diff --git a/docs/src/public/plot-get-started.png b/docs/src/public/plot-get-started.png new file mode 100644 index 0000000..acc8412 Binary files /dev/null and b/docs/src/public/plot-get-started.png differ diff --git a/docs/src/public/plot-globe.png b/docs/src/public/plot-globe.png new file mode 100644 index 0000000..6cdb8ca Binary files /dev/null and b/docs/src/public/plot-globe.png differ diff --git a/docs/src/public/plot-map.png b/docs/src/public/plot-map.png new file mode 100644 index 0000000..e31bfc7 Binary files /dev/null and b/docs/src/public/plot-map.png differ diff --git a/docs/src/public/plot-native.png b/docs/src/public/plot-native.png new file mode 100644 index 0000000..b073e26 Binary files /dev/null and b/docs/src/public/plot-native.png differ diff --git a/docs/src/public/rings-disks.png b/docs/src/public/rings-disks.png new file mode 100644 index 0000000..f710b19 Binary files /dev/null and b/docs/src/public/rings-disks.png differ diff --git a/ext/DGGSMakie.jl b/ext/DGGSMakie.jl index 4ea0192..de9abc5 100644 --- a/ext/DGGSMakie.jl +++ b/ext/DGGSMakie.jl @@ -150,6 +150,11 @@ function Makie.plot( # use colormap if only one layer is supplied if dggs isa DGGSArray || (dggs isa DGGSPyramid && args isa Tuple{Symbol}) filtered_data = filter(x -> !ismissing(x) && !isnan(x), data[]) + + if filtered_data isa Vector{Missing} + return fig + end + cb_limits = (minimum(filtered_data), maximum(filtered_data)) label = if dggs isa DGGSPyramid && length(args) == 1 diff --git a/ext/DGGSZarr.jl b/ext/DGGSZarr.jl index c97319a..024c310 100644 --- a/ext/DGGSZarr.jl +++ b/ext/DGGSZarr.jl @@ -12,17 +12,119 @@ using FillArrays # Currently, YAXArrays does not support saving the experimental nested DimTree type -function DGGS.save_dggs_pyramid(path::String, dggs_p::DGGSPyramid, args...; storetype=DirectoryStore) +""" + _local_ranges(r, i, chunk_size) + +Compute local indices within a chunk for writing data. +Handles boundary chunks that may be smaller than chunk_size. +""" +function _local_ranges(r, i, chunk_size) + ntuple(length(r)) do d + start_local = r[d].start - (i[d] - 1) * chunk_size[d] + end_local = r[d].stop - (i[d] - 1) * chunk_size[d] + start_local:end_local + end +end + +""" + _write_tile_array_to_zarr!(disk_array, tile_array, chunks) + +Write only present tiles from a TileArray to a Zarr array. +Missing tiles (chunks that are `missing`) are not written to disk. +""" +function _write_tile_array_to_zarr!(disk_array, tile_array, chunks) + if !isnothing(chunks) + disk_array = setchunks(disk_array, chunks) + end + + for (r, i) in zip(DGGS.ranges(tile_array), findall(!ismissing, tile_array.data)) + chunk_data = tile_array.data[i] + local_ranges = _local_ranges(r, i, tile_array.chunk_size) + disk_array[r...] .= chunk_data[local_ranges...] + end +end + +""" + save_dggs_array(file_path, dggs_array; chunks=nothing, kwargs...) + +Save a DGGSArray to Zarr format, optimized for TileArrays. +Missing tiles are not written to disk, leveraging Zarr's sparse chunk storage. +Uses `fill_value=nothing` so chunks can be filled with any value. +""" +function DGGS.save_dggs_array(file_path::String, dggs_array::DGGSArray; chunks=nothing, kwargs...) + yax_array = YAXArray(dggs_array) + if !isnothing(chunks) + yax_array = setchunks(yax_array, chunks) + end + + # Check if underlying data is a TileArray + if dggs_array.data isa TileArray + # Create skeleton with fill_value=nothing + ds = Dataset(; Dict(DD.name(dggs_array) => yax_array)...) + disk_ds = savedataset(ds; path=file_path, skeleton=true, driver=:zarr, fill_value=nothing, kwargs...) + + # Write only present tiles + _write_tile_array_to_zarr!(disk_ds[DD.name(dggs_array)], dggs_array.data, chunks) + else + # Fallback for non-TileArray data + ds = Dataset(; Dict(DD.name(dggs_array) => yax_array)...) + savedataset(ds; path=file_path, driver=:zarr, fill_value=nothing, kwargs...) + end + return file_path +end + +""" + save_dggs_dataset(file_path, dggs_ds; chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1), kwargs...) + +Save a DGGSDataset to Zarr format, optimized for TileArrays. +Missing tiles are not written to disk, leveraging Zarr's sparse chunk storage. +Uses `fill_value=nothing` so chunks can be filled with any value. +""" +function DGGS.save_dggs_dataset(file_path::String, dggs_ds::DGGSDataset; chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1), kwargs...) + yax_ds = Dataset(dggs_ds) + if !isnothing(chunks) + yax_ds = setchunks(yax_ds, chunks) + end + + # Check if any underlying data is a TileArray + has_tile_array = any(k -> getproperty(dggs_ds, k).data isa TileArray, keys(dggs_ds)) + + + if has_tile_array + # Create skeleton with fill_value=nothing + disk_ds = savedataset(yax_ds; path=file_path, skeleton=true, driver=:zarr, fill_value=nothing, kwargs...) + + # Write only present tiles for each variable + for key in keys(dggs_ds) + if getproperty(dggs_ds, key).data isa TileArray + _write_tile_array_to_zarr!(disk_ds[key], getproperty(dggs_ds, key).data, chunks) + end + end + else + # Fallback for non-TileArray data + savedataset(yax_ds; path=file_path, driver=:zarr, fill_value=nothing, kwargs...) + end + return file_path +end + +""" + save_dggs_pyramid(path, dggs_p; storetype=DirectoryStore, chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1), kwargs...) + +Save a DGGSPyramid to Zarr format, optimized for TileArrays. +Missing tiles are not written to disk, leveraging Zarr's sparse chunk storage. +Uses `fill_value=nothing` so chunks can be filled with any value. +""" +function DGGS.save_dggs_pyramid(path::String, dggs_p::DGGSPyramid, args...; storetype=DirectoryStore, chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1), kwargs...) pyramid_attrs = Dict( "dggs_bbox" => dggs_p.bbox, "dggs_dggsrs" => dggs_p.dggsrs, ) store = storetype(path, args...) group = zgroup(store; attrs=pyramid_attrs) + for key in keys(dggs_p.branches) dggs_ds = getproperty(dggs_p, key) |> x -> x isa DGGSArray ? DGGSDataset(x) : x - ds = Dataset(dggs_ds) - savedataset(ds; path="$(path)/$(key)", driver=:zarr) + DGGS.save_dggs_dataset("$(path)/$(key)", dggs_ds; chunks=chunks, kwargs...) end return path end @@ -84,7 +186,7 @@ function DGGS.init_global_dggs_dataset( for (key, geo_array) in pairs(geo_ds.cubes) is_spatial = x_dim_name in name(geo_array.axes) && y_dim_name in name(geo_array.axes) if is_spatial - spatial_dims = (Dim{:dggs_i}(0:2*2^resolution-1), Dim{:dggs_j}(0:2^resolution-1), Dim{:dggs_n}(0:4)) + spatial_dims = (Dim{:dggs_i}(0:(2*2^resolution-1)), Dim{:dggs_j}(0:(2^resolution-1)), Dim{:dggs_n}(0:4)) non_spatial_dims = filter(x -> !(name(x) in [x_dim_name, y_dim_name]), geo_array.axes) dims = (spatial_dims..., non_spatial_dims...) else diff --git a/src/arrays.jl b/src/arrays.jl index 2683e35..915f744 100644 --- a/src/arrays.jl +++ b/src/arrays.jl @@ -37,6 +37,43 @@ function get_dggs_bbox(cells) ) end +# Version that works with any iterable of cells (e.g., dictionary keys) +function get_dggs_bbox_cells(cells) + cell = first(cells) + resolution = cell.resolution + + i_min = cell.i + i_max = cell.i + j_min, j_max = cell.j, cell.j + n_min, n_max = cell.n, cell.n + + for cell in cells + if cell.i < i_min + i_min = cell.i + elseif cell.i > i_max + i_max = cell.i + end + + if cell.j < j_min + j_min = cell.j + elseif cell.j > j_max + j_max = cell.j + end + + if cell.n < n_min + n_min = cell.n + elseif cell.n > n_max + n_max = cell.n + end + end + + return ( + Dim{:dggs_i}(i_min:i_max), + Dim{:dggs_j}(j_min:j_max), + Dim{:dggs_n}(n_min:n_max) + ) +end + "Infere max possible geo extent" function get_geo_bbox(x::Union{DGGSArray,DGGSDataset}) i_min, i_max = dims(x, :dggs_i).val.data |> x -> (first(x), last(x)) @@ -76,130 +113,136 @@ function get_geo_bbox(geo_array::AbstractDimArray, crs::String; x_name=:X, y_nam end end -function to_dggs_array( - geo_array::AbstractDimArray, - cells, - cell_coords, - dggs_bbox, - geo_bbox::Extent, - agg_func::Function - ; - outtype=Float64, - backend=:array, - path=tempname() * ".dggs.zarr", - name=get_name(geo_array), - x_name=:X, - y_name=:Y, - kwargs... -) - resolution = first(cells).resolution - - # re-grid - res = mapCube( - # mapCube can't find axes of other AbstractDimArrays e.g. Raster - YAXArray(dims(geo_array), geo_array.data, metadata(geo_array)); - indims=InDims(dims(geo_array, x_name), dims(geo_array, y_name)), - outdims=OutDims( - dggs_bbox..., - outtype=outtype, - backend=backend, - path=path - ), kwargs...) do xout, xin - for ci in CartesianIndices(xout) - i, j, n = ci.I - try - cell = Cell(dggs_bbox[1][i], dggs_bbox[2][j], dggs_bbox[3][n], resolution) - cells = cell_coords[cell] - res = agg_func(view(xin, cells)) - xout[i, j, n] = res - catch - # fill gap by averaging available neighbors - xmin = clamp(i, 2, size(xout, 1) - 1) - 1 - ymin = clamp(j, 2, size(xout, 2) - 1) - 1 - res = filter(!ismissing, xout[xmin:xmin+2, ymin:ymin+2, n]) |> agg_func - xout[i, j, n] = res - end - end +function cells_to_coord_dict(cells::DimArray{Cell{Int64},2}) + cell_coords = Dict{Cell{Int64},Vector{CartesianIndex{2}}}() + for cell_idx in CartesianIndices(cells) + cell = cells[cell_idx] + current_cells = get!(() -> CartesianIndex{2}[], cell_coords, cell) + push!(current_cells, cell_idx) end + return cell_coords +end - return DGGSArray( - res.data, dims(res), refdims(res), name, metadata(geo_array), - resolution, "ISEA4D.Penta", geo_bbox - ) +# Fused version: directly builds the cell coordinate dictionary from dimensions +# without creating an intermediate cell array +function cells_to_coord_dict(x_dim, y_dim, resolution, crs) + trans = Proj.Transformation(crs, crs_isea; ctx=Proj.proj_context_create(), always_xy=true) + cell_coords = Dict{Cell{Int64},Vector{CartesianIndex{2}}}() + + # Pre-size dictionary based on expected number of unique cells + # At high resolution, many pixels will map to the same cell + expected_cells = min(length(x_dim) * length(y_dim), 2 * 2^resolution * 2^resolution * 5) + sizehint!(cell_coords, expected_cells) + + for (j_idx, y) in enumerate(y_dim) + for (i_idx, x) in enumerate(x_dim) + cell = to_cell(x, y, resolution, trans) + current_cells = get!(() -> CartesianIndex{2}[], cell_coords, cell) + push!(current_cells, CartesianIndex(i_idx, j_idx)) + end + end + return cell_coords end -"Fast iterative version only supporting mean" + function to_dggs_array( geo_array::AbstractDimArray, - cells, - dggs_bbox, + resolution::Integer, + cell_coords, geo_bbox::Extent ; - outtype=eltype(geo_array), - outtype_counts=UInt16, - outtype_sums=eltype(geo_array), - backend=:array, - path=tempname() * ".dggs.zarr", + agg_func::Function=mean, name=get_name(geo_array), - x_name=:X, - y_name=:Y, + out_eltype=Union{Missing,eltype(geo_array)}, + chunk_length=2^12, kwargs... ) - resolution = first(cells).resolution - - # re-grid - # mean = sum first, then divide by count - # no slow dict building and lookup needed + dggsrs = "ISEA4D.Penta" - counts = zeros(outtype_counts, length.(dggs_bbox)...) + # Identify non-spatial dimensions (e.g. time, band) + spatial_dim_names = (:dggs_i, :dggs_j, :dggs_n, :X, :Y) + non_spatial = filter(d -> !(DD.name(d) in spatial_dim_names), dims(geo_array)) + non_spatial_sizes = map(length, non_spatial) - if any(size(geo_array) .> typemax(outtype_counts)) - error("Input array too large for outtype_counts, consider using a larger integer type for counts") + # Create spatial dims + spatial_dims = (Dim{:dggs_i}(0:(2*2^resolution-1)), Dim{:dggs_j}(0:(2^resolution-1)), Dim{:dggs_n}(0:4)) + all_dims = if isempty(non_spatial) + spatial_dims + else + (spatial_dims..., non_spatial...) end - if outtype isa Union || outtype_sums isa Union - outtype = outtype.b - outtype_sums = outtype_sums.b + # Build data array with all dimensions + all_sizes = length.(all_dims) + spatial_chunk = (chunk_length, chunk_length, 1) + all_chunks = if isempty(non_spatial) + spatial_chunk + else + (spatial_chunk..., ntuple(i -> non_spatial_sizes[i], length(non_spatial))...) end + data = TileArray{out_eltype}(missing, all_sizes, all_chunks) - sums = mapCube( - # mapCube can't find axes of other AbstractDimArrays e.g. Raster - YAXArray(dims(geo_array), geo_array.data, metadata(geo_array)); - indims=InDims(dims(geo_array, x_name), dims(geo_array, y_name)), - outdims=OutDims( - dggs_bbox..., - outtype=outtype_sums, - backend=backend, - path=path - ), kwargs...) do xout, xin - for ci in CartesianIndices(xin) - ismissing(xin[ci]) && continue - isnan(xin[ci]) && continue - - cell = cells[ci] - i_pos, j_pos, n_pos = cell.i + 1 - dggs_bbox[1][1], cell.j + 1 - dggs_bbox[2][1], cell.n + 1 - dggs_bbox[3][1] - if ismissing(xout[i_pos, j_pos, n_pos]) - xout[i_pos, j_pos, n_pos] = xin[ci] - else - xout[i_pos, j_pos, n_pos] += xin[ci] + # Pre-sized reusable buffer + buf = Vector{eltype(geo_array)}(undef, 32) + + # Helper: aggregate spatial pixels for one slice into data + function _fill_spatial!(data, cell_coords, geo_slice, buf) + for (k, v) in cell_coords + try + buf_len = 0 + for idx in v + val = geo_slice[idx] + if val !== missing + buf_len += 1 + if buf_len > length(buf) + resize!(buf, length(buf) * 2) + end + @inbounds buf[buf_len] = val + end + end + buf_len == 0 && continue + res = agg_func(@view buf[1:buf_len]) + data[Int(k.i)+1, Int(k.j)+1, Int(k.n)+1] = res + catch end - counts[i_pos, j_pos, n_pos] += 1 end end - means = sums.data ./ counts - data = if outtype <: Integer || outtype <: Union{Missing,Integer} - round.(means) + if isempty(non_spatial) + # Pure spatial case — original fast path + _fill_spatial!(data, cell_coords, geo_array, buf) else - means + # Iterate over all non-spatial dimension index combinations + for extra_ci in CartesianIndices(non_spatial_sizes) + extra_idxs = Tuple(extra_ci) + # Build indexer: spatial dims get Colon(), non-spatial get specific index + spatial_names_set = Set(spatial_dim_names) + idxs = map(dims(geo_array)) do d + if DD.name(d) in spatial_names_set + Colon() + else + # find which position this dim has in non_spatial + pos = findfirst(nd -> DD.name(nd) == DD.name(d), non_spatial) + extra_idxs[pos] + end + end + geo_slice = view(geo_array, idxs...) + + # Build data indexer: spatial dims get Colon(), non-spatial get specific index + data_idxs = (Colon(), Colon(), Colon(), extra_idxs...) + data_view = view(data, data_idxs...) + + _fill_spatial!(data_view, cell_coords, geo_slice, buf) + end end return DGGSArray( - data, dims(sums), refdims(sums), name, metadata(geo_array), - resolution, "ISEA4D.Penta", geo_bbox + data, all_dims, (), name, metadata(geo_array), + resolution, dggsrs, geo_bbox ) end + function to_dggs_array( geo_array::AbstractDimArray, resolution::Integer, crs::String, agg_func::Function; x_name=:X, y_name=:Y, kwargs... @@ -214,22 +257,18 @@ function to_dggs_array( properties = metadata(geo_array) delete!(properties, "projection") - cells = to_cell_array(x_dim, y_dim, resolution, crs) - - # get pixels to aggregate for each cell - cell_coords = Dict{eltype(cells),Vector{CartesianIndex{2}}}() - for cell_idx in CartesianIndices(cells) - cell = cells[cell_idx] - current_cells = get!(() -> CartesianIndex{2}[], cell_coords, cell) - push!(current_cells, cell_idx) - end - - dggs_bbox = get_dggs_bbox(keys(cell_coords)) + cell_coords = cells_to_coord_dict(x_dim, y_dim, resolution, crs) geo_bbox = get_geo_bbox(geo_array, crs) dggs_array = to_dggs_array( - geo_array, cells, cell_coords, dggs_bbox, geo_bbox, agg_func; - x_name=x_name, y_name=y_name, kwargs... + geo_array::AbstractDimArray, + resolution, + cell_coords, + geo_bbox::Extent, + agg_func::Function + ; + name=get_name(geo_array), + kwargs... ) return dggs_array end @@ -244,11 +283,14 @@ function to_dggs_array(geo_array::AbstractDimArray, resolution::Integer, crs::St properties = metadata(geo_array) - cells = to_cell_array(x_dim, y_dim, resolution, crs) - dggs_bbox = get_dggs_bbox(cells) + cell_coords = cells_to_coord_dict(x_dim, y_dim, resolution, crs) + + # Compute dggs_bbox from cell_coords keys + cells = keys(cell_coords) + dggs_bbox = get_dggs_bbox_cells(cells) geo_bbox = get_geo_bbox(geo_array, crs) - dggs_array = to_dggs_array(geo_array, cells, dggs_bbox, geo_bbox; x_name=x_name, y_name=y_name, kwargs...) + dggs_array = to_dggs_array(geo_array, resolution, cell_coords, geo_bbox; x_name=x_name, y_name=y_name, kwargs...) return dggs_array end @@ -257,7 +299,7 @@ function to_geo_array(dggs_array::DGGSArray, cells::AbstractDimArray; backend=:a lat_dim = dims(cells, :Y) # dggs_array may only contain parts of the world, having only parts of the dimension - get_extent(i_dim) = dggs_array.dims[i_dim].val.data |> x -> (first(x), last(x)) + get_extent(i_dim) = dggs_array.dims[i_dim].val |> x -> (first(x), last(x)) i_min, i_max = get_extent(1) j_min, j_max = get_extent(2) n_min, n_max = get_extent(3) @@ -363,7 +405,7 @@ function DGGSArray(array::AbstractDimArray) end function DGGSArray(resolution; chunk_length=2^12) - spatial_dims = (Dim{:dggs_i}(0:2*2^resolution-1), Dim{:dggs_j}(0:2^resolution-1), Dim{:dggs_n}(0:4)) + spatial_dims = (Dim{:dggs_i}(0:(2*2^resolution-1)), Dim{:dggs_j}(0:(2^resolution-1)), Dim{:dggs_n}(0:4)) data = TileArray{Union{Missing,Float64}}(missing, length.(spatial_dims), (chunk_length, chunk_length, 1)) dggsrs = "ISEA4D.Penta" bbox = Extent(X=(-180, 180), Y=(-90, 90)) @@ -410,6 +452,8 @@ Base.getindex(a::DGGSArray, c::Cell) = YAXArray(a)[dggs_i=At(c.i), dggs_j=At(c.j # DGGSArrays are usually big. Like YAXArrays, avoid DiskArray to load everything in memory Base.getindex(a::DGGSArray; i...) = view(a; i...) +Base.setindex!(a::DGGSArray, val, c::Cell) = YAXArray(a)[dggs_i=At(c.i), dggs_j=At(c.j), dggs_n=At(c.n)] = val + # # IO:: Serialization of DGGS Arrays @@ -423,10 +467,12 @@ function open_dggs_array(file_path::String) return DGGSArray(arr) end -function save_dggs_array(file_path::String, dggs_array::DGGSArray; kwargs...) - ds = Dataset(; Dict(DD.name(dggs_array) => YAXArray(dggs_array))...) - savedataset(ds; path=file_path, kwargs...) -end +""" + save_dggs_array(file_path, dggs_array; kwargs...) + +Save a DGGSArray to disk. This is a stub function that is extended by the DGGSZarr extension. +""" +function save_dggs_array end # # Operations diff --git a/src/datasets.jl b/src/datasets.jl index ac5c696..4046b14 100644 --- a/src/datasets.jl +++ b/src/datasets.jl @@ -33,38 +33,23 @@ function Base.getproperty(ds::DGGSDataset, s::Symbol) end -function to_dggs_dataset(geo_ds::Dataset, resolution::Integer, crs::String, agg_func::Function; metadata=Dict(), x_name=:X, y_name=:Y, kwargs...) - cells = to_cell_array(dims(geo_ds, x_name), dims(geo_ds, y_name), resolution, crs) - - # get pixels to aggregate for each cell - cell_coords = Dict{eltype(cells),Vector{CartesianIndex{2}}}() - for cell_idx in CartesianIndices(cells) - cell = cells[cell_idx] - current_cells = get!(() -> CartesianIndex{2}[], cell_coords, cell) - push!(current_cells, cell_idx) - end - - dggs_bbox = get_dggs_bbox(keys(cell_coords)) - - dggs_arrays = [] - dggs_arrays_lock = ReentrantLock() - Threads.@threads for (name, geo_array) in collect(geo_ds.cubes) - dggs_array = to_dggs_array(geo_array, cells, cell_coords, dggs_bbox, agg_func; name=name, x_name=x_name, y_name=y_name, kwargs...) - @lock dggs_arrays_lock push!(dggs_arrays, dggs_array) - end - return DGGSDataset(dggs_arrays...; metadata=metadata) -end +function to_dggs_dataset(geo_ds::Dataset, resolution::Integer, crs::String; agg_func::Function=mean, metadata=Dict(), x_name=:X, y_name=:Y, kwargs...) + # Fused algorithm: directly build cell coordinate dictionary from dimensions + # without creating an intermediate cell array + cell_coords = cells_to_coord_dict(geo_ds[x_name], geo_ds[y_name], resolution, crs) -"Fast iterative version only supporting mean" -function to_dggs_dataset(geo_ds::Dataset, resolution::Integer, crs::String; x_name=:X, y_name=:Y, metadata=Dict(), kwargs...) - cells = to_cell_array(geo_ds.axes[x_name], geo_ds.axes[y_name], resolution, crs) - dggs_bbox = get_dggs_bbox(cells) - geo_bbox = get_geo_bbox(geo_ds.cubes |> values |> first, crs; x_name=x_name, y_name=y_name) + # Compute geo_bbox once for the entire dataset instead of per-cube + # All cubes in a dataset share the same spatial extent + first_cube = first(geo_ds.cubes)[2] + geo_bbox = get_geo_bbox(first_cube, crs; x_name=x_name, y_name=y_name) dggs_arrays = [] dggs_arrays_lock = ReentrantLock() Threads.@threads for (name, geo_array) in collect(geo_ds.cubes) - dggs_array = to_dggs_array(geo_array, cells, dggs_bbox, geo_bbox; name=name, x_name=x_name, y_name=y_name, kwargs...) + dggs_array = to_dggs_array( + geo_array, resolution, cell_coords, geo_bbox; + agg_func=agg_func, name=name, x_name=x_name, y_name=y_name, kwargs... + ) @lock dggs_arrays_lock push!(dggs_arrays, dggs_array) end return DGGSDataset(dggs_arrays...; metadata=metadata) @@ -122,8 +107,11 @@ end open_dggs_dataset(file_path::String; kwargs...) = file_path |> x -> open_dataset(x; kwargs...) |> cache |> DGGSDataset -function save_dggs_dataset(file_path::String, dggs_ds::DGGSDataset; kwargs...) - dggs_ds |> Dataset |> x -> savedataset(x; path=file_path, kwargs...) -end +""" + save_dggs_dataset(file_path, dggs_ds; kwargs...) + +Save a DGGSDataset to disk. This is a stub function that is extended by the DGGSZarr extension. +""" +function save_dggs_dataset end -init_global_dggs_dataset() = @error("Please load package Zarr") \ No newline at end of file +init_global_dggs_dataset() = @error("Please load package Zarr") diff --git a/src/pyramids.jl b/src/pyramids.jl index b55cbdb..62d02c7 100644 --- a/src/pyramids.jl +++ b/src/pyramids.jl @@ -1,3 +1,36 @@ +""" + coarsen(A::AbstractArray, factors::Tuple) + +Coarsen an array by aggregating blocks of elements. Each dimension is reduced +by the corresponding factor using mean aggregation. + +# Arguments +- `A::AbstractArray`: Input array to coarsen +- `factors::Tuple`: Tuple of coarsening factors, one per dimension. + Use `1` to keep a dimension unchanged. + +# Example +```julia +a = rand(64, 32, 10) +coarse_a = coarsen(a, (2, 2, 1)) # Result: 32×16×10 +``` +""" +function coarsen(A::AbstractArray, factors::Tuple) + # Build the reshaped dimensions: interleave (new_dim, factor) pairs + reshaped_dims = Int[] + for (s, f) in zip(size(A), factors) + push!(reshaped_dims, s ÷ f) + push!(reshaped_dims, f) + end + + reshaped = reshape(A, Tuple(reshaped_dims)) + + # Average over the factor dimensions (every even dimension: 2, 4, 6, ...) + reduce_dims = Tuple(2:2:length(reshaped_dims)) + result = mean(reshaped, dims=reduce_dims) + return dropdims(result, dims=reduce_dims) +end + function DGGSPyramid(data::AbstractDict{T,A}, dggsrs, bbox) where {T,A<:DGGSDataset} dimtree = DimTree() # add all res levels as branches @@ -67,52 +100,130 @@ end function aggregate_by_factor( xin::AbstractArray, xout::AbstractArray, - pyramid_agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean + agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean ) fac = ceil(Int, size(xin, 1) / size(xout, 1)) for j in axes(xout, 2) for i in axes(xout, 1) - xview = ((i-1)*fac+1):min(size(xin, 1), (i * fac)) - yview = ((j-1)*fac+1):min(size(xin, 2), (j * fac)) - xout[i, j] = pyramid_agg_func(view(xin, xview, yview)) + xview = ((i-1)*fac+1):min(size(xin, 1), (i*fac)) + yview = ((j-1)*fac+1):min(size(xin, 2), (j*fac)) + xout[i, j] = agg_func(view(xin, xview, yview)) end end end +""" + coarsen(dggs_array::DGGSArray{<:Any,<:Any,<:Any,<:Any,<:TileArray}; agg_func) + +Coarsen a DGGSArray backed by a TileArray by aggregating 2x2 blocks. +Missing tiles are skipped entirely, making this efficient for sparse arrays. +""" function coarsen( - dggs_array::DGGSArray; - pyramid_agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean + dggs_array::DGGSArray{<:Any,<:Any,<:Any,<:Any,<:TileArray}; + agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean ) + tile_array = dggs_array.data coarser_level = dggs_array.resolution - 1 - # analog to GeoTIFF overviews - coarser_dims = map((:dggs_i, :dggs_j)) do dim - dim_min, dim_max = dims(dggs_array, dim) |> extrema - dim_min = floor(dim_min / 2) |> Int - dim_max = floor(dim_max / 2) |> Int - Dim{dim}(dim_min:dim_max) - end + # Get dimension extents (0-based DGGS coordinates) + i_min, i_max = extrema(dims(dggs_array, :dggs_i)) + j_min, j_max = extrema(dims(dggs_array, :dggs_j)) + n_min, n_max = extrema(dims(dggs_array, :dggs_n)) + + # Compute coarser dimensions (0-based) + coarser_i_min = floor(Int, i_min / 2) + coarser_i_max = floor(Int, i_max / 2) + coarser_j_min = floor(Int, j_min / 2) + coarser_j_max = floor(Int, j_max / 2) + + # Output TileArray dimensions (1-based extent) + out_dims = ( + coarser_i_max - coarser_i_min + 1, + coarser_j_max - coarser_j_min + 1, + n_max - n_min + 1 + ) + + # Create output TileArray with same chunk size + out_eltype = eltype(tile_array) + out_tile_array = TileArray{out_eltype}(missing, out_dims, tile_array.chunk_size) + + cs = tile_array.chunk_size + + # Get present tile ranges (1-based internal indices) + present_ranges = ranges(tile_array) + + # Process each present tile + for tile_range in present_ranges + i_range, j_range, n_range = tile_range - coarser_arr = mapCube( - dggs_array; - indims=InDims(:dggs_i, :dggs_j), - outdims=OutDims(coarser_dims...) - ) do xout, xin - xout = aggregate_by_factor(xin, xout, pyramid_agg_func) + # Get chunk data using 1-based chunk indices + ci = div(first(i_range) - 1, cs[1]) + 1 + cj = div(first(j_range) - 1, cs[2]) + 1 + cn = div(first(n_range) - 1, cs[3]) + 1 + chunk_data = tile_array.data[ci, cj, cn] + chunk_data === missing && continue + + # Global 0-based coordinate of chunk start + global_i_start_0 = first(i_range) - 1 + i_min + global_j_start_0 = first(j_range) - 1 + j_min + + # Iterate over 2x2 blocks aligned to the global grid + for li in 1:length(i_range) + gi_0 = global_i_start_0 + li - 1 + gi_0 % 2 != 0 && continue # Skip non-aligned positions + + for lj in 1:length(j_range) + gj_0 = global_j_start_0 + lj - 1 + gj_0 % 2 != 0 && continue + + # Collect values from 2x2 block + values = out_eltype[] + for di in 0:1, dj in 0:1 + ni, nj = li + di, lj + dj + if ni <= length(i_range) && nj <= length(j_range) + for ln in 1:length(n_range) + val = chunk_data[ni, nj, ln] + val !== missing && push!(values, val) + end + end + end + + if !isempty(values) + agg_val = agg_func(values) + + # Write to output (convert to 1-based) + out_i_1 = div(gi_0, 2) - coarser_i_min + 1 + out_j_1 = div(gj_0, 2) - coarser_j_min + 1 + for n_idx in n_range + out_tile_array[out_i_1, out_j_1, n_idx] = agg_val + end + end + end + end end + # Build coarser DGGSArray directly (bypass YAXArray to preserve TileArray) + coarser_dims = ( + Dim{:dggs_i}(coarser_i_min:coarser_i_max), + Dim{:dggs_j}(coarser_j_min:coarser_j_max), + Dim{:dggs_n}(n_min:n_max) + ) + properties = Dict{String,Any}(metadata(dggs_array)) properties["dggs_dggsrs"] = dggs_array.dggsrs properties["dggs_resolution"] = coarser_level properties["dggs_bbox"] = dggs_array.bbox - coarser_dggs_arr = YAXArray(dims(coarser_arr), coarser_arr.data, properties) |> DGGSArray - coarser_dggs_arr = rebuild(coarser_dggs_arr; name=name(dggs_array)) + coarser_dggs_arr = DGGSArray( + out_tile_array, coarser_dims, (), name(dggs_array), properties, + coarser_level, dggs_array.dggsrs, dggs_array.bbox + ) return coarser_dggs_arr end + function coarsen(dggs_ds::DGGSDataset; kwargs...) coarser_arrays = [] for key in keys(dggs_ds) @@ -127,7 +238,7 @@ end function to_dggs_pyramid(dggs_ds::DGGSDataset; kwargs...) pyramid = DGGSDataset[] push!(pyramid, dggs_ds) - for resolution in dggs_ds.resolution-1:-1:1 + for resolution in (dggs_ds.resolution-1):-1:1 current_dggs_ds = pyramid[end] coarser_ds = coarsen(current_dggs_ds; kwargs...) push!(pyramid, coarser_ds) @@ -143,24 +254,15 @@ function to_dggs_pyramid(dggs_array::DGGSArray; kwargs...) return pyramid end -function to_dggs_pyramid( - geo_ds::YAXArrays.Dataset, resolution::Integer, crs::String, agg_func::Function; - kwargs... -) - dggs_ds = to_dggs_dataset(geo_ds, resolution, crs, agg_func; kwargs...) - dggs_pyramid = to_dggs_pyramid(dggs_ds) - return dggs_pyramid -end - function to_dggs_pyramid( geo_ds::YAXArrays.Dataset, resolution::Integer, crs::String; - pyramid_agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean, + agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean, kwargs... ) - dggs_ds = to_dggs_dataset(geo_ds, resolution, crs; kwargs...) - dggs_pyramid = to_dggs_pyramid(dggs_ds; pyramid_agg_func=pyramid_agg_func) + dggs_ds = to_dggs_dataset(geo_ds, resolution, crs; agg_func=agg_func, kwargs...) + dggs_pyramid = to_dggs_pyramid(dggs_ds; agg_func=agg_func) return dggs_pyramid end @@ -168,11 +270,11 @@ function to_dggs_pyramid( geo_array::YAXArrays.YAXArray, resolution::Integer, crs::String; - pyramid_agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean, + agg_func::Function=x -> filter(y -> !ismissing(y) && !isnan(y), x) |> mean, kwargs... ) - dggs_array = to_dggs_array(geo_array, resolution, crs; kwargs...) - dggs_pyramid = to_dggs_pyramid(dggs_array; pyramid_agg_func=pyramid_agg_func) + dggs_array = to_dggs_array(geo_array, resolution, crs; agg_func=agg_func, kwargs...) + dggs_pyramid = to_dggs_pyramid(dggs_array; agg_func=agg_func) return dggs_pyramid end diff --git a/src/tile-array.jl b/src/tile-array.jl index 21f1d23..725a855 100644 --- a/src/tile-array.jl +++ b/src/tile-array.jl @@ -68,7 +68,7 @@ function Base.iterate(A::TileArray, state=1) if state > length(A) return nothing end - I = ntuple(i -> div(state - 1, prod(A.dims[i+1:end])) % A.dims[i] + 1, length(A.dims)) + I = ntuple(i -> div(state - 1, prod(A.dims[(i+1):end])) % A.dims[i] + 1, length(A.dims)) return (getindex(A, I...), state + 1) end @@ -76,4 +76,384 @@ function Base.show(io::IO, ::MIME"text/plain", a::TileArray) print(io, join(a.dims, "x")) print(io, " ") print(io, typeof(a)) +end + +""" + ranges(A::TileArray) + +Return a vector of tuples of UnitRanges representing the index ranges of all +non-missing tiles in `A`. Each tuple corresponds to one present chunk and +contains the global index range along each dimension. + +# Example +```julia +a = DGGSArray(19) +a.data[1:4096, 1:4096, 1] .= 1 +a.data[1:4096, 1:4096, 2] .= 2 +ds = DGGSDataset(a) +tile_array = ds.a.data +ranges(tile_array) # returns e.g. [(1:4096, 1:4096, 1:1), (1:4096, 1:4096, 2:2)] +``` +""" +function ranges(A::TileArray{T,N}) where {T,N} + result = Tuple{Vararg{UnitRange{Int},N}}[] + for ck in CartesianIndices(size(A.data)) + if A.data[ck] !== missing + ranges_tuple = ntuple(N) do d + start_idx = (ck[d] - 1) * A.chunk_size[d] + 1 + end_idx = min(ck[d] * A.chunk_size[d], A.dims[d]) + start_idx:end_idx + end + push!(result, ranges_tuple) + end + end + return result +end + +""" + _tile_reduce_contribution(op, base_result, n::Int) + +Given `base_result = f(x)` and the fact that value `x` appears `n` times, +return the reduction of these `n` identical values: `reduce(op, fill(base_result, n))`. +""" +function _tile_reduce_contribution(op, base_result, n::Int) + n <= 0 && throw(DomainError(n, "n must be positive")) + n == 1 && return base_result + if op === + + return n * base_result + end + result = base_result + remaining = n - 1 + while remaining > 0 + result = op(result, base_result) + remaining -= 1 + end + result +end + +# Actual size of chunk `ck_d` (1-based) along dimension `d` of array `A`, +# accounting for boundary chunks that may be smaller than `chunk_size[d]`. +@inline _tile_actual_chunk_dim(A, d, ck_d::Int) = + min(A.chunk_size[d], A.dims[d] - (ck_d - 1) * A.chunk_size[d]) + +""" + Base.mapreduce(f, op, A::TileArray; dims=:, init) + +Efficiently compute `mapreduce(f, op, A)` over a `TileArray`. + +Missing tiles (which contain only `A.default`) are never materialized: +their contribution is computed analytically from the count of missing +elements and the value `f(A.default)`. + +Supported forms: +- `mapreduce(f, op, A::TileArray)` — full reduction to a scalar. +- `mapreduce(f, op, A::TileArray; dims=...)` — reduction along `dims`, + returning an array of shape `(d -> d in dims ? 1 : size(A, d))`. +- `init` optionally seeds the reduction; with `init=Val(:init)` a value can be + passed as the fourth positional argument following Base's convention. +""" +function Base.mapreduce(f, op, A::TileArray{T,N}; dims=:, init=nothing) where {T,N} + # When default is missing, missing tiles are "empty" and should be skipped + A.default === missing && return _mapreduce_tile_skipmissing(A, f, op, dims, init) + _mapreduce_tile(A, f, op, dims, init) +end + +# Internal: mapreduce for TileArray where A.default === missing +# Skip all missing tiles entirely and treat missing elements within chunks as absent +function _mapreduce_tile_skipmissing(A::TileArray{T,N}, f, op, dims, init) where {T,N} + chunk_dims_tuple = size(A.data) + + if dims === (:) + # -------------- Full reduction to scalar -------------- + acc = nothing + for ck in CartesianIndices(chunk_dims_tuple) + chunk = A.data[ck] + chunk === missing && continue # skip missing tiles + for v in chunk + v === missing && continue # skip missing elements within chunks + val = f(v) + acc = acc === nothing ? val : op(acc, val) + end + end + if acc === nothing + return init === nothing ? + throw(ArgumentError("reducing over an empty collection with no initial value is not allowed")) : + init + end + return init === nothing ? acc : op(init, acc) + end + + # -------------- Reduction along `dims` to an array -------------- + dims_tuple = dims isa Union{Tuple,AbstractVector} ? Tuple(dims) : (dims,) + out_dims = ntuple(d -> d in dims_tuple ? 1 : A.dims[d], N) + OutEltype = if init === nothing + typeof(first(_nonmiss_from_present_chunk(A))) + else + typeof(op(init, first(_nonmiss_from_present_chunk(A)))) + end + + present_acc = Array{Union{Nothing,OutEltype},N}(undef, out_dims...) + fill!(present_acc, nothing) + seen_present = fill(false, out_dims...) + + for ck in CartesianIndices(chunk_dims_tuple) + chunk = A.data[ck] + chunk === missing && continue + local_shape = ntuple(d -> _tile_actual_chunk_dim(A, d, ck[d]), N) + for li in CartesianIndices(local_shape) + v = chunk[li] + v === missing && continue + li_vals = Tuple(li) + gi = ntuple(d -> (ck[d] - 1) * A.chunk_size[d] + li_vals[d], N) + out_idx = ntuple(d -> d in dims_tuple ? 1 : gi[d], N) + val = f(v) + if !seen_present[out_idx...] + present_acc[out_idx...] = val + seen_present[out_idx...] = true + else + present_acc[out_idx...] = op(present_acc[out_idx...], val) + end + end + end + + acc = Array{OutEltype,N}(undef, out_dims...) + for out_idx in CartesianIndices(out_dims...) + p = present_acc[out_idx] + if init === nothing + if p === nothing + acc[out_idx] = zero(OutEltype) + else + acc[out_idx] = p + end + else + acc[out_idx] = p !== nothing ? op(init, p) : init + end + end + return acc +end + +# Handle `mapreduce` over `skipmissing(a)`. Missing tiles (chunks of only `missing`) +# are skipped entirely — we only fold over elements from present chunks. +function Base.mapreduce(f, op, itr::Base.SkipMissing{<:TileArray{T,N}}; dims=:, init=nothing) where {T,N} + A = itr.x + chunk_dims_tuple = size(A.data) + + if dims === (:) + acc = nothing + for ck in CartesianIndices(chunk_dims_tuple) + chunk = A.data[ck] + chunk === missing && continue # whole tile is missing — skip + for v in chunk + v === missing && continue + val = f(v) + acc = acc === nothing ? val : op(acc, val) + end + end + if acc === nothing + return init === nothing ? + throw(ArgumentError("reducing over an empty collection with no initial value is not allowed")) : + init + end + return init === nothing ? acc : op(init, acc) + end + + # Reduction along dims — delegate to generic TileArray reduction that skips missing values + dims_tuple = dims isa Union{Tuple,AbstractVector} ? Tuple(dims) : (dims,) + out_dims = ntuple(d -> d in dims_tuple ? 1 : A.dims[d], N) + OutEltype = if init === nothing + typeof(first(_nonmiss_from_present_chunk(A))) + else + typeof(op(init, first(_nonmiss_from_present_chunk(A)))) + end + miss_count = zeros(Int, out_dims...) # not used for skipmissing but kept for uniformity + present_acc = Array{Union{Nothing,OutEltype},N}(undef, out_dims...) + fill!(present_acc, nothing) + seen_present = fill(false, out_dims...) + + for ck in CartesianIndices(chunk_dims_tuple) + chunk = A.data[ck] + chunk === missing && continue + local_shape = ntuple(d -> _tile_actual_chunk_dim(A, d, ck[d]), N) + for li in CartesianIndices(local_shape) + v = chunk[li] + v === missing && continue + li_vals = Tuple(li) + gi = ntuple(d -> (ck[d] - 1) * A.chunk_size[d] + li_vals[d], N) + out_idx = ntuple(d -> d in dims_tuple ? 1 : gi[d], N) + val = f(v) + if !seen_present[out_idx...] + present_acc[out_idx...] = val + seen_present[out_idx...] = true + else + present_acc[out_idx...] = op(present_acc[out_idx...], val) + end + end + end + + acc = Array{OutEltype,N}(undef, out_dims...) + for out_idx in CartesianIndices(out_dims...) + p = present_acc[out_idx] + if init === nothing + if p === nothing + acc[out_idx] = zero(OutEltype) # shouldn't happen for skipmissing with data + else + acc[out_idx] = p + end + else + acc[out_idx] = p !== nothing ? op(init, p) : init + end + end + return acc +end + +# Helper: find first non-missing element from any present chunk (used only for type inference) +function _nonmiss_from_present_chunk(A::TileArray{T,N}) where {T,N} + for chunk in A.data + chunk === missing && continue + for v in chunk + v === missing || return Some(v) + end + end + # Fallback if all present chunks are missing (shouldn't happen in practice) + return Some(A.default) +end +Base.first(s::Some) = s.value + +# Internal entry point. The `kw...` splat lets us also accept the Base-style +# four-argument form `mapreduce(f, op, ::InitialValueOperator, A)` but here +# we only handle the keyword-based `init`. +function _mapreduce_tile(A::TileArray{T,N}, f, op, dims, init) where {T,N} + f_default = f(A.default) + chunk_dims_tuple = size(A.data) + skip_missing = f_default === missing + + if dims === (:) + # -------------- Full reduction to scalar -------------- + present_acc = nothing + n_missing_elems = 0 + for ck in CartesianIndices(chunk_dims_tuple) + chunk = A.data[ck] + if chunk === missing + n_missing_elems += prod(_tile_actual_chunk_dim(A, d, ck[d]) for d in 1:N) + else + if skip_missing + # Skip missing elements within the chunk + for v in chunk + v === missing && continue + val = f(v) + present_acc = present_acc === nothing ? val : op(present_acc, val) + end + else + chunk_val = mapreduce(f, op, chunk) + present_acc = present_acc === nothing ? chunk_val : op(present_acc, chunk_val) + end + end + end + result = if n_missing_elems > 0 && !skip_missing + mc = _tile_reduce_contribution(op, f_default, n_missing_elems) + present_acc === nothing ? mc : op(present_acc, mc) + else + present_acc + end + if result === nothing + if init === nothing + throw(ArgumentError( + "reducing over an empty collection with no initial value is not allowed")) + end + return init + end + return init === nothing ? result : op(init, result) + end + + # -------------- Reduction along `dims` to an array -------------- + dims_tuple = dims isa Union{Tuple,AbstractVector} ? Tuple(dims) : (dims,) + out_dims = ntuple(d -> d in dims_tuple ? 1 : A.dims[d], N) + + # Infer the output element type. + OutEltype = if init === nothing + typeof(f_default) + else + typeof(op(init, f_default)) + end + + # Accumulators: per output cell, store (i) count of `default` elements, + # (ii) reduction over values from present chunks. The latter is `nothing` + # until at least one present element folds into it. + miss_count = zeros(Int, out_dims...) + present_acc = Array{Union{Nothing,OutEltype},N}(undef, out_dims...) + fill!(present_acc, nothing) + + # Scatter pass: visit each chunk once. + # For a missing chunk covering non-reduced element indices `I_nr`, each + # output cell `out_idx` indexed by `I_nr` gets `n_miss = prod(actual_size + # along reduced dims)` elements of `A.default`. For a present chunk, we + # iterate every element and fold `f(chunk[I])` into the corresponding + # output cell. + # + # For missing chunks that span many non-reduced output cells, this is + # cheaper than the gather approach because we don't re-iterate the + # reduced-dim chunk range for each output cell. + for ck in CartesianIndices(size(A.data)) + chunk = A.data[ck] + # Reduced-dim actual sizes of this chunk: + n_miss = prod(_tile_actual_chunk_dim(A, d, ck[d]) for d in dims_tuple) + if chunk === missing + # Iterate over non-reduced element indices within this chunk. + # For each such index, the output-cell index has non-reduced coords + # equal to the global element index and reduced coords equal to 1. + non_red_ranges = ntuple(d -> d in dims_tuple ? (1:1) : + (1:_tile_actual_chunk_dim(A, d, ck[d])), N) + for li_nr in CartesianIndices(non_red_ranges) + li_vals = Tuple(li_nr) + out_idx = ntuple(d -> d in dims_tuple ? 1 : + (ck[d] - 1) * A.chunk_size[d] + li_vals[d], N) + @inbounds miss_count[out_idx...] += n_miss + end + else + # Present chunk: fold contributions into per-output-cell accumulator. + local_shape = ntuple(d -> _tile_actual_chunk_dim(A, d, ck[d]), N) + for li in CartesianIndices(local_shape) + li_vals = Tuple(li) + gi = ntuple(d -> (ck[d] - 1) * A.chunk_size[d] + li_vals[d], N) + out_idx = ntuple(d -> d in dims_tuple ? 1 : gi[d], N) + val = f(chunk[li]) + prev = present_acc[out_idx...] + @inbounds present_acc[out_idx...] = prev === nothing ? val : op(prev, val) + end + end + end + + # Final pass: merge miss_count and present_acc into the output array. + # Skip missing-tile contributions when f_default === missing (they're "empty" tiles). + skip_missing_default = f_default === missing + acc = Array{OutEltype,N}(undef, out_dims...) + for out_idx in CartesianIndices(out_dims) + m = miss_count[out_idx] + p = present_acc[out_idx] + if skip_missing_default + m = 0 # treat missing-tile elements as absent + end + if init === nothing + if m == 0 && p === nothing + throw(ArgumentError( + "reducing over an empty collection with no initial value is not allowed")) + elseif m == 0 + acc[out_idx] = p + elseif p === nothing + acc[out_idx] = _tile_reduce_contribution(op, f_default, m) + else + acc[out_idx] = op(_tile_reduce_contribution(op, f_default, m), p) + end + else + base = init + if m > 0 + base = op(base, _tile_reduce_contribution(op, f_default, m)) + end + if p !== nothing + base = op(base, p) + end + acc[out_idx] = base + end + end + return acc end \ No newline at end of file diff --git a/src/types.jl b/src/types.jl index d5b7601..59126b3 100644 --- a/src/types.jl +++ b/src/types.jl @@ -50,7 +50,7 @@ struct DGGSArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Na,Me} <: AbstractD end "Set of DGGSArrays with aligned and shared dimensions at the same resolution." -struct DGGSDataset{K,T,N,L,D<:Tuple,R<:Tuple,LD,M,LM} <: AbstractDimStack{K,T,N,L} +struct DGGSDataset{K,T,N,L,D<:Tuple,R<:Tuple,LD,M,LM} <: AbstractDimStack{K,T,N,L,D} # DimStack fields data::L dims::D diff --git a/test/runtests.jl b/test/runtests.jl index 8a261fe..d5a7e3b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -37,8 +37,8 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) @testset "DGGSArray" begin resolution = 3 - i_dim = Dim{:dggs_i}(0:2*2^resolution-1) - j_dim = Dim{:dggs_j}(0:2^resolution-1) + i_dim = Dim{:dggs_i}(0:(2*2^resolution-1)) + j_dim = Dim{:dggs_j}(0:(2^resolution-1)) n_dim = Dim{:dggs_n}(0:4) time_dim = Ti(1:10) dim_array = rand(i_dim, j_dim, n_dim, time_dim) @@ -76,6 +76,35 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) @test YAXArray(TileArray(0, (100, 100), (10, 10))) isa YAXArray end + @testset "TileArray mapreduce" begin + # Test mapreduce with non-missing default + a = TileArray(0, (10, 10), (5, 5)) + @test mapreduce(identity, +, a) == 0 + @test mapreduce(x -> x + 1, +, a) == 100 + a[1, 1] = 42 + @test mapreduce(identity, +, a) == 42 + @test maximum(a) == 42 + @test minimum(a) == 0 + + # Test mapreduce with missing default + b = TileArray{Union{Missing,Int}}(missing, (10, 10, 10), (5, 5, 5)) + b[1, 1, 1] = 2 + b[2, 1, 1] = 3 + @test maximum(b) == 3 + @test minimum(skipmissing(b)) == 2 + @test sum(skipmissing(b)) == 5 + @test maximum(skipmissing(b)) == 3 + @test b |> skipmissing |> maximum == 3 + + # Test mapreduce with dims + c = TileArray(1, (4, 6), (4, 6)) + c[1, 1] = 10 + c[2, 3] = 5 + dense = [c[i, j] for i in 1:4, j in 1:6] + @test mapreduce(identity, +, c; dims=1) == sum(dense; dims=1) + @test mapreduce(identity, +, c; dims=2) == sum(dense; dims=2) + end + @testset "Coordinate transformations" begin resolution = 20 geo_points = [(lon, lat) for lat in -90:5:90 for lon in -180:5:180] @@ -88,19 +117,19 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) @test sum(dists .< 10) / length(dists) >= 0.99 # cell ids must be in bounds - @test all(map(x -> x.i in 0:2*2^resolution-1, cell_ids)) - @test all(map(x -> x.j in 0:2^resolution-1, cell_ids)) + @test all(map(x -> x.i in 0:(2*2^resolution-1), cell_ids)) + @test all(map(x -> x.j in 0:(2^resolution-1), cell_ids)) @test all(map(x -> x.n in 0:4, cell_ids)) end @testset "Integer index" begin resolution = 5 - cells = [Cell(i, j, n, resolution) for n in 0:4 for j in 0:2^resolution-1 for i in 0:2*2^resolution-1] + cells = [Cell(i, j, n, resolution) for n in 0:4 for j in 0:(2^resolution-1) for i in 0:(2*2^resolution-1)] cells_int = Int64.(cells) cells2 = Cell.(cells_int, resolution) @test length(cells) == length(cells_int |> unique) - @test cells_int == 0:length(cells)-1 + @test cells_int == 0:(length(cells)-1) @test cells == cells2 end @@ -125,7 +154,7 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) @test dggs_array3 isa DGGSArray # other agg_func - dggs_array4 = to_dggs_array(geo_array3, 10, projection, median) + dggs_array4 = to_dggs_array(geo_array3, 10, projection; agg_func=median) @test dggs_array4 isa DGGSArray end @@ -160,8 +189,8 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) @testset "DGGSDataset" begin resolution = 3 - i_dim = Dim{:dggs_i}(0:2*2^resolution-1) - j_dim = Dim{:dggs_j}(0:2^resolution-1) + i_dim = Dim{:dggs_i}(0:(2*2^resolution-1)) + j_dim = Dim{:dggs_j}(0:(2^resolution-1)) n_dim = Dim{:dggs_n}(0:4) time_dim = Ti(1:10) dim_array = rand(i_dim, j_dim, n_dim, time_dim) @@ -195,22 +224,30 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) @test length(dggs_p.branches) == dggs_ds.resolution @test dggs_p.dggsrs == dggs_ds.dggsrs @test dggs_p.bbox == dggs_ds.bbox - @test dggs_p.dggs_s3 == dggs_p[3] - - # save and open pyramid - temp_dir = tempname() * ".dggs.zarr" - @info temp_dir - save_dggs_pyramid(temp_dir, dggs_p) - dggs_p2 = open_dggs_pyramid(temp_dir) - @test dggs_p.bbox == dggs_p2.bbox - @test dggs_p.dggsrs == dggs_p2.dggsrs - @test length(dggs_p.data) == length(dggs_p2.data) - @test all(keys(dggs_p.data) .== keys(dggs_p2.data)) - - # both layers must be present after save and open - @test name(dggs_p.dggs_s3.air_temperature) == name(dggs_p2.dggs_s3.air_temperature) - @test name(dggs_p.dggs_s3.precipitation) == name(dggs_p2.dggs_s3.precipitation) - rm(temp_dir, recursive=true) + @test dggs_p.dggs_s3.resolution == dggs_p[3].resolution + + @testset "all values are present" begin + data = collect(dggs_p[3].precipitation) + for n in 1:5 + @test length(data[:, :, n]) == length(filter(!ismissing, data[:, :, n])) + end + end + + @testset "save and open pyramid" begin + temp_dir = tempname() * ".dggs.zarr" + @info temp_dir + save_dggs_pyramid(temp_dir, dggs_p) + dggs_p2 = open_dggs_pyramid(temp_dir) + @test dggs_p.bbox == dggs_p2.bbox + @test dggs_p.dggsrs == dggs_p2.dggsrs + @test length(dggs_p.data) == length(dggs_p2.data) + @test all(keys(dggs_p.data) .== keys(dggs_p2.data)) + + # both layers must be present after save and open + @test name(dggs_p.dggs_s3.air_temperature) == name(dggs_p2.dggs_s3.air_temperature) + @test name(dggs_p.dggs_s3.precipitation) == name(dggs_p2.dggs_s3.precipitation) + rm(temp_dir, recursive=true) + end # pyramid from just one array dggs_p2 = to_dggs_pyramid(dggs_array)