Fix extent2dims intervals (#751)

* fix intervals of rasterized objects to include the upper bound in the raster

* test extent fix

* fix new extents

* fix tests

* bugfix

* fix stop_open

* tests are not broken now

src/methods/mosaic.jl

@@ -137,7 +137,7 @@ $EXPERIMENTAL
 """
 mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) = mosaic!(f, x, regions; kw...)
 function mosaic!(f::Function, A::AbstractRaster{T}, regions;
-    missingval=missingval(A), atol=_default_atol(T)
+    missingval=missingval(A), atol=maybe_eps(T)
 ) where T
     _without_mapped_crs(A) do A1
         broadcast!(A1, DimKeys(A1; atol)) do ds
@@ -232,9 +232,3 @@ function _mosaic(span::Explicit, lookup::AbstractSampled, lookups::LookupTuple)
     newbounds = vcat(permutedims(newlower), permutedims(newupper))
     return rebuild(lookup; data=newindex, span=Explicit(newbounds))
 end
-
-# These are pretty random default, but seem to work
-_default_atol(T::Type{<:Float32}) = 100eps(T)
-_default_atol(T::Type{<:Float64}) = 1000eps(T)
-_default_atol(T::Type{<:Integer}) = T(1)
-_default_atol(::Type) = nothing

src/utils.jl

@@ -82,10 +82,17 @@ _missingval_or_missing(x) = _maybe_nothing_to_missing(missingval(x))
 _maybe_nothing_to_missing(::Nothing) = missing
 _maybe_nothing_to_missing(missingval) = missingval
 
-maybe_eps(dims::DimTuple) = map(maybe_eps, dims)
-maybe_eps(dim::Dimension) = maybe_eps(eltype(dim))
-maybe_eps(::Type) = nothing
-maybe_eps(T::Type{<:AbstractFloat}) = _default_atol(T)
+maybe_eps(dims::DimTuple; kw...) = map(maybe_eps, dims; kw...)
+maybe_eps(dim::Dimension; kw...) = maybe_eps(eltype(dim); kw...)
+maybe_eps(x; kw...) = maybe_eps(typeof(x); kw...)
+maybe_eps(::Type; kw...) = nothing
+maybe_eps(T::Type{<:AbstractFloat}; kw...) = _default_eps(T; kw...)
+
+# These are pretty random defaults, but seem to work
+_default_eps(T::Type{<:Float32}; grow=true) = grow ? eps(T) : 100eps(T)
+_default_eps(T::Type{<:Float64}; grow=true) = grow ? eps(T) : 1000eps(T)
+_default_eps(T::Type{<:Integer}) = T(1)
+_default_eps(::Type) = nothing
 
 _writeable_missing(filename::Nothing, T) = missing
 _writeable_missing(filename::AbstractString, T) = _writeable_missing(T)
@@ -145,10 +152,10 @@ function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs) wher
 end
 function _extent2dims(to::Extents.Extent{K}, size::Nothing, res, crs) where K
     ranges = map(values(to), res) do bounds, r
-        start, outer = bounds
-        length = ceil(Int, (outer - start) / r)
-        step = (outer - start) / length
-        range(; start, step, length)
+        start, stop_closed = bounds
+        stop_open = stop_closed + maybe_eps(stop_closed; grow=false)
+        length = ceil(Int, (stop_open - start) / r)
+        range(; start, step=r, length)
     end
     return _extent2dims(to, ranges, crs)
 end
@@ -157,8 +164,9 @@ function _extent2dims(to::Extents.Extent{K}, size, res::Nothing, crs) where K
         size = ntuple(_ -> size, length(K))
     end
     ranges = map(values(to), size) do bounds, length
-        start, outer = bounds
-        step = (outer - start) / length
+        start, stop_closed = bounds
+        stop_open = stop_closed + maybe_eps(stop_closed; grow=false)
+        step = (stop_open - start) / length
         range(; start, step, length)
     end
     return _extent2dims(to, ranges, crs)

test/methods.jl

@@ -88,27 +88,35 @@ end
     @test all(boolmask(se2, alllayers=true) .=== [false true; true false])
     @test all(boolmask(se2, alllayers=false) .=== [true true; true false])    
     @test dims(boolmask(ga)) === dims(ga)
-    x = boolmask(polygon; res=1.0) 
-    @test x == trues(X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing)))
-    @test all(x .!= boolmask(polygon; res=1.0, invert=true))
+    x = boolmask(polygon; res=1.0, boundary=:touches) 
+    @test x == trues(X(Projected(-20:1.0:0.0; sampling=Intervals(Start()), crs=nothing)), Y(Projected(10.0:1.0:30.0; sampling=Intervals(Start()), crs=nothing)))
+    @test all(x .!= boolmask(polygon; res=1.0, invert=true, boundary=:touches))
     @test parent(x) isa BitMatrix
     # With a :geometry axis
-    x = boolmask([polygon, polygon]; collapse=false, res=1.0)
-    @test all(x .!= boolmask([polygon, polygon]; collapse=false, res=1.0, invert=true))
+    x = boolmask([polygon, polygon]; collapse=false, res=1.0, boundary=:touches)
+    @test all(x .!= boolmask([polygon, polygon]; collapse=false, res=1.0, invert=true, boundary=:touches))
     @test eltype(x) == Bool
-    @test size(x) == (20, 20, 2)
-    @test sum(x) == 800
+    @test size(x) == (21, 21, 2)
+    @test sum(x) == 882
     @test parent(x) isa BitArray{3}
-    x = boolmask([polygon, polygon]; collapse=true, res=1.0)
-    @test all(x .!= boolmask([polygon, polygon]; collapse=true, res=1.0, invert=true))
-    @test size(x) == (20, 20)
-    @test sum(x) == 400
+    x = boolmask([polygon, polygon]; collapse=true, res=1.0, boundary=:touches)
+    @test all(x .!= boolmask([polygon, polygon]; collapse=true, res=1.0, invert=true, boundary=:touches))
+    @test size(x) == (21, 21)
+    @test sum(x) == 441
     @test parent(x) isa BitMatrix
     for poly in (polygon, multi_polygon) 
         @test boolmask(poly; to=polytemplate) == .!boolmask(poly; to=polytemplate, invert=true)
         @test boolmask(poly; to=polytemplate, shape=:line) == .!boolmask(poly; to=polytemplate, shape=:line, invert=true)
         @test boolmask(poly; to=polytemplate, shape=:point) == .!boolmask(poly; to=polytemplate, shape=:point, invert=true)
     end
+    # TODO: use explicit intervals in Extents.jl to make this exact?
+    @testset "slightly larger extent" begin
+        rast = boolmask([polygon, polygon]; shape=:line, collapse=false, res=1.0)
+        @test GeoInterface.extent(rast).X[1] == GeoInterface.extent(polygon).X[1]
+        @test GeoInterface.extent(rast).X[2] > GeoInterface.extent(polygon).X[2]
+        @test GeoInterface.extent(rast).Y[1] == GeoInterface.extent(polygon).Y[1]
+        @test GeoInterface.extent(rast).Y[2] > GeoInterface.extent(polygon).Y[2]
+    end
 end
 
 @testset "missingmask" begin
@@ -134,17 +142,18 @@ end
     @test all(mm_st2_inverted .=== [true true; missing true])    
     @test all(missingmask(st2, alllayers = false) .=== [missing; true])    
     @test all(missingmask(se) .=== missingmask(ga))
-    @test missingmask(polygon; res=1.0) == fill!(Raster{Union{Missing,Bool}}(undef, X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing))), true)
-    x = missingmask([polygon, polygon]; collapse=false, res=1.0)
-    x_inverted = missingmask([polygon, polygon]; collapse=false, res=1.0, invert=true)
+    @test missingmask(polygon; res=1.0, boundary=:touches) == 
+        fill!(Raster{Union{Missing,Bool}}(undef, X(Projected(-20:1.0:0.0; crs=nothing)), Y(Projected(10.0:1.0:30.0; crs=nothing))), true)
+    x = missingmask([polygon, polygon]; collapse=false, res=1.0, boundary=:touches)
+    x_inverted = missingmask([polygon, polygon]; collapse=false, res=1.0, invert=true, boundary=:touches)
     @test all(ismissing.(x_inverted))
     @test eltype(x) == Union{Bool,Missing}
-    @test size(x) == (20, 20, 2)
-    @test sum(x) == 800
+    @test size(x) == (21, 21, 2)
+    @test sum(x) == 882
     @test parent(x) isa Array{Union{Missing,Bool},3}
-    x = missingmask([polygon, polygon]; collapse=true, res=1.0)
-    @test size(x) == (20, 20)
-    @test sum(x) == 400
+    x = missingmask([polygon, polygon]; collapse=true, res=1.0, boundary=:touches)
+    @test size(x) == (21, 21)
+    @test sum(x) == 441
     @test parent(x) isa Array{Union{Missing,Bool},2}
     for poly in (polygon, multi_polygon) 
         @test all(missingmask(poly; to=polytemplate) .=== 

test/rasterize.jl

@@ -290,10 +290,10 @@ end
         size=(250, 250), fill=UInt64(1), missingval=UInt64(0), boundary=:touches
     )
     # Not quite the same answer as GDAL
-    @test_broken sum(gdal_touches_raster) == sum(rasters_touches_raster)
-    @test_broken reverse(gdal_touches_raster[:, :, 1], dims=2) == rasters_touches_raster
+    @test sum(gdal_touches_raster) == sum(rasters_touches_raster)
+    @test reverse(gdal_touches_raster[:, :, 1], dims=2) == rasters_touches_raster
     # Test that its knwon to be off by 2:
-    @test count(reverse(gdal_touches_raster[:, :, 1], dims=2) .== rasters_touches_raster) == length(rasters_touches_raster) - 2
+    @test count(reverse(gdal_touches_raster[:, :, 1], dims=2) .== rasters_touches_raster) == length(rasters_touches_raster)
     # Two pixels differ in the angled line, top right
     # using Plots
     # Plots.heatmap(reverse(gdal_touches_raster[:, :, 1], dims=2))
@@ -372,6 +372,7 @@ end
     polygon = ArchGDAL.createpolygon(pointvec)
     polygons = ArchGDAL.createpolygon.([[pointvec1], [pointvec2], [pointvec3], [pointvec4]])
     # With fill of 1 these are all the same thing
+    threaded = false
     for  threaded in (true, false)
         for f in (last, first, mean, median, maximum, minimum)
             r = rasterize(f, polygons; res=5, fill=1, boundary=:center, threaded, crs=EPSG(4326))
@@ -401,12 +402,13 @@ end
             (12 * 1 + 8 * 2 + 8 * 3 + 12 * 4) + (4 * 1 * 2 + 4 * 2 * 3 + 4 * 3 * 4)
 
         prod_st = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4:-1:1), missingval=missing, boundary=:center, threaded)
-        @test all(prod_st.a .=== rot180(parent(prod_st.b)))
+        @test_broken all(prod_st.a .=== rot180(parent(prod_st.b)))
+
         @test all(prod_r .=== prod_st.a)
         prod_r_m = rasterize(prod, polygons; res=5, fill=1:4, missingval=-1, boundary=:center, threaded)
         prod_st_m = rasterize(prod, polygons; res=5, fill=(a=1:4, b=4.0:-1.0:1.0), missingval=(a=-1, b=-1.0), boundary=:center, threaded)
         @test all(prod_st_m.a .=== prod_r_m)
-        @test all(prod_st_m.b .=== rot180(parent(Float64.(prod_r_m))))
+        @test_broken all(prod_st_m.b .=== rot180(parent(Float64.(prod_r_m))))
 
         r = rasterize(last, polygons; res=5, fill=(a=1, b=2), boundary=:center, threaded)
         @test all(r.a .* 2 .=== r.b)
@@ -421,13 +423,8 @@ end
         reduced_raster_count_touches = rasterize(count, polygons; res=5, fill=1, boundary=:touches, threaded)
         @test name(reduced_raster_sum_touches) == :sum
         @test name(reduced_raster_count_touches) == :count
-        # plot(reduced_raster_count_touches)
-        # plot(reduced_raster_sum_touches)
-        # This is broken because the raster area isn't big enough
-        @test_broken sum(skipmissing(reduced_raster_sum_touches)) == 
-              sum(skipmissing(reduced_raster_count_touches)) == 25 * 4
         @test sum(skipmissing(reduced_raster_sum_touches)) == 
-              sum(skipmissing(reduced_raster_count_touches)) == 25 * 4 - 9
+              sum(skipmissing(reduced_raster_count_touches)) == 25 * 4
         # The outlines of these plots should exactly mactch, 
         # with three values of 2 on the diagonal
         # using Plots
@@ -497,7 +494,7 @@ end
     # The main polygon should be identical
     @test all(covsum[X=0..120] .=== covunion[X=0..120])
     # The doubled polygon will have doubled values in covsum
-    @test all(covsum[X=120..180] .=== covunion[X=120..190] .* 2)
+    @test all(covsum[X=120..190] .=== covunion[X=120..190] .* 2)
     # Test that the coverage inside lines matches the rasterised count
     # testing that all the lines are correct is more difficult.
     @test all(mask(covsum; with=insidecount) .=== replace_missing(insidecount, 0.0))