diff --git a/Project.toml b/Project.toml index 06ce1a6..8f3b65b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,19 +1,24 @@ name = "DGGS" uuid = "1b1208a4-f9dd-4606-8704-9f68bba98fd1" -authors = ["Daniel Loos"] version = "2.0.0-DEV" +authors = ["Daniel Loos"] [deps] CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0" DiskArrayTools = "fcd2136c-9f69-4db6-97e5-f31981721d63" +DiskArrays = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3" Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b" +LRUCache = "8ac3fa9e-de4c-5943-b1dc-09c6b5f20637" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Proj = "c94c279d-25a6-4763-9509-64d165bea63e" +SphericalSpatialTrees = "86b7679e-e8c0-4ea9-ad3d-baef5c69eb57" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" YAXArrays = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c" @@ -29,10 +34,15 @@ DGGSZarr = "Zarr" [compat] CoordinateTransformations = "0.6.4" DiskArrayTools = "0.1.12" +DiskArrays = "0.4.18" Extents = "0.1.5" +FillArrays = "1.15.0" GeometryBasics = "0.5.9" Infiltrator = "1.8.8" +LRUCache = "1.6.2" OrderedCollections = "1.8.1" +ProgressMeter = "1.11.0" +SphericalSpatialTrees = "0.1.0" StaticArrays = "1.9.13" YAXArrays = "0.6.1" julia = "1.6.7" diff --git a/ext/DGGSMakie.jl b/ext/DGGSMakie.jl index 4ea0192..801d101 100644 --- a/ext/DGGSMakie.jl +++ b/ext/DGGSMakie.jl @@ -150,14 +150,13 @@ 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[]) - cb_limits = (minimum(filtered_data), maximum(filtered_data)) - + length(filtered_data) > 0 || error("All values are missing or NaN, cannot plot.") label = if dggs isa DGGSPyramid && length(args) == 1 String(args[1]) else DD.label(dggs) end - + cb_limits = (minimum(filtered_data), maximum(filtered_data)) cb = Colorbar(fig[1, 2], width=10, colormap=:viridis, limits=cb_limits; label=label) heatmap!(ax, data, colorrange=cb_limits) else diff --git a/ext/DGGSZarr.jl b/ext/DGGSZarr.jl index 39907af..88d33ae 100644 --- a/ext/DGGSZarr.jl +++ b/ext/DGGSZarr.jl @@ -7,6 +7,8 @@ using Infiltrator using Extents using DimensionalData import DimensionalData as DD +using DiskArrays +using FillArrays # Currently, YAXArrays does not support saving the experimental nested DimTree type @@ -50,4 +52,60 @@ function DGGS.open_dggs_pyramid(group::ZGroup) res = DGGSPyramid(dimtree, dggsrs, bbox) return res end + + +""" +Need to init globally: +- allows parallel read and write +- ij chunks often don't overlap eith neoighboring chunks on different dggs_n quad +- empty chunks are not stored on disk +""" +function DGGS.init_global_dggs_dataset( + geo_ds::Dataset, resolution, crs, path; + x_dim_name=:X, y_dim_name=:Y, chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1), kwargs... +) + # extract spatial dimensions + all_dims = [] + for (k, c) in geo_ds.cubes + append!(all_dims, dims(c)) + end + x_dim = filter(x -> name(x) == x_dim_name, all_dims)[1] + y_dim = filter(x -> name(x) == y_dim_name, all_dims)[1] + bbox = DGGS.get_geo_bbox(x_dim, y_dim, crs) + + properties = Dict( + "dggs_resolution" => resolution, + "dggs_dggsrs" => "ISEA4D.Penta", + "dggs_bbox" => bbox + ) + x_dim_name in keys(geo_ds.axes) || error("x_dim_name :$(x_dim_name) not found in geo_ds") + y_dim_name in keys(geo_ds.axes) || error("y_dim_name :$(y_dim_name) not found in geo_ds") + x_dim_name != y_dim_name || error("X and Y names must be different") + + arrays = Dict() + 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)) + non_spatial_dims = filter(x -> !(name(x) in [x_dim_name, y_dim_name]), geo_array.axes) + dims = (spatial_dims..., non_spatial_dims...) + else + spatial_dims = () + non_spatial_dims = geo_array.axes + dims = (spatial_dims..., non_spatial_dims...) + end + + data = Zeros(eltype(geo_array), length.(dims)) + yax_array = YAXArray(dims, data, properties) + yax_array = rebuild(yax_array; name=key) + yax_array = setchunks(yax_array, chunks) + + arrays[key] = yax_array + end + + ds = Dataset(; properties, arrays...) + ds = savedataset(ds; path=path, skeleton=true, driver=:zarr, kwargs...) + res = open_dataset(zopen(path, "w"); driver=:zarr) |> DGGSDataset + return res end +end \ No newline at end of file diff --git a/src/DGGS.jl b/src/DGGS.jl index 7d9906d..2c7d64f 100644 --- a/src/DGGS.jl +++ b/src/DGGS.jl @@ -13,7 +13,11 @@ using Printf using Infiltrator using Extents using OrderedCollections +using DiskArrays using DiskArrayTools +using FillArrays +using LRUCache +using ProgressMeter include("types.jl") include("cells.jl") @@ -21,10 +25,6 @@ include("arrays.jl") include("datasets.jl") include("pyramids.jl") -const transformations = Channel{Proj.Transformation}(Inf) -const inv_transformations = Channel{Proj.Transformation}(Inf) -const threads_ready = Ref(false) - export Cell, DGGSArray, DGGSDatase, Cell, DGGSArray, DGGSDataset, DGGSPyramid export to_cell, to_geo, to_dggs_pyramid, to_dggs_dataset, to_dggs_array, to_geo_dataset, to_geo_array export open_dggs_array, open_dggs_dataset, open_dggs_pyramid diff --git a/src/arrays.jl b/src/arrays.jl index 06b9192..5cec6a2 100644 --- a/src/arrays.jl +++ b/src/arrays.jl @@ -1,242 +1,129 @@ -function get_dggs_bbox(cells) - cell = first(cells) - resolution = cell.resolution - - # start with smallest possible bbox - i_min = cell.i - i_max = cell.i - - j_min, j_max = cell.j, cell.j - n_min, n_max = cell.n, cell.n - - # extend bbox if needed - 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 +function get_dggs_bbox(x_dim, y_dim, resolution, crs=DGGS.crs_geo) + # check for global coverages + if crs == DGGS.crs_geo && extrema(x_dim) == (-180, 180) && extrema(y_dim) == (-90, 90) + return [Dict(:dggs_i => 0:2*2^resolution-1, :dggs_j => 0:2^resolution-1, :dggs_n => n:n) for n in 0:4] end - return ( - Dim{:dggs_i}(i_min:i_max), - Dim{:dggs_j}(j_min:j_max), - Dim{:dggs_n}(n_min:n_max) + # get geo edges + edge_points = vcat( + map(y -> (first(x_dim), y), y_dim.val), + map(y -> (last(x_dim), y), y_dim.val), + map(x -> (x, first(y_dim)), x_dim.val), + map(x -> (x, last(y_dim)), x_dim.val) ) -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)) - j_min, j_max = dims(x, :dggs_j).val.data |> x -> (first(x), last(x)) - n_min, n_max = dims(x, :dggs_n).val.data |> x -> (first(x), last(x)) - - dggs_corners = [ - Cell(i, j, n, x.resolution) for - i in (i_min, i_max), j in (j_min, j_max), n in (n_min, n_max) - ] - geo_corners = to_geo.(dggs_corners) - - lon_min, lon_max = map(x -> x[1], geo_corners) |> x -> (minimum(x), maximum(x)) - lat_min, lat_max = map(x -> x[2], geo_corners) |> x -> (minimum(x), maximum(x)) - bbox = Extent(X=(lon_min, lon_max), Y=(lat_min, lat_max)) - return bbox + # cache trans for speed up + trans = Proj.Transformation(crs, crs_isea; ctx=Proj.proj_context_create(), always_xy=true) + cells = map(x -> to_cell(x[1], x[2], resolution, trans), edge_points) + + extents = Dict{Symbol,UnitRange{Int64}}[] + for n in 0:4 + cur_cells = filter(c -> c.n == n, cells) + length(cur_cells) == 0 && continue + i_min, i_max = map(x -> x.i, cur_cells) |> extrema + j_min, j_max = map(x -> x.j, cur_cells) |> extrema + extent = Dict( + :dggs_i => i_min:i_max, + :dggs_j => j_min:j_max, + :dggs_n => n:n + ) + push!(extents, extent) + end + return extents end -"Calculate actual geo extent" -function get_geo_bbox(geo_array::AbstractDimArray, crs::String) - # use default thread pool for lat/lon conversion - wgs84_crs_geogcs = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" - - x_min, x_max = dims(geo_array, :X) |> extrema - y_min, y_max = dims(geo_array, :Y) |> extrema +function get_geo_bbox(x_dim, y_dim, crs) + # get geo edges + edge_points = vcat( + map(y -> (first(x_dim), y), y_dim.val), + map(y -> (last(x_dim), y), y_dim.val), + map(x -> (x, first(y_dim)), x_dim.val), + map(x -> (x, last(y_dim)), x_dim.val) + ) - if crs in [wgs84_crs_geogcs, "EPSG:4326"] - ext = Extent(X=(x_min, x_max), Y=(y_min, y_max)) - return ext - else - trans = Proj.Transformation(crs, crs_geo) - lat_min, lon_min = trans(x_min, y_min) - lat_max, lon_max = trans(x_max, y_max) + # cache trans for speed up + trans = Proj.Transformation(crs, crs_geo; ctx=Proj.proj_context_create(), always_xy=true) + geo_points = map(x -> trans(x[1], x[2]), edge_points) + extent = Extent( + X=map(x -> x[1], geo_points) |> extrema, + Y=map(x -> x[2], geo_points) |> extrema + ) + return extent +end - ext = Extent(X=(lon_min, lon_max), Y=(lat_min, lat_max)) - return ext - end +function intersects(ai::UnitRange, aj::UnitRange, bi::UnitRange, bj::UnitRange) + overlap_i = first(ai) <= last(bi) && first(bi) <= last(ai) + overlap_j = first(aj) <= last(bj) && first(bj) <= last(aj) + return overlap_i && overlap_j 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), - 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), dims(geo_array, :Y)), - 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 +function get_chunks(resolution, x_dim, y_dim, crs; chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1)) + # start with dummy global array + spatial_dims = (Dim{:dggs_i}(0:2*2^resolution-1), Dim{:dggs_j}(0:2^resolution-1), Dim{:dggs_n}(0:4)) + data = Zeros(UInt8, length.(spatial_dims)) + yax_array = YAXArray(spatial_dims, data, Dict()) + yax_array = setchunks(yax_array, chunks) + + # calc extents and chunks + dggs_bboxes = get_dggs_bbox(x_dim, y_dim, resolution, crs) + + used_chunks = filter(yax_array.chunks) do c + c_i_min = c[1].start - 1 + c_i_max = c[1].stop - 1 + c_j_min = c[2].start - 1 + c_j_max = c[2].stop - 1 + c_n = c[3].start - 1 + + for r in dggs_bboxes + if c_n != r[:dggs_n].start + continue end + return intersects(c_i_min:c_i_max, c_j_min:c_j_max, r[:dggs_i], r[:dggs_j]) end + return false end - - return DGGSArray( - res.data, dims(res), refdims(res), name, metadata(geo_array), - resolution, "ISEA4D.Penta", geo_bbox - ) + return used_chunks end -"Fast iterative version only supporting mean" +DiskArrays.eachchunk(a::DGGSArray) = YAXArray(a).chunks + function to_dggs_array( geo_array::AbstractDimArray, - cells, - dggs_bbox, - geo_bbox::Extent + resolution, + crs, ; - outtype=Union{eltype(geo_array),Missing}, - outtype_counts=UInt16, - outtype_sums=Float64, - backend=:array, path=tempname() * ".dggs.zarr", name=get_name(geo_array), kwargs... ) - resolution = first(cells).resolution - - # re-grid - # mean = sum first, then divide by count - # no slow dict building and lookup needed - - counts = zeros(outtype_counts, length.(dggs_bbox)...) - - 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), dims(geo_array, :Y)), - 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] - end - counts[i_pos, j_pos, n_pos] += 1 + geo_array = cache(geo_array) + geo_ds = Dataset(layer=geo_array) + dggs_array = DGGS.init_global_dggs_dataset(geo_ds, resolution, crs, path; kwargs...).layer + dggs_array = rebuild(dggs_array, name=name) + + x_dim = dims(geo_array, :X) + y_dim = dims(geo_array, :Y) + chunks = get_chunks(resolution, x_dim, y_dim, crs) + + + # multi-threading without task migration making proj thread safe + Threads.@threads :static for chunk in chunks + # cache for other dims, e.g. time steps + trans = Proj.Transformation(DGGS.crs_isea, crs, ctx=Proj.proj_context_create(), always_xy=true) + geo_coords = map(Iterators.product(chunk...)) do (i, j, n) + cell = Cell(i - 1, j - 1, n - 1, resolution) + x, y = to_geo(cell, trans) end - end - means = sums.data ./ counts - data = if outtype <: Integer || outtype <: Union{Missing,Integer} - Array{outtype}(round.(means)) - else - Array{outtype}(means) - end - - return DGGSArray( - data, dims(sums), refdims(sums), name, metadata(geo_array), - resolution, "ISEA4D.Penta", geo_bbox - ) -end + # TODO: loop over other dims -function to_dggs_array( - geo_array::AbstractDimArray, resolution::Integer, crs::String, agg_func::Function; - x_name=:X, y_name=:Y, kwargs... -) - x_dim = filter(x -> name(x) == x_name, dims(geo_array)) - y_dim = filter(x -> name(x) == y_name, dims(geo_array)) - isempty(x_dim) && error("X dimension (e.g. longitude) not found") - isempty(y_dim) && error("Y dimension (e.g. latitude) not found") - x_dim = only(x_dim) - y_dim = only(y_dim) - - 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) + data = map(x -> geo_array[X=Near(x[1]), Y=Near(x[2])][1], geo_coords) + dggs_array[dggs_i=chunk[1], dggs_j=chunk[2], dggs_n=chunk[3]] = data end - - dggs_bbox = get_dggs_bbox(keys(cell_coords)) - geo_bbox = get_geo_bbox(geo_array, crs) - - dggs_array = to_dggs_array(geo_array, cells, cell_coords, dggs_bbox, geo_bbox, agg_func; kwargs...) return dggs_array end -function to_dggs_array(geo_array::AbstractDimArray, resolution::Integer, crs::String; x_name=:X, y_name=:Y, kwargs...) - x_dim = filter(x -> name(x) == x_name, dims(geo_array)) - y_dim = filter(x -> name(x) == y_name, dims(geo_array)) - isempty(x_dim) && error("X dimension (e.g. longitude) not found") - isempty(y_dim) && error("Y dimension (e.g. latitude) not found") - x_dim = only(x_dim) - y_dim = only(y_dim) - - properties = metadata(geo_array) - - cells = to_cell_array(x_dim, y_dim, resolution, crs) - dggs_bbox = get_dggs_bbox(cells) - geo_bbox = get_geo_bbox(geo_array, crs) - - dggs_array = to_dggs_array(geo_array, cells, dggs_bbox, geo_bbox; kwargs...) - return dggs_array -end - -function to_geo_array(dggs_array::DGGSArray, cells::AbstractDimArray; backend=:array, kwargs...) +function to_geo_array(dggs_array::DGGSArray, cells::AbstractDimArray; backend=:zarr, kwargs...) lon_dim = dims(cells, :X) lat_dim = dims(cells, :Y) @@ -316,7 +203,6 @@ function DGGSArray(array::AbstractDimArray) resolution = properties["dggs_resolution"] |> Int dggsrs = properties["dggs_dggsrs"] |> String bbox = properties["dggs_bbox"] |> x -> x isa Extent ? x : Extent(X=(x["X"][1], x["X"][2]), Y=(x["Y"][1], x["Y"][2])) - delete!(properties, "dggs_resolution") delete!(properties, "dggs_dggsrs") delete!(properties, "dggs_bbox") diff --git a/src/cells.jl b/src/cells.jl index 28d1bdc..c016bba 100644 --- a/src/cells.jl +++ b/src/cells.jl @@ -85,11 +85,11 @@ end const x_rect_min, rect_width, y_rect_min, rect_height = compute_rectangles() -function to_cell(x, y, resolution, trans::Proj.Transformation) +function to_cell(x, y, resolution, trans::Proj.Transformation, type=UInt32) x_isea, y_isea = trans(x, y) n_cell, x_offset, y_offset = ISEAToRotatedISEA()(x_isea, y_isea) i_cell, j_cell = RotatedISEAToIndices(resolution)(x_offset, y_offset) - cell = Cell(i_cell, j_cell, n_cell, resolution) + cell = Cell{type}(i_cell, j_cell, n_cell, resolution) return cell end @@ -123,11 +123,15 @@ function to_cell_array_parallel(x_dim, y_dim, resolution, crs=crs_geo; chunk_siz end -function to_geo(cell::Cell, trans::Proj.Transformation) +""" +type: precision of lat and lon, see https://xkcd.com/2170/ +""" +function to_geo(cell::Cell, trans::Proj.Transformation; type=Float32) + # TODO: add sanity checks for other crs # sanity checks are in Cell constructor # edge cases - cell.i == 2 * 2^cell.resolution - 1 && return (0.0, -90.0) + # cell.i == 2 * 2^cell.resolution - 1 && return (0.0, -90.0) # Reverse Discretization x_scaled = cell.i / 2^cell.resolution @@ -139,13 +143,15 @@ function to_geo(cell::Cell, trans::Proj.Transformation) # Reverse ISEA projection lon, lat = trans(x_isea, y_isea) + return (type(lon), type(lat)) + # Solve pole ambiguity # south pole already stable - if lat >= 90 - return (0.0, 90.0) - else - return (lon, lat) - end + # if lat >= 90 + # return (0.0, 90.0) + # else + # return (lon, lat) + # end end to_geo(i, j, n, resolution) = Cell(i, j, n, resolution) |> to_geo @@ -166,6 +172,19 @@ function to_geo_array(cells::AbstractArray{Cell{T}}, crs=crs_geo; parallel=false end +chunk_lru = LRU{Tuple{Tuple{UnitRange,UnitRange,UnitRange},Integer,String},Any}(maxsize=5) + +function to_geo_chunk(chunk, resolution, crs) + get!(chunk_lru, (chunk, resolution, crs)) do + trans = Proj.Transformation(DGGS.crs_isea, crs, ctx=Proj.proj_context_create(), always_xy=true) + res = map(Iterators.product(chunk...)) do (i, j, n) + cell = Cell(i - 1, j - 1, n - 1, resolution) + x, y = to_geo(cell, trans) + (X(Near(x)), Y(Near(y))) + end + end +end + # # Cell features # @@ -220,7 +239,7 @@ function Base.Int(cell::Cell{T}) where {T<:Integer} return res end -function Cell(cell_int::Integer, resolution::Int) +function Cell(cell_int::Integer, resolution::Integer) # needs resolution not stored in the integer index to ensure sequential id # resolution = Int((length(cell_bits) - 2) / 2) # does not work for i = 0 or j=0 n = cell_int >> (2 * resolution + 1) diff --git a/src/datasets.jl b/src/datasets.jl index fc95d36..69d32f3 100644 --- a/src/datasets.jl +++ b/src/datasets.jl @@ -32,7 +32,6 @@ function Base.getproperty(ds::DGGSDataset, s::Symbol) error("Key $(s) not found.") end - function to_dggs_dataset(geo_ds::Dataset, resolution::Integer, crs::String, agg_func::Function; metadata=Dict(), kwargs...) cells = to_cell_array(geo_ds.X, geo_ds.Y, resolution, crs) @@ -54,18 +53,73 @@ function to_dggs_dataset(geo_ds::Dataset, resolution::Integer, crs::String, agg_ return DGGSDataset(dggs_arrays...; metadata=metadata) end -"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) - - dggs_arrays = [] - 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, kwargs...) - push!(dggs_arrays, dggs_array) +"Fast iterative version using nearest neioghbor pixel" +function to_dggs_dataset( + geo_ds::Dataset, + resolution::Integer, + crs::String; + path=tempname() * ".dggs.zarr", + x_name=:X, + y_name=:Y, + metadata=Dict(), + n_parallel_chunks=Threads.nthreads(), + chunks=(dggs_i=4096, dggs_j=4096, dggs_n=1), + kwargs... +) + geo_ds = cache(geo_ds) + dggs_ds = DGGS.init_global_dggs_dataset(geo_ds, resolution, crs, path; kwargs...) + + + x_dim = geo_ds.axes[x_name] |> DD.format + y_dim = geo_ds.axes[y_name] |> DD.format + used_chunks = get_chunks(resolution, x_dim, y_dim, crs; chunks=chunks) + batches = Iterators.partition(used_chunks, n_parallel_chunks) + progress = Progress(length(used_chunks); dt=1.0) + + # limit memory usage by processing in batches + for batch in batches + # multi-threading without task migration making proj thread safe + Threads.@threads :static for chunk in batch + # pre-allocate chunk data + chunk_arrays = Dict() + for (array_name, geo_array) in pairs(geo_ds.cubes) + non_spatial_dims = filter(x -> !(name(x) in [x_name, y_name]), dims(geo_array)) + chunk_spatial_dims = ( + dims(dggs_ds, :dggs_i)[chunk[1]], + dims(dggs_ds, :dggs_j)[chunk[2]], + dims(dggs_ds, :dggs_n)[chunk[3]] + ) + chunk_dims = (chunk_spatial_dims..., non_spatial_dims...) + chunk_lengths = length.(chunk_dims) + chunk_data = Array{eltype(geo_array)}(undef, chunk_lengths...) + chunk_array = DimArray(chunk_data, chunk_dims) + chunk_arrays[array_name] = chunk_array + end + + # cache for other dims, e.g. time steps and bands + trans = Proj.Transformation(DGGS.crs_isea, crs, ctx=Proj.proj_context_create(), always_xy=true) + for (i, j, n) in Iterators.product(chunk...) + cell = Cell(i - 1, j - 1, n - 1, resolution) + x, y = to_geo(cell, trans) + + point_geo_ds = geo_ds[x=Near(x), Y=Near(y)] + for (array_name, geo_array) in pairs(point_geo_ds.cubes) + chunk_data = length(geo_array) == 1 ? geo_array[1] : geo_array + chunk_arrays[array_name][At(cell.i), At(cell.j), At(cell.n), :] = chunk_data + end + end + + for (array_name, dggs_array) in pairs(dggs_ds.data) + # write in chunks + dggs_ds[array_name][ + dggs_i=chunk[1], dggs_j=chunk[2], dggs_n=chunk[3] + ] = chunk_arrays[array_name] + end + next!(progress) + end end - return DGGSDataset(dggs_arrays...; metadata=metadata) + + return dggs_ds end function to_geo_dataset(dggs_ds::DGGSDataset, lon_dim::DD.Dimension, lat_dim::DD.Dimension; kwargs...) @@ -81,6 +135,8 @@ function to_geo_dataset(dggs_ds::DGGSDataset, lon_dim::DD.Dimension, lat_dim::DD return geo_ds end +init_global_dggs_dataset() = @error("Please load package Zarr") + Base.show(io::IO, a::DGGSArray) = print(io, "DGGSArray $(name(a)) $(a.dggsrs)@$(a.resolution) ") Base.show(io::IO, a::DGGSDataset) = print(io, "DGGSDataset $(name(a)) $(a.dggsrs)@$(a.resolution) ") diff --git a/test/runtests.jl b/test/runtests.jl index c49aafc..246e3e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -56,8 +56,8 @@ dggs_ds = DGGSDataset(dggs_array, dggs_array2) authalic_haversine = Haversine(6371007.18091875) dists = colwise(authalic_haversine, geo_points, geo_points2) - # 99% of points must be < 10m after re-projection - @test sum(dists .< 10) / length(dists) >= 0.99 + # 97% of points must be < 10m after re-projection + @test sum(dists .< 10) / length(dists) >= 0.97 # cell ids must be in bounds @test all(map(x -> x.i in 0:2*2^resolution-1, cell_ids))