Various minor improvements

docs/src/man/operations.md

@@ -5,6 +5,7 @@ coarsen
 slidingwindow
 latitudes
 longitudes
+boundingbox
 clip
 mask
 rescale!

docs/src/man/overloads.md

@@ -30,6 +30,9 @@ Base.min
 -
 *
 /
+==
+isequal
+hash
 ```
 
 ## From `Broadcast`

src/SimpleSDMLayers.jl

@@ -15,7 +15,7 @@ include(joinpath("lib", "overloads.jl"))
 include(joinpath("lib", "generated.jl"))
 
 include(joinpath("lib", "basics.jl"))
-export latitudes, longitudes
+export latitudes, longitudes, boundingbox
 
 include(joinpath("lib", "iteration.jl"))
 
@@ -71,11 +71,11 @@ export clip
 
 function __init__()
     @require GBIF="ee291a33-5a6c-5552-a3c8-0f29a1181037" begin
-        @info "GBIF integration loaded"
+        @info "Loading GBIF support for SimpleSDMLayers.jl"
         include("integrations/GBIF.jl")
     end
     @require DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
-        @info "DataFrames integration loaded"
+        @info "Loading DataFrames support for SimpleSDMLayers.jl"
         include("integrations/DataFrames.jl")
     end
 

src/datasets/geotiff.jl

@@ -102,7 +102,7 @@ Write a single `layer` to a `file`, where the `nodata` field is set to an
 arbitrary value.
 """
 function geotiff(file::AbstractString, layer::SimpleSDMPredictor{T}; nodata::T=convert(T, -9999)) where {T <: Number}
-    array_t = _prepare_layer_for_burnin(layer)
+    array_t = _prepare_layer_for_burnin(layer, nodata)
     width, height = size(array_t)
 
     # Geotransform
@@ -138,9 +138,9 @@ function geotiff(file::AbstractString, layer::SimpleSDMPredictor{T}; nodata::T=c
     return file
 end
 
-function _prepare_layer_for_burnin(layer::T) where {T <: SimpleSDMLayer}
+function _prepare_layer_for_burnin(layer::SimpleSDMPredictor{T}, nodata::T) where {T <: Number}
     @assert eltype(layer) <: Number
-    array = replace(layer.grid, nothing => NaN)
+    array = replace(layer.grid, nothing => nodata)
     array = convert(Matrix{eltype(layer)}, array)
     dtype = eltype(array)
     array_t = reverse(permutedims(array, [2, 1]); dims=2)
@@ -156,7 +156,7 @@ Stores a series of `layers` in a `file`, where every layer in a band. See
 function geotiff(file::AbstractString, layers::Vector{SimpleSDMPredictor{T}}; nodata::T=convert(T, -9999)) where {T <: Number}
     bands = 1:length(layers)
     _layers_are_compatible(layers)
-    width, height = size(_prepare_layer_for_burnin(layers[1]))
+    width, height = size(_prepare_layer_for_burnin(layers[1], nodata))
 
     # Geotransform
     gt = zeros(Float64, 6)
@@ -179,7 +179,7 @@ function geotiff(file::AbstractString, layers::Vector{SimpleSDMPredictor{T}}; no
             band = ArchGDAL.getband(dataset, i)
 
             # Write data to band
-            ArchGDAL.write!(band, _prepare_layer_for_burnin(layers[i]))
+            ArchGDAL.write!(band, _prepare_layer_for_burnin(layers[i], nodata))
 
             # Write nodata and projection info
             ArchGDAL.setnodatavalue!(band, nodata)

src/lib/basics.jl

@@ -14,6 +14,14 @@ This returns the longitudes at the center of each cell in the grid.
 """
 longitudes(layer::T) where {T <: SimpleSDMLayer} = range(layer.left+stride(layer, 1), layer.right-stride(layer, 1); length=size(layer,2))
 
+"""
+    boundingbox(layer::T) where {T <: SimpleSDMLayer}
+
+Returns the bounding coordinates of a layer as `NamedTuple`.
+"""
+
+boundingbox(layer::T) where {T <: SimpleSDMLayer} =  (left=layer.left, right=layer.right, bottom=layer.bottom, top=layer.top)
+
 """
     _layers_are_compatible(l1::X, l2::Y) where {X <: SimpleSDMLayer, Y <: SimpleSDMLayer}
 

src/lib/overloads.jl

@@ -14,6 +14,9 @@ import Base.Broadcast: broadcast
 import Base: hcat
 import Base: vcat
 import Base: show
+import Base: ==
+import Base: isequal
+import Base: hash
 
 """
     Base.show(io::IO, ::MIME"text/plain", layer::T) where {T <: SimpleSDMLayer}
@@ -263,11 +266,20 @@ end
     Base.getindex(layer1::T1, layer2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
 
 Extract a layer based on a second layer. Note that the two layers must be
-*compatible*, which is to say they must have the same bounding box and grid
-size.
+*compatible*, which is to say they must have the same stride and the bounding
+coordinates of layer2 must be contained in layer1.
 """
 function Base.getindex(layer1::T1, layer2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
-    SimpleSDMLayers._layers_are_compatible(layer1, layer2)
+    iscompat = all(
+        [
+            layer2.left >= layer1.left,
+            layer2.right <= layer1.right,
+            layer2.bottom >= layer1.bottom,
+            layer2.top <= layer1.top,
+        ]
+    )
+    iscompat || throw(ArgumentError("layer2 has bounding coordinates that are not contained in layer1"))
+    stride(layer1) == stride(layer2) || throw(ArgumentError("The layers have different strides"))
     return layer1[left=layer2.left, right=layer2.right, bottom=layer2.bottom, top=layer2.top]
 end
 
@@ -399,11 +411,14 @@ function Base.hcat(l1::T, l2::T) where {T <: SimpleSDMLayer}
 end
 
 """
-    Base.replace!(layer::T, old_new::Pair...) where {T <: SimpleSDMResponse}
+    Base.replace!(layer::T, old_new::Pair...) where {T <: SimpleSDMLayer}
 
-Replaces the elements of `layer` according to a series of pairs. In place.
+Replaces the elements of `layer` according to a series of pairs. In place. Only
+possible for `SimpleSDMResponse` elements (which are mutable) and will throw an
+error if called on a `SimpleSDMPredictor` element (which is not mutable).
 """
-function Base.replace!(layer::T, old_new::Pair...) where {T <: SimpleSDMResponse}
+function Base.replace!(layer::T, old_new::Pair...) where {T <: SimpleSDMLayer}
+    layer isa SimpleSDMResponse || throw(ArgumentError("`SimpleSDMPredictor` elements are immutable. Convert to a `SimpleSDMResponse` first or call `replace!` directly on the grid element."))
     replace!(layer.grid, old_new...)
     return layer
 end
@@ -420,7 +435,7 @@ function Base.replace(layer::T, old_new::Pair...) where {T <: SimpleSDMResponse}
 end
 
 """
-    Base.replace!(layer::T, old_new::Pair...) where {T <: SimpleSDMResponse}
+    Base.replace(layer::T, old_new::Pair...) where {T <: SimpleSDMPredictor}
 
 Replaces the elements of `layer` according to a series of pairs. Copies the
 layer as a response before.
@@ -438,4 +453,48 @@ Returns the quantiles of `layer` at `p`, using `Statistics.quantile`.
 """
 function Statistics.quantile(layer::T, p) where {T <: SimpleSDMLayer}
     return quantile(collect(layer), p)
-end
+end
+
+"""
+    ==(layer1::SimpleSDMLayer, layer2::SimpleSDMLayer)
+
+Tests whether two `SimpleSDMLayer` elements are equal. The layers are equal if 
+all their fields (`grid`, `left`, `right`, `bottom`, `top`) are equal, as 
+verified with `==` (e.g., `layer1.grid == layer2.grid`).
+"""
+
+function Base.:(==)(layer1::SimpleSDMLayer, layer2::SimpleSDMLayer)
+    return all(
+        [
+            layer1.grid == layer2.grid,
+            layer1.left == layer2.left,
+            layer1.right == layer2.right,
+            layer1.bottom == layer2.bottom,
+            layer1.top == layer2.top,
+        ]
+    )
+end
+
+function Base.hash(layer::SimpleSDMLayer, h::UInt)
+    return hash((layer.grid, layer.left, layer.right, layer.bottom, layer.top), h)
+end
+
+"""
+    isequal(layer1::SimpleSDMLayer, layer2::SimpleSDMLayer)
+
+Tests whether two `SimpleSDMLayer` elements are equal. The layers are equal if 
+all their fields (`grid`, `left`, `right`, `bottom`, `top`) are equal, as 
+verified with `isequal` (e.g., `isequal(layer1.grid, layer2.grid)`).
+"""
+
+function Base.isequal(layer1::SimpleSDMLayer, layer2::SimpleSDMLayer)
+    return all(
+        [
+            isequal(layer1.grid, layer2.grid),
+            isequal(layer1.left, layer2.left),
+            isequal(layer1.right, layer2.right),
+            isequal(layer1.bottom, layer2.bottom),
+            isequal(layer1.top, layer2.top),
+        ]
+    )
+end

test/basics.jl

@@ -14,4 +14,10 @@ S = SimpleSDMPredictor(M, 0.2, 1.8, -1.0, 2.0)
 @test longitudes(S) == range(S.left+stride(S,1), S.right-stride(S,1); length=size(S,2))
 @test latitudes(S) == range(S.bottom+stride(S,2), S.top-stride(S,2); length=size(S,1))
 
+bbox = boundingbox(S)
+@test bbox.left == S.left
+@test bbox.right == S.right
+@test bbox.bottom == S.bottom
+@test bbox.top == S.top
+
 end

test/dataread.jl

@@ -3,11 +3,21 @@ using SimpleSDMLayers
 using Test
 
 l = SimpleSDMPredictor(WorldClim, BioClim, 1)
-f = tempname()
-geotiff(f, l)
-mp = geotiff(SimpleSDMResponse, f)
 
-@test typeof(mp) <: SimpleSDMResponse
-@test size(mp) == size(l)
+f1 = tempname()
+geotiff(f1, l)
+mp1 = geotiff(SimpleSDMResponse, f1)
+
+@test typeof(mp1) <: SimpleSDMResponse
+@test size(mp1) == size(l)
+@test mp1 == l
+
+f2 = tempname()
+geotiff(f2, l; nodata=-3.4f38)
+mp2 = geotiff(SimpleSDMPredictor, f2)
+
+@test typeof(mp2) <: SimpleSDMPredictor
+@test size(mp2) == size(l)
+@test mp2 == l
 
 end

test/overloads.jl

@@ -101,4 +101,47 @@ s2 = replace(s1, 1 => 2, 3 => 2, 9 => nothing)
 @test s2.grid[3,1] == 2
 @test s2.grid[1,3] == 7
 
+# ==, isequal, hash
+l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2; left = 0.0, right = 10.0, bottom = 0.0, top = 10.0)
+l3 = copy(l1)
+l4 = similar(l1)
+replace!(l4, nothing => NaN)
+l5 = SimpleSDMPredictor(replace(l1.grid, nothing => missing), l1)
+
+@test l1 == l1
+@test l1 === l1
+@test l2 != l1
+@test l3 == l1
+@test l3 !== l1
+
+@test l4 != l1
+@test l4 != l4
+@test !isequal(l4, l1)
+@test isequal(l4, l4)
+
+@test ismissing(l5 == l1)
+@test ismissing(l5 == l5)
+@test !isequal(l5, l1)
+@test isequal(l5, l5)
+
+@test hash(l1) == hash(l1)
+@test hash(l2) != hash(l1)
+@test hash(l3) == hash(l1)
+@test hash(l4) != hash(l1)
+@test hash(l4) == hash(l4)
+@test hash(l5) != hash(l1)
+@test hash(l5) == hash(l5)
+
+# getindex(layer1, layer2)
+l1, l2 = SimpleSDMPredictor(WorldClim, BioClim, 1:2; left = 0.0, right = 10.0, bottom = 0.0, top = 10.0)
+l3 = SimpleSDMPredictor(WorldClim, BioClim, 1; left = 5.0, right = 10.0, bottom = 5.0, top = 10.0)
+@test stride(l1) == stride(l3)
+
+l4 = l1[l2]
+@test l4 == l1
+
+l5 = l1[l3]
+@test l5 == l3
+@test_throws ArgumentError l3[l1]
+
 end

test/subsetting.jl

@@ -10,7 +10,7 @@ l1 = temp[coords]
 l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
 tempfile = tempname()
 geotiff(tempfile, l2)
-l3 = replace(geotiff(SimpleSDMPredictor, tempfile), NaN => nothing)
+l3 = geotiff(SimpleSDMPredictor, tempfile)
 
 @test size(l1) == size(l2)
 @test size(l1) == size(l3)