diff --git a/src/Rasters.jl b/src/Rasters.jl
index 94af83dfe..0e7f77416 100644
--- a/src/Rasters.jl
+++ b/src/Rasters.jl
@@ -131,6 +131,7 @@ include("methods/burning/line.jl")
 include("methods/burning/polygon.jl")
 include("methods/burning/extents.jl")
 include("methods/burning/utils.jl")
+include("methods/burning/reductions.jl")
 
 include("methods/mask.jl")
 include("methods/rasterize.jl")
diff --git a/src/methods/burning/allocs.jl b/src/methods/burning/allocs.jl
index 7956e266f..45ffe4a5c 100644
--- a/src/methods/burning/allocs.jl
+++ b/src/methods/burning/allocs.jl
@@ -11,30 +11,35 @@ function Allocs(buffer)
     return Allocs(buffer, edges, scratch, crossings)
 end
 
+function _burning_allocs(x::Nothing; 
+    nthreads=_nthreads(), 
+    threaded=true, 
+    burncheck_metadata=nothing,
+    kw...
+) 
+    threaded ? [Allocs(nothing) for _ in 1:nthreads] : [Allocs(nothing)]
+end
 function _burning_allocs(x; 
     nthreads=_nthreads(), 
     threaded=true, 
     burncheck_metadata=Metadata(),
     kw...
 ) 
+    dims = commondims(x, (XDim, YDim))
     if threaded
-        if isnothing(x)
-            [Allocs(nothing) for _ in 1:nthreads]
-        else
-            dims = commondims(x, DEFAULT_POINT_ORDER)
-            [Allocs(_init_bools(dims; metadata=deepcopy(burncheck_metadata))) for _ in 1:nthreads]
-        end
+        [Allocs(_init_bools(dims; metadata=deepcopy(burncheck_metadata))) for _ in 1:nthreads]
     else
-        if isnothing(x)
-            Allocs()
-        else
-            dims = commondims(x, DEFAULT_POINT_ORDER)
-            Allocs(_init_bools(dims; metadata=burncheck_metadata))
-        end
+        [Allocs(_init_bools(dims; metadata=burncheck_metadata))]
     end
 end
 
-_get_alloc(allocs::Vector{<:Allocs}) = _get_alloc(allocs[Threads.threadid()])
+function _get_alloc(allocs::Vector{<:Allocs}) 
+    if length(allocs) == 1
+        only(allocs)
+    else
+        allocs[Threads.threadid()]
+    end
+end
 _get_alloc(allocs::Allocs) = allocs
 
 # TODO include these in Allocs
diff --git a/src/methods/burning/line.jl b/src/methods/burning/line.jl
index 24f924f6c..290bcc395 100644
--- a/src/methods/burning/line.jl
+++ b/src/methods/burning/line.jl
@@ -1,3 +1,44 @@
+#= This mutable object is passed into line burning closures,
+as updating object fields is type-stable in closures but 
+updating variables is not
+i: number of valid pixels found while line burning so far
+prev: previous pixel index to avoid repeats at start/end of lines
+vals: vector for extracting values to
+acc: an optional accumulator for reductions
+=#
+mutable struct LineRefs{T}
+    i::Int
+    prev::CartesianIndex{2}
+    vals::Vector{T}
+    acc::T
+    function LineRefs{T}(; init::Union{T,Nothing}=nothing) where T
+        i = 1
+        prev = CartesianIndex((typemin(Int), typemin(Int)))
+        vals = Vector{T}()
+        if isnothing(init)
+            new{T}(i, prev, vals)
+        else
+            new{T}(i, prev, vals, init)
+        end
+    end
+end
+
+# Re-initialise line refs
+function _initialise!(lr::LineRefs{T}; 
+    vals=true, init::Union{T,Nothing}=nothing
+) where T
+    lr.i = 1
+    lr.prev = CartesianIndex((typemin(Int), typemin(Int)))
+    if !isnothing(init)
+        lr.acc = init
+    end
+    if vals 
+        lr.vals = Vector{T}()
+    end
+    return lr
+end
+
+
 # _burn_lines!
 # Fill a raster with `fill` where pixels touch lines in a geom
 # Usually `fill` is `true` of `false`
