diff --git a/Project.toml b/Project.toml
index a7ba6b713..02d79c98a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,6 +25,7 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
@@ -58,7 +59,7 @@ CommonDataModel = "0.2.3, 0.3"
 ConstructionBase = "1"
 CoordinateTransformations = "0.6.2"
 DataFrames = "1"
-DimensionalData = "0.28.2"
+DimensionalData = "0.29.4"
 DiskArrays = "0.3, 0.4"
 Extents = "0.1"
 FillArrays = "0.12, 0.13, 1"
@@ -106,8 +107,7 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "Statistics", "StatsBase", "Test", "ZarrDatasets"]
+test = ["Aqua", "ArchGDAL", "CFTime", "CoordinateTransformations", "DataFrames", "GeoDataFrames", "GeometryBasics", "GRIBDatasets", "NCDatasets", "Plots", "Proj", "RasterDataSources", "SafeTestsets", "Shapefile", "StableRNGs", "StatsBase", "Test", "ZarrDatasets"]
diff --git a/src/Rasters.jl b/src/Rasters.jl
index 57d02f577..94af83dfe 100644
--- a/src/Rasters.jl
+++ b/src/Rasters.jl
@@ -28,19 +28,21 @@ import Adapt,
        Mmap,
        RecipesBase,
        Reexport,
-       Setfield
+       Setfield,
+       Statistics
 
 Reexport.@reexport using DimensionalData, GeoFormatTypes
 
 using DimensionalData.Tables,
       DimensionalData.Lookups,
-      DimensionalData.Dimensions
+      DimensionalData.Dimensions,
       DimensionalData.Lookups.IntervalSets
 
 using DimensionalData: Name, NoName
 using .Dimensions: StandardIndices, DimTuple
 using .Lookups: LookupTuple
 
+using Statistics: mean
 using RecipesBase: @recipe, @series
 using Base: tail, @propagate_inbounds
 
diff --git a/src/array.jl b/src/array.jl
index 7b0a02d7d..92f5d52e7 100644
--- a/src/array.jl
+++ b/src/array.jl
@@ -105,20 +105,17 @@ function DD.modify(f, A::AbstractRaster)
     return rebuild(A, newdata)
 end
 
-function DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{AbstractRaster}}...)
-    DD.DimTable(DimStack(map(read, As...)))
-end
+DD.DimTable(As::Tuple{<:AbstractRaster,Vararg{AbstractRaster}}) =
+    DD.DimTable(DimStack(map(read, As)))
 
 # DiskArrays methods
 
 DiskArrays.eachchunk(A::AbstractRaster) = DiskArrays.eachchunk(parent(A))
 DiskArrays.haschunks(A::AbstractRaster) = DiskArrays.haschunks(parent(A))
-function DA.readblock!(A::AbstractRaster, dst, r::AbstractUnitRange...)
+DA.readblock!(A::AbstractRaster, dst, r::AbstractUnitRange...) =
     DA.readblock!(parent(A), dst, r...)
-end
-function DA.writeblock!(A::AbstractRaster, src, r::AbstractUnitRange...)
+DA.writeblock!(A::AbstractRaster, src, r::AbstractUnitRange...) =
     DA.writeblock!(parent(A), src, r...)
-end
 
 # Base methods
 
diff --git a/src/lookup.jl b/src/lookup.jl
index 0c0c05fa1..07d1a68e0 100644
--- a/src/lookup.jl
+++ b/src/lookup.jl
@@ -89,18 +89,30 @@ GeoInterface.crs(lookup::Projected) = lookup.crs
 mappedcrs(lookup::Projected) = lookup.mappedcrs
 dim(lookup::Projected) = lookup.dim
 
-function LA.selectindices(l::Projected, sel::Contains)
+@inline function LA.selectindices(l::Projected, sel::LA.Selector; kw...)
     selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
-    LA.contains(l, rebuild(sel; val=selval))
+    LA._selectindices(l, rebuild(sel; val=selval); kw...)
 end
