diff --git a/Project.toml b/Project.toml
index 7c6985c79..fd88e8331 100644
--- a/Project.toml
+++ b/Project.toml
@@ -60,7 +60,7 @@ ConstructionBase = "1"
 CoordinateTransformations = "0.6.2"
 DataFrames = "1"
 DimensionalData = "0.29.4"
-DiskArrays = "0.3, 0.4"
+DiskArrays = "0.4"
 Extents = "0.1"
 FillArrays = "0.12, 0.13, 1"
 Flatten = "0.4"
diff --git a/src/methods/burning/allocs.jl b/src/methods/burning/allocs.jl
index 5d916cf76..7956e266f 100644
--- a/src/methods/burning/allocs.jl
+++ b/src/methods/burning/allocs.jl
@@ -11,22 +11,36 @@ function Allocs(buffer)
     return Allocs(buffer, edges, scratch, crossings)
 end
 
-function _burning_allocs(x; nthreads=_nthreads(), threaded=true, kw...) 
-    dims = commondims(x, DEFAULT_POINT_ORDER)
+function _burning_allocs(x; 
+    nthreads=_nthreads(), 
+    threaded=true, 
+    burncheck_metadata=Metadata(),
+    kw...
+) 
     if threaded
-        [Allocs(_init_bools(dims; metadata=Metadata())) for _ in 1:nthreads]
+        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
     else
-        Allocs(_init_bools(dims; metadata=Metadata()))
+        if isnothing(x)
+            Allocs()
+        else
+            dims = commondims(x, DEFAULT_POINT_ORDER)
+            Allocs(_init_bools(dims; metadata=burncheck_metadata))
+        end
     end
 end
 
 _get_alloc(allocs::Vector{<:Allocs}) = _get_alloc(allocs[Threads.threadid()])
 _get_alloc(allocs::Allocs) = allocs
 
-
 # TODO include these in Allocs
 _alloc_burnchecks(n::Int) = fill(false, n)
 _alloc_burnchecks(x::AbstractArray) = _alloc_burnchecks(length(x))
+
 function _set_burnchecks(burnchecks, metadata::Metadata{<:Any,<:Dict}, verbose)
     metadata["missed_geometries"] = .!burnchecks
     verbose && _burncheck_info(burnchecks)
diff --git a/src/methods/burning/array_init.jl b/src/methods/burning/array_init.jl
index 10ddd5937..5f130b992 100644
--- a/src/methods/burning/array_init.jl
+++ b/src/methods/burning/array_init.jl
@@ -1,14 +1,19 @@
 
-# Like `create` but without disk writes, mostly for Bool/Union{Missing,Boo},
+# Like `create` but without disk writes, mostly for Bool or Union{Missing,Bool},
 # and uses `similar` where possible
 # TODO merge this with `create` somehow
 _init_bools(to; kw...) = _init_bools(to, BitArray; kw...)
 _init_bools(to, T::Type; kw...) = _init_bools(to, T, nothing; kw...)
-_init_bools(to::AbstractRasterSeries, T::Type, data; kw...) = _init_bools(first(to), T, data; kw...)
-_init_bools(to::AbstractRasterStack, T::Type, data; kw...) = _init_bools(first(to), T, data; kw...)
-_init_bools(to::AbstractRaster, T::Type, data; kw...) = _init_bools(to, dims(to), T, data; kw...)
-_init_bools(to::Extents.Extent, T::Type, data; kw...) = _init_bools(to, _extent2dims(to; kw...), T, data; kw...)
-_init_bools(to::DimTuple, T::Type, data; kw...) = _init_bools(to, to, T, data; kw...)
+_init_bools(to::AbstractRasterSeries, T::Type, data; kw...) =
+    _init_bools(dims(first(to)), T, data; kw...)
+_init_bools(to::AbstractRasterStack, T::Type, data; kw...) =
+    _init_bools(dims(to), dims(to), T, data; kw...)
+_init_bools(to::AbstractRaster, T::Type, data; kw...) =
+    _init_bools(to, dims(to), T, data; kw...)
+_init_bools(to::Extents.Extent, T::Type, data; kw...) =
+    _init_bools(to, _extent2dims(to; kw...), T, data; kw...)
+_init_bools(to::DimTuple, T::Type, data; kw...) =
+    _init_bools(to, to, T, data; kw...)
 function _init_bools(to::Nothing, T::Type, data; geometrycolumn=nothing,kw...)
     # Get the extent of the geometries
     ext = _extent(data; geometrycolumn, kw...)
@@ -38,7 +43,7 @@ function _alloc_bools(to, dims::DimTuple, ::Type{BitArray}; missingval::Bool=fal
 end
 function _alloc_bools(to, dims::DimTuple, ::Type{<:Array{T}}; missingval=false, metadata=NoMetadata(), kw...) where T
     # Use an `Array`
-    data = fill!(Raster{T}(undef, dims), missingval) 
+    data = fill!(Raster{T}(undef, dims), missingval)
     return rebuild(data; missingval, metadata)
 end
 
@@ -53,7 +58,9 @@ function _prepare_for_burning(B, locus=Center())
 end
 
 # Convert to Array if its not one already
-_lookup_as_array(d::Dimension) = parent(lookup(d)) isa Array ? d : modify(Array, d) 
+_lookup_as_array(x) = setdims(x, _lookup_as_array(dims(x)))
+_lookup_as_array(dims::Tuple) = map(_lookup_as_array, dims)
+_lookup_as_array(d::Dimension) = parent(lookup(d)) isa Array ? d : modify(Array, d)
 
 function _forward_ordered(B)
     reduce(dims(B); init=B) do A, d
diff --git a/src/methods/burning/geometry.jl b/src/methods/burning/geometry.jl
index fa61c1f24..2ee27211b 100644
--- a/src/methods/burning/geometry.jl
+++ b/src/methods/burning/geometry.jl
@@ -22,10 +22,11 @@ function _burn_geometry!(B::AbstractRaster, trait::Nothing, data;
     lock=Threads.SpinLock(), 
     verbose=true, 
     progress=true, 
-    threaded=true,
+    threaded=false,
     fill=true,
-    allocs=_burning_allocs(B; threaded), 
     geometrycolumn=nothing,
+    burncheck_metadata=Metadata(),
+    allocs=_burning_allocs(B; threaded, burncheck_metadata), 
     kw...
 )::Bool
     geoms = _get_geometries(data, geometrycolumn)
@@ -36,8 +37,8 @@ function _burn_geometry!(B::AbstractRaster, trait::Nothing, data;
             geom = _getgeom(geoms, i)
             ismissing(geom) && return nothing
             a = _get_alloc(allocs)
-            B1 = a.buffer
-            burnchecks[i] = _burn_geometry!(B1, geom; fill, allocs=a, lock, kw...)
+            buffer = a.buffer
+            burnchecks[i] = _burn_geometry!(buffer, geom; fill, allocs=a, lock, kw...)
             return nothing
         end
         if fill
@@ -99,14 +100,17 @@ function _burn_geometry!(B::AbstractRaster, ::GI.AbstractGeometryTrait, geom;
             return false
         end
         # Take a view of the geometry extent
-        B1 = view(B, Touches(geomextent))
-        buf1 = _init_bools(B1; missingval=false)
+        V = view(B, Touches(geomextent))
+        buffer = _init_bools(V; missingval=false)
         # Burn the polygon into the buffer
-        allocs = isnothing(allocs) ? Allocs(B) : allocs
-        hasburned = _burn_polygon!(buf1, geom; shape, geomextent, allocs, boundary, kw...)
-        @inbounds for i in eachindex(B1)
-            if buf1[i]
-                B1[i] = fill
+        # We allocate a new bitarray for the view for performance
+        # and always fill with `true`. 
+        allocs = isnothing(allocs) ? Allocs(nothing) : allocs
+        hasburned = _burn_polygon!(buffer, geom; shape, geomextent, allocs, boundary, kw...)
+        # We then transfer burned `true` values to B via the V view
+        @inbounds for i in eachindex(V)
+            if buffer[i]
+                V[i] = fill
             end
         end
     else
diff --git a/src/methods/burning/line.jl b/src/methods/burning/line.jl
index 9155479bb..c5d973dbe 100644
--- a/src/methods/burning/line.jl
+++ b/src/methods/burning/line.jl
@@ -1,33 +1,53 @@
 # _burn_lines!
 # Fill a raster with `fill` where pixels touch lines in a geom
-# Separated for a type stability function barrier
-function _burn_lines!(B::AbstractRaster, geom; fill=true, kw...)
-    _check_intervals(B)
+# Usually `fill` is `true` of `false`
+function _burn_lines!(
+    B::AbstractRaster, geom; fill=true, verbose=false, kw...
+)
+    _check_intervals(B, verbose)
     B1 = _prepare_for_burning(B)
-    return _burn_lines!(B1, geom, fill)
+
+    xdim, ydim = dims(B, DEFAULT_POINT_ORDER)
+    regular = map((xdim, ydim)) do d
+        # @assert (parent(lookup(d)) isa Array)
+        lookup(d) isa AbstractSampled && span(d) isa Regular
+    end
+    msg = """
+        Can only fill lines where dimensions have `Regular` lookups.
+        Consider using `boundary=:center`, reprojecting the crs,
+        or make an issue in Rasters.jl on github if you need this to work.
+        """
+    all(regular) || throw(ArgumentError(msg))
+
+    # Set indices of B as `fill` when a cell is found to burn.
+    _burn_lines!(identity, dims(B), geom) do D
+        @inbounds B[D] = fill
+    end
 end
 
-_burn_lines!(B, geom, fill) =
-    _burn_lines!(B, GI.geomtrait(geom), geom, fill)
-function _burn_lines!(B::AbstractArray, ::Union{GI.MultiLineStringTrait}, geom, fill)
+_burn_lines!(f::F, c::C, dims::Tuple, geom) where {F<:Function,C<:Function} =
+    _burn_lines!(f, c, dims, GI.geomtrait(geom), geom)
+function _burn_lines!(
+    f::F, c::C, dims::Tuple, ::Union{GI.MultiLineStringTrait}, geom
+) where {F<:Function,C<:Function}
     n_on_line = 0
     for linestring in GI.getlinestring(geom)
-        n_on_line += _burn_lines!(B, linestring, fill)
+        n_on_line += _burn_lines!(f, c, dims, linestring)
     end
     return n_on_line
 end
 function _burn_lines!(
-    B::AbstractArray, ::Union{GI.MultiPolygonTrait,GI.PolygonTrait}, geom, fill
-)
+    f::F, c::C, dims::Tuple, ::Union{GI.MultiPolygonTrait,GI.PolygonTrait}, geom
+) where {F<:Function,C<:Function}
     n_on_line = 0
     for ring in GI.getring(geom)
-        n_on_line += _burn_lines!(B, ring, fill)
+        n_on_line += _burn_lines!(f, c, dims, ring)
     end
     return n_on_line
 end
 function _burn_lines!(
-    B::AbstractArray, ::GI.AbstractCurveTrait, linestring, fill
-)
+    f::F, c::C, dims::Tuple, ::GI.AbstractCurveTrait, linestring
+) where {F<:Function,C<:Function}
     isfirst = true
     local firstpoint, laststop
     n_on_line = 0
@@ -46,19 +66,19 @@ function _burn_lines!(
             stop=(x=GI.x(point), y=GI.y(point)),
         )
         laststop = line.stop
-        n_on_line += _burn_line!(B, line, fill)
+        n_on_line += _burn_line!(f, c, dims, line)
     end
     return n_on_line
 end
 function _burn_lines!(
-    B::AbstractArray, t::GI.LineTrait, line, fill
-)
+    f::F, c::C, dims::Tuple, t::GI.LineTrait, line
+) where {F<:Function,C<:Function}
     p1, p2 = GI.getpoint(t, line)
     line1 = (
         start=(x=GI.x(p1), y=GI.y(p1)),
         stop=(x=GI.x(p2), y=GI.y(p2)),
     )
