diff --git a/ext/RastersArchGDALExt/reproject.jl b/ext/RastersArchGDALExt/reproject.jl
index 384ef332d..e6046d999 100644
--- a/ext/RastersArchGDALExt/reproject.jl
+++ b/ext/RastersArchGDALExt/reproject.jl
@@ -4,6 +4,7 @@ function Rasters._reproject(source::GeoFormat, target::GeoFormat, ::XDim, vals::
     push!(paired_vals, (first(vals), one(first(vals))))
     reps = AG.reproject(paired_vals, source, target; order=:trad)
     at_one = pop!(reps)
+    # TODO is this test sufficient to reject the reprojection?
     first(reps)[1] == at_one[1] || _reproject_crs_error(source, target)
     return map(r -> r[1], reps)
 end
@@ -16,4 +17,5 @@ function Rasters._reproject(source::GeoFormat, target::GeoFormat, dim::YDim, val
     return map(r -> r[2], reps)
 end
 
-_reproject_crs_error(source, target) = throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target")) 
+_reproject_crs_error(source, target) = 
+    throw(ArgumentError("Cannot reproject from: \n $source \nto: \n $target, coordinate reference systems are not aligned with on all axes. You may need to use `resample` instead")) 
diff --git a/ext/RastersArchGDALExt/resample.jl b/ext/RastersArchGDALExt/resample.jl
index f7fa32edc..6bb610ae3 100644
--- a/ext/RastersArchGDALExt/resample.jl
+++ b/ext/RastersArchGDALExt/resample.jl
@@ -2,10 +2,16 @@
 	resample(x; kw...)
     resample(xs...; to=first(xs), kw...)
 
-`resample` uses `ArchGDAL.gdalwarp` to resample a [`Raster`](@ref) or
-[`RasterStack`](@ref) to a new `resolution` and optionally new `crs`,
+`resample` uses `warp` (which uses GDALs `gdalwarp`) to resample a [`Raster`](@ref)
+or [`RasterStack`](@ref) to a new `resolution` and optionally new `crs`,
 or to snap to the bounds, resolution and crs of the object `to`.
 
+Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
+are iteratively resampled with GDAL and joined back int a single array.
+
+If projections can be converted for each axis independently, it may 
+be faster and more accurate to use [`reproject`](@ref).
+
 # Arguments
 
 - `x`: the object/s to resample.
@@ -39,8 +45,8 @@ $FILENAME_KEYWORD
 $SUFFIX_KEYWORD
 
 Note:
-- GDAL may cause some unexpected changes in the data, such as returning a reversed Y dimension or
-  changing the `crs` type from `EPSG` to `WellKnownText` (it will represent the same CRS).
+- GDAL may cause some unexpected changes in the raster, such as changing the `crs`
+    type from `EPSG` to `WellKnownText` (it will represent the same CRS).
 
 # Example
 
diff --git a/ext/RastersNCDatasetsExt/ncdatasets_source.jl b/ext/RastersNCDatasetsExt/ncdatasets_source.jl
index 8c4acfe0f..45793467b 100644
--- a/ext/RastersNCDatasetsExt/ncdatasets_source.jl
+++ b/ext/RastersNCDatasetsExt/ncdatasets_source.jl
@@ -201,7 +201,7 @@ function DD.layermetadata(ds::NCD.Dataset)
         end
         md
     end
-    NamedTuple{map(Symbol, keys)}(dimtypes)
+    return NamedTuple{map(Symbol, keys)}(dimtypes)
 end
 
 RA.missingval(var::NCD.CFVariable{T}) where T = missing isa T ? missing : nothing
@@ -209,40 +209,22 @@ RA.missingval(var::NCD.Dataset) = missing
 
 function RA.layerkeys(ds::NCD.Dataset)
     dimkeys = _dimkeys(ds)
-    toremove = if "bnds" in dimkeys
-        dimkeys = setdiff(dimkeys, ("bnds",))
-        boundskeys = String[]
-        for k in dimkeys
-            var = NCD.variable(ds, k)
-            if haskey(var.attrib, "bounds")
-                push!(boundskeys, var.attrib["bounds"])
-            end
-        end
-        union(dimkeys, boundskeys)::Vector{String}
-    else
-        dimkeys::Vector{String}
-    end
-    nondim = setdiff(keys(ds), toremove)
-    grid_mapping = String[]
-    for k in nondim
-        var = NCD.variable(ds, k)
-        if haskey(var.attrib, "grid_mapping")
-            push!(grid_mapping, var.attrib["grid_mapping"])
-        end
+    return filter(keys(ds)) do key
+        key in dimkeys && return false
+        endswith(lowercase(key), r"bnds|bounds") return false
+        !haskey(NCD.variable(ds, k).attrib, "grid_mapping")
     end
-    nondim = setdiff(nondim, grid_mapping)
 end
 
 function RA.FileStack{NCDsource}(ds::NCD.Dataset, filename::AbstractString; write=false, keys)
     keys = map(Symbol, keys isa Nothing ? RA.layerkeys(ds) : keys) |> Tuple
-    type_size_ec_hc = map(keys) do key
-        var = ds[string(key)]
-        Union{Missing,eltype(var)}, size(var), _ncd_eachchunk(var), _ncd_haschunks(var)
+    vars = map(keys) do key
+        ds[string(key)]
     end
-    layertypes = map(x->x[1], type_size_ec_hc)
-    layersizes = map(x->x[2], type_size_ec_hc)
-    eachchunk = map(x->x[3], type_size_ec_hc)
-    haschunks = map(x->x[4], type_size_ec_hc)
+    layertypes = map(v -> Union{Missing,eltype(v)}, vars)
+    layersizes = map(size, vars)
+    eachchunk = map(_ncd_eachchunk, vars)
+    haschunks = map(_ncd_haschunks, vars)
     return RA.FileStack{NCDsource,keys}(filename, layertypes, layersizes, eachchunk, haschunks, write)
 end
 function RA.OpenStack(fs::RA.FileStack{NCDsource,K}) where K
diff --git a/src/methods/crop_extend.jl b/src/methods/crop_extend.jl
index b153c2a64..911b90243 100644
--- a/src/methods/crop_extend.jl
+++ b/src/methods/crop_extend.jl
@@ -13,6 +13,7 @@ to match the size of the object `to`, or smallest of any dimensions that are sha
     area of all `xs` is used.
 - `touches`: `true` or `false`. Whether to use `Touches` wraper on the object extent.
    When lines need to be included in e.g. zonal statistics, `true` should be used.
+- `boundary`: `:inside`, `:center` or `:touches`
 
 As `crop` is lazy, `filename` and `suffix` keywords are not used.
 
@@ -100,11 +101,12 @@ function _crop_to(x, to::DimTuple; kw...)
     end
     return _crop_to(x, Extents.extent(sampled); kw...)
 end
-function _crop_to(x, to::Extents.Extent; touches=false)
+function _crop_to(x, to::Extents.Extent; touches=false, boundary=:inside)
     # Take a view over the bounds
     _without_mapped_crs(x) do x1
         if touches 
             view(x1, Touches(to))
+            e
         else
             view(x1, to)
         end
@@ -124,7 +126,7 @@ covered by all `xs`, or by the keyword argument `to`.
 - `to`: the Raster or dims to extend to. If no `to` keyword is passed, the largest
     shared area of all `xs` is used.
 - `touches`: `true` or `false`. Whether to use `Touches` wraper on the object extent.
-   When lines need to be included in e.g. zonal statistics, `true` shoudle be used.
+   When lines need to be included in e.g. zonal statistics, `true` should be used.
 $FILENAME_KEYWORD
 $SUFFIX_KEYWORD
 
@@ -215,6 +217,7 @@ function _extend_to(x::RasterStackOrArray, extent::Extents.Extent{K}; kw...) whe
     shareddims = dims(x, dims(extent))
     bnds = map(val, dims(extent, shareddims))
     newdims = map(shareddims, bnds) do d, b
+        span(dim) isa Regular || throw(ArgumentError("Can only `extend` rasters with `Regular` lookups"))
         l = lookup(d)
         # Use ranges for math because they have TwicePrecision magic
         # Define a range down to the lowest value,
diff --git a/src/methods/reproject.jl b/src/methods/reproject.jl
index 75f6df34b..c9f1a10be 100644
--- a/src/methods/reproject.jl
+++ b/src/methods/reproject.jl
@@ -1,17 +1,26 @@
 
 """
-    reproject(target::GeoFormat, x)
+    reproject(obj; crs)
 
-Reproject the dimensions of `x` to a different crs.
+Reproject the lookups of `obj` to a different crs. 
 
-# Arguments
+This is a lossless operation for the raster data, as only the 
+lookup values change. This is only possible when the axes of source
+and destination projections are alligned: the change is usually from
+a [`Regular`](@ref) and an [`Irregular`](@ref) lookup spans.
 
-- `target`: any crs in a GeoFormatTypes.jl wrapper, e.g. `EPSG`, `WellKnownText`, `ProjString`.
-- `x`: a `Dimension`, `Tuple` of `Dimension`, `Raster` or `RasterStack`.
+For converting between projections that are rotated, 
+skewed or warped in any way, `resample` use [`resample`](@ref).
 
 Dimensions without an `AbstractProjected` lookup (such as a `Ti` dimension)
 are silently returned without modification.
+
+# Arguments
+
+- `obj`: a `LookupArray`, `Dimension`, `Tuple` of `Dimension`, `Raster` or `RasterStack`.
+$CRS_KEYWORD
 """
+reproject(x; crs::GeoFormat) = reproject(crs, x)
 reproject(target::GeoFormat, x) = rebuild(x; dims=reproject(target, dims(x)))
 reproject(target::GeoFormat, dims::Tuple) = map(d -> reproject(target, d), dims)
 reproject(target::GeoFormat, l::LookupArray) = l
diff --git a/test/reproject.jl b/test/reproject.jl
index eb84cfecf..01cb07723 100644
--- a/test/reproject.jl
+++ b/test/reproject.jl
@@ -14,7 +14,6 @@ using Rasters: reproject, convertlookup
     @test reproject(wktcea, EPSG(4326), Y(), -3.658789324855012e6) ≈ -30.0
     @test reproject(projcea, wkt4326, Y(), -3.658789324855012e6) ≈ -30.0
     @test reproject(cea, proj4326, Y(), [-3.658789324855012e6]) ≈ [-30.0]
-    
 
     @test reproject(proj4326, cea, X(), 180.0) ≈ 1.7367530445161372e7
     @test reproject(cea, EPSG(4326), X(), 1.7367530445161372e7) ≈ 180.0
@@ -23,15 +22,23 @@ using Rasters: reproject, convertlookup
 
     x = X(Projected(0.0:1.0:360.0; crs=EPSG(4326), dim=X(), order=ForwardOrdered(), span=Regular(1.0), sampling=Intervals(Start())))
     y = Y(Projected(-80.0:1.0:80.0; crs=EPSG(4326), dim=Y(), order=ReverseOrdered(), span=Regular(1.0), sampling=Intervals(Start())))
-    x1, y1 = reproject(projcea, (x, y))
+
+    x1, y1 = reproject((x, y); crs=projcea)
     @test span(x1) isa Irregular
     @test span(y1) isa Irregular
-    x2, y2 = reproject(EPSG(4326), (x, y))
+    x2, y2 = reproject((x, y); crs=EPSG(4326))
     @test all(x .== x2)
     @test all(y .== y2)
     @test span(x2) == Regular(1.0)
     @test span(y2) == Regular(1.0)
 
+    raster = rand(x, y)
+    reprojected_raster = reproject(raster; crs=projcea)
+    # Dims have changed
+    @test dims(reprojected_raster) == (x1, y1)
+    # Raster has not changed
+    @test reprojected_raster == raster
+
     @test_throws ArgumentError reproject(cea, EPSG(32618), Y(), [-3.658789324855012e6])
     @test_throws ArgumentError reproject(cea, EPSG(32618), X(), [-3.658789324855012e6])
 end