-function LA.selectindices(l::Projected, sel::At)
+@inline function LA.selectindices(l::Projected, sel::LA.Selector{<:AbstractVector}; kw...)
     selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
-    LA.at(l, rebuild(sel; val=selval))
+    _selectvec(l, rebuild(sel; val=selval)sel; kw...)
 end
-function LA.selectindices(l::Projected, sel::Between)
+@inline function LA.selectindices(l::Projected, sel::LA.IntSelector{<:Tuple}; kw...)
+    selval = reproject(mappedcrs(l), crs(l), dim(l), val(sel))
+    _selecttuple(l, rebuild(sel; val=selval); kw...)
+end
+@inline LA.selectindices(l::Projected{<:Tuple}, sel::LA.IntSelector{<:Tuple}; kw...) = LA._selectindices(l, sel; kw...)
+@inline LA.selectindices(l::Projected{<:Tuple}, sel::LA.IntSelector{<:Tuple{<:Tuple,<:Tuple}}; kw...) = 
+    _selecttuple(l, sel; kw...)
+
+function LA.selectindices(l::Projected, sel::Between{<:Tuple})
     selval = map(v -> reproject(mappedcrs(l), crs(l), dim(l), v), val(sel))
     LA.between(l, rebuild(sel; val=selval))
 end