-    return _burn_line!(B, line1, fill)
+    return _burn_line!(f, c, dims, line1)
 end
 
 # _burn_line!
@@ -66,26 +86,23 @@ end
 # Line-burning algorithm
 # Burns a single line into a raster with value where pixels touch a line
 #
+# Function `f` does the actual work when passed a Dimension Tuple of a pixel to burn, 
+# and `c` is an initialisation callback that is passed the maximyum 
+# number of times `f` will be called. It may be called less than that.
+#
 # TODO: generalise to Irregular spans?
-function _burn_line!(A::AbstractRaster, line, fill)
+function _burn_line!(f::Function, c::Function, dims::Tuple, line::NamedTuple)
+    xdim, ydim = dims
+    di = DimIndices(dims)
 
-    xdim, ydim = dims(A, DEFAULT_POINT_ORDER)
-    regular = map((xdim, ydim)) do d
-        @assert (parent(lookup(d)) isa Array)
-        lookup(d) isa AbstractSampled && span(d) isa Regular
-    end
-    msg = """
-        Can only fill lines where dimensions have `Regular` lookups.
-        Consider using `boundary=:center`, reprojecting the crs,
-        or make an issue in Rasters.jl on github if you need this to work.
-        """
-    all(regular) || throw(ArgumentError(msg))
+    @assert xdim isa XDim
+    @assert ydim isa YDim
 
     @assert order(xdim) == order(ydim) == Lookups.ForwardOrdered()
     @assert locus(xdim) == locus(ydim) == Lookups.Center()
 
-    raster_x_step = abs(step(span(A, X)))
-    raster_y_step = abs(step(span(A, Y)))
+    raster_x_step = abs(step(span(xdim)))
+    raster_y_step = abs(step(span(ydim)))
     raster_x_offset = @inbounds xdim[1] - raster_x_step / 2 # Shift from center to start of pixel
     raster_y_offset = @inbounds ydim[1] - raster_y_step / 2
 
@@ -117,15 +134,18 @@ function _burn_line!(A::AbstractRaster, line, fill)
     # Int starting points for the line. +1 converts to julia indexing
     j, i = trunc(Int, relstart.x) + 1, trunc(Int, relstart.y) + 1 # Int
 
-    # For arbitrary dimension indexing
-    dimconstructors = map(DD.basetypeof, (xdim, ydim))
+    n_on_line = 0
+    
+    # Pass of number of runs of `f` to callback `c`
+    # This can help with e.g. allocating a vector
+    c(manhattan_distance + 1)
 
     if manhattan_distance == 0
-        D = map((d, o) -> d(o), dimconstructors, (j, i))
-        if checkbounds(Bool, A, D...)
-            @inbounds A[D...] = fill
+        D = map(rebuild, dims, (j, i))
+        if checkbounds(Bool, di, D...)
+            f(D)
+            n_on_line += 1
         end
-        n_on_line = 1
         return n_on_line
     end
 
@@ -153,18 +173,17 @@ function _burn_line!(A::AbstractRaster, line, fill)
     n_on_line = 0
     countx = county = 0
 
-
     # Int steps to move allong the line
     step_j = signbit(diff_x) * -2 + 1
     step_i = signbit(diff_y) * -2 + 1
 
     # Travel one grid cell at a time. Start at zero for the current cell
     for _ in 0:manhattan_distance
-        D = map((d, o) -> d(o), dimconstructors, (j, i))
-        if checkbounds(Bool, A, D...)
-            @inbounds A[D...] = fill
+        D = map(rebuild, dims, (j, i))
+        if checkbounds(Bool, di, D...)
+            f(D)
+            n_on_line += 1
         end
-
         # Only move in either X or Y coordinates, not both.
         if abs(max_x) < abs(max_y)
             max_x += delta_x
@@ -176,12 +195,6 @@ function _burn_line!(A::AbstractRaster, line, fill)
             county +=1
         end
     end
-    return n_on_line
-end
-function _burn_line!(A::AbstractRaster, line, fill, order::Tuple{Vararg{Dimension}})
-    msg = """"
-        Converting a `:line` geometry to raster is currently only implemented for 2d lines.
-        Make a Rasters.jl github issue if you need this for more dimensions.
-        """
-    throw(ArgumentError(msg))
+
+    return n_on_line 
 end
diff --git a/src/methods/burning/polygon.jl b/src/methods/burning/polygon.jl
index f6d4512e4..1ce6c48a9 100644
--- a/src/methods/burning/polygon.jl
+++ b/src/methods/burning/polygon.jl
@@ -24,14 +24,14 @@ function _burn_polygon!(B::AbstractDimArray, trait, geom;
     # Lines
     n_on_line = 0
     if boundary !== :center
-        _check_intervals(B, boundary)
+        _check_intervals(B, boundary, verbose)
         if boundary === :touches 
-            if _check_intervals(B, boundary)
+            if _check_intervals(B, boundary, verbose)
                 # Add line pixels
                 n_on_line = _burn_lines!(B, geom; fill)::Int
             end
         elseif boundary === :inside 
-            if _check_intervals(B, boundary)
+            if _check_intervals(B, boundary, verbose)
                 # Remove line pixels
                 n_on_line = _burn_lines!(B, geom; fill=!fill)::Int
             end
@@ -135,9 +135,9 @@ end
 
 const INTERVALS_INFO = "makes more sense on `Intervals` than `Points` and will have more correct results. You can construct dimensions with a `X(values; sampling=Intervals(Center()))` to achieve this"
 
-@noinline _check_intervals(B) = 
-    _chki(B) ? true : (@info "burning lines $INTERVALS_INFO"; false)
-@noinline _check_intervals(B, boundary) =
-    _chki(B) ? true : (@info "`boundary=:$boundary` $INTERVALS_INFO"; false)
+@noinline _check_intervals(B, verbose::Bool) = 
+    _chki(B) ? true : (verbose && @info "burning lines $INTERVALS_INFO"; false)
+@noinline _check_intervals(B, boundary, verbose::Bool) =
+    _chki(B) ? true : (verbose && @info "`boundary=:$boundary` $INTERVALS_INFO"; false)
 
 _chki(B) = all(map(s -> s isa Intervals, sampling(dims(B)))) 
diff --git a/src/methods/extract.jl b/src/methods/extract.jl
index ef41b04e7..ca2ba6752 100644
--- a/src/methods/extract.jl
+++ b/src/methods/extract.jl
@@ -1,9 +1,85 @@
+# Object to hold extract keywords
+struct Extractor{T,P,K,N<:NamedTuple{K},G,I,S,F}
+    names::N
+    geometry::G
+    index::I
+    skipmissing::S
+    flatten::F
+end
+Base.@constprop :aggressive @inline function Extractor(x, data;
+    name::NTuple{<:Any,Symbol},
+    skipmissing,
+    flatten,
+    geometry,
+    index,
+    kw...
+)
+    nt = NamedTuple{name}(name)
+    g = _booltype(geometry)
+    i = _booltype(index)
+    sm = _booltype(skipmissing)
+    f = _booltype(flatten)
+    Extractor(x, data, nt, g, i, sm, f)
+end
+function Extractor(x, data, nt::N, g::G, i::I, sm::S, f::F) where {N<:NamedTuple{K},G,I,S,F} where K
+    P = _proptype(x; skipmissing=sm, names=nt)
+    T = _geom_rowtype(x, data; geometry=g, index=i, names=nt, skipmissing=sm)
+    Extractor{T,P,K,N,G,I,S,F}(nt, g, i, sm, f)
+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)
+    lr.i = 1
+    lr.prev = CartesianIndex((typemin(Int), typemin(Int)))
+end
+
+_geom_rowtype(A, geom; kw...) =
+    _geom_rowtype(A, GI.trait(geom), geom; kw...)
+_geom_rowtype(A, ::Nothing, geoms; kw...) =
+    _geom_rowtype(A, first(geoms); kw...)
+_geom_rowtype(A, ::GI.AbstractGeometryTrait, geom; kw...) =
+    _geom_rowtype(A, first(GI.getpoint(geom)); kw...)
+function _geom_rowtype(A, ::GI.AbstractPointTrait, p; kw...)
+    tuplepoint = if GI.is3d(p)
+        (GI.x(p), GI.y(p), GI.z(p))
+    else
+        (GI.x(p), GI.y(p))
+    end
+    return _rowtype(A, tuplepoint; kw...)
+end
+
+# skipinvalid: can G and I be missing. skipmissing: can nametypes be missing
+function _proptype(x;
+    skipmissing, names::NamedTuple{K}, kw...
+) where K
+    NamedTuple{K,Tuple{_nametypes(x, names, skipmissing)...}}
+end
+
 """
-    extract(x, data; kw...)
+    extract(x, geometries; kw...)
+
+Extracts the value of `Raster` or `RasterStack` for the passed in geometries,
+returning an `Vector{NamedTuple}` with properties for `:geometry` and `Raster`
+or `RasterStack` layer values.
 
-Extracts the value of `Raster` or `RasterStack` at given points, returning
-an iterable of `NamedTuple` with properties for `:geometry` and raster or
-stack layer values.
+For lines, linestrings and linear rings points are extracted for each pixel that
+the line touches.
+
+For polygons, all cells witih centers covered by the polygon are returned.
 
 Note that if objects have more dimensions than the length of the point tuples,
 sliced arrays or stacks will be returned instead of single values.
@@ -15,13 +91,22 @@ $DATA_ARGUMENT
 
 # Keywords
 
-- `geometry`: include `:geometry` in returned `NamedTuple`, `true` by default.
-- `index`: include `:index` of the `CartesianIndex` in returned `NamedTuple`, `false` by default.
-- `name`: a `Symbol` or `Tuple` of `Symbol` corresponding to layer/s of a `RasterStack` to extract. All layers by default.
+- `geometry`: include a `:geometry` field in rows, which will be a
+    tuple point. Either the original point for points or the pixel
+    center point for line and polygon extract. `true` by default.
+- `index`: include `:index` field of extracted points in rows, `false` by default.
+- `name`: a `Symbol` or `Tuple` of `Symbol` corresponding to layer/s
+    of a `RasterStack` to extract. All layers are extracted by default.
 - `skipmissing`: skip missing points automatically.
+- `flatten`: flatten extracted points from multiple 
+    geometries into a single vector. `true` by default. 
+    Unmixed point geometries are always flattened.
+    Flattening is slow and single threaded, `flatten=false` may be a 
+    large performance improvement in combination with `threaded=true`.
 - `atol`: a tolerance for floating point lookup values for when the `Lookup`
     contains `Points`. `atol` is ignored for `Intervals`.
--$GEOMETRYCOLUMN_KEYWORD
+$BOUNDARY_KEYWORD
+$GEOMETRYCOLUMN_KEYWORD
 
 # Example
 
@@ -50,206 +135,349 @@ extract(st, obs; skipmissing=true)
  (geometry = (0.32, 40.24), bio1 = 16.321388f0, bio3 = 41.659454f0, bio5 = 30.029825f0, bio7 = 25.544561f0, bio12 = 480.0f0)
 ```
 
-Note: passing in arrays, geometry collections or feature collections 
+Note: passing in arrays, geometry collections or feature collections
 containing a mix of points and other geometries has undefined results.
 """
 function extract end
