🩹 SimpleSDMLayers handling of clip with coordinates on cell boundaries (

#153)

* 🩹 SimpleSDMLayers version bump

* ✅ test for clipping with natural coordinates (#152)

* 🐛 no use of stride in clip / expand keyword

* 🐛 Update the stitching code to get the right bounding box (#155)

* 🐛 Stitch returns the wrong bounding box in the end
Fixes #154

* 🟩 test the stitch patch

* 🟩 some more tests for the right coordinate

* 🟩 removed a test where the layer type changes

* ➖ drop dependencies for SimpleSDMLayers test

* 🟩 tests with the expand keyword

* GBIF throws 404 (not 410) for deleted resources (#156)

* 🤢 GBIF throws 404 (not 410) for deleted resources

* 🐛 throws an error when the occurrence is not available

* ✅ 404 error on single occurrence is handled

* 🟥 remove the subsetting test (handled by clip)

* 🐛 fix clipping with natural coordinates

* 🐛 fix approximation in displayed coordinate

* 🐛 tests when clipping exactly on the limit

* ✅ add test comparing clip with gdalwarp

* 🐛 missing end for test module

---------

Co-authored-by: Gabriel Dansereau <gabriel.dansereau@umontreal.ca>

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.1"
+version = "0.9.2"
 
 [deps]
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"

SimpleSDMLayers/src/lib/clip.jl

@@ -1,62 +1,142 @@
 """
-    clip(layer::T, p1::Point, p2::Point) where {T <: SimpleSDMLayer}
+    clip(layer::T, p1::Point, p2::Point; expand=Symbol[]) where {T <: SimpleSDMLayer}
 
 Return a raster by defining a bounding box through two points. The order of the
 points (in terms of bottom/top/left/right) is not really important, as the
 correct coordinates will be extracted.
+
+**This method is the one that other versions of `clip` ultimately end up
+calling**.
+
+This operation takes an additional keyword argument `expand`, which is a vector
+of `Symbol`s, the symbols being any combination of `:left`, `:right`, `:bottom`,
+and `:top`. This argument specifies, in the case where the coordinate to clip is
+at the limit between two cells in a raster, whether the outer neighboring
+row/column should be included. This defaults to `expand=Symbol[]`, *i.e.* the
+outer neighboring cell will *not* be included.
 """
-function clip(layer::T, p1::Point, p2::Point) where {T <: SimpleSDMLayer}
-    latextr = extrema([p1[2], p2[2]])
-    lonextr = extrema([p1[1], p2[1]])
-    pmin = _point_to_cartesian(layer, Point(minimum(lonextr), minimum(latextr)); side=:bottomleft)
-    pmax = _point_to_cartesian(layer, Point(maximum(lonextr), maximum(latextr)); side=:topright)
-    R = T <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
-    left = longitudes(layer)[last(pmin.I)]-stride(layer, 1)
-    right = longitudes(layer)[last(pmax.I)]+stride(layer, 1)
-    bottom = latitudes(layer)[first(pmin.I)]-stride(layer, 2)
-    top = latitudes(layer)[first(pmax.I)]+stride(layer, 2)
-    return R(
-        layer.grid[pmin:pmax],
-        left = isapprox(left, round(left)) ? round(left) : left,
-        right = isapprox(right, round(right)) ? round(right) : right,
-        bottom = isapprox(bottom, round(bottom)) ? round(bottom) : bottom,
-        top = isapprox(top, round(top)) ? round(top) : top
+function clip(
+    layer::T,
+    p1::Point,
+    p2::Point;
+    expand = Symbol[],
+) where {T <: SimpleSDMLayer}
+    # Find the new bounding box for the clipped layer
+    bottom, top = extrema([p1[2], p2[2]])
+    left, right = extrema([p1[1], p2[1]])
+
+    # Get the boundaries between pixels
+    lon_boundaries = LinRange(layer.left, layer.right, size(layer, 2) + 1)
+    lat_boundaries = LinRange(layer.bottom, layer.top, size(layer, 1) + 1)
+
+    # Left to cell coordinate
+    left_cell = findlast(isapprox.(lon_boundaries, left) .|| lon_boundaries .<= left)
+    if left == lon_boundaries[left_cell]
+        if :left in expand
+            left_cell = max(1, left_cell - 1)
+        end
+    end
+
+    # Right to cell coordinate
+    right_cell = findfirst(isapprox.(lon_boundaries, right) .||lon_boundaries .>= right)
+    if right == lon_boundaries[right_cell]
+        if :right in expand
+            right_cell = min(size(layer, 2) + 1, right_cell + 1)
+        end
+    end
+
+    # Bottom to cell coordinate
+    bottom_cell = findlast(isapprox.(lat_boundaries, bottom) .|| lat_boundaries .<= bottom)
+    if bottom == lat_boundaries[bottom_cell]
+        if :bottom in expand
+            bottom_cell = max(1, bottom_cell - 1)
+        end
+    end
+
+    # Top to cell coordinate
+    top_cell = findfirst(isapprox.(lat_boundaries, top) .|| lat_boundaries .>= top)
+    if top == lat_boundaries[top_cell]
+        if :top in expand
+            top_cell = min(size(layer, 1) + 1, top_cell + 1)
+        end
+    end
+
+    # New grid
+    newgrid = layer.grid[bottom_cell:(top_cell - 1), left_cell:(right_cell - 1)]
+
+    # New return type
+    RT = typeof(layer) <: SimpleSDMResponse ? SimpleSDMResponse : SimpleSDMPredictor
+
+    # Return with the correct bounding box
+    return RT(
+        newgrid;
+        left = isapprox(left, lon_boundaries[left_cell]) ? left : lon_boundaries[left_cell],
+        right = isapprox(right, lon_boundaries[right_cell]) ? right : lon_boundaries[right_cell],
+        bottom = isapprox(bottom, lat_boundaries[bottom_cell]) ? bottom : lat_boundaries[bottom_cell],
+        top = isapprox(top, lat_boundaries[top_cell]) ? top : lat_boundaries[top_cell],
     )
 end
 
 """
-    clip(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer}
+    clip(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing, kwargs...) where {T <: SimpleSDMLayer}
 
-Clips a raster by giving the (optional) limites `left`, `right`, `bottom`, and `top`.
+Clips a raster by giving the (optional) limits `left`, `right`, `bottom`, and
+`top`.
+
+This operation takes an additional keyword argument `expand`, which is a vector
+of `Symbol`s, the symbols being any combination of `:left`, `:right`, `:bottom`,
+and `:top`. This argument specifies, in the case where the coordinate to clip is
+at the limit between two cells in a raster, whether the outer neighboring
+row/column should be included. This defaults to `expand=Symbol[]`, *i.e.* the
+outer neighboring cell will *not* be included.
 """
-function clip(layer::T; left=nothing, right=nothing, top=nothing, bottom=nothing) where {T <: SimpleSDMLayer}
+function clip(
+    layer::T;
+    left = nothing,
+    right = nothing,
+    top = nothing,
+    bottom = nothing,
+    kwargs...,
+) where {T <: SimpleSDMLayer}
     p1 = Point(
         isnothing(left) ? layer.left : left,
-        isnothing(bottom) ? layer.bottom : bottom
+        isnothing(bottom) ? layer.bottom : bottom,
     )
     p2 = Point(
         isnothing(right) ? layer.right : right,
-        isnothing(top) ? layer.top : top
+        isnothing(top) ? layer.top : top,
     )
-    return clip(layer, p1, p2)
+    return clip(layer, p1, p2; kwargs...)
 end
 
 """
-    clip(origin::T1, destination::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
+    clip(origin::T1, destination::T2; kwargs...) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
 
 Clips a layer by another layer, *i.e.* subsets the first layer so that it has
-the dimensions of the second layer. This operation applies a very small
-displacement to the limits (`5eps()`) to ensure that if the coordinate to cut at
-falls exactly on a cell boundary, the correct cell will be used in the return
-layer.
+the dimensions of the second layer.
+
+This operation takes an additional keyword argument `expand`, which is a vector
+of `Symbol`s, the symbols being any combination of `:left`, `:right`, `:bottom`,
+and `:top`. This argument specifies, in the case where the coordinate to clip is
+at the limit between two cells in a raster, whether the outer neighboring
+row/column should be included. This defaults to `expand=Symbol[]`, *i.e.* the
+outer neighboring cell will *not* be included.
 """
-function clip(origin::T1, destination::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
-    _d = 5eps()
+function clip(
+    origin::T1,
+    destination::T2;
+    kwargs...,
+) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
     err = false
     destination.right > origin.right && (err = true)
     destination.left < origin.left && (err = true)
     destination.bottom < origin.bottom && (err = true)
     destination.top > origin.top && (err = true)
     err && throw(ArgumentError("The two layers are not compatible"))
-    return clip(origin, Point(destination.left+_d, destination.bottom+_d), Point(destination.right-_d, destination.top-_d))
+    return clip(
+        origin,
+        Point(destination.left, destination.top),
+        Point(destination.right, destination.bottom);
+        kwargs...,
+    )
 end

SimpleSDMLayers/src/operations/tiling.jl

@@ -73,10 +73,10 @@ function stitch(tiles::Matrix{T}) where {T <: SimpleSDMLayer}
         end
     end
     # Bounding box of the tile
-    left = tiles[1].left
-    right = tiles[end].right
-    bottom = tiles[end].bottom
-    top = tiles[end].top
+    left = minimum([t.left for t in tiles])
+    right = maximum([t.right for t in tiles])
+    bottom = minimum([t.bottom for t in tiles])
+    top = maximum([t.top for t in tiles])
     # Size of the raster
     rows = size.(tiles[:,1], 1)
     cols = size.(tiles[1,:], 2)

SimpleSDMLayers/test/Project.toml

@@ -1,9 +1,3 @@
 [deps]
-GBIF = "ee291a33-5a6c-5552-a3c8-0f29a1181037"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[compat]
-GBIF = "0.2, 0.3"

SimpleSDMLayers/test/operations/clip.jl

@@ -5,19 +5,66 @@ using Test
 M = rand(Bool, (10, 10))
 S = SimpleSDMPredictor(M, 0.0, 1.0, 0.0, 1.0)
 
-cl1 = clip(S; left=0.2, right=0.6, bottom=0.5, top=1.0)
+cl1 = clip(S; left = 0.2, right = 0.6, bottom = 0.5, top = 1.0)
 @test typeof(cl1) == typeof(S)
 @test cl1.top ≈ 1.0
 @test cl1.bottom ≈ 0.5
 @test cl1.right ≈ 0.6
 @test cl1.left ≈ 0.2
-@test clip(S; left=0.19).left <= 0.2
+@test clip(S; left = 0.19).left <= 0.2
 
-cl2 = clip(S; left=0.2, bottom=0.5)
+cl2 = clip(S; left = 0.2, bottom = 0.5)
 @test typeof(cl2) == typeof(S)
 @test cl2.top ≈ 1.0
 @test cl2.bottom ≈ 0.5
 @test cl2.right ≈ 1.0
 @test cl2.left ≈ 0.2
 
+# Test the expand keyword to include points that fall on the limit between two cells
+M2 = rand(Bool, (1080, 2160))
+S2 = SimpleSDMPredictor(M2, -180.0, 180.0, -90.0, 90.0)
+
+# Test that expand works, including at the edges of the layer
+@test clip(S2; left = 0.0).left == 0.0
+@test clip(S2; left = 0.0, expand = [:left]).left < 0.0
+@test clip(S2; left = -180.0).left == -180.0
+@test clip(S2; left = -180.0, expand = [:left]).left == -180.0
+
+# Test clipping with natural coordinates exactly on the limit
+exact_lats = (S2.bottom+1.0):1.0:(S2.top-1.0)
+exact_lons = (S2.left+1.0):1.0:(S2.right-1.0)
+@test all([isequal(l, clip(S2; top=l).top) for l in exact_lats])
+@test all([isequal(l, clip(S2; bottom=l).bottom) for l in exact_lats])
+@test all([isequal(l, clip(S2; left=l).left) for l in exact_lons])
+@test all([isequal(l, clip(S2; right=l).right) for l in exact_lons])
+
+# Test clipping with non-natural coordinates
+@test all([isequal(l, clip(S; top=l).top) for l in 0.1:0.1:0.9])
+@test all([isequal(l, clip(S; bottom=l).bottom) for l in 0.1:0.1:0.9])
+@test all([isequal(l, clip(S; left=l).left) for l in 0.1:0.1:0.9])
+@test all([isequal(l, clip(S; right=l).right) for l in 0.1:0.1:0.9])
+
+# Test that clip is consistent with GDAL.jl
+# This test requires to read and write the data, which is handled by
+# SpeciesDistributionToolkit, not SimpleSDMLayers. We'll leave the code
+# commented out here and test only with the values that should be returned when
+# using GDAL.
+#=
+using SpeciesDistributionToolkit
+using GDAL
+F2 = "$(tempname()).tif"
+F3 = "$(tempname()).tif"
+SpeciesDistributionToolkit.save(F2, convert(Float32, S2));
+run(`$(GDAL.gdalwarp_path()) $F2 $F3 -te -180.0 -90.0 180.0 8.0 -overwrite`)
+S3 = SpeciesDistributionToolkit._read_geotiff(F3, SimpleSDMPredictor);
+cl3 = clip(convert(Float32, S2); top=8.0);
+@test cl3 == S3
+@test size(cl3) == size(S3)
+@test boundingbox(cl3) == boundingbox(S3)
+=#
+
+cl3 = clip(convert(Float32, S2); top=8.0);
+@test size(cl3) == (588, 2160)
+@test boundingbox(cl3) == (left = -180.0, right = 180.0, bottom = -90.0, top = 8.0)
+
 end

SimpleSDMLayers/test/runtests.jl

@@ -16,7 +16,6 @@ tests = [
     "rescale" => "operations/rescale.jl",
     "mosaic" => "operations/mosaic.jl",
     "coarsen" => "operations/coarsen.jl",
-    "subsetting" => "operations/subsetting.jl",
     "generated" => "generated.jl",
 ]
 

test/edgecases/01_stitch_wrong_bb.jl

@@ -0,0 +1,43 @@
+module TestThatStitchIsNotInAMood
+
+# The stitching sometimes gave the wrong bounding box after using `map`, so this
+# checks this doesn't happen anymore
+
+using SpeciesDistributionToolkit
+using Test
+
+# Get some data
+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
+])
+
+for tilesize in [(3, 4), (8, 8), (4, 3), (9, 5), (2, 2), (12, 1), (3, 6)]
+
+    # make the tiles
+    tiles = tile(trees, tilesize)
+    @test eltype(tiles) == typeof(trees)
+
+    # stitch the tiles as they are
+    tiled_trees = stitch(tiles)
+    @test typeof(tiled_trees) == typeof(trees)
+    @test size(tiled_trees) == size(trees)
+    @test tiled_trees.grid == trees.grid
+    @test tiled_trees.left == trees.left
+    @test tiled_trees.right == trees.right
+    @test tiled_trees.bottom == trees.bottom
+    @test tiled_trees.top == trees.top
+
+    # stitch the tiles after a transform
+    more_trees  = stitch(map(x -> 2x, tiles))
+    @test size(more_trees) == size(trees)
+    @test more_trees.left == trees.left
+    @test more_trees.right == trees.right
+    @test more_trees.bottom == trees.bottom
+    @test more_trees.top == trees.top
+
+end
+
+end

test/runtests.jl

@@ -5,6 +5,7 @@ global anyerrors = false
 
 tests = [
     "read/write layers" => "01_integration_read.jl",
+    "EDGE: stitch bounding box" => "edgecases/01_stitch_wrong_bb.jl",
 ]
 
 for test in tests