+function LA.selectindices(l::Projected, sel::T) where T<:DD.IntervalSets.Interval
+    left, right = map(v -> reproject(mappedcrs(l), crs(l), dim(l), v), (sel.left, sel.right))
+    LA.between(l, T(left, right))
+end
 
 """
     Mapped <: AbstractProjected
diff --git a/src/methods/burning/utils.jl b/src/methods/burning/utils.jl
index d51af3c8f..0e87d3e73 100644
--- a/src/methods/burning/utils.jl
+++ b/src/methods/burning/utils.jl
@@ -4,16 +4,16 @@ _geomindices(::GI.FeatureCollectionTrait, geoms) = 1:GI.nfeature(geoms)
 _geomindices(::GI.FeatureTrait, geoms) = _geomindices(GI.geometry(geoms))
 _geomindices(::GI.AbstractGeometryTrait, geoms) = 1:GI.ngeom(geoms)
 
-_getgeom(geoms, i::Integer) = _getgeom(GI.trait(geoms), geoms, i)
-_getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i))
-_getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i)
-_getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i)
-_getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached")
-_getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index?
+_getgeom(geoms, i::Integer) = __getgeom(GI.trait(geoms), geoms, i)
+__getgeom(::GI.FeatureCollectionTrait, geoms, i::Integer) = GI.geometry(GI.getfeature(geoms, i))
+__getgeom(::GI.FeatureTrait, geoms, i::Integer) = GI.getgeom(GI.geometry(geoms), i)
+__getgeom(::GI.AbstractGeometryTrait, geoms, i::Integer) = GI.getgeom(geom, i)
+__getgeom(::GI.PointTrait, geom, i::Integer) = error("PointTrait should not be reached")
+__getgeom(::Nothing, geoms, i::Integer) = geoms[i] # Otherwise we can probably just index?
 
-_getgeom(geoms) = _getgeom(GI.trait(geoms), geoms)
-_getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms))
-_getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms))
-_getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms)
-_getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached")
-_getgeom(::Nothing, geoms) = geoms
+_getgeom(geoms) = __getgeom(GI.trait(geoms), geoms)
+__getgeom(::GI.FeatureCollectionTrait, geoms) = (GI.geometry(f) for f in GI.getfeature(geoms))
+__getgeom(::GI.FeatureTrait, geoms) = GI.getgeom(GI.geometry(geoms))
+__getgeom(::GI.AbstractGeometryTrait, geoms) = GI.getgeom(geoms)
+__getgeom(::GI.PointTrait, geom) = error("PointTrait should not be reached")
+__getgeom(::Nothing, geoms) = geoms
diff --git a/src/methods/classify.jl b/src/methods/classify.jl
index 22a587944..be046ecb9 100644
--- a/src/methods/classify.jl
+++ b/src/methods/classify.jl
@@ -152,7 +152,11 @@ function classify!(A::AbstractRaster, pairs;
     end
     return rebuild(out; missingval=missingval)
 end
-function classify!(xs::RasterSeriesOrStack, pairs...; kw...)
+function classify!(xs::RasterStack, pairs...; kw...)
+    maplayers(x -> classify!(x, pairs...; kw...),  xs)
+    return xs
+end
+function classify!(xs::RasterSeries, pairs...; kw...)
     map(x -> classify!(x, pairs...; kw...),  xs)
     return xs
 end
diff --git a/src/methods/mask.jl b/src/methods/mask.jl
index 55b0f3732..1b19b44cb 100644
--- a/src/methods/mask.jl
+++ b/src/methods/mask.jl
@@ -170,7 +170,7 @@ function _mask!(x::RasterStackOrArray, geom; kw...)
 end
 # Array mask
 function _mask!(st::RasterStack, with::AbstractRaster; kw...)
-    map(A -> mask!(A; with, kw...), st)
+    maplayers(A -> mask!(A; with, kw...), st)
     return st
 end
 
@@ -431,7 +431,9 @@ end
 
 _false_to_missing(b::Bool) = (b ? true : missing)::Union{Missing,Bool}
 
-function _mask_multilayer(layers::Union{<:AbstractRasterStack,<:AbstractRasterSeries}, to; 
+_mask_multilayer(st::AbstractRasterStack, to; kw...) =
+    _mask_multilayer(layers(st), to; kw...)
+function _mask_multilayer(layers::Union{<:NamedTuple,<:AbstractRasterSeries}, to; 
     _dest_presentval,
     _dest_missingval, 
     missingval=nothing, 
diff --git a/src/methods/mosaic.jl b/src/methods/mosaic.jl
index 17f13c6e3..0c913d173 100644
--- a/src/methods/mosaic.jl
+++ b/src/methods/mosaic.jl
@@ -74,7 +74,7 @@ function _mosaic(f::Function, ::AbstractRaster, regions;
     l1 = first(regions)
     A = create(filename, T, dims; name=name(l1), missingval, metadata=metadata(l1))
     open(A; write=true) do a
-        mosaic!(f, a, regions; missingval, kw...)
+        _mosaic!(f, a, regions; missingval, kw...)
     end
     return A
 end
@@ -135,8 +135,10 @@ nothing
 
 $EXPERIMENTAL
 """