-@inline function extract(x::RasterStackOrArray, data;
-   names=_names(x), name=names, skipmissing=false, geometry=true, index=false, geometrycolumn=nothing, kw...
+Base.@constprop :aggressive @inline function extract(x::RasterStackOrArray, data;
+    geometrycolumn=nothing,
+    names::NTuple{<:Any,Symbol}=_names(x),
+    name::NTuple{<:Any,Symbol}=names,
+    skipmissing=false,
+    flatten=true,
+    geometry=true,
+    index=false,
+    kw...
 )
     n = DD._astuple(name)
-    _extract(x, data;
-         dims=DD.dims(x, DEFAULT_POINT_ORDER),
-         names=NamedTuple{n}(n),
-         # These keywords are converted to _True/_False for type stability later on
-         # The @inline above helps constant propagation of the Bools
-         geometry=_booltype(geometry), 
-         index=_booltype(index), 
-         skipmissing=_booltype(skipmissing), 
-         geometrycolumn,
-         kw...
-    )
+    g, g1 = if GI.isgeometry(data)
+        data, data
+    elseif GI.isfeature(data)
+        g = GI.geometry(data)
+        g, g
+    else
+        gs = _get_geometries(data, geometrycolumn)
+        gs, first(Base.skipmissing(gs))
+    end
+    open(x) do xo
+        xp = _prepare_for_burning(xo)
+        e = Extractor(xp, g1; name, skipmissing, flatten, geometry, index)
+        return _extract(xp, g, e; kw...)
+    end
 end
 
-function _extract(A::RasterStackOrArray, geom::Missing, names, kw...)
-    T = _rowtype(A, geom; names, kw...)
-    [_maybe_add_fields(T, map(_ -> missing, names), missing, missing)]
+# TODO use a GeometryOpsCore method like `applyreduce` here?
+function _extract(A::RasterStackOrArray, geom::Missing, e; kw...)
+    return if istrue(e.skipmissing)
+        T[]
+    else
+        T[_maybe_add_fields(e, map(_ -> missing, e.names), missing, missing)]
+    end
 end
-function _extract(A::RasterStackOrArray, geom; names, kw...)
-    _extract(A, GI.geomtrait(geom), geom; names, kw...)
+function _extract(A::RasterStackOrArray, geom, e; kw...)
+    _extract(A, GI.geomtrait(geom), geom, e; kw...)
 end
-function _extract(A::RasterStackOrArray, ::Nothing, data; 
-    names, skipmissing, geometrycolumn, kw...
-)
-    geoms = _get_geometries(data, geometrycolumn)
-    T = _rowtype(A, eltype(geoms); names, skipmissing, kw...)
+function _extract(A::RasterStackOrArray, ::Nothing, geoms::AbstractArray, e::Extractor{T};
+    threaded=false, progress=true, kw...
+) where T
     # Handle empty / all missing cases
     (length(geoms) > 0 && any(!ismissing, geoms)) || return T[]
-    
-    geom1 = first(Base.skipmissing(geoms))
+
+    geom1 = first(skipmissing(geoms))
     trait1 = GI.trait(geom1)
-    # We need to split out points from other geoms
-    # TODO this will fail with mixed point/geom vectors
+    # We split out points from other geoms for performance
     if trait1 isa GI.PointTrait
-    rows = Vector{T}(undef, length(geoms))
-        if istrue(skipmissing)
-            j = 1
-            for i in eachindex(geoms)
-                g = geoms[i]
-                ismissing(g) && continue
-                e = _extract_point(T, A, g, skipmissing; names, kw...)
-                if !ismissing(e) 
-                    rows[j] = e
-                    j += 1
+        allpoints = true
+        i = 1
+        rows = Vector{T}(undef, length(geoms))
+        if istrue(e.skipmissing)
+            if threaded
+                nonmissing = Vector{Bool}(undef, length(geoms))
+                _run(1:length(geoms), threaded, progress, "Extracting points...") do i
+                    g = geoms[i]
+                    geomtrait(g) isa GI.PointTrait || return nothing
+                    nonmissing[i] = _extract_point!(rows, A, g, e, i; kw...)::T
+                    return nothing
                 end
-                nothing
+                # This could use less memory if we reuse `rows` and shorten it
+                rows = rows[nonmissing]
+            else
+                j = Ref(1)
+                # For non-threaded be more memory efficient
+                _run(1:length(geoms), threaded, progress, "Extracting points...") do _
+                    g = geoms[j[]]
+                    geomtrait(g) isa GI.PointTrait || return nothing
+                    j[] += _extract_point!(rows, A, g, e, i; kw...)::T
+                    return nothing
+                end
+                deleteat!(rows, j[]:length(rows))
             end
-            deleteat!(rows, j:length(rows))
         else
-            for i in eachindex(geoms)
+            _run(1:length(geoms), threaded, progress, "Extracting points...") do i
                 g = geoms[i]
-                rows[i] = _extract_point(T, A, g, skipmissing; names, kw...)::T
-                nothing
+                geomtrait(g) isa GI.PointTrait || return nothing
+                _extract_point!(rows, A, g, e, i; kw...)::T
+                return nothing
             end
         end
-        return rows
+        # If we found a non-point geometry,
+        # ignore these rows and start again generically
+        allpoints && return rows
+    end
+
+    
+        row_vecs = Vector{Vector{T}}(undef, length(geoms))
+    if threaded
+        thread_line_refs = [LineRefs{T}() for _ in 1:Threads.nthreads()]
+        _run(1:length(geoms), threaded, progress, "Extracting geometries...") do i
+            line_refs = thread_line_refs[Threads.threadid()]
+            row_vecs[i] = _extract(A, geoms[i], e; line_refs, kw...)
+        end
     else
-        # This will be a list of vectors, we need to flatten it into one
-        rows = Iterators.flatten(_extract(A, g; names, skipmissing, kw...) for g in geoms)
-        if istrue(skipmissing) 
-            return collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) 
-        else
-            return collect(rows)
+        line_refs = LineRefs{T}()
+        _run(1:length(geoms), threaded, progress, "Extracting geometries...") do i
+            row_vecs[i] = _extract(A, geoms[i], e; line_refs, kw...)
         end
     end
+
+    # TODO this is a bit slow and only on one thread
+    if istrue(e.flatten)
+        n = sum(map(length, row_vecs))
+        out_rows = Vector{T}(undef, n)
+        k::Int = 1
+        for row_vec in row_vecs
+            for row in row_vec
+                out_rows[k] = row
+                k += 1
+            end
+        end
+        return out_rows
+    else
+        return row_vecs
+    end
 end
-function _extract(A::RasterStackOrArray, ::GI.AbstractFeatureTrait, feature; kw...)
-    _extract(A, GI.geometry(feature); kw...)
-end
-function _extract(A::RasterStackOrArray, ::GI.FeatureCollectionTrait, fc; kw...)
-    # Fall back to the Array/iterator method for feature collections
-    _extract(A, [GI.geometry(f) for f in GI.getfeature(fc)]; kw...)
-end
-function _extract(A::RasterStackOrArray, ::GI.AbstractMultiPointTrait, geom; 
-    skipmissing, kw...
-)
-    T = _rowtype(A, GI.getpoint(geom, 1); names, skipmissing, kw...)
-    rows = (_extract_point(T, A, p, skipmissing; kw...) for p in GI.getpoint(geom))
-    return skipmissing isa _True ? collect(_skip_missing_rows(rows, _missingval_or_missing(A), names)) : collect(rows)
+function _extract(
+    A::RasterStackOrArray, ::GI.AbstractMultiPointTrait, geom, e::Extractor{T}; kw...
+)::Vector{T} where T
+    n = GI.npoint(geom)
+    rows = _init_rows(e, n)
+    for p in GI.getpoint(geom)
+        i += _extract_point!(rows, A, e, p, i; kw...)
+    end
+    if istrue(skipmissing)
+        # Remove excees rows where missing
+        deleteat!(rows, i:length(rows))
+    end
+    return rows
 end
