fix zonal for missing values in output (#442)

* fix zonal for missing values

* fix test

Author: rafaqz

src/methods/zonal.jl

@@ -127,7 +127,7 @@ end
 
 function _alloc_zonal(f, x, geom, n; kw...)
     z1 = _zonal(f, x, geom; kw...)
-    zs = Vector{typeof(z1)}(undef, n)
+    zs = Vector{Union{Missing,typeof(z1)}}(undef, n)
     zs[1] = z1
     return zs
 end

test/methods.jl

@@ -144,10 +144,13 @@ 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)]
+    polygon_empty = ArchGDAL.createpolygon(pointvec_empty)
     zonal(sum, a; of=polygon) ==
         zonal(sum, a; of=[polygon, polygon])[1] ==
+        zonal(sum, a; of=[polygon, polygon_empty])[1] ==
         zonal(sum, a; of=(geometry=polygon, x=:a, y=:b)) ==
-        zonal(sum, a; of=[(geometry=polygon, x=:a, y=:b)])[1]
+        zonal(sum, a; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
         zonal(sum, a; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
         sum(skipmissing(mask(a; with=polygon)))
     @test zonal(sum, a; of=a) ==