mask and mask! now remove missing values in the source raster (#507)

* Update skipmissing.jl Fixed bug when calling skipmissing on a Raster of floats for which missingval = nothing. * Wrote test cases for skipmissing * Additional skipmissing test cases * Implemented additional _missing methods. * Implemented additional _missing methods. * mask and mask! now remove missing values in the source raster * skipmissing does not drop nothing * Update src/skipmissing.jl * Update src/skipmissing.jl * missingval must be convertable in both mask and mask * mask now allows missingval to be either a convertable value or missing * Test cases for mask * Test cases for mask! * Test cases for mask! * mask! test cases * mask! test cases * mask! test cases * mask! test cases --------- Co-authored-by: Rafael Schouten <[email protected]>
rafaqz · Oct 8, 2023 · 96a56ed · 96a56ed
1 parent ff90b8f
commit 96a56ed
Show file tree

Hide file tree

Showing 2 changed files with 60 additions and 5 deletions.
diff --git a/src/methods/mask.jl b/src/methods/mask.jl
@@ -76,6 +76,7 @@ end
 function _mask(A::AbstractRaster, with::AbstractRaster;
     filename=nothing, suffix=nothing, missingval=_missingval_or_missing(A), kw...
 )
+    missingval = ismissing(missingval) ? missing : convert(eltype(A), missingval)
     A1 = create(filename, A; suffix, missingval)
     open(A1; write=true) do a
         # The values array will be be written to A1 in `mask!`
@@ -168,20 +169,21 @@ function _mask!(st::RasterStack, with::AbstractRaster; kw...)
     map(A -> mask!(A; with, kw...), st)
     return st
 end
-function _mask!(A::AbstractRaster, with::AbstractRaster;
-    missingval=missingval(A), values=A
-)
+
+function _mask!(A::AbstractRaster, with::AbstractRaster; missingval=missingval(A), values=A)
     missingval isa Nothing && _nomissingerror()
     missingval = convert(eltype(A), missingval)
 
     broadcast_dims!(A, values, with) do x, w
-        if isequal(w, Rasters.missingval(with))
+        if ismissing(w) || ismissing(x)
+            return missingval
+        elseif isequal(w, Rasters.missingval(with)) || isequal(x, Rasters.missingval(values))
             missingval
         else
             convert(eltype(A), x)
         end
     end
-    return A
+    return rebuild(A, missingval=missingval)
 end
 
 _nomissingerror() = throw(ArgumentError("Array has no `missingval`. Pass a `missingval` keyword compatible with the type, or use `rebuild(A; missingval=somemissingval)` to set it."))

diff --git a/test/methods.jl b/test/methods.jl
@@ -142,6 +142,59 @@ end
     end
 end
 
+@testset "mask_replace_missing" begin
+    # Floating point rasters
+    a = Raster([1.0 0.0; 1.0 1.0], dims=(X, Y), missingval=0)
+    b = Raster([1.0 1.0; 1.0 0.0], dims=(X, Y), missingval=0)
+
+    # Integer rasters
+    c = Raster([1 0; 1 1], dims=(X, Y), missingval=0)
+    d = Raster([1 1; 1 0], dims=(X, Y), missingval=0)
+
+    # Test that missingval is replaced in source mask (Floats)
+    @test isequal(mask(a, with=b, missingval=3.14), [1.0 3.14; 1.0 3.14]) # Test missingval = 3.14
+    @test isequal(mask(a, with=b, missingval=missing), [1.0 missing; 1.0 missing]) # Test missingval = missing
+    @test isequal(mask(a, with=b, missingval=NaN), [1.0 NaN; 1.0 NaN]) # Test missingval = NaN
+    @test isequal(mask(a, with=b, missingval=NaN32), [1.0 NaN; 1.0 NaN]) # Test convert NaN32 to NaN
+    @test isequal(mask(a, with=b, missingval=Inf), [1.0 Inf; 1.0 Inf]) # Test missingval = Inf
+    @test_throws MethodError mask(a, with=b, missingval=nothing)
+
+    # Test that missingval is replaced in source mask (Ints)
+    @test isequal(mask(c, with=d, missingval=missing), [1 missing; 1 missing]) # Test missingval = missing
+    @test isequal(mask(c, with=d, missingval=-1.0), [1 -1; 1 -1])
+    @test_throws MethodError mask(c, with=d, missingval=nothing)
+    @test_throws InexactError mask(c, with=d, missingval=NaN)
+    @test_throws InexactError mask(c, with=d, missingval=3.14)
+    @test_throws InexactError mask(c, with=d, missingval=Inf)
+
+    # Test Type Stability
+    @test eltype(mask(a, with=b, missingval=0)) == Float64
+    @test eltype(mask(a, with=b, missingval=-1)) == Float64
+    @test eltype(mask(a, with=b, missingval=Inf32)) == Float64
+    @test eltype(mask(Float32.(a), with=b, missingval=Inf)) == Float32
+    @test eltype(mask(Float32.(a), with=b, missingval=NaN)) == Float32
+    @test eltype(mask(Float32.(a), with=b, missingval=0.0)) == Float32
+    @test eltype(mask(Float32.(a), with=b, missingval=0)) == Float32
+    @test eltype(mask(Float32.(a), with=b, missingval=-1)) == Float32
+    @test eltype(mask(c, with=d, missingval=-1.0)) == Int64
+    @test eltype(mask(c, with=d, missingval=0.0f0)) == Int64
+    @test eltype(mask(c, with=Float64.(d), missingval=-1.0)) == Int64
+    @test eltype(mask(c, with=Float64.(d), missingval=0.0f0)) == Int64
+
+    # Test mask!
+    @test_throws MethodError mask!(a, with=b, missingval=missing)
+    @test isequal(mask!(deepcopy(a), with=b, missingval=3.14), [1.0 3.14; 1.0 3.14]) # Test missingval = 3.14
+    @test isequal(mask!(deepcopy(a), with=b, missingval=NaN), [1.0 NaN; 1.0 NaN]) # Test missingval = NaN
+    @test isequal(mask!(deepcopy(a), with=b, missingval=NaN32), [1.0 NaN; 1.0 NaN]) # Test convert NaN32 to NaN
+    @test isequal(mask!(deepcopy(a), with=b, missingval=Inf), [1.0 Inf; 1.0 Inf]) # Test missingval = Inf
+    @test isequal(mask(deepcopy(c), with=d, missingval=-1.0), [1 -1; 1 -1])
+    @test_throws MethodError mask!(deepcopy(a), with=b, missingval=missing)
+    @test_throws ArgumentError mask!(deepcopy(c), with=d, missingval=nothing)
+    @test_throws InexactError mask!(deepcopy(c), with=d, missingval=NaN)
+    @test_throws InexactError mask!(deepcopy(c), with=d, missingval=3.14)
+    @test_throws InexactError mask!(deepcopy(c), with=d, missingval=Inf)
+end
+
 @testset "zonal" begin
     a = Raster((1:26) * (1:31)', (X(-20:5), Y(0:30)))
     pointvec_empty = [(-100.0, 0.0), (-100.0, 0.0), (-100.0, 0.0), (-100.0, 0.0), (-100.0, 0.0)]