-function _extract(A::RasterStackOrArray, ::GI.PointTrait, geom; 
-    skipmissing, kw...
-)
-    T = _rowtype(A, geom; names, skipmissing, kw...)
-    _extract_point(T, A, geom, skipmissing; kw...)
+function _extract(A::RasterStackOrArray, ::GI.PointTrait, p, e::Extractor{T}; kw...) where T
+    rows = _init_rows(e, 1)
+    _extract_point!(rows, A, e, p, 1; kw...)
 end
-function _extract(A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom;
-    names, skipmissing, kw...
-)
-    # Make a raster mask of the geometry
-    template = view(A, Touches(GI.extent(geom)))
-    ods = otherdims(A, DEFAULT_POINT_ORDER)
-    if length(ods) > 0
-        template = view(template, map(d -> rebuild(d, firstindex(d)), ods)) 
+@noinline function _extract(
+    A::RasterStackOrArray, trait::GI.AbstractLineStringTrait, geom, e::Extractor{T};
+    line_refs=LineRefs{T}(),
+    kw...
+)::Vector{T} where T
+    _initialise!(line_refs)
+    # Subset/offst is not really needed for line buring
+    # But without it we can get different fp errors
+    # to polygons and eng up with lines in different
+    # places when they are right on the line.
+    dims, offset = _template(A, geom)
+    dp = DimPoints(A)
+    function _length_callback(n)
+        rows = line_refs.rows
+        resize!(rows, n + line_refs.i - 1)
     end
-    p = first(GI.getpoint(geom))
-    tuplepoint = if GI.is3d(geom) 
-        GI.x(p), GI.y(p), GI.z(p)
-    else
-        GI.x(p), GI.y(p)
+
+    _burn_lines!(_length_callback, dims, geom) do D
+        I = CartesianIndex(map(val, D))
+        # Avoid duplicates from adjacent line segments
+        line_refs.prev == I && return nothing
+        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, A, dp, offset, I, line_refs.i)
+        return nothing
     end
-    T = _rowtype(A, tuplepoint; names, skipmissing, kw...)
-    B = boolmask(geom; to=template, kw...)
-    offset = CartesianIndex(map(x -> first(x) - 1, parentindices(parent(template))))
+    rows = line_refs.rows
+    deleteat!(rows, line_refs.i:length(rows))
+    return rows
+end
+@noinline function _extract(
+    A::RasterStackOrArray, t::GI.AbstractGeometryTrait, geom, e::Extractor{T}; kw...
+)::Vector{T} where T
+    # Make a raster mask of the geometry
+    dims, offset = _template(A, geom)
+    B = boolmask(geom; to=dims, burncheck_metadata=NoMetadata(), kw...)
+    n = count(B)
+    dp = DimPoints(A)
+    i = 1
+    rows = _init_rows(e, n)
     # Add a row for each pixel that is `true` in the mask
-    rows = (_missing_or_fields(T, A, Tuple(I + offset), names, skipmissing) for I in CartesianIndices(B) if B[I])
-    # Maybe skip missing rows
-    if istrue(skipmissing)
-        return collect(Base.skipmissing(rows)) 
+    for I in CartesianIndices(B)
+        B[I] || continue
+        i += _maybe_set_row!(rows, e.skipmissing, e, A, dp, offset, I, i)
+    end
+    # Cleanup
+    deleteat!(rows, i:length(rows))
+    return rows
+end
+
+function _template(x, geom)
+    ods = otherdims(x, DEFAULT_POINT_ORDER)
+    # Build a dummy raster in case this is a stack
+    # views are just easier on an array than dims
+    t1 = Raster(FillArrays.Zeros(size(x)), dims(x))
+    t2 = if length(ods) > 0
+        view(t1, map(d -> rebuild(d, firstindex(d)), ods))
     else
-        return collect(rows)
+        t1
     end
+    I = dims2indices(dims(t2), Touches(GI.extent(geom)))
+    t3 = view(t2, I...)
+    offset = CartesianIndex(map(i -> first(i) - 1, I))
+    return dims(t3), offset
 end
-_extract(::Type{T}, A::RasterStackOrArray, trait::GI.PointTrait, point; skipmissing, kw...) where T =
-    _extract_point(T, A, point, skipmissing; kw...)
 
-function _missing_or_fields(::Type{T}, A, I, names, skipmissing) where T
-    props = _prop_nt(A, I, names)
-    if istrue(skipmissing) && _ismissingval(A, props) 
-        missing
+Base.@assume_effects :foldable function _maybe_set_row!(
+    rows, skipmissing::_True, e, A, dp, offset, I, i;
+)
+    props = _prop_nt(e, A, I, skipmissing)
+    return if ismissing(props)
+        0
     else
-        geom = DimPoints(A)[I...]
-        _maybe_add_fields(T, props, geom, I)
+        _maybe_set_row!(rows, _False(), e, A, dp, offset, I, i; props)
     end
 end
+Base.@assume_effects :foldable function _maybe_set_row!(
+    rows::Vector{T}, skipmissing::_False, e, A, dp, offset, I, i;
+    props=_prop_nt(e, A, I, skipmissing)
+) where T
+    Ioff = I + offset
+    geom = dp[Ioff]
+    rows[i] = _maybe_add_fields(T, props, geom, Tuple(Ioff))
+    return 1
+end
 
-_ismissingval(A::Union{Raster,RasterStack}, props) = 
+_init_rows(e::Extractor{T}, n=0) where T = Vector{T}(undef, n)
+
+Base.@assume_effects :foldable _ismissingval(A::Union{Raster,RasterStack}, props)::Bool =
     _ismissingval(missingval(A), props)
-_ismissingval(A::Union{Raster,RasterStack}, props::NamedTuple) = 
+Base.@assume_effects :foldable _ismissingval(A::Union{Raster,RasterStack}, props::NamedTuple)::Bool =
     _ismissingval(missingval(A), props)
-_ismissingval(mvs::NamedTuple, props::NamedTuple{K}) where K = 
-    any(k -> ismissing(props[k]) || props[k] === mvs[k], K)
-_ismissingval(mv, props::NamedTuple) = any(x -> ismissing(x) || x === mv, props)
-_ismissingval(mv, prop) = (mv === prop)
+Base.@assume_effects :foldable _ismissingval(mvs::NamedTuple, props::NamedTuple)::Bool =
+    any(DD.unrolled_map((p, mv) -> ismissing(p) || p === mv, props, mvs))
+Base.@assume_effects :foldable _ismissingval(mv, props::NamedTuple)::Bool =
+    any(DD.unrolled_map(p -> ismissing(p) || p === mv, props))
+Base.@assume_effects :foldable _ismissingval(mv, prop)::Bool = (mv === prop)
 
-@inline _prop_nt(st::AbstractRasterStack, I, ::NamedTuple{K}) where K = st[I...][K]
-@inline _prop_nt(A::AbstractRaster, I, ::NamedTuple{K}) where K = NamedTuple{K}((A[I...],))
+# We always return NamedTuple
+Base.@assume_effects :foldable function _prop_nt(
+    ::Extractor{<:Any,P,K}, st::AbstractRasterStack, I, sm::_False
+)::P where {P,K}
+    P(st[K][I])
+end
+Base.@assume_effects :foldable function _prop_nt(
+    ::Extractor{<:Any,P,K}, st::AbstractRasterStack{K}, I, sm::_False
+)::P where {P,K}
+    P(st[I])
+end
+Base.@assume_effects :foldable function _prop_nt(
+    ::Extractor{<:Any,P,K}, A::AbstractRaster, I, sm::_False
+)::P where {P,K}
+    P(NamedTuple{K}((A[I],)))
+end
+Base.@assume_effects :foldable function _prop_nt(
+    ::Extractor{<:Any,P,K}, st::AbstractRasterStack, I, sm::_True
+)::Union{P,Missing} where {P,K}
+    x = st[K][I]
+    _ismissingval(A, x) ? missing : x::P
+end
+Base.@assume_effects :foldable function _prop_nt(
+    ::Extractor{<:Any,P,K}, st::AbstractRasterStack{K}, I, sm::_True
+)::Union{P,Missing} where {P,K}
+    x = st[I]
+    _ismissingval(st, x) ? missing : x::P
+end
+Base.@assume_effects :foldable function _prop_nt(
+    ::Extractor{<:Any,P,K}, A::AbstractRaster, I, sm::_True
+)::Union{P,Missing} where {P,K}
+    x = A[I]
+    _ismissingval(A, x) ? missing : NamedTuple{K}((x,))::P
+end
 
 # Extract a single point
-@inline function _extract_point(::Type, x::RasterStackOrArray, point::Missing, skipmissing::_True; kw...)
-    return missing
+# Missing point to remove
+@inline function _extract_point!(
+    rows::Vector{T}, x::RasterStackOrArray, skipmissing::_True, point::Missing, i;
+    kw...
+) where T
+    return false
 end