diff --git a/src/methods/burning/reductions.jl b/src/methods/burning/reductions.jl
new file mode 100644
index 000000000..00ad19dac
--- /dev/null
+++ b/src/methods/burning/reductions.jl
@@ -0,0 +1,39 @@
+# _TakeFirst: `first` with a missing value check
+# TODO: clarify why we need this  for `first` and not `last`
+struct _TakeFirst{MV}
+    missingval::MV
+end
+(tf::_TakeFirst)(a, b) = a === tf.missingval ? b : a
+_take_last(a, b) = b
+
+# Get the operator for a reduction. 
+# This can be used for online stats style optimizations
+_reduce_op(::typeof(sum)) = Base.add_sum
+_reduce_op(::typeof(prod)) = Base.mul_prod
+_reduce_op(::typeof(minimum)) = min
+_reduce_op(::typeof(maximum)) = max
+_reduce_op(::typeof(last)) = _take_last
+_reduce_op(f, missingval) = _reduce_op(f)
+_reduce_op(::typeof(first), missingval) = _TakeFirst(missingval)
+_reduce_op(x) = nothing
+
+# Can we split the reducion and combine later
+_is_op_threadsafe(::typeof(sum)) = true
+_is_op_threadsafe(::typeof(prod)) = true
+_is_op_threadsafe(::typeof(minimum)) = true
+_is_op_threadsafe(::typeof(maximum)) = true
+_is_op_threadsafe(f) = false
+
+# Get an initial value for the reduction
+# TODO: clarify the difference between returning missing and zero
+_reduce_init(reducer, st::AbstractRasterStack, missingval) = map(A -> _reduce_init(reducer, A, missingval), st)
+_reduce_init(reducer, ::AbstractRaster{T}, missingval) where T = _reduce_init(reducer, T, missingval)
+_reduce_init(reducer, nt::NamedTuple, missingval) = map(x -> _reduce_init(reducer, x, missingval), nt)
+_reduce_init(f, x, missingval) = _reduce_init(f, typeof(x), missingval)
+_reduce_init(::Nothing, x::Type{T}, missingval) where T = zero(T)
+_reduce_init(f::Function, ::Type{T}, missingval) where T = zero(f((zero(nonmissingtype(T)), zero(nonmissingtype(T)))))
+_reduce_init(::typeof(sum), ::Type{T}, missingval) where T = zero(nonmissingtype(T))
+_reduce_init(::typeof(prod), ::Type{T}, missingval) where T = oneunit(nonmissingtype(T))
+_reduce_init(::typeof(minimum), ::Type{T}, missingval) where T = typemax(nonmissingtype(T))
+_reduce_init(::typeof(maximum), ::Type{T}, missingval) where T = typemin(nonmissingtype(T))
+_reduce_init(::typeof(last), ::Type{T}, missingval) where T = _maybe_nothing_to_missing(missingval)
\ No newline at end of file
diff --git a/src/methods/extract.jl b/src/methods/extract.jl
index ce1ee4b93..b9c2e78e3 100644
--- a/src/methods/extract.jl
+++ b/src/methods/extract.jl
@@ -34,25 +34,6 @@ end
 
 Base.eltype(::Extractor{T}) where T = T
 
-# This mutable object is passed into closures as 
-# fields are type-stable in clusores but variables are not
-mutable struct LineRefs{T}
-    i::Int
-    prev::CartesianIndex{2}
-    rows::Vector{T}
-    function LineRefs{T}() where T
-        i = 1
-        prev = CartesianIndex((typemin(Int), typemin(Int)))
-        new{T}(i, prev, Vector{T}())
-    end
-end
-
-function _initialise!(lr::LineRefs{T}) where T
-    lr.i = 1
-    lr.prev = CartesianIndex((typemin(Int), typemin(Int)))
-    lr.rows = Vector{T}()
-end
-
 _geom_rowtype(A, geom; kw...) =
     _geom_rowtype(A, GI.trait(geom), geom; kw...)
 _geom_rowtype(A, ::Nothing, geoms; kw...) =
@@ -297,9 +278,7 @@ end
     # places when they are right on the line.
     dims, offset = _template(A, geom)
     dp = DimPoints(A)