-mosaic!(f::Function, x::RasterStackOrArray, regions::RasterStackOrArray...; kw...) = mosaic!(f, x, regions; kw...)
-function mosaic!(f::Function, A::AbstractRaster{T}, regions;
+mosaic!(f::Function, dest::RasterStackOrArray, regions::RasterStackOrArray...; kw...) =
+    _mosaic!(f, dest, regions; kw...)
+
+function _mosaic!(f::Function, A::AbstractRaster{T}, regions::Union{Tuple,AbstractArray};
     missingval=missingval(A), atol=maybe_eps(T)
 ) where T
     _without_mapped_crs(A) do A1
@@ -169,10 +171,11 @@ function mosaic!(f::Function, A::AbstractRaster{T}, regions;
     end
     return A
 end
-function mosaic!(f::Function, st::AbstractRasterStack, regions::Tuple; kw...)
-    map(st, regions...) do A, r...
-        mosaic!(f, A, r; kw...)
+function _mosaic!(f::Function, st::AbstractRasterStack, regions::Union{Tuple,AbstractArray}; kw...)
+    map(values(st), map(values, regions)...) do A, r...
+        mosaic!(f, A, r...; kw...)
     end
+    return st
 end
 
 _mosaic(alldims::Tuple{<:DimTuple,Vararg{DimTuple}}) = map(_mosaic, alldims...)
diff --git a/src/methods/rasterize.jl b/src/methods/rasterize.jl
index 0b1640b0e..db0194165 100644
--- a/src/methods/rasterize.jl
+++ b/src/methods/rasterize.jl
@@ -420,7 +420,7 @@ function rasterize(reducer::typeof(count), data; fill=nothing, init=nothing, kw.
 end
 # `mean` is sum ./ count. This is actually optimal with threading, 
 # as its means order is irrelevant so its threadsafe.
-function rasterize(reducer::typeof(DD.Statistics.mean), data; fill, kw...)
+function rasterize(reducer::typeof(Statistics.mean), data; fill, kw...)
     sums = rasterize(sum, data; kw..., fill)
     counts = rasterize(count, data; kw..., fill=nothing)
     rebuild(sums ./ counts; name=:mean)
@@ -831,6 +831,7 @@ end
 # We get 64 Bool values to a regular `Int` meaning this doesn't scale too
 # badly for large tables of geometries. 64k geometries and a 1000 * 1000
 # raster needs 1GB of memory just for the `BitArray`.
+# TODO combine these theyre nearly the same
 function _reduce_bitarray!(f, st::AbstractRasterStack, geoms, fill::NamedTuple, r::Rasterizer, allocs)
     (; lock, shape, boundary, verbose, progress, threaded) = r
     # Define mask dimensions, the same size as the spatial dims of x
@@ -840,7 +841,7 @@ function _reduce_bitarray!(f, st::AbstractRasterStack, geoms, fill::NamedTuple,
     # Use a generator over the array axis in case the iterator has no length
     geom_axis = axes(masks, Dim{:geometry}())
     fill = map(itr -> [v for (_, v) in zip(geom_axis, itr)], fill)
-    T = NamedTuple{keys(st),Tuple{map(eltype, st)...}}
+    T = eltype(st)
     range = axes(st, Y())
     _run(range, threaded, progress, "Reducing...") do y
         _reduce_bitarray_loop(f, st, T, fill, masks, y)
diff --git a/src/methods/trim.jl b/src/methods/trim.jl
index 21830f8ef..588eb7dcf 100644
--- a/src/methods/trim.jl
+++ b/src/methods/trim.jl
@@ -131,5 +131,5 @@ end
 
 # Broadcast over the array and tracker to mark axis indices as being missing or not
 _update!(tr::AxisTrackers, A::AbstractRaster) = tr .= A .!== missingval(A)
-_update!(tr::AxisTrackers, st::AbstractRasterStack) = map(A -> tr .= A .!== missingval(A), st)
+_update!(tr::AxisTrackers, st::AbstractRasterStack) = maplayers(A -> tr .= A .!== missingval(A), st)
 
diff --git a/src/methods/zonal.jl b/src/methods/zonal.jl
index 476386688..1263a0acf 100644
--- a/src/methods/zonal.jl
+++ b/src/methods/zonal.jl
@@ -84,7 +84,7 @@ _zonal(f, x::Raster, of::Extents.Extent; skipmissing=true) =
 function _zonal(f, x::RasterStack, ext::Extents.Extent; skipmissing=true)
     cropped = crop(x; to=ext, touches=true)
     prod(size(cropped)) > 0 || return missing
-    return map(cropped) do A
+    return maplayers(cropped) do A
         _maybe_skipmissing_call(f, A, skipmissing)
     end
 end
@@ -107,9 +107,9 @@ function _zonal(f, st::AbstractRasterStack, ::GI.AbstractGeometryTrait, geom;
     skipmissing=true, kw...
 )
     cropped = crop(st; to=geom, touches=true)
-    prod(size(cropped)) > 0 || return map(_ -> missing, st)
+    prod(size(cropped)) > 0 || return map(_ -> missing, layerdims(st))
     masked = mask(cropped; with=geom, kw...)
-    return map(masked) do A
+    return maplayers(masked) do A
         prod(size(A)) > 0 || return missing
         _maybe_skipmissing_call(f, A, skipmissing)
     end
diff --git a/src/read.jl b/src/read.jl
index eb6c7f622..af01f16ef 100644
--- a/src/read.jl
+++ b/src/read.jl
@@ -50,6 +50,10 @@ function Base.read!(src::AbstractRasterSeries, dst::AbstractRasterSeries)
     map(read!, src, dst)
     return dst
 end
+Base.read!(::AbstractRaster{<:Union{AbstractString,NamedTuple},1}, ::AbstractRasterSeries) =
+    error("Cant read Raster to RasterSeries")
+Base.read!(::AbstractRaster, ::AbstractRasterSeries) =
+    error("Cant read Raster to RasterSeries")
 
 # Filename methods
 function Base.read!(filename::AbstractString, dst::AbstractRaster)
diff --git a/src/series.jl b/src/series.jl
index 1b23caa6f..b4070f4a9 100644
--- a/src/series.jl
+++ b/src/series.jl
@@ -111,10 +111,21 @@ When loading a series from a single `String` path:
 Others:
 - `refdims`: existing reference dimension/s, normally not required.
 """
-struct RasterSeries{T,N,D,R,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
+struct RasterSeries{T<:Union{AbstractRaster,AbstractRasterStack},N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N}} <: AbstractRasterSeries{T,N,D,A}
     data::A
     dims::D
     refdims::R
+    # Needed for ambiguity with DD
+    function RasterSeries{T,N,D,R,A}(
+        data::AbstractArray, dims::Tuple, refdims::Tuple
+    ) where {T,N,D,R,A}
+        new{T,N,D,R,A}(data, dims, refdims)
+    end
+end
+function RasterSeries(
+    data::A, dims::D, refdims::R
+) where {A<:AbstractArray{T,N},D<:Tuple,R<:Tuple} where {T,N}
+    RasterSeries{T,N,D,R,A}(data, dims, refdims)
 end
 function RasterSeries(data::AbstractArray{<:Union{AbstractRasterStack,AbstractRaster}}, dims;
     refdims=()
@@ -205,6 +216,17 @@ end
 
 @inline function DD.rebuild(
     A::RasterSeries, data, dims::Tuple, refdims=(), name=nothing, metadata=nothing,
+)
+    # if `data` is not an AbstractArray of Raster or RasterStacks, return a Raster
+    Raster(data, dims; refdims, name, metadata)
+end
+@inline function DD.rebuild(
+    A::RasterSeries, 
+    data::AbstractArray{<:Union{AbstractRaster,AbstractRasterStack}}, 
+    dims::Tuple,
+    refdims=(),
+    name=nothing,
+    metadata=nothing,
 )
     RasterSeries(data, dims, refdims)
 end
@@ -212,7 +234,7 @@ end
     A::RasterSeries;
     data=parent(A), dims=dims(A), refdims=refdims(A), name=nothing, metadata=nothing,
 )
-    RasterSeries(data, dims, refdims)
+    rebuild(A, data, dims, refdims)
 end
 
 function Base.map(f, series::RasterSeries)
diff --git a/src/stack.jl b/src/stack.jl
index cfdd62b74..796793e73 100644
--- a/src/stack.jl
+++ b/src/stack.jl
@@ -31,8 +31,10 @@ abstract type AbstractRasterStack{K,T,N,L} <: AbstractDimStack{K,T,N,L} end
 missingval(stack::AbstractRasterStack) = getfield(stack, :missingval)
 missingval(s::AbstractRasterStack, name::Symbol) = _singlemissingval(missingval(s), name)
 
-filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = map(s -> filename(s), stack)
-filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = filename(parent(stack))
+filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:NamedTuple}) = 
+    map(s -> filename(s), layers(stack))
+filename(stack::AbstractRasterStack{<:Any,<:Any,<:Any,<:Union{FileStack,OpenStack}}) = 
+    filename(parent(stack))
 
 isdisk(st::AbstractRasterStack) = any(isdisk, layers(st))
 
@@ -83,12 +85,15 @@ function DD.rebuild(s::AbstractRasterStack;
 end
 
 function DD.rebuild_from_arrays(
-    s::AbstractRasterStack{<:Union{FileStack{<:Any,Keys},OpenStack{<:Any,Keys}}}, das::Tuple{Vararg{AbstractDimArray}}; kw...
+    s::AbstractRasterStack{<:Union{FileStack{<:Any,Keys},OpenStack{<:Any,Keys}}}, 
+    das::Tuple{Vararg{AbstractDimArray}}; 
+    kw...
 ) where Keys
     DD.rebuild_from_arrays(s, NamedTuple{Keys}(das); kw...)
 end
 function DD.rebuild_from_arrays(
-    s::AbstractRasterStack, das::NamedTuple{<:Any,<:Tuple{Vararg{AbstractDimArray}}};
+    s::AbstractRasterStack, 
+    das::NamedTuple{<:Any,<:Tuple{Vararg{AbstractDimArray}}};
     refdims=refdims(s),
     metadata=DD.metadata(s),
     dims=nothing,
@@ -112,7 +117,7 @@ end
 
 DD.name(s::AbstractRasterStack) = keys(s)
 Base.names(s::AbstractRasterStack) = keys(s)
-Base.copy(stack::AbstractRasterStack) = map(copy, stack)
+Base.copy(stack::AbstractRasterStack) = maplayers(copy, stack)
 
 #### Stack getindex ####
 # Different to DimensionalData as we construct a Raster
@@ -205,17 +210,17 @@ function RasterStack(
     data::Union{FileStack,OpenStack,NamedTuple};
     dims::Tuple,
     refdims::Tuple=(),
-    layerdims,
+    layerdims::NamedTuple,
     metadata=nokw,
     layermetadata=nokw,
     missingval=nokw,
     kw...
 )
-    # Handle values that musbe be `NamedTuple`
+    # Handle values that must be be `NamedTuple`
     layermetadata = if layermetadata isa NamedTuple
         layermetadata
-    elseif layermetadata isa Union{Nothing,NoMetadata}
-        map(_ -> NoMetadata(), layers)
+    elseif layermetadata isa Union{NoKW,Nothing,NoMetadata}
+        map(_ -> NoMetadata(), layerdims)
     else
         throw(ArgumentError("$layermetadata is not a valid input for `layermetadata`. Try a `NamedTuple` of `Dict`, `MetaData` or `NoMetadata`"))
     end
@@ -392,10 +397,10 @@ function RasterStack(filename::AbstractString;
             l_st = _layer_stack(filename; source, name, lazy, group, replace_missing, kw...)
 
             # Maybe split the stack into separate arrays to remove extra dims.
-            if !isnokw(name)
-                map(identity, l_st)
-            else
+            if isnokw(name)
                 l_st
+            else
+                maplayers(identity, l_st)
             end
         else
             # With bands actings as layers
@@ -415,12 +420,14 @@ function RasterStack(filename::AbstractString;
     end
 end
 
+# TODO test this properly
 function DD.modify(f, s::AbstractRasterStack{<:FileStack{<:Any,K}}) where K
-    open(s) do o
+    data = open(s) do o
         map(K) do k
-            Array(parent(ost)[k])
+            f(parent(ost)[k])
         end
     end
+    rebuild(s; data)
 end
 
 # Open a single file stack
diff --git a/src/utils.jl b/src/utils.jl
index e069755c3..18aead165 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -115,26 +115,29 @@ _without_mapped_crs(f, x, ::Nothing) = f(x)
 function _without_mapped_crs(f, dims::DimTuple, mappedcrs::GeoFormat)
     dims1 = setmappedcrs(dims, nothing)
     x = f(dims1)
-    if x isa DimTuple
-        x = setmappedcrs(x, mappedcrs)
+    return if x isa DimTuple
+        setmappedcrs(x, mappedcrs)
+    else
+        x
     end
-    return x
 end
 function _without_mapped_crs(f, A::AbstractRaster, mappedcrs::GeoFormat)
     A = setmappedcrs(A, nothing)
     x = f(A)
-    if x isa AbstractRaster
-        x = setmappedcrs(x, mappedcrs)
+    return if x isa AbstractRaster
+        setmappedcrs(x, mappedcrs)
+    else
+        x
     end
-    return x
 end
 function _without_mapped_crs(f, st::AbstractRasterStack, mappedcrs::GeoFormat)
-    st1 = map(A -> setmappedcrs(A, nothing), st)
+    st1 = maplayers(A -> setmappedcrs(A, nothing), st)
     x = f(st1)
-    if x isa AbstractRasterStack
-        x = map(A -> setmappedcrs(A, mappedcrs), x)
+    return if x isa AbstractRasterStack
+        setmappedcrs(x, mappedcrs)
+    else
+        x
     end
-    return x
 end
 
 function _extent2dims(to; size=nothing, res=nothing, crs=nothing, kw...) 
diff --git a/test/methods.jl b/test/methods.jl
index 324edaa8f..ead6fd173 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -59,7 +59,7 @@ gaMi = replace_missing(ga)
     @test all(isequal.(dNaN, [NaN32 7.0f0; 2.0f0 NaN32]))
     rm("test.tif")
     stNaN = replace_missing(st, NaN32; filename="teststack.tif")
-    @test all(map(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
+    @test all(maplayers(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
         all(x .=== y)
     end)
     rm("teststack_a.tif")
@@ -303,7 +303,7 @@ end
         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] ==
-        map(sum ∘ skipmissing, mask(st; with=polygon))  
+        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)) == 
@@ -755,4 +755,4 @@ test = rebuild(ga; name = :test)
 
     @test_throws "strictly positive" Rasters.sample(StableRNG(123), test, 3, skipmissing = true, replace = false)
     @test_throws "Cannot draw" Rasters.sample(StableRNG(123), test, 5, replace = false)
-end
\ No newline at end of file
+end
diff --git a/test/rasterize.jl b/test/rasterize.jl
index d975396d8..6d91f38f8 100644
--- a/test/rasterize.jl
+++ b/test/rasterize.jl
@@ -198,7 +198,7 @@ end
             @test eltype(rst) == @NamedTuple{fill1::Union{Missing,Int64}, fill2::Union{Missing,Float32}}
             @test keys(rst) == (:fill1, :fill2)
             @test dims(rst) == dims(A)
-            @test map(sum ∘ skipmissing, rst) === (fill1=15, fill2=30.0f0)
+            @test maplayers(sum ∘ skipmissing, rst) === (fill1=15, fill2=30.0f0)
         end
         @testset "Single value fill makes an array (ignoring table vals)" begin
             ra = rasterize(sum, data; to=A, fill=0x03, missingval=0x00)
@@ -216,7 +216,7 @@ end
             @test parent(rst.val2) isa Array{Union{Missing,Float32},2}
             @test keys(rst) == (:val1, :val2)
             @test dims(rst) == dims(A)
-            @test map(sum ∘ skipmissing, rst) === (val1=4, val2=8.0f0)
+            @test maplayers(sum ∘ skipmissing, rst) === (val1=4, val2=8.0f0)
         end
         @testset "Symbol fill makes an array" begin
             ra = rasterize(feature; to=A, fill=:val1)
@@ -235,7 +235,7 @@ end
                 @test parent(rst.val1) isa Array{Union{Missing,Int},2}
                 @test parent(rst.val2) isa Array{Union{Missing,Float32},2}
                 @test keys(rst) == (:val1, :val2)
-                @test map(sum ∘ skipmissing, rst) === (val1=15, val2=30.0f0)
+                @test maplayers(sum ∘ skipmissing, rst) === (val1=15, val2=30.0f0)
                 @test_throws ArgumentError rasterize(data; to=A, fill=(:val1, :not_a_column), threaded)
             end
             @testset "Symbol fill makes an array" begin
diff --git a/test/runtests.jl b/test/runtests.jl
index 589bf472d..e4337545c 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,16 +1,15 @@
 using Rasters, Test, Aqua, SafeTestsets
 
-if VERSION >= v"1.9.0"
-    # Aqua.test_ambiguities([Rasters, Base, Core])
+@testset "Aqua" begin
+    Aqua.test_ambiguities(Rasters)
     Aqua.test_unbound_args(Rasters)
     Aqua.test_stale_deps(Rasters)
     Aqua.test_undefined_exports(Rasters)
     Aqua.test_project_extras(Rasters)
     # Aqua.test_deps_compat(Rasters) # This breaks GDAL downstream tests
-    # Aqua.test_project_toml_formatting(Rasters) # This seems to change between versions for extensions
-    @time @safetestset "extensions" begin include("extensions.jl") end
 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
diff --git a/test/sources/gdal.jl b/test/sources/gdal.jl
index efa636d6d..4ff6a813d 100644
--- a/test/sources/gdal.jl
+++ b/test/sources/gdal.jl
@@ -624,20 +624,20 @@ end
 
     @testset "methods" begin
         @testset "mean" begin
-            means = map(A -> mean(parent(A); dims=2), gdalstack)
-            @test map((a, b) -> all(a .== b), mean(gdalstack; dims=Y), means) |> all
+            means = maplayers(A -> mean(parent(A); dims=2), gdalstack)
+            @test maplayers((a, b) -> all(a .== b), mean(gdalstack; dims=Y), means) |> all
         end
         @testset "trim, crop, extend" begin
             mv = zero(eltype(gdalstack[:a]))
             st = read(replace_missing(gdalstack, mv))
-            st = map(A -> (view(A, X(1:100)) .= mv; A), st)
+            st = maplayers(A -> (view(A, X(1:100)) .= mv; A), st)
             trimmed = trim(st)
             @test size(trimmed) == (414, 514)
             cropped = crop(st; to=trimmed)
             @test size(cropped) == (414, 514)
-            @test map((c, t) -> all(collect(c .=== t)), cropped, trimmed) |> all
+            @test all(maplayers((c, t) -> all(collect(c .=== t)), cropped, trimmed))
             extended = extend(read(cropped); to=st)
-            @test all(map((s, e) -> all(s .=== e), st, extended))
+            @test all(maplayers((s, e) -> all(s .=== e), st, extended))
         end
         @testset "mask and mask!" begin
             st = read(gdalstack)
@@ -859,7 +859,7 @@ end
         read!([(a=gdalpath, b=gdalpath), (a=gdalpath, b=gdalpath)], ser2)
         read!(ser1, ser3)
         @test map(ser1, ser2, ser3) do st1, st2, st3
-            map(st1, st2, st3) do A1, A2, A3
+            maplayers(st1, st2, st3) do A1, A2, A3
                 (A2 .=== A2 .=== A3) |> all
             end |> all
         end |> all
diff --git a/test/sources/ncdatasets.jl b/test/sources/ncdatasets.jl
index 2fa3a24ea..225fe7fcd 100644
--- a/test/sources/ncdatasets.jl
+++ b/test/sources/ncdatasets.jl
@@ -461,7 +461,7 @@ end
         cp(ncmulti, ncmulti_custom, force=true)
         @time ncstack_custom = RasterStack(ncmulti_custom, source=Rasters.NCDsource())
         @test ncstack_custom isa RasterStack
-        @test map(read(ncstack_custom), read(ncstack)) do a, b
+        @test maplayers(read(ncstack_custom), read(ncstack)) do a, b
             all(a .=== b)
         end |> all
     end