-@inline function _extract_point(::Type{T}, x::RasterStackOrArray, point::Missing, skipmissing::_False;
+# Missing point to keep
+@inline function _extract_point!(
+    rows::Vector{T}, x::RasterStackOrArray, skipmissing::_False, point::Missing, i;
     names, kw...
 ) where T
     props = map(_ -> missing, names)
     geom = missing
     I = missing
-    return _maybe_add_fields(T, props, geom, I)
+    rows[i] = _maybe_add_fields(T, props, geom, I)
+    return true
 end
-@inline function _extract_point(::Type{T}, x::RasterStackOrArray, point, skipmissing::_True;
-    dims, names::NamedTuple{K}, atol=nothing, kw...
+# Normal point with missing / out of bounds data removed
+@inline function _extract_point!(
+    rows::Vector{T}, x::RasterStackOrArray, skipmissing::_True, point, i;
+    dims=DD.dims(x, DEFAULT_POINT_ORDER),
+    names::NamedTuple{K},
+    atol=nothing,
+    kw...
 ) where {T,K}
     # Get the actual dimensions available in the object
     coords = map(DD.commondims(x, dims)) do d
         _dimcoord(d, point)
     end
     # If any are coordinates missing, also return missing for everything
-    if any(map(ismissing, coords))
-        return missing
-    else
+    if !any(map(ismissing, coords))
         selector_dims = map(dims, coords) do d, c
             _at_or_contains(d, c, atol)
         end
         selectors = map(val, DD.dims(selector_dims, DD.dims(x)))
         # Check the selector is in bounds / actually selectable
         I = DD.selectindices(DD.dims(x, dims), selectors; err=_False())::Union{Nothing,Tuple{Int,Int}}
-        if isnothing(I)
-            return missing
-        else
+        if !isnothing(I)
             D = map(rebuild, selector_dims, I)
-            if x isa Raster 
+            if x isa Raster
                 prop = x[D]
-                _ismissingval(missingval(x), prop) && return missing
+                _ismissingval(missingval(x), prop) && return false
                 props = NamedTuple{K,Tuple{eltype(x)}}((prop,))
             else
                 props = x[D][K]
-                _ismissingval(missingval(x), props) && return missing
+                _ismissingval(missingval(x), props) && return false
             end
-            return _maybe_add_fields(T, props, point, I)
+            rows[i] = _maybe_add_fields(T, props, point, I)
+            return true
         end
     end
+    return false
 end
-@inline function _extract_point(::Type{T}, x::RasterStackOrArray, point, skipmissing::_False;
-    dims, names::NamedTuple{K}, atol=nothing, kw...
+# Normal point with missing / out of bounds data kept with `missing` fields
+@inline function _extract_point!(
+    row::Vector{T}, x::RasterStackOrArray, skipmissing::_False, point, i;
+    dims=DD.dims(x, DEFAULT_POINT_ORDER),
+    names::NamedTuple{K},
+    atol=nothing,
+    kw...
 ) where {T,K}
     # Get the actual dimensions available in the object
     coords = map(DD.commondims(x, dims)) do d
         _dimcoord(d, point)
     end
     # If any are coordinates missing, also return missing for everything
-    if any(map(ismissing, coords))
-        return _maybe_add_fields(T, map(_ -> missing, names), missing, missing)
+    rows[i] = if any(map(ismissing, coords))
+        _maybe_add_fields(T, map(_ -> missing, names), missing, missing)
     else
         selector_dims = map(dims, coords) do d, c
             _at_or_contains(d, c, atol)
@@ -258,40 +486,28 @@ end
         # Check the selector is in bounds / actually selectable
         I = DD.selectindices(DD.dims(x, dims), selectors; err=_False())::Union{Nothing,Tuple{Int,Int}}
         if isnothing(I)
-            return _maybe_add_fields(T, map(_ -> missing, names), point, missing)
+            _maybe_add_fields(T, map(_ -> missing, names), point, missing)
         else
             D = map(rebuild, selector_dims, I)
-            props = if x isa Raster 
+            props = if x isa Raster
                 NamedTuple{K,Tuple{eltype(x)}}((x[D],))
             else
                 NamedTuple(x[D])[K]
             end
-            return _maybe_add_fields(T, props, point, I)
+            _maybe_add_fields(T, props, point, I)
         end
     end
+    return 1
 end
 
 # Maybe add optional fields
-Base.@assume_effects :total function _maybe_add_fields(::Type{T}, props::NamedTuple, point, I)::T where {T<:NamedTuple{K}} where K
+# It is critically important for performance that this is type stable
+Base.@assume_effects :total function _maybe_add_fields(
+    ::Type{T}, props::NamedTuple, point, I
+)::T where {T<:NamedTuple{K}} where K
     if :geometry in K
         :index in K ? merge((; geometry=point, index=I), props) : merge((; geometry=point), props)
     else
         :index in K ? merge((; index=I), props) : props
     end |> T
 end
-
-@inline _skip_missing_rows(rows, ::Missing, names) = 
-    Iterators.filter(row -> !any(ismissing, row), rows)
-@inline _skip_missing_rows(rows, missingval, names) = 
-    Iterators.filter(row -> !any(x -> ismissing(x) || x === missingval, row), rows)
-@inline function _skip_missing_rows(rows, missingval::NamedTuple, names::NamedTuple{K}) where K
-    # first check if all fields are equal - if so just call with the first value
-    if Base.allequal(missingval) == 1
-        return _skip_missing_rows(rows, first(missingval), names)
-    else
-        Iterators.filter(rows) do row
-            # rows may or may not contain a :geometry field, so map over keys instead
-            !any(key -> ismissing(row[key]) || row[key] === missingval[key], K)
-        end
-    end
-end
diff --git a/src/methods/mask.jl b/src/methods/mask.jl
index 1b19b44cb..8c8f5901c 100644
--- a/src/methods/mask.jl
+++ b/src/methods/mask.jl
@@ -28,7 +28,8 @@ $SUFFIX_KEYWORD
 
 These can be used when `with` is a GeoInterface.jl compatible object:
 
-$SHAPE_KEYWORDS
+$SHAPE_KEYWORD
+$BOUNDARY_KEYWORD
 $GEOMETRYCOLUMN_KEYWORD
 
 # Example
@@ -336,7 +337,8 @@ function boolmask!(dest::AbstractRaster, data;
     lock=nothing, 
     progress=true, 
     threaded=false, 
-    allocs=_burning_allocs(dest; threaded), 
+    burncheck_metadata=Metadata(),
+    allocs=_burning_allocs(dest; threaded, burncheck_metadata), 
     geometrycolumn=nothing,
     kw...
 )
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index db0194165..a2cd328d8 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -312,15 +312,10 @@ const RASTERIZE_KEYWORDS = """
 - `op`: A reducing function that accepts two values and returns one, like `min` to `minimum`.
     For common methods this will be assigned for you, or is not required. But you can use it
     instead of a `reducer` as it will usually be faster.
-- `shape`: force `data` to be treated as `:polygon`, `:line` or `:point`, where possible
-    Points can't be treated as lines or polygons, and lines may not work as polygons, but
-    an attempt will be made.
-- `geometrycolumn`: `Symbol` to manually select the column the geometries are in
-    when `data` is a Tables.jl compatible table, or a tuple of `Symbol` for columns of
-    point coordinates.
-- `progress`: show a progress bar, `true` by default, `false` to hide..
-- `verbose`: print information and warnings when there are problems with the rasterisation.
-    `true` by default.
+$GEOM_KEYWORDS
+$GEOMETRYCOLUMN_KEYWORD
+$PROGRESS_KEYWORD
+$VERBOSE_KEYWORD
 $THREADED_KEYWORD
 - `threadsafe`: specify that custom `reducer` and/or `op` functions are thread-safe, 
     in that the order of operation or blocking does not matter. For example, 
@@ -356,8 +351,6 @@ $RASTERIZE_ARGUMENTS
 
 These are detected automatically from `data` where possible.
 
-$GEOMETRYCOLUMN_KEYWORD
-$GEOM_KEYWORDS
 $RASTERIZE_KEYWORDS
 $FILENAME_KEYWORD
 $SUFFIX_KEYWORD
diff --git a/src/methods/shared_docstrings.jl b/src/methods/shared_docstrings.jl
index 443808750..f2a690488 100644
--- a/src/methods/shared_docstrings.jl
+++ b/src/methods/shared_docstrings.jl
@@ -18,15 +18,15 @@ const CRS_KEYWORD = """
 - `crs`: a `crs` which will be attached to the resulting raster when `to` not passed
    or is an `Extent`. Otherwise the crs from `to` is used.
 """
-
-const SHAPE_KEYWORDS = """
-- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point` geometries.
-    using points or lines as polygons may have unexpected results.
+const BOUNDARY_KEYWORD = """
 - `boundary`: for polygons, include pixels where the `:center` is inside the polygon,
     where the polygon `:touches` the pixel, or that are completely `:inside` the polygon.
     The default is `:center`.
 """
-
+const SHAPE_KEYWORD = """
+- `shape`: Force `data` to be treated as `:polygon`, `:line` or `:point` geometries.
+    using points or lines as polygons may have unexpected results.
+"""
 const THREADED_KEYWORD = """
 - `threaded`: run operations in parallel, `false` by default. In some circumstances `threaded` 
     can give large speedups over single-threaded operation. This can be true for complicated 
@@ -41,7 +41,8 @@ $TO_KEYWORD
 $RES_KEYWORD
 $SIZE_KEYWORD
 $CRS_KEYWORD
-$SHAPE_KEYWORDS
+$BOUNDARY_KEYWORD
+$SHAPE_KEYWORD
 """
 
 const DATA_ARGUMENT = """
diff --git a/src/utils.jl b/src/utils.jl
index 18aead165..795ad63e2 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -252,7 +252,7 @@ _warn_disk() = @warn "Disk-based objects may be very slow here. User `read` firs
 
 _filenotfound_error(filename) = throw(ArgumentError("file \"$filename\" not found"))
 
-_progress(args...; kw...) = ProgressMeter.Progress(args...; color=:blue, barlen=50, kw...)
+_progress(args...; kw...) = ProgressMeter.Progress(args...; dt=0.1, color=:blue, barlen=50, kw...)
 
 # Function barrier for splatted vector broadcast
 @noinline _do_broadcast!(f, x, args...) = broadcast!(f, x, args...)
@@ -371,7 +371,11 @@ _names(A::AbstractRaster) = (Symbol(name(A)),)
 _names(A::AbstractRasterStack) = keys(A)
 
 using DimensionalData.Lookups: _True, _False
+
+_booltype(::_True) = _True()
+_booltype(::_False) = _False()
 _booltype(x) = x ? _True() : _False()
+
 istrue(::_True) = true
 istrue(::_False) = false
 
@@ -379,7 +383,7 @@ istrue(::_False) = false
 _rowtype(x, g, args...; kw...) = _rowtype(x, typeof(g), args...; kw...)
 function _rowtype(
     x, ::Type{G}, i::Type{I} = typeof(size(x)); 
-    geometry, index, skipmissing, skipinvalid = skipmissing, names, kw...
+    geometry, index, skipmissing, skipinvalid=skipmissing, names, kw...
 ) where {G, I}
     _G = istrue(skipinvalid) ? nonmissingtype(G) : G
     _I = istrue(skipinvalid) ? I : Union{Missing, I}
@@ -388,7 +392,6 @@ function _rowtype(
     NamedTuple{keys,types}
 end
 
-
 function _rowtypes(
     x, ::Type{G}, ::Type{I}, geometry::_True, index::_True, skipmissing, names::NamedTuple{Names}
 ) where {G,I,Names}
@@ -410,8 +413,10 @@ function _rowtypes(
     Tuple{_nametypes(x, names, skipmissing)...}
 end
 
-@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, skipmissing::_True) where {T,Names} = (nonmissingtype(T),)
-@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, skipmissing::_False) where {T,Names} = (Union{Missing,T},)
+@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, sm::_True) where {T,Names} = 
+    (nonmissingtype(T),)
+@inline _nametypes(::Raster{T}, ::NamedTuple{Names}, sm::_False) where {T,Names} = 
+    (Union{Missing,T},)
 # This only compiles away when generated
 @generated function _nametypes(
     ::RasterStack{<:Any,T}, ::NamedTuple{PropNames}, skipmissing::_True
diff --git a/test/extract.jl b/test/extract.jl
new file mode 100644
index 000000000..c9668d03a
--- /dev/null
+++ b/test/extract.jl
@@ -0,0 +1,455 @@
+using Rasters, Test, DataFrames, Extents
+import GeoInterface as GI
+using Rasters.Lookups, Rasters.Dimensions 
+
+include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
+
+dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2))
+rast = Raster(Union{Int,Missing}[1 2; 3 4], dimz; name=:test, missingval=missing)
+rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=5)
+rast_m = Raster([1 2; 3 missing], dimz; name=:test, missingval=missing)
+st = RasterStack(rast, rast2)
+
+points = [missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3), (10.0, 0.2)]
+poly = GI.Polygon([[(8.0, 0.0), (11.0, 0.0), (11.0, 0.4), (8.0, 0.0)]])
+linestring = GI.LineString([(8.0, 0.0), (9.5, 0.0), (10.0, 0.4)])
+line = GI.Line([(8.0, 0.0), (12.0, 0.4)])
+table = (geometry=points, foo=zeros(4))
+
+@testset "Points" begin
+    @testset "From Raster" begin
+        # Tuple points
+        ex = extract(rast, points)
+        T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}}
+        @test eltype(ex) == T
+        @test all(ex .=== T[
+            (geometry = missing, test = missing)
+            (geometry = (9.0, 0.1), test=1)
+            (geometry = (9.0, 0.2), test=2)
+            (geometry = (10.0, 0.3), test=missing)
+            (geometry = (10.0, 0.2), test=4)
+        ])
+        ex = extract(rast_m, points; skipmissing=true)
+        T = @NamedTuple{geometry::Tuple{Float64, Float64}, test::Int64}
+        @test eltype(ex) == T
+        @test all(ex .=== T[(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)])
+        ex = extract(rast_m, points; skipmissing=true, geometry=false)
+        T = @NamedTuple{test::Int64}
+        @test eltype(ex) == T
+        @test all(ex .=== T[(test = 1,), (test = 2,)])
+        @test all(extract(rast_m, points; skipmissing=true, geometry=false, index=true) .=== [
+            (index = (1, 1), test = 1,)
+            (index = (1, 2), test = 2,)
+        ])
+        # NamedTuple (reversed) points - tests a Table that iterates over points
+        T = @NamedTuple{geometry::Union{@NamedTuple{Y::Float64,X::Float64}},test::Union{Missing,Int64}}
+        @test all(extract(rast, [(Y=0.1, X=9.0), (Y=0.2, X=10.0), (Y=0.3, X=10.0)]) .=== T[
+            (geometry = (Y = 0.1, X = 9.0), test = 1)
+            (geometry = (Y = 0.2, X = 10.0), test = 4)
+            (geometry = (Y = 0.3, X = 10.0), test = missing)
+        ])
+        # Vector points
+        @test extract(rast, [[9.0, 0.1], [10.0, 0.2]] == [
+            (geometry = [9.0, 0.1], test = 1)
+            (geometry = [10.0, 0.2], test = 4)
+        ]
+    end
+
+    @testset "From RasterStack" begin
+        T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64},test2::Union{Missing,Int64}}
+        @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]) .=== T[
+            (geometry = missing, test = missing, test2 = missing)
+            (geometry = (9.0, 0.1), test = 1, test2 = 5)
+            (geometry = (10.0, 0.2), test = 4, test2 = 8)
+            (geometry = (10.0, 0.3), test = missing, test2 = missing)
+        ])
+        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true) == [
+            (geometry = (10.0, 0.2), test = 4, test2 = 8)
+        ]
+        @test extract(st2, [missing, (2, 2), (2, 1)]; skipmissing=true) == [
+            (geometry = (2, 1), a = 7.0, b = 2.0)
+        ]
+        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false) == [
+            (test = 4, test2 = 8)
+        ]
+        T = @NamedTuple{index::Union{Missing, Tuple{Int,Int}}, test::Union{Missing, Int64}, test2::Union{Missing, Int64}}
+        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false, index=true) == T[
+            (index = (2, 2), test = 4, test2 = 8)
+        ]
+        # Subset with `names`
+        T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test2::Union{Missing, Int64}}
+        @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,)) .=== T[
+            (geometry = missing, test2 = missing)
+            (geometry = (9.0, 0.1), test2 = 5)
+            (geometry = (10.0, 0.2), test2 = 8)
+            (geometry = (10.0, 0.3), test2 = missing)
+        ])
+        # Subset with `names` and `skipmissing` with mixed missingvals
+        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,), skipmissing=true) == [
+            (geometry = (10.0, 0.2), test2 = 8)
+        ]
+        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test,), skipmissing=true) == [
+            (geometry = (9.0, 0.1), test = 1)
+            (geometry = (10.0, 0.2), test = 4)
+        ]
+    end
+
+    @testset "Errors" begin
+        # Missing coord errors
+        @test_throws ArgumentError extract(rast, [(0.0, missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)])
+        @test_throws ArgumentError extract(rast, [(9.0, 0.1), (0.0, missing), (9.0, 0.2), (10.0, 0.3)])
+        @test_throws ArgumentError extract(rast, [(X=0.0, Y=missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)])
+    end
+end
+
+@testset "Polygons" begin
+    # Extract a polygon
+    T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},test::Union{Missing,Int64}}
+    @test all(extract(rast_m, poly) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+    # Test all the keyword combinations
+    @test all(extract(rast_m, poly) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+    T = @NamedTuple{test::Union{Missing,Int64}}
+    @test all(extract(rast_m, poly; geometry=false) .=== T[
+        (test = 1,)
+        (test = 3,)
+        (test = missing,)
+    ])
+    T = @NamedTuple{index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}}
+    @test all(extract(rast_m, poly; geometry=false, index=true) .=== T[
+        (index = (1, 1), test = 1)
+        (index = (2, 1), test = 3)
+        (index = (2, 2), test = missing)
+    ])
+    T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}}
+    @test all(extract(rast_m, poly; index=true) .=== T[
+         (geometry = (9.0, 0.1), index = (1, 1), test = 1)
+         (geometry = (10.0, 0.1), index = (2, 1), test = 3)
+         (geometry = (10.0, 0.2), index = (2, 2), test = missing)
+    ])
+    @test extract(rast_m, poly; skipmissing=true) == [
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+    ]                                                         
+    @test extract(rast_m, poly; skipmissing=true, geometry=false) == [
+        (test = 1,)
+        (test = 3,)
+    ]                                                         
+    @test extract(rast_m, poly; skipmissing=true, geometry=false, index=true) == [
+        (index = (1, 1), test = 1)
+        (index = (2, 1), test = 3)
+    ]                                                         
+    @test extract(rast_m, poly; skipmissing=true, index=true) == [
+        (geometry = (9.0, 0.1), index = (1, 1), test = 1)
+        (geometry = (10.0, 0.1), index = (2, 1), test = 3)
+    ]          
+    @test extract(rast2, poly; skipmissing=true) == [
+        (geometry = (10.0, 0.1), test2 = 7)
+        (geometry = (10.0, 0.2), test2 = 8)
+    ]                                               
+    T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},test::Union{Missing,Int64}}
+    @test all(extract(rast_m, poly) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+
+    @testset "Vector of polygons" begin
+        ex = extract(rast_m, [poly, poly, poly])
+        @test eltype(ex) == T
+        @test all(ex .=== T[
+            (geometry = (9.0, 0.1), test = 1)
+            (geometry = (10.0, 0.1), test = 3)
+            (geometry = (10.0, 0.2), test = missing)
+            (geometry = (9.0, 0.1), test = 1)
+            (geometry = (10.0, 0.1), test = 3)
+            (geometry = (10.0, 0.2), test = missing)
+        ])
+    end
+end
+
+@testset "Extract a linestring" begin
+    T = @NamedTuple{geometry::Tuple{Float64,Float64},test::Union{Missing,Int64}}
+    Tsm = @NamedTuple{geometry::Tuple{Float64,Float64},test::Int64}
+    linestrings = [linestring, linestring, linestring]
+    fc = GI.FeatureCollection(map(GI.Feature, linestrings))
+
+    # Single LineString
+    @test extract(rast, linestring) isa Vector{T}
+    @test all(extract(rast_m, linestring) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+    @test all(extract(rast_m, linestring; skipmissing=true) .=== Tsm[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+    ])
+
+    # Multiple linstrings
+    # Test all combinations of skipmissing, flatten and threaded
+    @test extract(rast_m, linestrings; skipmissing=true) isa Vector{Tsm}
+    @test extract(rast_m, fc; skipmissing=true) == 
+          extract(rast_m, fc; skipmissing=true, threaded=true) == 
+          extract(rast_m, linestrings; skipmissing=true) == 
+          extract(rast_m, linestrings; skipmissing=true, threaded=true) == Tsm[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+    ]
+
+    @test extract(rast_m, linestrings; skipmissing=false) isa Vector{T}
+    @test all(
+              extract(rast_m, fc; skipmissing=false) .=== 
+              extract(rast_m, fc; skipmissing=false, threaded=true) .===
+              extract(rast_m, linestrings; skipmissing=false) .=== 
+              extract(rast_m, linestrings; skipmissing=false, threaded=true) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (10.0, 0.1), test = 3)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+
+    @test extract(rast_m, linestrings; skipmissing=true, flatten=false) isa Vector{Vector{Tsm}}
+    @test extract(rast_m, linestrings; skipmissing=true, flatten=false) == 
+        extract(rast_m, linestrings; skipmissing=true, flatten=false, threaded=true) == Vector{Tsm}[
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3)],
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3)],
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3)],
+    ]
+
+    # Nested Vector holding missing needs special handling to check equality
+    @test extract(rast_m, linestrings; skipmissing=false, flatten=false) isa Vector{Vector{T}}
+    ref = Vector{T}[
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3), (geometry = (10.0, 0.2), test = missing)] ,
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3), (geometry = (10.0, 0.2), test = missing)],
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (10.0, 0.1), test = 3), (geometry = (10.0, 0.2), test = missing)],
+    ]
+    matching(a, b) = all(map(===, a, b))
+    @test all(map(matching, extract(rast_m, linestrings; skipmissing=false, flatten=false, threaded=false), ref))
+    @test all(map(matching, extract(rast_m, linestrings; skipmissing=false, flatten=false, threaded=true), ref))
+end
+
+@testset "Extract a line" begin
+    T = @NamedTuple{geometry::Tuple{Float64,Float64},test::Union{Missing,Int64}}
+    Tsm = @NamedTuple{geometry::Tuple{Float64,Float64},test::Int64}
+    lines = [line, line]
+    fc = GI.FeatureCollection(map(GI.Feature, lines))
+
+    # Single LineString
+    @test extract(rast, line) isa Vector{T}
+    @test all(extract(rast_m, line) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+    @test all(extract(rast_m, line; skipmissing=true) .=== Tsm[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+    ])
+
+    # Multiple linstrings
+    # Test all combinations of skipmissing, flatten and threaded
+    @test extract(rast_m, lines; skipmissing=true) isa Vector{Tsm}
+    @test extract(rast_m, fc; skipmissing=true) == 
+          extract(rast_m, fc; skipmissing=true, threaded=true) == 
+          extract(rast_m, lines; skipmissing=true) == 
+          extract(rast_m, lines; skipmissing=true, threaded=true) == Tsm[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+    ]
+
+    @test extract(rast_m, lines; skipmissing=false) isa Vector{T}
+    @test all(
+              extract(rast_m, fc; skipmissing=false) .=== 
+              extract(rast_m, fc; skipmissing=false, threaded=true) .===
+              extract(rast_m, lines; skipmissing=false) .=== 
+              extract(rast_m, lines; skipmissing=false, threaded=true) .=== T[
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+        (geometry = (10.0, 0.2), test = missing)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+        (geometry = (10.0, 0.2), test = missing)
+    ])
+
+    @test extract(rast_m, lines; skipmissing=true, flatten=false) isa Vector{Vector{Tsm}}
+    @test extract(rast_m, lines; skipmissing=true, flatten=false) == 
+        extract(rast_m, lines; skipmissing=true, flatten=false, threaded=true) == Vector{Tsm}[
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)],
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)],
+    ]
+
+    # Nested Vector holding missing needs special handling to check equality
+    @test extract(rast_m, lines; skipmissing=false, flatten=false) isa Vector{Vector{T}}
+    ref = Vector{T}[
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2), (geometry = (10.0, 0.2), test = missing)] ,
+        [(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2), (geometry = (10.0, 0.2), test = missing)],
+    ]
+    matching(a, b) = all(map(===, a, b))
+    @test all(map(matching, extract(rast_m, lines; skipmissing=false, flatten=false, threaded=false), ref))
+    @test all(map(matching, extract(rast_m, lines; skipmissing=false, flatten=false, threaded=true), ref))
+end
+
+@testset "Table" begin
+    T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test::Union{Missing, Int64}}
+    @test all(extract(rast, table) .=== T[
+        (geometry = missing, test = missing)
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+        (geometry = (10.0, 0.3), test = missing)
+        (geometry = (10.0, 0.2), test = 4)
+    ])
+    @test extract(rast, table; skipmissing=true) == [
+        (geometry = (9.0, 0.1), test = 1)
+        (geometry = (9.0, 0.2), test = 2)
+        (geometry = (10.0, 0.2), test = 4)
+    ]
+    @test extract(rast, table; skipmissing=true, geometry=false) == [
+        (test = 1,)
+        (test = 2,)
+        (test = 4,)
+    ]
+    @test extract(rast, table; skipmissing=true, geometry=false, index=true) == [
+        (index = (1, 1), test = 1,)
+        (index = (1, 2), test = 2,)
+        (index = (2, 2), test = 4,)
+    ]
+    @test_throws ArgumentError extract(rast, (foo = zeros(4),))
+end
+
+@testset "Empty geoms" begin
+    @test extract(rast, []) == NamedTuple{(:geometry, :test),Tuple{Missing,Missing}}[]
+    @test extract(rast, []; geometry=false) == NamedTuple{(:test,),Tuple{Missing}}[]
+end
+
+#=
+Some benchmark and plotting code
+
+using BenchmarkTools
+using ProfileView
+using Chairmarks
+using Cthulhu
+using Rasters
+using DataFrames
+import GeoInterface as GI
+using Rasters.Lookups, Rasters.Dimensions 
+
+dimz = (X(5.0:0.1:15.0; sampling=Intervals(Start())), Y(-0.1:0.01:0.5; sampling=Intervals(Start())))
+rast_m = Raster(rand([missing, 1, 2], dimz); name=:test, missingval=missing)
+rast = Raster(rand(Int, dimz); name=:test2, missingval=5)
+st = RasterStack((a=rast, b=rast, c=Float32.(rast)))
+ks = ntuple(x -> gensym(), 100)
+st100 = RasterStack(NamedTuple{ks}(ntuple(x -> rast, length(ks))))
+
+poly = GI.Polygon([[(8.0, 0.0), (11.0, 0.0), (11.0, 0.4), (8.0, 0.0)]])
+linestring = GI.LineString([(8.0, 0.0), (11.0, 0.0), (11.0, 0.4), (8.0, 0.0)])
+line = GI.Line([(8.0, 0.0), (12.0, 4.0)])
+polys = [poly for i in 1:10000]
+lines = [line for i in 1:100000]
+linestrings = [linestring for i in 1:100000]
+
+extract(rast, linestring)
+extract(rast, polys; threaded=false)
+@time extract(rast, polys; geometry=false, threaded=true, flatten=true) |> length
+@time extract(rast, lines; geometry=false, threaded=true, flatten=true) |> length
+@time extract(st, polys; geometry=false, threaded=true, flatten=true) |> length
+@time extract(st, lines; geometry=false, threaded=false, flatten=true) |> length
+
+@time extract(rast_m, lines; geometry=false, threaded=true, flatten=false) |> length
+@time extract(rast, lines; geometry=false, threaded=true, flatten=false) |> length
+@time extract(rast_m, lines; skipmissing=true, geometry=false, threaded=true, flatten=false) |> length
+@time extract(rast,   lines; skipmissing=true, geometry=false, threaded=true, flatten=false) |> length
+@time extract(st,     lines; skipmissing=true, geometry=false, threaded=true, flatten=false) |> length
+@time extract(st100, lines; geometry=false, threaded=false, flatten=false) |> length
+@time extract(st100, lines; geometry=false, threaded=true, flatten=false) |> length
+@time extract(rast, lines; geometry=false, threaded=true, flatten=true) |> length
+@time extract(rast, lines; threaded=true, flatten=true, progress=false) |> length
+@time extract(rast, lines; threaded=false, flatten=false, progress=false) |> length
+@benchmark extract(rast, lines; threaded=false, flatten=true)
+@benchmark extract(rast, lines; threaded=true, flatten=true)
+@benchmark extract(rast, lines; threaded=false, flatten=false)
+@benchmark extract(rast, lines; threaded=true, flatten=false)
+@time extract(rast, polys; flatten=false);
+extract(rast, polys; flatten=true);
+extract(st100, linestring)
+extract(st100, poly)
+
+@profview extract(rast, lines)
+@profview extract(rast, lines; threaded=true)
+@profview extract(rast, linestrings)
+@profview extract(rast_m, linestrings)
+@profview extract(rast, polys)
+@profview extract(rast, polys; flatten=false)
+@profview extract(st, linestring)
+@profview extract(st, poly)
+
+# Profile running one function a lot. 
+# People will do this
+f(rast, ls, n; skipmissing=true) = for _ in 1:n extract(rast, ls; skipmissing) end
+@profview f(rast, polies, 10; skipmissing=false)
+@profview f(rast, poly, 10000; skipmissing=false)
+@profview f(rast_m, poly, 10000; skipmissing=false)
+@profview f(rast, linestring, 10000; skipmissing=false)
+@profview f(rast_m, linestring, 10000; skipmissing=false)
+@profview f(st, linestring, 10000; skipmissing=false)
+@profview f(rast, linestring, 10000)
+@profview f(st, linestring, 100000; skipmissing=false)
+@profview f(st, linestring, 100000; skipmissing=true)
+@profview f(st100, linestring, 1000)
+@profview f(st100, line, 1000)
+@profview f(st100, poly, 10000; skipmissing=true)
+
+@profview f(st100, linestring, 1000; skipmissing=true)
+@profview f(st100, linestring, 1000; skipmissing=false)
+@profview f(st100, linestring, 1000; names=skipmissing=false)
+keys(st100)[1:50]
+
+
+
+# Demo plots
+using GLMakie
+using GeoInterfaceMakie
+
+Makie.plot(rast; colormap=:Reds)
+Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens)
+Makie.plot!(polys; alpha=0.2)
+ps = getindex.(extract(rast, poly; index=true, flatten=true), :geometry)
+Makie.scatter!(ps; color=:yellow)
+
+Makie.plot(rast; colormap=:Reds)
+Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens)
+Makie.plot!(poly; alpha=0.7)
+ps = getindex.(extract(rast, poly; index=true, flatten=true, boundary=:touches), :geometry);
+Makie.scatter!(ps; color=:pink)
+
+Makie.plot(rast; colormap=:Reds)
+Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens)
+Makie.plot!(polys; alpha=0.7)
+ps = getindex.(extract(rast, poly; index=true, flatten=true, boundary=:inside), :geometry)
+Makie.scatter!(ps; color=:pink)
+
+Makie.plot(rast; colormap=:Reds)
+Makie.plot!(rast[Touches(GI.extent(first(polys)))] ; colormap=:Greens)
+Makie.plot!(linestring; color=:violet, linewidth=5)
+ps = getindex.(extract(rast, linestring; index=true), :geometry);
+Makie.scatter!(ps; color=:yellow)
+
+=#
diff --git a/test/methods.jl b/test/methods.jl
index 81823e524..c2c6b1978 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -370,192 +370,6 @@ end
               [(9.0, 0.1) (9.0, 0.2); (10.0, 0.1) (10.0, 0.2)])
 end
 
-createpoint(args...) = ArchGDAL.createpoint(args...)
-
-@testset "extract" begin
-    dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2))
-    rast = Raster(Union{Int,Missing}[1 2; 3 4], dimz; name=:test, missingval=missing)
-    rast2 = Raster([5 6; 7 8], dimz; name=:test2, missingval=5)
-    rast_m = Raster([1 2; 3 missing], dimz; name=:test, missingval=missing)
-    mypoints = [missing, (9.0, 0.1), (9.0, 0.2), (10.0, 0.3), (10.0, 0.2)]
-    table = (geometry=mypoints, foo=zeros(4))
-    st = RasterStack(rast, rast2)
-    @testset "from Raster" begin
-        # Tuple points
-        ex = extract(rast, mypoints)
-        T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64}}
-        @test eltype(ex) == T
-        @test all(ex .=== T[
-            (geometry = missing, test = missing)
-            (geometry = (9.0, 0.1), test=1)
-            (geometry = (9.0, 0.2), test=2)
-            (geometry = (10.0, 0.3), test=missing)
-            (geometry = (10.0, 0.2), test=4)
-        ])
-        ex = extract(rast_m, mypoints; skipmissing=true)
-        T = @NamedTuple{geometry::Tuple{Float64, Float64}, test::Int64}
-        @test eltype(ex) == T
-        @test all(ex .=== T[(geometry = (9.0, 0.1), test = 1), (geometry = (9.0, 0.2), test = 2)])
-        ex = extract(rast_m, mypoints; skipmissing=true, geometry=false)
-        T = @NamedTuple{test::Int64}
-        @test eltype(ex) == T
-        @test all(ex .=== T[(test = 1,), (test = 2,)])
-        @test all(extract(rast_m, mypoints; skipmissing=true, geometry=false, index=true) .=== [
-            (index = (1, 1), test = 1,)
-            (index = (1, 2), test = 2,)
-        ])
-        # NamedTuple (reversed) points - tests a Table that iterates over points
-        T = @NamedTuple{geometry::Union{@NamedTuple{Y::Float64,X::Float64}},test::Union{Missing,Int64}}
-        @test all(extract(rast, [(Y=0.1, X=9.0), (Y=0.2, X=10.0), (Y=0.3, X=10.0)]) .=== T[
-            (geometry = (Y = 0.1, X = 9.0), test = 1)
-            (geometry = (Y = 0.2, X = 10.0), test = 4)
-            (geometry = (Y = 0.3, X = 10.0), test = missing)
-        ])
-        # Vector points
-        @test all(extract(rast, [[9.0, 0.1], [10.0, 0.2]]) .== [
-            (geometry = [9.0, 0.1], test = 1)
-            (geometry = [10.0, 0.2], test = 4)
-        ])
-        # Extract a polygon
-        p = ArchGDAL.createpolygon([[[8.0, 0.0], [11.0, 0.0], [11.0, 0.4], [8.0, 0.0]]])
-        T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},test::Union{Missing,Int64}}
-        @test all(extract(rast_m, p) .=== T[
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (10.0, 0.1), test = 3)
-            (geometry = (10.0, 0.2), test = missing)
-        ])
-        # Extract a vector of polygons
-        ex = extract(rast_m, [p, p])
-        @test eltype(ex) == T
-        @test all(ex .=== T[
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (10.0, 0.1), test = 3)
-            (geometry = (10.0, 0.2), test = missing)
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (10.0, 0.1), test = 3)
-            (geometry = (10.0, 0.2), test = missing)
-        ])
-        # Test all the keyword combinations
-        @test all(extract(rast_m, p) .=== T[
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (10.0, 0.1), test = 3)
-            (geometry = (10.0, 0.2), test = missing)
-        ])
-        T = @NamedTuple{test::Union{Missing,Int64}}
-        @test all(extract(rast_m, p; geometry=false) .=== T[
-            (test = 1,)
-            (test = 3,)
-            (test = missing,)
-        ])
-        T = @NamedTuple{index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}}
-        @test all(extract(rast_m, p; geometry=false, index=true) .=== T[
-            (index = (1, 1), test = 1)
-            (index = (2, 1), test = 3)
-            (index = (2, 2), test = missing)
-        ])
-        T = @NamedTuple{geometry::Union{Tuple{Float64,Float64}},index::Union{Missing,Tuple{Int,Int}},test::Union{Missing,Int64}}
-        @test all(extract(rast_m, p; index=true) .=== T[
-             (geometry = (9.0, 0.1), index = (1, 1), test = 1)
-             (geometry = (10.0, 0.1), index = (2, 1), test = 3)
-             (geometry = (10.0, 0.2), index = (2, 2), test = missing)
-        ])
-        @test extract(rast_m, p; skipmissing=true) == [
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (10.0, 0.1), test = 3)
-        ]                                                         
-        @test extract(rast_m, p; skipmissing=true, geometry=false) == [
-            (test = 1,)
-            (test = 3,)
-        ]                                                         
-        @test extract(rast_m, p; skipmissing=true, geometry=false, index=true) == [
-            (index = (1, 1), test = 1)
-            (index = (2, 1), test = 3)
-        ]                                                         
-        @test extract(rast_m, p; skipmissing=true, index=true) == [
-            (geometry = (9.0, 0.1), index = (1, 1), test = 1)
-            (geometry = (10.0, 0.1), index = (2, 1), test = 3)
-        ]          
-        @test extract(rast2, p; skipmissing=true) == [
-            (geometry = (10.0, 0.1), test2 = 7)
-            (geometry = (10.0, 0.2), test2 = 8)
-        ]                                               
-        # Empty geoms
-        @test extract(rast, []) == NamedTuple{(:geometry, :test),Tuple{Missing,Missing}}[]
-        @test extract(rast, []; geometry=false) == NamedTuple{(:test,),Tuple{Missing}}[]
-        # Missing coord errors
-        @test_throws ArgumentError extract(rast, [(0.0, missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)])
-        @test_throws ArgumentError extract(rast, [(9.0, 0.1), (0.0, missing), (9.0, 0.2), (10.0, 0.3)])
-        @test_throws ArgumentError extract(rast, [(X=0.0, Y=missing), (9.0, 0.1), (9.0, 0.2), (10.0, 0.3)])
-    end
-
-    @testset "with table" begin
-        T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test::Union{Missing, Int64}}
-        @test all(extract(rast, table) .=== T[
-            (geometry = missing, test = missing)
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (9.0, 0.2), test = 2)
-            (geometry = (10.0, 0.3), test = missing)
-            (geometry = (10.0, 0.2), test = 4)
-        ])
-        @test extract(rast, table; skipmissing=true) == [
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (9.0, 0.2), test = 2)
-            (geometry = (10.0, 0.2), test = 4)
-        ]
-        @test extract(rast, table; skipmissing=true, geometry=false) == [
-            (test = 1,)
-            (test = 2,)
-            (test = 4,)
-        ]
-        @test extract(rast, table; skipmissing=true, geometry=false, index=true) == [
-            (index = (1, 1), test = 1,)
-            (index = (1, 2), test = 2,)
-            (index = (2, 2), test = 4,)
-        ]
-
-        @test_throws ArgumentError extract(rast, (foo = zeros(4),))
-    end
-
-    @testset "from stack" begin
-        T = @NamedTuple{geometry::Union{Missing,Tuple{Float64,Float64}},test::Union{Missing,Int64},test2::Union{Missing,Int64}}
-        @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]) .=== T[
-            (geometry = missing, test = missing, test2 = missing)
-            (geometry = (9.0, 0.1), test = 1, test2 = 5)
-            (geometry = (10.0, 0.2), test = 4, test2 = 8)
-            (geometry = (10.0, 0.3), test = missing, test2 = missing)
-        ])
-        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true) == [
-            (geometry = (10.0, 0.2), test = 4, test2 = 8)
-        ]
-        @test extract(st2, [missing, (2, 2), (2, 1)]; skipmissing=true) == [
-            (geometry = (2, 1), a = 7.0, b = 2.0)
-        ]
-        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false) == [
-            (test = 4, test2 = 8)
-        ]
-        T = @NamedTuple{index::Union{Missing, Tuple{Int,Int}}, test::Union{Missing, Int64}, test2::Union{Missing, Int64}}
-        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; skipmissing=true, geometry=false, index=true) == T[
-            (index = (2, 2), test = 4, test2 = 8)
-        ]
-        # Subset with `names`
-        T = @NamedTuple{geometry::Union{Missing, Tuple{Float64, Float64}}, test2::Union{Missing, Int64}}
-        @test all(extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,)) .=== T[
-            (geometry = missing, test2 = missing)
-            (geometry = (9.0, 0.1), test2 = 5)
-            (geometry = (10.0, 0.2), test2 = 8)
-            (geometry = (10.0, 0.3), test2 = missing)
-        ])
-        # Subset with `names` and `skipmissing` with mixed missingvals
-        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test2,), skipmissing=true) == [
-            (geometry = (10.0, 0.2), test2 = 8)
-        ]
-        @test extract(st, [missing, (9.0, 0.1), (10.0, 0.2), (10.0, 0.3)]; name=(:test,), skipmissing=true) == [
-            (geometry = (9.0, 0.1), test = 1)
-            (geometry = (10.0, 0.2), test = 4)
-        ]
-    end
-end
-
 @testset "trim, crop, extend" begin
     A = [missing missing missing
          missing 2.0     0.5
diff --git a/test/runtests.jl b/test/runtests.jl
index e4337545c..93b2c9ec6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,15 +10,16 @@ using Rasters, Test, Aqua, SafeTestsets
 end
 
 @time @safetestset "extensions" begin include("extensions.jl") end
-@time @safetestset "methods" begin include("methods.jl") end
 @time @safetestset "array" begin include("array.jl") end
 @time @safetestset "stack" begin include("stack.jl") end
 @time @safetestset "series" begin include("series.jl") end
 @time @safetestset "utils" begin include("utils.jl") end
 @time @safetestset "set" begin include("set.jl") end
+@time @safetestset "adapt" begin include("adapt.jl") end
+@time @safetestset "methods" begin include("methods.jl") end
 @time @safetestset "aggregate" begin include("aggregate.jl") end
 @time @safetestset "rasterize" begin include("rasterize.jl") end
-@time @safetestset "adapt" begin include("adapt.jl") end
+@time @safetestset "extract" begin include("extract.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