Add DGGS chunked resampling

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"

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

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

src/DGGS.jl

@@ -13,18 +13,18 @@ using Printf
 using Infiltrator
 using Extents
 using OrderedCollections
+using DiskArrays
 using DiskArrayTools
+using FillArrays
+using LRUCache
+using ProgressMeter
 
 include("types.jl")
 include("cells.jl")
 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