diff --git a/SimpleSDMLayers/Project.toml b/SimpleSDMLayers/Project.toml index f868dad1a..0cff01476 100644 --- a/SimpleSDMLayers/Project.toml +++ b/SimpleSDMLayers/Project.toml @@ -1,7 +1,7 @@ name = "SimpleSDMLayers" uuid = "2c645270-77db-11e9-22c3-0f302a89c64c" authors = ["Timothée Poisot ", "Gabriel Dansereau "] -version = "0.9.1" +version = "0.9.2" [deps] Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" diff --git a/SimpleSDMLayers/src/lib/clip.jl b/SimpleSDMLayers/src/lib/clip.jl index 79f14656e..9785209ee 100644 --- a/SimpleSDMLayers/src/lib/clip.jl +++ b/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 diff --git a/SimpleSDMLayers/src/operations/tiling.jl b/SimpleSDMLayers/src/operations/tiling.jl index 8ccbd8b9c..bc96915a7 100644 --- a/SimpleSDMLayers/src/operations/tiling.jl +++ b/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) diff --git a/SimpleSDMLayers/test/Project.toml b/SimpleSDMLayers/test/Project.toml index 265081820..644bf61f5 100644 --- a/SimpleSDMLayers/test/Project.toml +++ b/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" diff --git a/SimpleSDMLayers/test/operations/clip.jl b/SimpleSDMLayers/test/operations/clip.jl index 44c8e8e0d..c8b242d20 100644 --- a/SimpleSDMLayers/test/operations/clip.jl +++ b/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 \ No newline at end of file diff --git a/SimpleSDMLayers/test/operations/subsetting.jl b/SimpleSDMLayers/test/operations/subsetting.jl deleted file mode 100644 index a7307737c..000000000 --- a/SimpleSDMLayers/test/operations/subsetting.jl +++ /dev/null @@ -1,22 +0,0 @@ -module SSLTestSubsetting -using SimpleSDMLayers -using Test - -temp = SimpleSDMPredictor(rand(1080, 2160)) -@test size(temp) == (1080, 2160) - -coords = (left=-145.0, right=-50.0, bottom=20.0, top=75.0) -l1 = clip(temp; coords...) - -@test l1.left == coords.left -@test l1.right == coords.right -@test l1.bottom == coords.bottom -@test l1.top == coords.top - -@test stride(l1) == stride(temp) -@test longitudes(l1)[1] ≈ coords.left atol = 0.1 -@test longitudes(l1)[end] ≈ coords.right atol = 0.1 -@test latitudes(l1)[1] ≈ coords.bottom atol = 0.1 -@test latitudes(l1)[end] ≈ coords.top atol = 0.1 - -end diff --git a/SimpleSDMLayers/test/runtests.jl b/SimpleSDMLayers/test/runtests.jl index 8020fa3bb..132b2a00a 100644 --- a/SimpleSDMLayers/test/runtests.jl +++ b/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", ] diff --git a/test/edgecases/01_stitch_wrong_bb.jl b/test/edgecases/01_stitch_wrong_bb.jl new file mode 100644 index 000000000..aba5b9ae0 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index cb006db78..bb13cade4 100644 --- a/test/runtests.jl +++ b/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