test and warn when rasterizing empty things

src/methods/extract.jl

@@ -84,7 +84,7 @@ _prop_nt(A::AbstractRaster, I, names::NamedTuple{K}) where K = NamedTuple{K}((A[
 
 function _extract(x::RasterStackOrArray, ::GI.PointTrait, point; dims, names, atol=nothing)
     # Get the actual dimensions available in the object
-    coords = map(DD.dims(x)) do d
+    coords = map(DD.commondims(x, dims)) do d
         _dimcoord(d, point)
     end
 

src/methods/rasterize.jl

@@ -460,6 +460,9 @@ function alloc_rasterize(f, r::RasterCreator;
     metadata=r.metadata,
     suffix=r.suffix,
 )
+    if prod(size(r.to)) == 0  
+        throw(ArgumentError("Destination array is is empty, with size $(size(r.to))). Rasterization is not possible"))
+    end
     A = create(r.filename, eltype, r.to; name, missingval, metadata, suffix)
     # TODO f should apply to the file when it is initially created
     # instead of reopening but we need a `create(f, filename, ...)` method
@@ -536,6 +539,10 @@ function rasterize!(reducer::typeof(count), x::RasterStackOrArray, data; fill=no
     rasterize!(x::RasterStackOrArray, data; kw..., reducer=nothing, op=nothing, fill=_count_fill, init=0)
 end
 function rasterize!(x::RasterStackOrArray, data; threaded=true, kw...)
+    if prod(size(x)) == 0  
+        @warn "Destination is empty, rasterization skipped"
+        return x
+    end
     r = Rasterizer(data; eltype=eltype(x), threaded, kw...)
     allocs = r.shape == :points ? nothing : _burning_allocs(dims(x); threaded)
     return _rasterize!(x, r; allocs)

src/polygon_ops.jl

@@ -130,15 +130,15 @@ Base.getindex(edges::Edges, I...) = getindex(parent(edges), I...)
 Base.setindex(edges::Edges, x, I...) = setindex!(parent(edges), x, I...)
 
 @noinline function _to_edges!(edges, geom, dims, edge_count)
-    GI.npoint(geom) > 0 || return edge_count
-    xlookup, ylookup = lookup(dims, (X(), Y())) 
-    (length(xlookup) > 0 && length(ylookup) > 0) || return edge_count
-
     # Dummy Initialisation
     local firstpos = prevpos = nextpos = Position((0.0, 0.0), 0)
     isfirst = true
     local max_ylen = 0
 
+    GI.npoint(geom) > 0 || return edge_count, max_ylen
+    xlookup, ylookup = lookup(dims, (X(), Y())) 
+    (length(xlookup) > 0 && length(ylookup) > 0) || return edge_count, max_ylen
+
     # Raster properties
     starts = (Float64(first(xlookup)), Float64(first(ylookup)))
     steps = (Float64(Base.step(xlookup)), Float64(Base.step(ylookup)))
@@ -262,6 +262,7 @@ function _burn_geometry!(B::AbstractRaster, ::GI.AbstractGeometryTrait, geom;
     shape=nothing, verbose=true, boundary=:center, allocs=nothing, kw...
 )::Bool
     hasburned = false
+    GI.npoint(geom) > 0 || return false
     # Use the specified shape or detect it
     shape = shape isa Symbol ? shape : _geom_shape(geom)
     if shape === :point

test/rasterize.jl

@@ -413,6 +413,18 @@ end
     @test sum(A3) == 1200
 end
 
+@testset "Rasterize an empty polygon" begin
+    empty_polygon = ArchGDAL.createpolygon(Tuple{Float64,Float64}[])
+    rast = rasterize(empty_polygon; to=A1, fill=1, missingval=0)
+    @test sum(rast) == 0
+end
+
+@testset "Cant rasterize to a new empty array" begin
+    # Its difficult to make a new empty raster here (we possibly could in future?)
+    @test_throws ArgumentError rasterize(polygon; to=A1[1:0, 1:0], fill=1, missingval=0)
+    # But we just warn in `rasterize!` because the result is still correct
+    @test_warn "Destination is empty" rasterize!(A1[1:0, 1:0], polygon; to=A1[1:0, 1:0], fill=1, missingval=0)
+end
 
 @testset "coverage" begin
     @time covsum = coverage(sum, shphandle.shapes; res=1, scale=10)