✨ Split a layer into tiles (#148)

* ✨ Split a layer into tiles
Fixes #147

* 🧑‍🏫 doc text

* 🟩 tests for tiling

* ✨ stitching

SimpleSDMLayers/Project.toml

@@ -1,7 +1,7 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <timothee.poisot@umontreal.ca>", "Gabriel Dansereau <gabriel.dansereau@umontreal.ca>"]
-version = "0.9.0"
+version = "0.9.1"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"

SimpleSDMLayers/src/SimpleSDMLayers.jl

@@ -40,6 +40,7 @@ include("operations/sliding.jl")
 include("operations/mask.jl")
 include("operations/rescale.jl")
 include("operations/mosaic.jl")
-export coarsen, slidingwindow, mask, rescale!, rescale, mosaic
+include("operations/tiling.jl")
+export coarsen, slidingwindow, mask, rescale!, rescale, mosaic, tile, tile!, stitch
 
 end # module

SimpleSDMLayers/src/operations/tiling.jl

@@ -0,0 +1,99 @@
+"""
+    tile!(tiles::Matrix{T}, layer::T, s::Tuple{Integer, Integer} = (5, 5)) where {T <: SimpleSDMLayer}
+
+Split a layer into a matrix of tiles (with size `s`), where the tiles are as
+evenly sized as possible, and have strictly matching boundaries.
+"""
+function tile!(
+    tiles::Matrix{T},
+    layer::T,
+    s::Tuple{Integer, Integer} = (5, 5),
+) where {T <: SimpleSDMLayer}
+    @assert s[1] < size(layer, 1)
+    @assert s[2] < size(layer, 2)
+    @assert size(tiles) == s
+    # Return type
+    RT = T <: SimpleSDMPredictor ? SimpleSDMPredictor : SimpleSDMResponse
+    # Boundaries from the bounding box
+    l = layer.left
+    r = layer.right
+    b = layer.bottom
+    t = layer.top
+    l_r_points = LinRange(l, r, s[2] + 1)
+    t_b_points = LinRange(t, b, s[1] + 1)
+    # Split the coordinates
+    dim_1_cutoff = floor.(Int, LinRange(1, size(layer, 1)+1, s[1] + 1))
+    dim_2_cutoff = floor.(Int, LinRange(1, size(layer, 2)+1, s[2] + 1))
+    # Add the layers to the tiles matrix
+    for i in eachindex(dim_1_cutoff)[2:end]
+        start1, stop1 = dim_1_cutoff[(i - 1):i]
+        for j in eachindex(dim_2_cutoff)[2:end]
+            start2, stop2 = dim_2_cutoff[(j - 1):j]
+            tiles[(i - 1), (j - 1)] =
+                RT(
+                    layer.grid[start1:(stop1-1), start2:(stop2-1)],
+                    l_r_points[j - 1],
+                    l_r_points[j],
+                    t_b_points[i],
+                    t_b_points[i - 1],
+                )
+        end
+    end
+    return tiles
+end
+
+"""
+    tile(layer::T, s::Tuple{Integer, Integer} = (5, 5)) where {T <: SimpleSDMLayer}
+
+Split a layer into a matrix of tiles (with size `s`), where the tiles are as
+evenly sized as possible, and have strictly matching boundaries. This function
+will allocate the return matrix.
+"""
+function tile(layer::T, s::Tuple{Integer, Integer} = (5, 5)) where {T <: SimpleSDMLayer}
+    # Create the tiles
+    tiles = Matrix{T}(undef, s...)
+    tile!(tiles, layer, s)
+    return tiles
+end
+
+"""
+    stitch(tiles::Matrix{T})
+
+Returns a layer from a series of layer tiles, as produced by `tile` or `tile!`.
+"""
+function stitch(tiles::Matrix{T}) where {T <: SimpleSDMLayer}
+    for i in axes(tiles, 1)[1:(end - 1)]
+        for j in axes(tiles, 2)[1:(end - 1)]
+            @assert tiles[i, j].left == tiles[(i + 1), j].left
+            @assert tiles[i, j].right == tiles[(i + 1), j].right
+            @assert tiles[i, j].bottom == tiles[(i + 1), j].top
+            @assert tiles[i, j].bottom == tiles[i, (j + 1)].bottom
+            @assert tiles[i, j].top == tiles[i, (j + 1)].top
+            @assert tiles[i, j].right == tiles[i, (j + 1)].left
+        end
+    end
+    # Bounding box of the tile
+    left = tiles[1].left
+    right = tiles[end].right
+    bottom = tiles[end].bottom
+    top = tiles[end].top
+    # Size of the raster
+    rows = size.(tiles[:,1], 1)
+    cols = size.(tiles[1,:], 2)
+    # Create a grid to store info of the correct type
+    RT = eltype(tiles) <: SimpleSDMPredictor ? SimpleSDMPredictor : SimpleSDMResponse
+    IT = SimpleSDMLayers._inner_type(tiles[1])
+    grid = convert(Matrix{Union{Nothing,IT}}, fill(nothing, sum(rows), sum(cols)))
+    # Fill
+    for i in axes(tiles, 1)
+        for j in axes(tiles, 2)
+            start_i = i == 1 ? 1 : sum(rows[1:(i-1)]) + 1
+            start_j = j == 1 ? 1 : sum(cols[1:(j-1)]) + 1
+            stop_i = sum(rows[1:i])
+            stop_j = sum(cols[1:j])
+            grid[start_i:stop_i,start_j:stop_j] .= tiles[i,j].grid
+        end
+    end
+    # Return
+    return RT(grid, left, right, bottom, top)
+end

SimpleSDMLayers/test/operations/tile.jl

@@ -0,0 +1,33 @@
+module SSLTestTile
+using SimpleSDMLayers
+using Test
+
+M = rand(Bool, (13, 19))
+S = SimpleSDMPredictor(M, -2rand(), 10rand(), -3rand(), 4rand())
+
+tiles = tile(S)
+@test eltype(tiles) == typeof(S)
+@test tiles[1, 1].left == S.left
+@test tiles[1, 1].top == S.top
+@test tiles[end, 1].bottom == S.bottom
+@test tiles[end, 1].left == S.left
+@test tiles[1, end].top == S.top
+@test tiles[1, end].right == S.right
+@test tiles[end, end].bottom == S.bottom
+@test tiles[end, end].right == S.right
+
+@test tiles[1, 1].right == tiles[1, 2].left
+@test tiles[1, 1].bottom == tiles[2, 1].top
+@test tiles[1, 1].left == tiles[2, 1].left
+@test tiles[1, 1].right == tiles[2, 1].right
+
+U = stitch(S)
+@test typeof(U) == typeof(S)
+@test size(U) == size(S)
+@test U.grid == S.grid
+@test U.left == S.left
+@test U.right == S.right
+@test U.bottom == S.bottom
+@test U.top == S.top
+
+end

docs/src/vignettes/layers/05_tiling.jl

@@ -0,0 +1,49 @@
+# # Splitting layers in tiles
+
+# In this vignette, we will illustrate the use of the `tile` function, to
+# rapidly split a layer into multiple layers that can be stitched back together
+# using `hcat` and `vcat`. This is useful if you want to perform operations on
+# raster that can easily be made parallel, and do not require the full raster.
+
+using SpeciesDistributionToolkit
+using CairoMakie
+
+# To illustrate the tiling, we will grab the tree cover as given in the
+# *EarthEnv* dataset, for a small spatial extent.
+
+dataprovider = RasterData(EarthEnv, LandCover)
+spatial_extent = (left = -80.00, bottom = 43.19, right = -70.94, top = 46.93)
+trees = sum([
+    SimpleSDMPredictor(dataprovider; layer = i, full = true, spatial_extent...) for
+    i in 1:4
+])
+
+# Splitting the data into tiles can be done by calling the `tile` function (or
+# `tile!`, if you want to overwrite an existing matrix of layers). This will
+# return a matrix with the same type as the layer given as its first argument.
+# The second argument (the size of the matrix) can be omitted, and will default
+# to `(5, 5)`.
+
+tiles = tile(trees, (8, 8))
+
+# This can now be plotted:
+
+tile_plot = heatmap(
+    tiles[1, 2];
+    colormap = :Greens,
+    figure = (; resolution = (800, 350)),
+    axis = (;
+        aspect = DataAspect(),
+        xlabel = "Latitude",
+        ylabel = "Longitude",
+        title = "Relative tree cover",
+    ),
+)
+Colorbar(tile_plot.figure[:, end + 1], tile_plot.plot; height = Relative(0.5))
+current_figure()
+
+# The inverse operation to `tile` is `stitch`, which (assuming you have not
+# moved tiles around!) will reconstruct a layer. For example, let's say we want
+# to double the quantity of trees, but we want to do so using `map`.
+
+more_trees  = stitch(map(x -> 2x, tiles))