kind of broken for points raster inputs

rafaqz · Dec 10, 2024 · 07efdaf · 07efdaf
1 parent 41ec62d
commit 07efdaf
Show file tree

Hide file tree

Showing 3 changed files with 63 additions and 29 deletions.
diff --git a/src/array.jl b/src/array.jl
@@ -294,13 +294,12 @@ function Raster(filename::AbstractString, dims::Tuple{<:Dimension,<:Dimension,Va
 )::Raster
     Raster(filename; dims, kw...)
 end
-Raster(ext::Extents.Extent; kw...) = Raster{Float64}(ext::Extents.Extent; kw...)
-function Raster{T}(ext::Extents.Extent; 
-    size=nothing, res=nothing, crs=nothing, mappedcrs=nothing, sampling=Points(),
-    kw...
+Raster(::UndefInitializer, ext::Extents.Extent; kw...) =
+    Raster{Float64}(undef, ext::Extents.Extent; kw...)
+function Raster{T}(::UndefInitializer, ext::Extents.Extent; 
+    size=nothing, res=nothing, crs=nothing, mappedcrs=nothing, sampling=Points(), closed=false, kw...
 ) where T
-@show sampling
-    dims = _extent2dims(ext; size, res, crs, mappedcrs, sampling)
+    dims = _extent2dims(ext; size, res, crs, mappedcrs, sampling, closed)
     Raster{T}(undef, dims; kw...)
 end
 function Raster(filename::AbstractString; 

diff --git a/src/utils.jl b/src/utils.jl
@@ -141,38 +141,64 @@ function _without_mapped_crs(f, st::AbstractRasterStack, mappedcrs::GeoFormat)
 end
 
 function _extent2dims(to; 
-    size=nothing, res=nothing, crs=nothing, sampling=Intervals(Start()), kw...
+    size=nothing, res=nothing, crs=nothing, sampling=Intervals(Start()), closed=true, kw...
 ) 
-    _extent2dims(to, size, res, crs, sampling)
+    _extent2dims(to, size, res, crs, sampling, closed)
 end
-function _extent2dims(to::Extents.Extent, size::Nothing, res::Nothing, crs, sampling)
+function _extent2dims(to::Extents.Extent, size::Nothing, res::Nothing, crs, sampling, closed)
     isnothing(res) && throw(ArgumentError("Pass either `size` or `res` keywords or a `Tuple` of `Dimension`s for `to`."))
 end
-function _extent2dims(to::Extents.Extent, size, res, crs, sampling)
+function _extent2dims(to::Extents.Extent, size, res, crs, sampling, closed)
     isnothing(res) || _size_and_res_error()
 end
-function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs, sampling) where K
+function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs, sampling, closed) where K
     tuple_res = ntuple(_ -> res, length(K))
-    _extent2dims(to, size, tuple_res, crs, sampling)
-end
-function _extent2dims(to::Extents.Extent{K}, size::Nothing, res, crs, sampling) where K
-    ranges = map(values(to), res) do bounds, r
-        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)
+    _extent2dims(to, size, tuple_res, crs, sampling, closed)
+end
+function _extent2dims(to::Extents.Extent{K}, size::Nothing, res, crs, sampling, closed) where K
+    ranges = if sampling isa Intervals 
+        map(values(to), res) do bounds, r
+            start, stop_closed = bounds
+            # Always add an eps buffer ignoring `closed`, it will 
+            # be trimmed in `range` anyway as the step size is fixed
+            stop_open = stop_closed + maybe_eps(stop_closed; grow=false) 
+            length = ceil(Int, (stop_open - start) / r)
+            range(; start, step=r, length)
+        end
+    else
+        map(values(to), res) do bounds, r
+            # No added eps for Points, as the Raster has closed interval extents
+            # TODO handle all of this in Extents.jl
+            start, stop = bounds
+            length = ceil(Int, (stop - start) / r) + 1
+            rnge = range(; start, step=r, length)
+            @show rnge r length start stop
+            rnge
+        end
     end
     return _extent2dims(to, ranges, crs, sampling)
 end
-function _extent2dims(to::Extents.Extent{K}, size, res::Nothing, crs, sampling) where K
+function _extent2dims(to::Extents.Extent{K}, size, res::Nothing, crs, sampling, closed) where K
     if size isa Int
         size = ntuple(_ -> size, length(K))
     end
-    ranges = map(values(to), size) do bounds, length
-        start, stop_closed = bounds
-        stop_open = stop_closed + maybe_eps(stop_closed; grow=false)
-        step = (stop_open - start) / length
-        range(; start, step, length)
+    ranges = if sampling isa Intervals
+        # Add an eps buffer for closed interval `Extent` with `Intervals`, 
+        # so that the open upper interval bounds of a Raster always include the 
+        # extent bounds and e.g. boundary points will be rasterized into the raster
+        map(values(to), size) do bounds, length
+            start, stop_closed = bounds
+            stop_open = closed ? stop_closed + maybe_eps(stop_closed; grow=false) : stop_closed
+            step = (stop_open - start) / length
+            range(; start, step, length)
+        end
+    else
+        # No added eps for Points, as the Raster has closed interval extents
+        map(values(to), size) do bounds, length
+            start, stop = bounds
+            step = (stop - start) / (length - 1)
+            range(; start, step, length)
+        end
     end
     return _extent2dims(to, ranges, crs, sampling)
 end

diff --git a/test/array.jl b/test/array.jl
@@ -79,20 +79,29 @@ end
     @test A6 isa Raster{Bool,2}
     @test A6 == [false false; false false] 
 
-    A7 = Raster(extent(A1); res=1.0)
+    A7 = Raster(undef, extent(A1); size=(2, 2))
     @test A7 isa Raster{Float64,2}
     @test size(A7) == (2, 2)
-    @test dims(A8) == dims(A1)
+    @test dims(A7) == dims(A1)
 
-    A8 = Raster{Int}(extent(A1); res=1.0)
+    A8 = Raster{Int}(undef, extent(A1); res=1.0)
     @test A8 isa Raster{Int,2}
     @test size(A8) == (2, 2)
     @test dims(A8) == dims(A1)
 
-    A9 = Raster{Bool}(extent(A1); res=1.0, sampling=Intervals(Center()))
+    A9 = Raster{Bool}(undef, extent(A1); res=1.0, sampling=Intervals(Center()), closed=false)
     @test A9 isa Raster{Bool,2}
     @test size(A9) == (2, 2)
+    # Has to the same lookups to include the bounds in the extent
+    @test collect.(dims(A9)) == collect.(dims(A1))
     @test sampling(A9, X) == Intervals(Center())
+
+    A10 = Raster{Bool}(undef, extent(A1); size=(2, 2), sampling=Intervals(Center()))
+    @test A10 isa Raster{Bool,2}
+    @test size(A10) == (2, 2)
+    # Has to the same lookups to include the bounds in the extent
+    @test collect.(dims(A9)) == collect.(dims(A1))
+    @test sampling(A10, X) == Intervals(Center())
 end