Fix Upscale and Aggregate (#60)

* Fix Upscale

* Add '_fitdims' helper

* Fix Aggregate

* Add more tests

* Update compat

* Apply suggestions

Project.toml

@@ -39,7 +39,7 @@ GeoStatsModels = "0.6"
 GeoStatsProcesses = "0.8"
 GeoTables = "1.21"
 LinearAlgebra = "1.9"
-Meshes = "0.50 - 0.52"
+Meshes = "0.52.2"
 OhMyThreads = "0.5, 0.6, 0.7"
 PrecompileTools = "1.2"
 Random = "1.9"

src/aggregate.jl

@@ -61,7 +61,7 @@ function apply(transform::Aggregate, geotable::AbstractGeoTable)
 end
 
 function _aggregate(sdom, tdom, cols, vars, aggfun)
-  if sdom isa Grid && tdom isa Grid && extrema(sdom) == extrema(tdom)
+  if sdom isa Grid && tdom isa Grid && extrema(sdom) == extrema(tdom) && all(iszero, size(sdom) .% size(tdom))
     # we have two grids overlaid, and can rely on
     # tiled iteration for efficient aggregation
     _gridagg(sdom, tdom, cols, vars, aggfun)
@@ -73,7 +73,7 @@ end
 
 function _gridagg(sdom, tdom, cols, vars, aggfun)
   # determine tile size for tiled iteration
-  tilesize = ceil.(Int, size(sdom) ./ size(tdom))
+  tilesize = size(sdom) .÷ size(tdom)
   if any(<(1), tilesize)
     throw(ArgumentError("cannot aggregate a coarse grid over a fine grid"))
   end

src/upscale.jl

@@ -29,8 +29,33 @@ Upscale(factors::Int...) = Upscale(factors)
 isrevertible(::Type{<:Upscale}) = false
 
 function apply(transform::Upscale, geotable::AbstractGeoTable)
-  grid = domain(geotable)
-  tgrid = coarsen(grid, RegularCoarsening(transform.factors))
-  newgeotable = geotable |> Aggregate(tgrid)
+  gtb = _adjustunits(geotable)
+  tab = values(gtb)
+  grid = domain(gtb)
+  cols = Tables.columns(tab)
+  vars = Tables.columnnames(cols)
+
+  # upscale the grid
+  factors = _fitdims(transform.factors, paramdim(grid))
+  tgrid = coarsen(grid, RegularCoarsening(factors))
+
+  # aggregate columns
+  pairs = map(vars) do var
+    svals = Tables.getcolumn(cols, var)
+    aggfun = _defaultagg(svals)
+    array = reshape(svals, size(grid))
+    titer = TileIterator(axes(array), factors)
+    tvals = tmap(titer) do sinds
+      aggfun(array[sinds...])
+    end |> vec
+    var => tvals
+  end
+
+  # construct new table
+  newtab = (; pairs...) |> Tables.materializer(tab)
+
+  # new spatial data
+  newgeotable = georef(newtab, tgrid)
+
   newgeotable, nothing
 end

src/utils.jl

@@ -2,6 +2,9 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
+# fit tuple `dims` to a given length `D` by repeating the last dimension.
+_fitdims(dims::Dims{N}, D) where {N} = ntuple(i -> i ≤ N ? dims[i] : last(dims), D)
+
 #-------------
 # AGGREGATION
 #-------------

test/upscale.jl

@@ -6,6 +6,8 @@
   gtb = georef((a=rand(Float64, 400), b=rand(Int, 400)), grid)
   ngtb = gtb |> Upscale(2, 2)
   @test domain(ngtb) == tgrid
+  @test length(ngtb.a) == nelements(domain(ngtb))
+  @test length(ngtb.b) == nelements(domain(ngtb))
   @test ngtb[(1, 1), :a] == mean(gtb[(1:2, 1:2), :a])
   @test ngtb[(1, 10), :a] == mean(gtb[(1:2, 19:20), :a])
   @test ngtb[(10, 1), :a] == mean(gtb[(19:20, 1:2), :a])
@@ -20,6 +22,8 @@
   gtb = georef((a=rand(Float64, 400), b=rand(Int, 400)), rgrid)
   ngtb = gtb |> Upscale(2, 2)
   @test domain(ngtb) == trgrid
+  @test length(ngtb.a) == nelements(domain(ngtb))
+  @test length(ngtb.b) == nelements(domain(ngtb))
   @test ngtb[(1, 1), :a] == mean(gtb[(1:2, 1:2), :a])
   @test ngtb[(1, 10), :a] == mean(gtb[(1:2, 19:20), :a])
   @test ngtb[(10, 1), :a] == mean(gtb[(19:20, 1:2), :a])
@@ -34,6 +38,8 @@
   gtb = georef((a=rand(Float64, 400), b=rand(Int, 400)), sgrid)
   ngtb = gtb |> Upscale(2, 2)
   @test domain(ngtb) == tsgrid
+  @test length(ngtb.a) == nelements(domain(ngtb))
+  @test length(ngtb.b) == nelements(domain(ngtb))
   @test ngtb[(1, 1), :a] == mean(gtb[(1:2, 1:2), :a])
   @test ngtb[(1, 10), :a] == mean(gtb[(1:2, 19:20), :a])
   @test ngtb[(10, 1), :a] == mean(gtb[(19:20, 1:2), :a])
@@ -47,6 +53,8 @@
   gtb = georef((a=rand(Float64, 400), b=rand(Int, 400)), grid)
   ngtb = gtb |> Upscale(2, 4)
   @test domain(ngtb) == tgrid
+  @test length(ngtb.a) == nelements(domain(ngtb))
+  @test length(ngtb.b) == nelements(domain(ngtb))
   @test ngtb[(1, 1), :a] == mean(gtb[(1:2, 1:4), :a])
   @test ngtb[(1, 5), :a] == mean(gtb[(1:2, 17:20), :a])
   @test ngtb[(10, 1), :a] == mean(gtb[(19:20, 1:4), :a])
@@ -55,4 +63,36 @@
   @test ngtb[(1, 5), :b] == first(gtb[(1:2, 17:20), :b])
   @test ngtb[(10, 1), :b] == first(gtb[(19:20, 1:4), :b])
   @test ngtb[(10, 5), :b] == first(gtb[(19:20, 17:20), :b])
+
+  # non-multiple dimensions
+  grid = CartesianGrid((0.0, 0.0), (13.0, 17.0), dims=(13, 17))
+  gtb = georef((a=rand(Float64, nelements(grid)), b=rand(Int, nelements(grid))), grid)
+  ngtb = gtb |> Upscale(5, 3)
+  @test size(domain(ngtb)) == (3, 6)
+  @test length(ngtb.a) == nelements(domain(ngtb))
+  @test length(ngtb.b) == nelements(domain(ngtb))
+  @test ngtb[(1, 1), :a] == mean(gtb[(1:5, 1:3), :a])
+  @test ngtb[(1, 6), :a] == mean(gtb[(1:5, 16:17), :a])
+  @test ngtb[(3, 1), :a] == mean(gtb[(11:13, 1:3), :a])
+  @test ngtb[(3, 6), :a] == mean(gtb[(11:13, 16:17), :a])
+  @test ngtb[(1, 1), :b] == first(gtb[(1:5, 1:3), :b])
+  @test ngtb[(1, 6), :b] == first(gtb[(1:5, 16:17), :b])
+  @test ngtb[(3, 1), :b] == first(gtb[(11:13, 1:3), :b])
+  @test ngtb[(3, 6), :b] == first(gtb[(11:13, 16:17), :b])
+
+  # large grid
+  grid = CartesianGrid((0.0, 0.0), (16200.0, 8100.0), dims=(16200, 8100))
+  gtb = georef((a=rand(Float64, nelements(grid)), b=rand(Int, nelements(grid))), grid)
+  ngtb = gtb |> Upscale(80, 40)
+  @test size(domain(ngtb)) == (203, 203)
+  @test length(ngtb.a) == nelements(domain(ngtb))
+  @test length(ngtb.b) == nelements(domain(ngtb))
+  @test ngtb[(1, 1), :a] ≈ mean(gtb[(1:80, 1:40), :a])
+  @test ngtb[(1, 203), :a] ≈ mean(gtb[(1:80, 8081:8100), :a])
+  @test ngtb[(203, 1), :a] ≈ mean(gtb[(16161:16200, 1:40), :a])
+  @test ngtb[(203, 203), :a] ≈ mean(gtb[(16161:16200, 8081:8100), :a])
+  @test ngtb[(1, 1), :b] == first(gtb[(1:80, 1:40), :b])
+  @test ngtb[(1, 203), :b] == first(gtb[(1:80, 8081:8100), :b])
+  @test ngtb[(203, 1), :b] == first(gtb[(16161:16200, 1:40), :b])
+  @test ngtb[(203, 203), :b] == first(gtb[(16161:16200, 8081:8100), :b])
 end