-    function _length_callback(n)
-        resize!(line_refs.rows, n)
-    end
+    _length_callback(n) = resize!(line_refs.vals, n)
 
     _burn_lines!(_length_callback, dims, geom) do D
         I = CartesianIndex(map(val, D))
@@ -308,11 +287,11 @@ end
         line_refs.prev = I
         # Make sure we only hit this pixel once
         # D is always inbounds
-        line_refs.i += _maybe_set_row!(line_refs.rows, e.skipmissing, e, id, A, dp, offset, I, line_refs.i)
+        line_refs.i += _maybe_set_row!(line_refs.vals, e.skipmissing, e, id, A, dp, offset, I, line_refs.i)
         return nothing
     end
-    deleteat!(line_refs.rows, line_refs.i:length(line_refs.rows))
-    return line_refs.rows
+    deleteat!(line_refs.vals, line_refs.i:length(line_refs.vals))
+    return line_refs.vals
 end
 @noinline function _extract(
     A::RasterStackOrArray, e::Extractor{T}, id::Int, ::GI.AbstractGeometryTrait, geom; kw...
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index 674b3b965..99b8b7422 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -1,37 +1,3 @@
-struct _TakeFirst{MV}
-    missingval::MV
-end
-(tf::_TakeFirst)(a, b) = a === tf.missingval ? b : a
-_take_last(a, b) = b
-
-_reduce_op(::typeof(sum)) = Base.add_sum
-_reduce_op(::typeof(prod)) = Base.mul_prod
-_reduce_op(::typeof(minimum)) = min
-_reduce_op(::typeof(maximum)) = max
-_reduce_op(::typeof(last)) = _take_last
-_reduce_op(f, missingval) = _reduce_op(f)
-_reduce_op(::typeof(first), missingval) = _TakeFirst(missingval)
-_reduce_op(x) = nothing
-
-_is_op_threadsafe(::typeof(sum)) = true
-_is_op_threadsafe(::typeof(prod)) = true
-_is_op_threadsafe(::typeof(minimum)) = true
-_is_op_threadsafe(::typeof(maximum)) = true
-_is_op_threadsafe(f) = false
-
-_reduce_init(reducer, st::AbstractRasterStack, missingval) = map(A -> _reduce_init(reducer, A, missingval), st)
-_reduce_init(reducer, ::AbstractRaster{T}, missingval) where T = _reduce_init(reducer, T, missingval)
-_reduce_init(reducer, nt::NamedTuple, missingval) = map(x -> _reduce_init(reducer, x, missingval), nt)
-_reduce_init(f, x, missingval) = _reduce_init(f, typeof(x), missingval)
-
-_reduce_init(::Nothing, x::Type{T}, missingval) where T = zero(T)
-_reduce_init(f::Function, ::Type{T}, missingval) where T = zero(f((zero(nonmissingtype(T)), zero(nonmissingtype(T)))))
-_reduce_init(::typeof(sum), ::Type{T}, missingval) where T = zero(nonmissingtype(T))
-_reduce_init(::typeof(prod), ::Type{T}, missingval) where T = oneunit(nonmissingtype(T))
-_reduce_init(::typeof(minimum), ::Type{T}, missingval) where T = typemax(nonmissingtype(T))
-_reduce_init(::typeof(maximum), ::Type{T}, missingval) where T = typemin(nonmissingtype(T))
-_reduce_init(::typeof(last), ::Type{T}, missingval) where T = _maybe_nothing_to_missing(missingval)
-
 struct FillChooser{F,I,M}
     fill::F
     init::I
diff --git a/src/methods/zonal.jl b/src/methods/zonal.jl
index 97214be2a..d8c41b4f6 100644
--- a/src/methods/zonal.jl
+++ b/src/methods/zonal.jl
@@ -72,17 +72,22 @@ insertcols!(january_stats, 1, :country => first.(split.(countries.ADMIN, r"[^A-Z
                                                   3 columns and 243 rows omitted
 ```
 """
-zonal(f, x::RasterStackOrArray; of, skipmissing=true, kw...) = 
-    _zonal(f, _prepare_for_burning(x), of; skipmissing, kw...)
+Base.@constprop :aggressive @inline function zonal(f::F, x::RasterStackOrArray; 
+    of, threaded=true, progress=true, skipmissing=true, kw...
+) where F
+    _zonal(f, _prepare_for_burning(x), of; 
+        threaded, skipmissing, progress, kw...
+    )
+end
 
-_zonal(f, x::RasterStackOrArray, of::RasterStackOrArray; kw...) = 
-    _zonal(f, x, Extents.extent(of); kw...)
-_zonal(f, x::RasterStackOrArray, of::DimTuple; kw...) = 
+_zonal(f::F, x::RasterStackOrArray, of::RasterStackOrArray; kw...) where F = 
+    _zonal(f::F, x, Extents.extent(of); kw...)
+_zonal(f::F, x::RasterStackOrArray, of::DimTuple; kw...) where F= 
     _zonal(f, x, Extents.extent(of); kw...)
 # We don't need to `mask` with an extent, it's square so `crop` will do enough.
-_zonal(f, x::Raster, of::Extents.Extent; skipmissing) =
+_zonal(f::F, x::Raster, of::Extents.Extent; skipmissing) where F =
     _maybe_skipmissing_call(f, crop(x; to=of, touches=true), skipmissing)
-function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing)
+function _zonal(f::F, x::RasterStack, ext::Extents.Extent; skipmissing) where F
     cropped = crop(x; to=ext, touches=true)
     prod(size(cropped)) > 0 || return missing
     return maplayers(cropped) do A
@@ -90,39 +95,51 @@ function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing)
     end
 end
 # Otherwise of is a geom, table or vector
-_zonal(f, x::RasterStackOrArray, of; kw...) = _zonal(f, x, GI.trait(of), of; kw...)
+_zonal(f::F, x::RasterStackOrArray, of; kw...) where F =
+    _zonal(f, x, GI.trait(of), of; kw...)
 
 _zonal(f, x, ::GI.AbstractFeatureCollectionTrait, fc; kw...) =
     _zonal(f, x, nothing, fc; kw...)
 _zonal(f, x::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...) =
     _zonal(f, x, GI.geometry(feature); kw...)
-function _zonal(f, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; 
-    skipmissing, kw...
-)
+function _zonal(f::F, x::AbstractRaster, ::GI.AbstractGeometryTrait, geom; 
+    line_refs=nothing, skipmissing, threaded, kw...
+) where F
     cropped = crop(x; to=geom, touches=true)
     prod(size(cropped)) > 0 || return missing
     # Fast path for lines
-    if skipmissing && trait isa GI.AbstractLineStringTrait
+    if istrue(skipmissing) && GI.trait(geom) isa GI.AbstractLineStringTrait
+        if isnothing(line_refs) 
+            T = Base.promote_op(f, eltype(x))
+            init = _reduce_init(f, T, missingval(x))
+            line_refs = LineRefs{T}(; init)
+        end
         # Read the lines directly and calculate stats on the fly
-        return _zonal_lines(f, _reduce_op(f), st, geom; kw...)
+        return _zonal_lines(f, _reduce_op(f), cropped, geom; line_refs, kw...)
     else
         # Mask the cropped raster 
-        # TODO use mask! on re-used buffer where possible
-        masked = mask(cropped; with=geom, kw...)
-        # Calculate stats on whole raster or `skipmissing(raster)``    
-        return _maybe_skipmissing_call(f, masked, skipmissing)
+        if istrue(skipmissing)
+            # `boolmask` allocates a lot less than `mask` so use it where we can
+            bitmask = boolmask(geom; to=cropped, kw...)
+            # Call `f` on an iterator of unmasked and non-missing values
+            f(v for (m, v) in zip(bitmask, cropped) if m && !_ismissingval(x, v))
+        else
+            masked = mask(cropped; with=geom, kw...)
+            f(masked)
+        end
     end
 end
-function _zonal(f, st::AbstractRasterStack, trait::GI.AbstractGeometryTrait, geom; 
-    skipmissing=true, kw...
-)
-    cropped = crop(st; to=geom, touches=true)
-    prod(size(cropped)) > 0 || return map(_ -> missing, layerdims(st))
-    # Fast path for lines
-    if skipmissing && trait isa GI.AbstractLineStringTrait
-        # Read the lines directly and calculate stats on the fly
-        return _zonal_lines(f, _reduce_op(f), cropped, geom; kw...)
+function _zonal(f::F, st::AbstractRasterStack, trait::GI.AbstractGeometryTrait, geom; 
+    skipmissing, threaded, kw...
+) where F
+    if istrue(skipmissing) && trait isa GI.AbstractLineStringTrait
+        # Fast path for lines
+        maplayers(cropped) do layer
+            _zonal(f, layers, trait, geom; skipmissing, kw...)
+        end
     else
+        cropped = crop(st; to=geom, touches=true)
+        prod(size(cropped)) > 0 || return map(_ -> missing, layerdims(st))
         # Mask the cropped stack 
         # TODO use mask! on re-used buffer where possible
         masked = mask(cropped; with=geom, kw...)
@@ -133,41 +150,113 @@ function _zonal(f, st::AbstractRasterStack, trait::GI.AbstractGeometryTrait, geo
         end
     end
 end
-function _zonal(f, x::RasterStackOrArray, ::Nothing, data; 
+function _zonal(f::F, x::RasterStackOrArray, ::Nothing, data; 
     geometrycolumn=nothing, kw...
-)
+) where F
     geoms = _get_geometries(data, geometrycolumn)
     _zonal(f, x, nothing, geoms; kw...) 
 end
-function _zonal(f, x::RasterStackOrArray, ::Nothing, geoms::AbstractArray; 
-    progress=true, threaded=true, geometrycolumn=nothing, kw...
-)
+function _zonal(f::F, x::RasterStackOrArray, ::Nothing, geoms::AbstractArray; 
+    progress, threaded, skipmissing, kw...
+) where F
     n = length(geoms)
     n == 0 && return [] # TODO make this type stable
     # Allocate vector by running _zonal on the first geoms
-    # until a non-missing value is returns
-    zs, start_index = _alloc_zonal(f, x, geoms, n; kw...)
+    # until a non-missing value is returned
+    zs, start_index = _alloc_zonal(f, x, geoms, n; skipmissing, kw...)
     # If all geometries are missing, return a vector of missings
     start_index == n + 1 && return zs
-    # Otherwise run _zonal for the rest of the geometries, possibly threaded
-    _run(start_index:n, threaded, progress, "Applying $f to each geometry...") do i
-        zs[i] = _zonal(f, x, geoms[i]; kw...)
+    # If skipmising and 2d for AbstractLineStringTrait we can use optimisations 
+    if x isa AbstractRaster && ndims(x) == 2 && istrue(skipmissing) && all(x -> ismissing(x) || GI.trait(x) isa GI.LineStringTrait, geoms)
+        # TODO make this work for RasterStack
+        T = Base.promote_op(f, eltype(x))
+        init = _reduce_init(f, T, missingval(x))
+        if istrue(threaded)
+            thread_line_refs = [LineRefs{T}(; init) for _ in 1:Threads.nthreads()]
+            _run(start_index:n, threaded, progress, "Applying $f to each geometry...") do i
+                line_refs = thread_line_refs[Threads.threadid()]
+                zs[i] = _zonal(f, x, geoms[i]; line_refs, skipmissing, kw...)
+            end
+        else
+            line_refs = LineRefs{T}(; init)
+            _run(start_index:n, threaded, progress, "Applying $f to each geometry...") do i
+                zs[i] = _zonal(f, x, geoms[i]; line_refs, skipmissing, kw...)
+            end
+        end
+    else
+        _run(start_index:n, threaded, progress, "Applying $f to each geometry...") do i
+            zs[i] = _zonal(f, x, geoms[i]; skipmissing, kw...)
+        end
     end
     return zs
 end
 
-function _zonal_lines(f, op, r::RasterStackOrArray, geom; 
-    buffer, kw...
-)
-    burn_lines(rast, geom) do i
-        buffer[i] = rast[i]
+# Generic reductions of a vector of values for the line
+function _zonal_lines(f::F, op::Nothing, A::AbstractRaster{<:Any,2}, geom; 
+    bitmask=_init_bools(A),
+    line_refs, 
+    kw...
+) where F
+        _initialise!(line_refs; vals=false)
+    # Define a callback to resize the row when 
+    # the line burner knowns how long it is
+    function _length_callback(n) 
+        resize!(line_refs.vals, line_refs.i + n - 1)
+    end
+    # Burn lines and update row values in the closure
+    _burn_lines!(_length_callback, dims(A), geom) do D
+        I = CartesianIndex(map(val, D))
+        # If we've already included this pixel, skip
+        bitmask[I] && return nothing
+        bitmask[I] = true
+        i = line_refs.i
+        v = A[D]
+        # We are skipping missing values
+        # _ismissingval(A, v) && return nothing
+        # Set the row at `i`` to the raster value at `I``
+        line_refs.vals[i] = v
+        # Increment row position
+        line_refs.i = i + 1
+        return nothing
     end
+    # Finally resize! again to the final lenght
+    resize!(line_refs.vals, line_refs.i - 1)
+    # And apply the reduction and return
+    return f(line_refs.vals)
 end
-function _zonal_lines(f, op::Nothing, st::RasterStackOrArray, geom; kw...)
-    # burn and count
+# Iterative non-allocating reductions for known functions
+function _zonal_lines(f::F, op::O, A::AbstractRaster{<:Any,2}, geom;
+    bitmask=_init_bools(A),
+    line_refs,
+    kw...
+) where {F,O}
+    function _length_callback(n) 
+        _initialise!(line_refs; vals=false)
+    end
+    _burn_lines!(identity, dims(A), geom) do D
+        I = CartesianIndex(map(val, D))
+        # If we've already included this pixel, skip
+        # If we've already included this pixel, skip
+        bitmask[I] && return nothing
+        bitmask[I] = true
+        # Retrive the value at I
+        v = A[I]
+        # Check v is not missing
+        _ismissingval(A, v) && return nothing
+        # If acc is not missing run op, otherwise v
+        newv = if _ismissingval(A, line_refs.acc)
+            v
+        else
+            op(line_refs.acc, v) 
+        end
+        line_refs.acc = newv
+        return nothing
+    end
+    # Return the accululated value
+    return line_refs.acc
 end
 
-function _alloc_zonal(f, x, geoms, n; kw...)
+function _alloc_zonal(f::F, x, geoms, n; kw...) where F
     # Find first non-missing entry and count number of missing entries
     n_missing::Int = 0
     z1 = missing
@@ -179,12 +268,15 @@ function _alloc_zonal(f, x, geoms, n; kw...)
     end
     zs = Vector{Union{Missing,typeof(z1)}}(undef, n)
     zs[1:n_missing] .= missing
-    # Exit early when all elements are missing
     if n_missing == n
+        # All elements are missing
         return zs, n_missing + 1
+    else
+        # Set the first non-missing value to z1
+        zs[n_missing + 1] = z1
+        # Start after the first non-missing value
+        return zs, n_missing + 2
     end
-    zs[n_missing + 1] = z1
-    return zs, n_missing + 1
 end
 
-_maybe_skipmissing_call(f, A, sm) = sm ? f(skipmissing(A)) : f(A)
+_maybe_skipmissing_call(f::F, A, sm) where F = sm ? f(skipmissing(A)) : f(A)
diff --git a/src/utils.jl b/src/utils.jl
index 125f925e8..9161d3752 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -274,9 +274,9 @@ const URLREGEX = r"^[a-zA-Z][a-zA-Z\d+\-.]*:"
 _isurl(str::AbstractString) = !occursin(WINDOWSREGEX, str) && occursin(URLREGEX, str)
 
 # Run `f` threaded or not, w
-function _run(f, range::OrdinalRange, threaded::Bool, progress::Bool, desc::String) 
-    p = progress ? _progress(length(range); desc) : nothing
-    if threaded
+function _run(f, range::OrdinalRange, threaded, progress, desc::String) 
+    p = istrue(progress) ? _progress(length(range); desc) : nothing
+    if istrue(threaded)
         Threads.@threads :static for i in range
             f(i)
             isnothing(p) || ProgressMeter.next!(p)
@@ -379,6 +379,7 @@ _booltype(::_True) = _True()
 _booltype(::_False) = _False()
 _booltype(x) = x ? _True() : _False()
 
+istrue(b::Bool) = true
 istrue(::_True) = true
 istrue(::_False) = false
 istrue(::Type{_True}) = true
diff --git a/test/methods.jl b/test/methods.jl
index 2975ea286..6aa860ca8 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -286,59 +286,6 @@ end
     @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)]
-    polygon_empty = ArchGDAL.createpolygon(pointvec_empty)
-    @test 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] ==
-        sum(skipmissing(mask(a; with=polygon)))
-    @test zonal(sum, a; of=a) == 
-        zonal(sum, a; of=dims(a)) ==
-        zonal(sum, a; of=Extents.extent(a)) == 
-        sum(a)
-
-    b = a .* 2
-    c = cat(a .* 3, a .* 3; dims=:newdim)
-    st = RasterStack((; a, b, c))
-    @test zonal(sum, st; of=polygon) == zonal(sum, st; of=[polygon])[1] ==
-        zonal(sum, st; of=(geometry=polygon, x=:a, y=:b)) ==
-        zonal(sum, st; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
-        zonal(sum, st; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
-        maplayers(sum ∘ skipmissing, mask(st; with=polygon))  
-    @test zonal(sum, st; of=st) == 
-        zonal(sum, st; of=dims(st)) == 
-        zonal(sum, st; of=Extents.extent(st)) == 
-        sum(st)
-
-    @testset "skipmissing" begin
-        a = Array{Union{Missing,Int}}(undef, 26, 31)
-        a .= (1:26) * (1:31)'
-        a[1:10, 3:10] .= missing
-        rast = Raster(a, (X(-20:5), Y(0:30)))
-        @test zonal(sum, rast; of=polygon, skipmissing=false) === missing
-        @test zonal(sum, rast; of=polygon, skipmissing=true) isa Int
-        @test !zonal(x -> x isa Raster, rast; of=polygon, skipmissing=true)
-        @test zonal(x -> x isa Raster, rast; of=polygon, skipmissing=false)
-    end
-end
-
-@testset "zonal return missing" begin
-    a = Raster((1:26) * (1:31)', (X(-20:5), Y(0:30)))
-    out_bounds_pointvec = [(-40.0, -40.0), (-40.0, -35.0), (-35.0, -35.0), (-35.0, -40.0)]
-    out_bounds_polygon = ArchGDAL.createpolygon(out_bounds_pointvec)
-    @test ismissing(zonal(sum, a; of=[polygon, out_bounds_polygon, polygon])[2]) &&
-        ismissing(zonal(sum, a; of=[out_bounds_polygon, polygon])[1]) &&
-        ismissing(zonal(sum, a; of=(geometry=out_bounds_polygon, x=:a, y=:b))) &&
-        ismissing(zonal(sum, a; of=[(geometry=out_bounds_polygon, x=:a, y=:b)])[1])
-    @test zonal(sum, a; of=[out_bounds_polygon, out_bounds_polygon, polygon])[3] == 
-        sum(skipmissing(mask(a; with=polygon)))
-end
-
 @testset "classify" begin
     A1 = [missing 1; 2 3]
     ga1 = Raster(A1, (X, Y); missingval=missing)
diff --git a/test/runtests.jl b/test/runtests.jl
index 93b2c9ec6..e940db9c3 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -20,6 +20,7 @@ end
 @time @safetestset "aggregate" begin include("aggregate.jl") end
 @time @safetestset "rasterize" begin include("rasterize.jl") end
 @time @safetestset "extract" begin include("extract.jl") end
+@time @safetestset "zonal" begin include("zonal.jl") end
 @time @safetestset "reproject" begin include("reproject.jl") end
 @time @safetestset "warp" begin include("warp.jl") end
 @time @safetestset "cellarea" begin include("cellarea.jl") end
diff --git a/test/zonal.jl b/test/zonal.jl
new file mode 100644
index 000000000..3f048c705
--- /dev/null
+++ b/test/zonal.jl
@@ -0,0 +1,89 @@
+using Rasters, Extents, ArchGDAL, Test
+import GeoInterface as GI
+
+a = Raster((1:26) * (1:31)', (X(-20:5), Y(0:30)))
+b = a .* 2
+c = cat(a .* 3, a .* 3; dims=:newdim)
+st = RasterStack((; a, b, c))
+a_m = let 
+    a = Array{Union{Missing,Int}}(undef, 26, 31)
+    a .= (1:26) * (1:31)'
+    a[1:10, 3:10] .= missing
+    Raster(a, (X(-20:5), Y(0:30)))
+end
+
+pointvec = [(-20.0, 30.0), (-20.0, 10.0), (0.0, 10.0), (0.0, 30.0), (-20.0, 30.0)]
+pointvec_empty = [(-100.0, 0.0), (-100.0, 0.0), (-100.0, 0.0), (-100.0, 0.0), (-100.0, 0.0)]
+pointvec_out_of_bounds = [(-40.0, -40.0), (-40.0, -35.0), (-35.0, -35.0), (-35.0, -40.0)]
+linestring = GI.LineString(pointvec)
+polygon = GI.Polygon([pointvec])
+polygon_empty = GI.Polygon([pointvec_empty])
+polygon_out_of_bounds = GI.Polygon([pointvec_out_of_bounds])
+
+@testset "polygons" begin
+    @test 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] ==
+        sum(skipmissing(mask(a; with=polygon)))
+    @test zonal(sum, a; of=a) == 
+        zonal(sum, a; of=dims(a)) ==
+        zonal(sum, a; of=Extents.extent(a)) == 
+        sum(a)
+
+    @test zonal(sum, st; of=polygon) == zonal(sum, st; of=[polygon])[1] ==
+        zonal(sum, st; of=(geometry=polygon, x=:a, y=:b)) ==
+        zonal(sum, st; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
+        zonal(sum, st; of=[(geometry=polygon, x=:a, y=:b)])[1] ==
+        maplayers(sum ∘ skipmissing, mask(st; with=polygon))  
+    @test zonal(sum, st; of=st) == 
+        zonal(sum, st; of=dims(st)) == 
+        zonal(sum, st; of=Extents.extent(st)) == 
+        sum(st)
+
+    @testset "skipmissing" begin
+        @test zonal(sum, a_m; of=polygon, skipmissing=false) === missing
+        @test zonal(sum, a_m; of=polygon, skipmissing=true) isa Int
+        @test !zonal(x -> x isa Raster, a_m; of=polygon, skipmissing=true)
+        @test zonal(x -> x isa Raster, a_m; of=polygon, skipmissing=false)
+    end
+end
+
+
+@testset "Lines" begin
+    @test 
+end
+
+using BenchmarkTools
+@btime zonal(sum, a; of=linestring, skipmissing=true)
+@btime zonal(x -> sum(x), a; of=linestring, skipmissing=true) 
+@btime zonal(sum ∘ skipmissing, a; of=linestring, skipmissing=false)
+using ProfileView
+f(a, n) = for _ in 1:n zonal(sum ∘ skipmissing, a; of=linestring, skipmissing=false) end
+g(a, n) = for _ in 1:n zonal(sum ∘ skipmissing, a; of=linestring, skipmissing=true) end
+using Cthulhu
+Cthulhu.@descend zonal(sum, a; of=linestring, skipmissing=true)
+ProfileView.descend_clicked()
+VSCodeServer.@profview 
+f(a, 1000)
+ProfileView.@profview f(a, 100000)
+ProfileView.@profview g(a, 100000)
+zonal(a; of=[linestring, linestring], skipmissing=true) do A
+    @show A
+    sum(A)
+end
+zonal(a; of=[linestring, linestring], skipmissing=false) do A
+    @show collect(skipmissing(A))
+    sum(skipmissing(A))
+end zonal(sum, a; of=linestring, skipmissing=true)
+
+@testset "return missing" begin
+    @test ismissing(zonal(sum, a; of=[polygon, polygon_out_of_bounds, polygon])[2]) &&
+        ismissing(zonal(sum, a; of=[polygon_out_of_bounds, polygon])[1]) &&
+        ismissing(zonal(sum, a; of=(geometry=polygon_out_of_bounds, x=:a, y=:b))) &&
+        ismissing(zonal(sum, a; of=[(geometry=polygon_out_of_bounds, x=:a, y=:b)])[1])
+    @test zonal(sum, a; of=[polygon_out_of_bounds, polygon_out_of_bounds, polygon])[3] == 
+        sum(skipmissing(mask(a; with=polygon)))
+end
\ No newline at end of file