From f45f22dc50226457c97c66231902aee31442f702 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 28 Oct 2023 18:12:20 +0200
Subject: [PATCH 01/13] add resample test file

---
 docs/src/index.md                     |   2 +-
 ext/RastersArchGDALExt/cellsize.jl    |  23 ---
 ext/RastersArchGDALExt/gdal_source.jl | 257 ++++++++++++--------------
 ext/RastersArchGDALExt/resample.jl    |  97 +---------
 ext/RastersArchGDALExt/warp.jl        |  98 ++++------
 ext/RastersMakieExt/plotrecipes.jl    |  23 ---
 src/array.jl                          |   1 +
 src/extensions.jl                     | 164 ++++++++++++++++
 src/plotrecipes.jl                    |  23 +++
 src/utils.jl                          |  10 +-
 test/methods.jl                       | 122 ------------
 test/runtests.jl                      |   1 +
 test/sources/gdal.jl                  |  11 +-
 test/sources/ncdatasets.jl            |  19 +-
 14 files changed, 363 insertions(+), 488 deletions(-)

diff --git a/docs/src/index.md b/docs/src/index.md
index 38640e65..0017ba13 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -82,4 +82,4 @@ Other functionality in extensions:
     3. Paste the complete stack trace of the error it produces, right to the bottom, into the issue template. Then we can be sure we reproduced the same problem.
 
     Good issues are really appreciated, but they do take just a little extra effort with Rasters.jl because of this need for files.
-```
\ No newline at end of file
+```
diff --git a/ext/RastersArchGDALExt/cellsize.jl b/ext/RastersArchGDALExt/cellsize.jl
index 833008d3..ab80a876 100644
--- a/ext/RastersArchGDALExt/cellsize.jl
+++ b/ext/RastersArchGDALExt/cellsize.jl
@@ -1,26 +1,3 @@
-"""
-    cellsize(x)
-
-Gives the approximate size of each cell in square km. 
-This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.
-
-# Arguments
-
-- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions. 
-
-## Example
-
-```jldoctest
-using Rasters, Rasters.LookupArrays
-dimz = X(Projected(90.0:0.1:120; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))),
-       Y(Projected(0.0:0.1:50; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
-
-cs = cellsize(dimz)
-
-```
-$EXPERIMENTAL
-"""
-
 function _great_circle_bearing(lon1::AbstractFloat, lat1::AbstractFloat, lon2::AbstractFloat, lat2::AbstractFloat)
     dLong = lon1 - lon2
 
diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index 54b21a06..ee88ad84 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -1,10 +1,6 @@
 const AG = ArchGDAL
 
-const GDAL_X_ORDER = ForwardOrdered()
-const GDAL_Y_ORDER = ReverseOrdered()
-
-const GDAL_X_LOCUS = Start()
-const GDAL_Y_LOCUS = Start()
+const GDAL_LOCUS = Start()
 
 # drivers supporting the gdal Create() method to directly write to disk
 const GDAL_DRIVERS_SUPPORTING_CREATE = ("GTiff", "HDF4", "KEA", "netCDF", "PCIDSK", "Zarr", "MEM"#=...=#) 
@@ -62,81 +58,38 @@ Write a `Raster` to file using GDAL.
 Returns `filename`.
 """
 function Base.write(
-    filename::AbstractString, ::Type{GDALsource}, A::AbstractRaster{T,2}; 
-    force=false, verbose=true, kw...
-) where T
-    RA.check_can_write(filename, force)
-    all(hasdim(A, (X, Y))) || error("Array must have Y and X dims")
-
-    correctedA = _maybe_correct_to_write(A)
-    nbands = 1
-    _write(filename, correctedA, nbands; kw...)
-end
-function Base.write(
-    filename::AbstractString, ::Type{GDALsource}, A::AbstractRaster{T,3}; 
+    filename::AbstractString, ::Type{GDALsource}, A::AbstractRaster{T}; 
     force=false, verbose=true, kw...
 ) where T
     RA.check_can_write(filename, force)
-    all(hasdim(A, (X, Y))) || error("Array must have Y and X dims")
-    hasdim(A, Band()) || error("Must have a `Band` dimension to write a 3-dimensional array")
-
-    correctedA = _maybe_correct_to_write(A)
-    nbands = size(correctedA, Band())
-    _write(filename, correctedA, nbands; kw...)
+    A1 = _maybe_correct_to_write(A)
+    _create_with_driver(filename, dims(A1), eltype(A1), missingval(A1); _block_template=A1, kw...) do dataset
+        _maybe_warn_south_up(A, verbose, "Writing South-up. Use `reverse(raster; dims=Y)` first to write conventional North-up")
+        open(A1; write=true) do O
+            AG.RasterDataset(dataset) .= parent(O)
+        end
+    end
+    return filename
 end
 
 function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
-    missingval=nothing, metadata=nothing, name=nothing, keys=(name,),
-    driver="", lazy=true, options=Dict{String,String}(), _block_template=nothing
+    missingval=nothing, metadata=nothing, name=nothing, lazy=true, verbose=true, kw...
 )
-    driver = _check_driver(filename, driver)
-    if !(keys isa Nothing || keys isa Symbol) && length(keys) > 1
-        throw(ArgumentError("GDAL cant write more than one layer per file, but keys $keys have $(length(keys))"))
-    end
-    x, y = map(DD.dims(dims, (XDim, YDim)), (GDAL_X_ORDER, GDAL_Y_ORDER)) do d, o
-        reorder(lookup(d) isa NoLookup ? set(d, Sampled) : d, o)
-    end
     T = Missings.nonmissingtype(T)
-
-    if ismissing(missingval)
-        missingval = RA._writeable_missing(T)
-    end
-
-    if hasdim(dims, Band)
-        b = DD.dims(dims, Band)
-        nbands = length(b)
-        newdims = (x, y, b)
-    else
-        nbands = 1
-        newdims = (x, y)
+    missingval = ismissing(missingval) ? RA._writeable_missing(T) : missingval
+    _create_with_driver(filename, dims, T, missingval; kw...) do _
+        _maybe_warn_south_up(A, verbose, "Creating South-up. Use `reverse(ydim)` first to write conventional North-up")
+        nothing
     end
 
-    kw = (width=length(x), height=length(y), nbands=nbands, dtype=T)
-    options_vec = _process_options(driver, options; _block_template)
-    gdaldriver = driver isa String ? AG.getdriver(driver) : driver
-    if driver in GDAL_DRIVERS_SUPPORTING_CREATE
-        AG.create(filename; driver=gdaldriver, options=options_vec, kw...) do ds
-            _set_dataset_properties!(ds, newdims, missingval)
-        end
-    else
-        tif_options_vec = _process_options("GTiff", Dict{String,String}(); _block_template)
-        # Create a tif and copy it to `filename`, as ArchGDAL.create
-        # does not support direct creation of ASCII etc. rasters
-        ArchGDAL.create(tempname() * ".tif"; driver=AG.getdriver("GTiff"), options=tif_options_vec, kw...) do ds
-            _set_dataset_properties!(ds, newdims, missingval)
-            target_ds = AG.copy(ds; filename=filename, driver=gdaldriver, options=options_vec)
-            AG.destroy(target_ds)
-        end
-    end
-    if hasdim(dims, Band)
-        return Raster(filename; source=GDALsource, lazy)
-    else
-        return view(Raster(filename; source=GDALsource, lazy), Band(1))
-    end
+    return Raster(filename; source=GDALsource, name, lazy, dropband=!hasdim(dims, Band))
 end
 
+_maybe_warn_south_up(A, verbose, msg) = 
+    verbose && lookup(A, Y, msg) isa AbstractProjected && order(A, Y, msg) isa ForwardOrdered && @warn msg
+
 function RA._open(f, ::Type{GDALsource}, filename::AbstractString; write=false, kw...)
-    # Check the file actually exists because GDALs error is unhelpful
+    # Check the file actually exists because the GDAL error is unhelpful
     if !isfile(filename)
         # Allow gdal virtual file systems
         # the respective string is prepended to the data source,
@@ -168,11 +121,13 @@ RA._open(f, ::Type{GDALsource}, ds::AG.RasterDataset; kw...) = RA.cleanreturn(f(
 
 # We allow passing in crs and mappedcrs manually
 function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
-    gt = try
+    gt_dims = try
         AG.getgeotransform(raster)
     catch
         GDAL_EMPTY_TRANSFORM
     end
+    # @show gt_dims
+    gt = gt_dims
     xsize, ysize = size(raster)
     nbands = AG.nraster(raster)
     bandnames = _bandnames(raster, nbands)
@@ -201,7 +156,7 @@ function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
         xindex = LinRange(xmin, xmax, xsize)
         xorder = xstep > 0 ? ForwardOrdered() : ReverseOrdered()
 
-        ystep = gt[GDAL_NS_RES] # A negative number
+        ystep = gt[GDAL_NS_RES] # Usually a negative number
         if ystep > 0 
             ymax = gt[GDAL_TOPLEFT_Y]
             ymin = gt[GDAL_TOPLEFT_Y] + ystep * (ysize - 1)
@@ -218,7 +173,7 @@ function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
             Points(), Points()
         else
             # GeoTiff uses the "pixelCorner" convention
-            Intervals(GDAL_X_LOCUS), Intervals(GDAL_Y_LOCUS)
+            Intervals(GDAL_LOCUS), Intervals(GDAL_LOCUS)
         end
 
         xlookup = Projected(xindex;
@@ -330,33 +285,11 @@ function AG.Dataset(f::Function, A::AbstractRaster; kw...)
         f(rds.ds)
     end
 end
-function AG.RasterDataset(f::Function, A::AbstractRaster; 
-    filename=nothing, driver="",
-)
-    driver = _check_driver(filename, driver)
-    all(hasdim(A, (X, Y))) || throw(ArgumentError("`AbstractRaster` must have both an `X` and `Y` to be converted to an ArchGDAL `Dataset`"))
-    if ndims(A) === 3
-        thirddim = otherdims(A, (X, Y))[1]
-        thirddim isa Band || throw(ArgumentError("ArchGDAL can't handle $(basetypeof(thirddim)) dims - only X, Y, and Band"))
-    elseif ndims(A) > 3
-        throw(ArgumentError("ArchGDAL can only accept 2 or 3 dimensional arrays"))
-    end
-
-    A_m = RA._maybe_use_type_missingval(A, GDALsource)
-    A_p = _maybe_permute_to_gdal(A_m)
-    kw = (;
-        width=length(DD.dims(A_p, X)),
-        height=length(DD.dims(A_p, Y)),
-        nbands=hasdim(A_p, Band) ? size(A_p, Band()) : 1,
-        dtype=eltype(A_p),
-    )
-    if filename == nothing 
-        filename = ""
-    end
-    _create_with_driver(filename, driver, kw; _block_template=A_p) do dataset
-        _set_dataset_properties!(dataset, dims(A_p), missingval(A_p))
+function AG.RasterDataset(f::Function, A::AbstractRaster; filename="", kw...)
+    A1 = _maybe_correct_to_write(A)
+    return _create_with_driver(filename, dims(A1), eltype(A1), missingval(A1); _block_template=A1, kw...) do dataset
         rds = AG.RasterDataset(dataset)
-        open(A_p) do a
+        open(A1) do a
             rds .= parent(a)
         end
         f(rds)
@@ -389,32 +322,19 @@ _missingval_from_gdal(T, x) = x
 _maybe_correct_to_write(A) = _maybe_correct_to_write(lookup(A, X()), A)
 _maybe_correct_to_write(lookup, A) = A
 function _maybe_correct_to_write(lookup::Union{AbstractSampled,NoLookup}, A)
-    _maybe_permute_to_gdal(A) |>
-        a -> RA.noindex_to_sampled(a) |>
-        a -> reorder(a, (X(GDAL_X_ORDER), Y(GDAL_Y_ORDER)))
-end
-
-# Write a Raster to disk using GDAL
-function _write(filename, A::AbstractRaster, nbands; driver="", kw...)
-    A = RA._maybe_use_type_missingval(A, GDALsource)
-    create_kw = (width=size(A, X()), height=size(A, Y()), nbands=nbands, dtype=eltype(A))
-    _create_with_driver(filename, driver, create_kw; _block_template=A, kw...) do dataset
-        _set_dataset_properties!(dataset, A)
-        rds = AG.RasterDataset(dataset)
-        open(A; write=true) do O
-            rds .= parent(O)
-        end
-    end
-
-    return filename
+    RA._maybe_use_type_missingval(A, GDALsource) |> _maybe_permute_to_gdal
 end
 
 _check_driver(filename::Nothing, driver) = "MEM"
 function _check_driver(filename::AbstractString, driver)
-    if isempty(driver) 
-        driver = AG.extensiondriver(filename)
-        if driver == "COG"
-            driver = "GTiff"
+    if isempty(driver)
+        if isempty(filename)
+            driver = "MEM"
+        else
+            driver = AG.extensiondriver(filename)
+            if driver == "COG"
+                driver = "GTiff"
+            end
         end
     end
     return driver
@@ -422,32 +342,56 @@ end
 
 # Handle creating a dataset with any driver, 
 # applying the function `f` to the created dataset
-function _create_with_driver(f, filename, driver, create_kw;
-    options=Dict{String,String}(), _block_template=nothing
+function _create_with_driver(f, filename, dims, T, missingval;
+    options=Dict{String,String}(), driver="", _block_template=nothing, kw...
 )
+    _gdal_validate(dims)
+
+    x, y = map(DD.dims(dims, (XDim, YDim))) do d
+        maybeshiftlocus(Start(), RA.nolookup_to_sampled(d))
+    end
+    newdims = hasdim(dims, Band()) ? (x, y, DD.dims(dims, Band)) : (x, y) 
+    nbands = hasdim(dims, Band) ? length(DD.dims(dims, Band())) : 1
+
     driver = _check_driver(filename, driver)
     options_vec = _process_options(driver, options; _block_template)
     gdaldriver = driver isa String ? AG.getdriver(driver) : driver
+
+    create_kw = (; width=length(x), height=length(y), nbands, dtype=T,)
+    filename = isnothing(filename) ? "" : filename
+
     if AG.shortname(gdaldriver) in GDAL_DRIVERS_SUPPORTING_CREATE
-        AG.create(filename; driver=gdaldriver, create_kw..., options=options_vec) do dataset
+        AG.create(filename; driver=gdaldriver, options=options_vec, create_kw...) do dataset
+            _set_dataset_properties!(dataset, newdims, missingval)
             f(dataset)
         end
     else
-        # Create a memory object and copy it to disk, as ArchGDAL.create
+        # Create a tif and copy it to `filename`, as ArchGDAL.create
         # does not support direct creation of ASCII etc. rasters
-        ArchGDAL.create(""; driver=AG.getdriver("MEM"), create_kw...) do dataset
-            result = f(dataset)
-            # This `copy` copies _to disk_
-            AG.copy(dataset; filename=filename, driver=gdaldriver, options=options_vec) |> AG.destroy
-            result
+        tif_options_vec = _process_options("GTiff", Dict{String,String}(); _block_template)
+        tif_driver = AG.getdriver("GTiff")
+        tif_name = tempname() * ".tif"
+        AG.create(tif_name; driver=tif_driver, options=tif_options_vec, create_kw...) do dataset
+            _set_dataset_properties!(dataset, newdims, missingval)
+            f(dataset)
+            target_ds = AG.copy(dataset; filename=filename, driver=gdaldriver, options=options_vec)
+            AG.destroy(target_ds)
         end
     end
 end
 
+@noinline function _gdal_validate(dims)
+    all(hasdim(dims, (XDim, YDim))) || throw(ArgumentError("`Raster` must have both an `X` and `Y` to be converted to an ArchGDAL `Dataset`"))
+    if length(dims) === 3
+        otherdim = otherdims(dims, (XDim, YDim))[1]
+        otherdim isa Band || throw(ArgumentError("ArchGDAL can't handle $(basetypeof(thirddim)) dims - only X, Y, and Band"))
+    elseif !(length(dims) in (2, 3))
+        throw(ArgumentError("ArchGDAL can only accept 2 or 3 dimensional arrays"))
+    end
+end
+
 # Convert a Dict of options to a Vector{String} for GDAL
-function _process_options(driver::String, options::Dict;
-    _block_template=nothing
-)
+function _process_options(driver::String, options::Dict; _block_template=nothing)
     options_str = Dict(string(k)=>string(v) for (k,v) in options)
     # Get the GDAL driver object
     gdaldriver = AG.getdriver(driver)
@@ -523,19 +467,37 @@ end
 _set_dataset_properties!(ds::AG.Dataset, A) = 
     _set_dataset_properties!(ds, dims(A), missingval(A))
 function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval)
+    # We cant write mixed Points/Intervals, so default to Intervals if mixed
+    xy = DD.dims(dims, (X, Y))
+    if any(x -> x isa Intervals, sampling.(xy)) && any(x -> x isa Points, sampling.(xy))
+        dims = set(dims, X => Intervals, Y => Intervals)
+    end
     # Convert the dimensions to `Projected` if they are `Converted`
-    # This allows saving NetCDF to Tiff
+    # This allows saving NetCDF to Tiff.
+    x = convertlookup(Projected, DD.dims(dims, X))
+    y = convertlookup(Projected, DD.dims(dims, Y))
+
     # Set the index loci to the start of the cell for the lat and lon dimensions.
     # NetCDF or other formats use the center of the interval, so they need conversion.
-    x = DD.maybeshiftlocus(GDAL_X_LOCUS, convertlookup(Projected, DD.dims(dims, X)))
-    y = DD.maybeshiftlocus(GDAL_Y_LOCUS, convertlookup(Projected, DD.dims(dims, Y)))
-    # Convert crs to WKT if it exists
+    x = DD.maybeshiftlocus(GDAL_LOCUS, x)
+    y = DD.maybeshiftlocus(GDAL_LOCUS, y)
+
+    # Set GDAL AREA_OR_POINT metadata
+    area_or_point = sampling(x) isa Points ? "Point" : "Area" 
+    AG.GDAL.gdalsetmetadataitem(dataset, "AREA_OR_POINT", area_or_point, "")
+
+    # Set crs if it exists, converting crs to WKT
     if !isnothing(crs(x))
         AG.setproj!(dataset, convert(String, convert(WellKnownText, crs(x))))
     end
-    # Get the geotransform from the updated lat/lon dims and write
-    AG.setgeotransform!(dataset, RA.dims2geotransform(x, y))
 
+    # Set the geotransform from the updated lookups
+    gt = RA.dims2geotransform(x, y)
+    # @show gt
+    AG.setgeotransform!(dataset, gt)
+
+    # Set the missing value/nodataval. This is a little complicated
+    # because gdal has separate method for 64 bit integers
     if !isnothing(missingval)
         bands = hasdim(dims, Band) ? axes(DD.dims(dims, Band), 1) : 1
         for i in bands
@@ -552,7 +514,7 @@ function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval)
 
     # Write band labels if they are not Integers.
     if hasdim(dims, Band)
-        bandlookup = DD.lookup(dims, Band)
+        bandlookup = DD.lookup(DD.dims(dims, Band))
         if !(eltype(bandlookup) <: Integer)
             for i in eachindex(bandlookup)
                 AG.getband(dataset, i) do band
@@ -585,9 +547,9 @@ _maybe_permute_to_gdal(A, dims::Tuple) = A
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = permutedims(A, dims)
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim}) = permutedims(A, dims)
 
-_maybe_permute_from_gdal(A, dims::Tuple) = permutedims(A, dims)
-_maybe_permute_from_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = A
-_maybe_permute_from_gdal(A, dims::Tuple{<:XDim,<:YDim}) = A
+_maybe_restore_from_gdal(A, dims::Tuple) = reorder(permutedims(A, dims), dims) 
+_maybe_restore_from_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = reorder(A, dims)
+_maybe_restore_from_gdal(A, dims::Tuple{<:XDim,<:YDim}) = reorder(A, dims)
 
 #= Geotranforms ########################################################################
 
@@ -608,13 +570,22 @@ const USING_COORDINATETRANSFORMATIONS_MESSAGE =
     "Run `using CoordinateTransformations` to load affine transformed rasters"
 
 function RA.dims2geotransform(x::XDim, y::YDim)
+    sampling(x) == sampling(y) || throw(ArgumentError("Sampling of x and y must match to generate a geotransform"))
+    (order(x) isa Ordered && order(y) isa Ordered) || throw(ArgumentError("GDAL can only write ordered lookups"))
+    (span(x) isa Regular && span(y) isa Regular) || throw(ArgumentError("GDAL can only write regular lookups"))
     gt = zeros(6)
-    gt[GDAL_TOPLEFT_X] = first(x)
-    gt[GDAL_WE_RES] = step(x)
     gt[GDAL_ROT1] = zero(eltype(gt))
-    gt[GDAL_TOPLEFT_Y] = first(y) - step(y)
     gt[GDAL_ROT2] = zero(eltype(gt))
+    gt[GDAL_WE_RES] = step(x)
     gt[GDAL_NS_RES] = step(y)
+    if sampling(x) isa Points
+        gt[GDAL_TOPLEFT_X] = first(x)
+        gt[GDAL_TOPLEFT_Y] = first(y)
+    else
+        sampling(x) isa Intervals{Start} || throw(ArgumentError("GDAL can only write intervals with Start locus"))
+        gt[GDAL_TOPLEFT_X] = order(x) isa ReverseOrdered ? first(x) - step(x) : first(x)
+        gt[GDAL_TOPLEFT_Y] = order(y) isa ReverseOrdered ? first(y) - step(y) : first(y)
+    end
     return gt
 end
 RA.geotransform2affine(gt) = error(USING_COORDINATETRANSFORMATIONS_MESSAGE)
diff --git a/ext/RastersArchGDALExt/resample.jl b/ext/RastersArchGDALExt/resample.jl
index 67b723c9..da4aa01f 100644
--- a/ext/RastersArchGDALExt/resample.jl
+++ b/ext/RastersArchGDALExt/resample.jl
@@ -1,81 +1,3 @@
-"""
-	resample(x; kw...)
-    resample(xs...; to=first(xs), kw...)
-
-`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 into 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.
-
-# Keywords
-
-- `to`: a `Raster`, `RasterStack`, `Tuple` of `Dimension` or `Extents.Extent`.
-    If no `to` object is provided the extent will be calculated from `x`,
-$RES_KEYWORD
-$SIZE_KEYWORD
-$CRS_KEYWORD
-- `method`: A `Symbol` or `String` specifying the method to use for resampling.
-    From the docs for [`gdalwarp`](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r):
-    * `:near`: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
-    * `:bilinear`: bilinear resampling.
-    * `:cubic`: cubic resampling.
-    * `:cubicspline`: cubic spline resampling.
-    * `:lanczos`: Lanczos windowed sinc resampling.
-    * `:average`: average resampling, computes the weighted average of all non-NODATA contributing pixels.
-        rms root mean square / quadratic mean of all non-NODATA contributing pixels (GDAL >= 3.3)
-    * `:mode`: mode resampling, selects the value which appears most often of all the sampled points.
-    * `:max`: maximum resampling, selects the maximum value from all non-NODATA contributing pixels.
-    * `:min`: minimum resampling, selects the minimum value from all non-NODATA contributing pixels.
-    * `:med`: median resampling, selects the median value of all non-NODATA contributing pixels.
-    * `:q1`: first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels.
-    * `:q3`: third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels.
-    * `:sum`: compute the weighted sum of all non-NODATA contributing pixels (since GDAL 3.1)
-
-    Where NODATA values are set to `missingval`.
-$FILENAME_KEYWORD
-$SUFFIX_KEYWORD
-
-Note:
-- 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
-
-Resample a WorldClim layer to match an EarthEnv layer:
-
-```jldoctest
-using Rasters, RasterDataSources, ArchGDAL, Plots
-A = Raster(WorldClim{Climate}, :prec; month=1)
-B = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
-
-a = plot(A)
-b = plot(resample(A; to=B))
-
-savefig(a, "docs/build/resample_example_before.png");
-savefig(b, "docs/build/resample_example_after.png"); nothing
-# output
-```
-
-### Before `resample`:
-
-![before resample](resample_example_before.png)
-
-### After `resample`:
-
-![after resample](resample_example_after.png)
-
-$EXPERIMENTAL
-"""
-function resample end
 resample(x, res; kw...) = resample(x; res, kw...)
 resample(xs::RasterStackOrArray...; kw...) = resample(xs; kw...)
 function resample(ser::AbstractRasterSeries, args...; kw...)
@@ -84,7 +6,7 @@ end
 function resample(xs::Union{Tuple,NamedTuple}; to=first(xs), kw...)
     map(x -> resample(x; to, kw...), xs)
 end
-function resample(x::RasterStackOrArray; 
+function resample(A::RasterStackOrArray; 
     to=nothing, res=nothing, crs=nothing, size=nothing, method=:near, kw...
 )
     (isnothing(size) || isnothing(res)) || _size_and_res_error()
@@ -105,12 +27,6 @@ function resample(x::RasterStackOrArray;
         end
     else
         all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` must have both `XDim` and `YDim` dimensions to resize with GDAL"))
-        if sampling(to, XDim) isa Points
-            to = set(to, dims(to, XDim) => Intervals(Start()))
-        end
-        if sampling(to, YDim) isa Points
-            to = set(to, dims(to, YDim) => Intervals(Start()))
-        end
 
         # Set res from `to` if it was not already set
         if isnothing(res) && isnothing(size)
@@ -127,7 +43,7 @@ function resample(x::RasterStackOrArray;
             nothing
         else
             # get crs from `to` or `x` if none was passed in
-            isnothing(Rasters.crs(to)) ? Rasters.crs(x) : Rasters.crs(to)
+            isnothing(Rasters.crs(to)) ? Rasters.crs(A) : Rasters.crs(to)
         end
     else
         crs
@@ -135,7 +51,7 @@ function resample(x::RasterStackOrArray;
     if !isnothing(crs)
         wkt = convert(String, convert(WellKnownText, crs))
         flags[:t_srs] = wkt
-        if isnothing(Rasters.crs(x))
+        if isnothing(Rasters.crs(A))
             @warn "You have set a crs to resample to, but the object does not have crs so GDAL will assume it is already in the target crs. Use `newraster = setcrs(raster, somecrs)` to fix this."
         end
     end
@@ -151,6 +67,7 @@ function resample(x::RasterStackOrArray;
         else
             throw(ArgumentError("`res` must be a `Real`, or a 2 `Tuple` of `Real` or `Dimension`s wrapping `Real`. Got $res"))
         end
+
         flags[:tr] = [yres, xres]
     end
 
@@ -161,7 +78,7 @@ function resample(x::RasterStackOrArray;
         elseif size isa Tuple{<:Dimension{Int},<:Dimension{Int}}
             map(val, dims(size, (YDim, XDim)))
         elseif size isa Tuple{Int,Int}
-            reverse(size)
+            dimnum(A, XDim) > dimnum(A, YDim) ? size : reverse(size)
         else
             throw(ArgumentError("`size` must be a `Int`, or a 2 `Tuple` of `Int` or `Dimension`s wrapping `Int`. Got $size"))
         end
@@ -169,11 +86,11 @@ function resample(x::RasterStackOrArray;
     end
 
     # resample with `warp`
-    resampled = warp(x, flags; kw...)
+    resampled = warp(A, flags; kw...)
 
     # Return crs to the original type, from GDAL it will always be WellKnownText
     if isnothing(crs)
-        return setcrs(resampled, Rasters.crs(x))
+        return setcrs(resampled, Rasters.crs(A))
     else
         return setcrs(resampled, crs)
     end
diff --git a/ext/RastersArchGDALExt/warp.jl b/ext/RastersArchGDALExt/warp.jl
index 1559cf47..8433d92b 100644
--- a/ext/RastersArchGDALExt/warp.jl
+++ b/ext/RastersArchGDALExt/warp.jl
@@ -1,59 +1,3 @@
-"""
-    warp(A::AbstractRaster, flags::Dict; kw...)
-
-Gives access to the GDALs `gdalwarp` method given a `Dict` of 
-`flag => value` arguments that can be converted to strings, or vectors
-where multiple space-separated arguments are required.
-
-Arrays with additional dimensions not handled by GDAL (other than `X`, `Y`, `Band`)
-are sliced, warped, and then combined to match the original array dimensions. 
-These slices will *not* be written to disk and loaded lazily at this stage -
-you will need to do that manually if required.
-
-See [the gdalwarp docs](https://gdal.org/programs/gdalwarp.html) for a list of arguments.
-
-# Keywords
-
-$FILENAME_KEYWORD
-$SUFFIX_KEYWORD
-
-Any additional keywords are passed to `ArchGDAL.Dataset`.
-
-## Example
-
-This simply resamples the array with the `:tr` (output file resolution) and `:r`
-flags, giving us a pixelated version:
-
-```jldoctest
-using Rasters, RasterDataSources, Plots
-A = Raster(WorldClim{Climate}, :prec; month=1)
-a = plot(A)
-
-flags = Dict(
-    :tr => [2.0, 2.0],
-    :r => :near,
-)
-b = plot(warp(A, flags))
-
-savefig(a, "docs/build/warp_example_before.png");
-savefig(b, "docs/build/warp_example_after.png"); nothing
-
-# output
-
-```
-
-### Before `warp`:
-
-![before warp](warp_example_before.png)
-
-### After `warp`:
-
-![after warp](warp_example_after.png)
-
-In practise, prefer [`resample`](@ref) for this. But `warp` may be more flexible.
-
-$EXPERIMENTAL
-"""
 function warp(A::AbstractRaster, flags::Dict; filename=nothing, kw...)
     odims = otherdims(A, (X, Y, Band))
     if length(odims) > 0
@@ -71,23 +15,26 @@ function warp(st::AbstractRasterStack, flags::Dict; filename=nothing, suffix=key
 end
 
 function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw...)
+    A1 = _set_gdalwarp_sampling(A)
     filename = RA._maybe_add_suffix(filename, suffix)
     flagvect = reduce([flags...]; init=String[]) do acc, (key, val)
         append!(acc, String[_asflag(key), _stringvect(val)...])
     end
     tempfile = isnothing(filename) ? nothing : tempname() * ".tif"
     warp_kw = isnothing(filename) || filename == "/vsimem/tmp" ? () : (; dest=filename)
-    warped = AG.Dataset(A; filename=tempfile, kw...) do dataset
+    out = AG.Dataset(A1; filename=tempfile, kw...) do dataset
+        rds = Raster(dataset)
         AG.gdalwarp([dataset], flagvect; warp_kw...) do warped
             # Read the raster lazily, dropping Band if there is none in `A`
             raster = Raster(warped; lazy=true, dropband=!hasdim(A, Band()))
-            # Either read the MEM dataset, or get the filename as a FileArray
-            # And permute the dimensions back to what they were in A
-            p_raster = _maybe_permute_from_gdal(raster, dims(A))
-            # Either read the MEM dataset to an Array, or keep the raster lazy
-            return isnothing(filename) ? read(p_raster) : p_raster
+            # Either read the MEM dataset to an Array, or keep a filename base raster lazy
+            return isnothing(filename) ? read(raster) : raster
         end
     end
+    # And permute the dimensions back to what they were in A
+    out1 = _maybe_restore_from_gdal(out, dims(A))
+    out2 = _reset_gdalwarp_sampling(out1, A)
+    return out2
 end
 
 _asflag(x) = string(x)[1] == '-' ? x : string("-", x)
@@ -96,3 +43,30 @@ _stringvect(x::AbstractVector) = Vector(string.(x))
 _stringvect(x::Tuple) = [map(string, x)...]
 _stringvect(x) = [string(x)]
 
+function _set_gdalwarp_sampling(A)
+    x = if sampling(A, X) isa Points
+        DD.maybeshiftlocus(Start(), set(convertlookup(Projected, dims(A, X)), Intervals(Center())))
+    else
+        DD.maybeshiftlocus(Start(), convertlookup(Projected, dims(A, X)))
+    end
+    y = if sampling(A, Y) isa Points
+        DD.maybeshiftlocus(Start(), set(convertlookup(Projected, dims(A, Y)), Intervals(Center())))
+    else
+        DD.maybeshiftlocus(Start(), convertlookup(Projected, dims(A, Y)))
+    end
+    return set(A, X => x, Y=> y)
+end
+
+function _reset_gdalwarp_sampling(A, template)
+    x = if sampling(template, X) isa Points
+        set(DD.maybeshiftlocus(Center(), lookup(A, X)), Points())
+    else
+        DD.maybeshiftlocus(locus(template, X), lookup(A, X))
+    end
+    y = if sampling(template, Y) isa Points
+        set(DD.maybeshiftlocus(Center(), lookup(A, Y)), Points())
+    else
+        DD.maybeshiftlocus(locus(template, Y), lookup(A, Y))
+    end
+    return set(A, X => x, Y => y)
+end
diff --git a/ext/RastersMakieExt/plotrecipes.jl b/ext/RastersMakieExt/plotrecipes.jl
index 299701d4..58cdf755 100644
--- a/ext/RastersMakieExt/plotrecipes.jl
+++ b/ext/RastersMakieExt/plotrecipes.jl
@@ -32,29 +32,6 @@ lift_layer(rs::RasterStack, ind::Symbol) = getproperty(rs, ind)
 lift_layer(s::RasterSeries, inds...) = getindex(s, inds...)
 
 # The all-inclusive plotting function for a 2D raster
-"""
-    Rasters.rplot([position::GridPosition], raster; kw...)
-
-`raster` may be a `Raster` (of 2 or 3 dimensions) or a `RasterStack` whose underlying rasters are 2 dimensional, or 3-dimensional with a singleton (length-1) third dimension.
-
-## Keywords
-
-- `plottype = Makie.Heatmap`: The type of plot. Can be any Makie plot type which accepts a `Raster`; in practice, `Heatmap`, `Contour`, `Contourf` and `Surface` are the best bets.
-- `axistype = Makie.Axis`: The type of axis. This can be an `Axis`, `Axis3`, `LScene`, or even a `GeoAxis` from GeoMakie.jl.
-- `X = XDim`: The X dimension of the raster.
-- `Y = YDim`: The Y dimension of the raster.
-- `Z = YDim`: The Y dimension of the raster.
-- `draw_colorbar = true`: Whether to draw a colorbar for the axis or not.
-- `colorbar_position = Makie.Right()`: Indicates which side of the axis the colorbar should be placed on.  Can be `Makie.Top()`, `Makie.Bottom()`, `Makie.Left()`, or `Makie.Right()`.
-- `colorbar_padding = Makie.automatic`: The amound of padding between the colorbar and its axis.  If `automatic`, then this is set to the width of the colorbar.
-- `title = Makie.automatic`: The titles of each plot. If `automatic`, these are set to the name of the band.
-- `xlabel = Makie.automatic`: The x-label for the axis. If `automatic`, set to the dimension name of the X-dimension of the raster.
-- `ylabel = Makie.automatic`: The y-label for the axis. If `automatic`, set to the dimension name of the Y-dimension of the raster.
-- `colorbarlabel = ""`: Usually nothing, but here if you need it. Sets the label on the colorbar.
-- `colormap = nothing`: The colormap for the heatmap. This can be set to a vector of colormaps (symbols, strings, `cgrad`s) if plotting a 3D raster or RasterStack.
-- `colorrange = Makie.automatic`: The colormap for the heatmap.  This can be set to a vector of `(low, high)` if plotting a 3D raster or RasterStack.
-- `nan_color = :transparent`: The color which `NaN` values should take. Default to transparent.
-"""
 function Rasters.rplot(position::GridPosition, raster::Union{AbstractRaster{T,2,<:Tuple{D1,D2}},Observable{<:AbstractRaster{T,2,<:Tuple{D1,D2}}}};
     plottype=Makie.Heatmap,
     axistype=Makie.Axis,
diff --git a/src/array.jl b/src/array.jl
index 8a56372a..13e31c39 100644
--- a/src/array.jl
+++ b/src/array.jl
@@ -34,6 +34,7 @@ function filename(A::AbstractRaster)
     if length(arrays) == 0
         nothing
     else
+        # TODO its not clear which array supplies the filename if there are multiple in a broadcast
         filename(first(arrays))
     end
 end
diff --git a/src/extensions.jl b/src/extensions.jl
index 0fdf95af..c18601b2 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -13,8 +13,172 @@ end
 
 
 # stubs that need ArchGDAL
+
+"""
+	resample(x; kw...)
+    resample(xs...; to=first(xs), kw...)
+
+`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 into a single array.
+
+If projections can be converted for each axis independently, it may 
+be faster and more accurate to use [`reproject`](@ref).
+
+Run `using ArchGDAL` to make this method available.
+
+# Arguments
+
+- `x`: the object/s to resample.
+
+# Keywords
+
+- `to`: a `Raster`, `RasterStack`, `Tuple` of `Dimension` or `Extents.Extent`.
+    If no `to` object is provided the extent will be calculated from `x`,
+$RES_KEYWORD
+$SIZE_KEYWORD
+$CRS_KEYWORD
+- `method`: A `Symbol` or `String` specifying the method to use for resampling.
+    From the docs for [`gdalwarp`](https://gdal.org/programs/gdalwarp.html#cmdoption-gdalwarp-r):
+    * `:near`: nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
+    * `:bilinear`: bilinear resampling.
+    * `:cubic`: cubic resampling.
+    * `:cubicspline`: cubic spline resampling.
+    * `:lanczos`: Lanczos windowed sinc resampling.
+    * `:average`: average resampling, computes the weighted average of all non-NODATA contributing pixels.
+        rms root mean square / quadratic mean of all non-NODATA contributing pixels (GDAL >= 3.3)
+    * `:mode`: mode resampling, selects the value which appears most often of all the sampled points.
+    * `:max`: maximum resampling, selects the maximum value from all non-NODATA contributing pixels.
+    * `:min`: minimum resampling, selects the minimum value from all non-NODATA contributing pixels.
+    * `:med`: median resampling, selects the median value of all non-NODATA contributing pixels.
+    * `:q1`: first quartile resampling, selects the first quartile value of all non-NODATA contributing pixels.
+    * `:q3`: third quartile resampling, selects the third quartile value of all non-NODATA contributing pixels.
+    * `:sum`: compute the weighted sum of all non-NODATA contributing pixels (since GDAL 3.1)
+
+    Where NODATA values are set to `missingval`.
+$FILENAME_KEYWORD
+$SUFFIX_KEYWORD
+
+Note:
+- 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
+
+Resample a WorldClim layer to match an EarthEnv layer:
+
+```jldoctest
+using Rasters, RasterDataSources, ArchGDAL, Plots
+A = Raster(WorldClim{Climate}, :prec; month=1)
+B = Raster(EarthEnv{HabitatHeterogeneity}, :evenness)
+
+a = plot(A)
+b = plot(resample(A; to=B))
+
+savefig(a, "docs/build/resample_example_before.png");
+savefig(b, "docs/build/resample_example_after.png"); nothing
+# output
+```
+
+### Before `resample`:
+
+![before resample](resample_example_before.png)
+
+### After `resample`:
+
+![after resample](resample_example_after.png)
+
+$EXPERIMENTAL
+"""
 resample(args...; kw...) = throw_extension_error(resample, "ArchGDAL", :RastersArchGDALExt, args)
+
+"""
+    warp(A::AbstractRaster, flags::Dict; kw...)
+
+Gives access to the GDALs `gdalwarp` method given a `Dict` of 
+`flag => value` arguments that can be converted to strings, or vectors
+where multiple space-separated arguments are required.
+
+Arrays with additional dimensions not handled by GDAL (other than `X`, `Y`, `Band`)
+are sliced, warped, and then combined to match the original array dimensions. 
+These slices will *not* be written to disk and loaded lazily at this stage -
+you will need to do that manually if required.
+
+See [the gdalwarp docs](https://gdal.org/programs/gdalwarp.html) for a list of arguments.
+
+Run `using ArchGDAL` to make this method available.
+
+# Keywords
+
+$FILENAME_KEYWORD
+$SUFFIX_KEYWORD
+
+Any additional keywords are passed to `ArchGDAL.Dataset`.
+
+## Example
+
+This simply resamples the array with the `:tr` (output file resolution) and `:r`
+flags, giving us a pixelated version:
+
+```jldoctest
+using Rasters, RasterDataSources, Plots
+A = Raster(WorldClim{Climate}, :prec; month=1)
+a = plot(A)
+
+flags = Dict(
+    :tr => [2.0, 2.0],
+    :r => :near,
+)
+b = plot(warp(A, flags))
+
+savefig(a, "docs/build/warp_example_before.png");
+savefig(b, "docs/build/warp_example_after.png"); nothing
+
+# output
+
+```
+
+### Before `warp`:
+
+![before warp](warp_example_before.png)
+
+### After `warp`:
+
+![after warp](warp_example_after.png)
+
+In practise, prefer [`resample`](@ref) for this. But `warp` may be more flexible.
+
+$EXPERIMENTAL
+"""
 warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALExt, args)
+
+"""
+    cellsize(x)
+
+Gives the approximate size of each cell in square km. 
+This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.
+
+Run `using ArchGDAL` to make this method available.
+
+# Arguments
+
+- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions. 
+
+## Example
+
+```jldoctest
+using Rasters, Rasters.LookupArrays
+dimz = X(Projected(90.0:0.1:120; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))),
+       Y(Projected(0.0:0.1:50; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
+
+cs = cellsize(dimz)
+
+```
+$EXPERIMENTAL
+"""
 cellsize(args...; kw...) = throw_extension_error(cellsize, "ArchGDAL", :RastersArchGDALExt, args)
 
 # Other shared stubs
diff --git a/src/plotrecipes.jl b/src/plotrecipes.jl
index 29198fa7..38384014 100644
--- a/src/plotrecipes.jl
+++ b/src/plotrecipes.jl
@@ -236,6 +236,29 @@ end
 # Makie.jl recipes
 # initial definitions of `rplot`, to get around the extension package availability question
 
+"""
+    Rasters.rplot([position::GridPosition], raster; kw...)
+
+`raster` may be a `Raster` (of 2 or 3 dimensions) or a `RasterStack` whose underlying rasters are 2 dimensional, or 3-dimensional with a singleton (length-1) third dimension.
+
+## Keywords
+
+- `plottype = Makie.Heatmap`: The type of plot. Can be any Makie plot type which accepts a `Raster`; in practice, `Heatmap`, `Contour`, `Contourf` and `Surface` are the best bets.
+- `axistype = Makie.Axis`: The type of axis. This can be an `Axis`, `Axis3`, `LScene`, or even a `GeoAxis` from GeoMakie.jl.
+- `X = XDim`: The X dimension of the raster.
+- `Y = YDim`: The Y dimension of the raster.
+- `Z = YDim`: The Y dimension of the raster.
+- `draw_colorbar = true`: Whether to draw a colorbar for the axis or not.
+- `colorbar_position = Makie.Right()`: Indicates which side of the axis the colorbar should be placed on.  Can be `Makie.Top()`, `Makie.Bottom()`, `Makie.Left()`, or `Makie.Right()`.
+- `colorbar_padding = Makie.automatic`: The amound of padding between the colorbar and its axis.  If `automatic`, then this is set to the width of the colorbar.
+- `title = Makie.automatic`: The titles of each plot. If `automatic`, these are set to the name of the band.
+- `xlabel = Makie.automatic`: The x-label for the axis. If `automatic`, set to the dimension name of the X-dimension of the raster.
+- `ylabel = Makie.automatic`: The y-label for the axis. If `automatic`, set to the dimension name of the Y-dimension of the raster.
+- `colorbarlabel = ""`: Usually nothing, but here if you need it. Sets the label on the colorbar.
+- `colormap = nothing`: The colormap for the heatmap. This can be set to a vector of colormaps (symbols, strings, `cgrad`s) if plotting a 3D raster or RasterStack.
+- `colorrange = Makie.automatic`: The colormap for the heatmap.  This can be set to a vector of `(low, high)` if plotting a 3D raster or RasterStack.
+- `nan_color = :transparent`: The color which `NaN` values should take. Default to transparent.
+"""
 function rplot(args...)
     @error("Please load `Makie.jl` and then call this function. If Makie is loaded, then you can't call `rplot` with no arguments!")
 end
diff --git a/src/utils.jl b/src/utils.jl
index f180331e..5554a303 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -17,12 +17,10 @@ function _cleankey(name::Union{Symbol,AbstractString,Name,NoName}, i=1)
     end
 end
 
-noindex_to_sampled(A) = rebuild(A; dims=noindex_to_sampled(dims(A)))
-function noindex_to_sampled(dims::DimTuple)
-    map(dims) do d
-        lookup(d) isa NoLookup ? set(d, Sampled) : d
-    end
-end
+nolookup_to_sampled(A) = rebuild(A; dims=nolookup_to_sampled(dims(A)))
+nolookup_to_sampled(dims::DimTuple) = map(nolookup_to_sampled, dims)
+nolookup_to_sampled(d::Dimension) =
+    lookup(d) isa NoLookup ? set(d, Sampled(; sampling=Points())) : d
 
 function _maybe_use_type_missingval(A::AbstractRaster{T}, source::Type{<:Source}) where T
     if ismissing(missingval(A))
diff --git a/test/methods.jl b/test/methods.jl
index 59779f37..52fdc6ac 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -514,125 +514,3 @@ end
                    0.0     0.0     missing    
                    missing missing missing], dims=3))
 end
-
-@testset "resample" begin
-    raster_path = maybedownload("https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif")
-
-    output_res = 0.0027
-    output_crs = EPSG(4326)
-    method = "near"
-
-    ## Resample cea.tif manually with ArchGDAL
-    wkt = convert(String, convert(WellKnownText, output_crs))
-    AG_output = ArchGDAL.read(raster_path) do dataset
-        ArchGDAL.gdalwarp([dataset], ["-t_srs", "$(wkt)",
-                                "-tr", "$(output_res)", "$(output_res)",
-                                "-r", "$(method)"]) do warped
-            ArchGDAL.read(ArchGDAL.getband(warped, 1))
-        end
-    end
-
-    # Resample cea.tif using resample
-    cea = Raster(raster_path; missingval=0x00)
-    raster_output = 
-    resample(cea; res=output_res, crs=output_crs, method)
-    disk_output = resample(cea; res=output_res, crs=output_crs, method, filename="resample.tif")
-
-    cea_permuted = permutedims(Raster(raster_path), (Y, X))
-    permuted_output = resample(cea_permuted, output_res; crs=output_crs, method)
-
-    # Compare ArchGDAL, resample and permuted resample 
-    @test AG_output ==
-        raster_output[Band(1)] ==
-        disk_output[Band(1)] ==
-        permutedims(permuted_output, (X, Y))
-    @test abs(step(dims(raster_output, Y))) ≈
-        abs(step(dims(raster_output, X))) ≈ 
-        abs(step(dims(disk_output, X))) ≈ 
-        abs(step(dims(permuted_output, X))) ≈ output_res
-
-    rm("resample.tif")
-
-    @testset "snapped size and dim index match" begin
-        snaptarget = raster_output
-        snapped = resample(cea; to=snaptarget)
-        disk_snapped = resample(cea; to=snaptarget, filename="raster.tif")
-        @test size(snapped) == size(disk_snapped) == size(snaptarget)
-        @test isapprox(index(snaptarget, Y), index(snapped, Y))
-        @test isapprox(index(snaptarget, X), index(snapped, X))
-        @test isapprox(index(snaptarget, Y), index(disk_snapped, Y))
-        @test isapprox(index(snaptarget, X), index(disk_snapped, X))
-    end
-
-    @testset "`method` only does nothing" begin
-        resampled = resample(cea; method)
-        @test crs(cea) == crs(resampled)
-        @test cea == resampled
-        # There is some floating point error here after Rasters -> GDAL -> Rasterss...
-        # Should we correct it by detecting almost identical extent and using the original?
-        @test_broken extent(cea) = extent(resampled)
-    end
-
-    @testset "only `res` kw changes the array size predictably" begin
-        res = step(span(cea, X)) / 2
-        resampled = resample(cea; res)
-        @test crs(cea) == crs(resampled)
-        @test size(dims(resampled, (X, Y))) == size(dims(cea, (X, Y))) .* 2
-        # GDAL fp error see above
-        @test_broken extent(cea) = extent(resampled)
-        resampled = resample(cea; res=(res, 2res))
-        @test size(dims(resampled, (X, Y))) == (size(cea, X) .* 2, size(cea, Y))
-        resampled = resample(cea; res=(X(2res), Y(res)))
-        @test size(dims(resampled, (X, Y))) == (size(cea, X), size(cea, Y) * 2)
-    end
-
-    @testset "only `size` kw sets the size" begin
-        res = step(span(cea, X)) / 2
-        resampled = resample(cea; size=(100, 200))
-        @test crs(cea) == crs(resampled)
-        @test size(dims(resampled, (X, Y))) == size(resampled[:, :, 1]) == (100, 200)
-        resampled = resample(cea; size=(X(99), Y(111)))
-        @test crs(cea) == crs(resampled)
-        @test size(dims(resampled, (X, Y))) == size(resampled[:, :, 1]) == (99, 111)
-        resampled = resample(cea; size=100)
-        @test size(dims(resampled, (X, Y))) == size(resampled[:, :, 1]) == (100, 100)
-    end
-
-    @testset "Extent `to` can resize arbitrarily" begin
-        to = Extent(X=(-10000.0, 2000.0), Y=(4.2e6, 4.3e6))
-        resampled = resample(cea; to)
-        @test all(map((bs...,) -> all(map(≈, bs...)), extent(resampled), to))
-        @test crs(cea) == crs(resampled)
-        # Most of the area is extended, and filled with missingval
-        @test sum(x -> x === missingval(resampled), resampled) > length(resampled) / 2
-    end
-
-    @testset "Geometries work as `to`" begin
-        geom = GeoInterface.Polygon([[(0.0, 4e6), (-1e5, 4.4e6), (-1e5, 4.2e6)]])
-        resampled = resample(cea; to=geom)
-        @test map(extent(resampled[Band(1)]), GeoInterface.extent(geom)) do bs1, bs2 
-            map(≈, bs2, bs2) |> all 
-        end |> all
-        @test crs(cea) == crs(resampled)
-        # Most of the area is extended, and filled with missingval
-        @test sum(x -> x === missingval(resampled), resampled) > length(resampled) / 2
-    end
-
-    @testset "only `crs` kw changes the array size" begin
-        resampled = resample(cea; crs=EPSG(3857), method)
-        @test size(dims(resampled, (X, Y))) !== size(dims(cea, (X, Y)))
-        @test crs(resampled) == EPSG(3857)
-    end
-
-    @testset "no existing crs warns" begin
-        nocrs = setcrs(cea, nothing)
-        @test crs(nocrs) == nothing
-        @test_warn "does not have crs" resample(nocrs; crs=output_crs, method)
-    end
-
-    @testset "resample eltype propagates" begin
-        r = Raster(rand(UInt8, X(1:10), Y(1:10)))
-        r1 = resample(r; to=r)
-        @test eltype(r1) == UInt8
-    end
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index 56e7c4b5..4c801816 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,6 +23,7 @@ end
 @time @safetestset "adapt" begin include("adapt.jl") end
 @time @safetestset "reproject" begin include("reproject.jl") end
 @time @safetestset "warp" begin include("warp.jl") end
+@time @safetestset "resample" begin include("resample.jl") end
 
 
 # Only test SMAP locally for now, also RasterDataSources because CI dowloads keep breaking
diff --git a/test/sources/gdal.jl b/test/sources/gdal.jl
index 635e4d4b..15ac536b 100644
--- a/test/sources/gdal.jl
+++ b/test/sources/gdal.jl
@@ -76,7 +76,7 @@ gdalpath = maybedownload(url)
         named = set(A, Band => string.(Ref("layer_"), dims(A, Band)))
         tempfile = tempname() * ".tif"
         write(tempfile, named)
-        @test parent(dims(Raster(tempfile), Band)) == ["layer_1", "layer_2"]
+        @test lookup(Raster(tempfile), Band) == ["layer_1", "layer_2"]
         @test keys(RasterStack(tempfile; layersfrom=Band)) == (:layer_1, :layer_2)
         rm(tempfile)
     end
@@ -265,6 +265,7 @@ gdalpath = maybedownload(url)
 
         @testset "2d" begin
             filename = tempname() * ".asc"
+            gdalarray
             @time write(filename, gdalarray)
             saved1 = Raster(filename);
             @test all(saved1 .== gdalarray)
@@ -351,13 +352,13 @@ gdalpath = maybedownload(url)
             ga = Raster(rand(100, 200), (X, Y))
             write(filename, ga)
             saved = Raster(filename)
-            @test all(reverse(saved[Band(1)]; dims=Y) .=== ga)
-            @test saved[1, end] == saved[At(1.0), At(1.0)]
-            @test saved[100, 1] == saved[At(100), At(200)]
+            @test all(saved[Band(1)] .=== ga)
+            @test saved[At(1.0), At(1.0)] == saved[1, 1]
+            @test saved[end, end] == saved[At(100), At(200)]
             filename2 = tempname() * ".tif"
             ga2 = Raster(rand(100, 200), (X(Sampled(101:200)), Y(Sampled(1:200))))
             write(filename2, ga2)
-            @test reverse(Raster(filename2)[Band(1)]; dims=Y) == ga2
+            @test Raster(filename2) == ga2
         end
 
         @testset "to netcdf" begin
diff --git a/test/sources/ncdatasets.jl b/test/sources/ncdatasets.jl
index b3c8c479..20be5240 100644
--- a/test/sources/ncdatasets.jl
+++ b/test/sources/ncdatasets.jl
@@ -1,7 +1,7 @@
 using Rasters, DimensionalData, Test, Statistics, Dates, CFTime, Plots
 using Rasters.LookupArrays, Rasters.Dimensions
 import ArchGDAL, NCDatasets
-using Rasters: FileArray, FileStack, NCDsource, crs
+using Rasters: FileArray, FileStack, NCDsource, crs, bounds
 testdir = realpath(joinpath(dirname(pathof(Rasters)), "../test"))
 include(joinpath(testdir, "test_utils.jl"))
 
@@ -43,8 +43,6 @@ end
 
 @testset "Raster" begin
     @time ncarray = Raster(ncsingle)
-    plot(ncarray)
-
     @time lazyarray = Raster(ncsingle; lazy=true);
     @time eagerarray = Raster(ncsingle; lazy=false);
     @test_throws ArgumentError Raster("notafile.nc")
@@ -211,9 +209,7 @@ end
         a = ncarray[X(At(21.0)), Y(Between(50, 52)), Ti(Near(DateTime360Day(2002, 12)))]
         @test bounds(a) == ((50.0, 52.0),)
         x = ncarray[X(Near(150)), Y(Near(30)), Ti(1)]
-        size(ncarray)
         @test x isa Float32
-        lookup(ncarray)
         dimz = X(Between(-0.0, 360)), Y(Between(-90, 90)), 
                Ti(Between(DateTime360Day(2001, 1, 1), DateTime360Day(2003, 01, 02)))
         @test size(ncarray[dimz...]) == (180, 170, 24)
@@ -308,9 +304,9 @@ end
             @test convert(ProjString, crs(gdalarray)) == convert(ProjString, EPSG(4326))
             @test bounds(gdalarray) == bounds(nccleaned)
             # Tiff locus = Start, Netcdf locus = Center
-            @test reverse(index(gdalarray, Y)) .+ 0.5 ≈ index(nccleaned, Y)
+            @test index(gdalarray, Y) .+ 0.5 ≈ index(nccleaned, Y)
             @test index(gdalarray, X) .+ 1.0  ≈ index(nccleaned, X)
-            @test reverse(Raster(gdalarray); dims=Y()) ≈ nccleaned
+            @test gdalarray ≈ nccleaned
         end
         @testset "to grd" begin
             nccleaned = replace_missing(ncarray[Ti(1)], -9999.0)
@@ -318,19 +314,16 @@ end
             grdarray = Raster("testgrd.gri");
             @test crs(grdarray) == convert(ProjString, EPSG(4326))
             @test bounds(grdarray) == bounds(nccleaned)
-            @test reverse(index(grdarray, Y)) ≈ index(nccleaned, Y) .- 0.5
+            @test index(grdarray, Y) ≈ reverse(index(nccleaned, Y)) .- 0.5
             @test index(grdarray, X) ≈ index(nccleaned, X) .- 1.0
-            @test Raster(grdarray) ≈ reverse(nccleaned; dims=Y)
+            @test reverse(grdarray; dims=Y) ≈ nccleaned
             rm("testgrd.gri")
             rm("testgrd.grd")
         end
     end
 
     @testset "no missing value" begin
-        write("nomissing.nc", 
-              boolmask(ncarray)
-              .* 1
-             )
+        write("nomissing.nc", boolmask(ncarray) .* 1)
         nomissing = Raster("nomissing.nc")
         @test missingval(nomissing) == nothing
         rm("nomissing.nc")

From 1293d5dd69bdba17fa1067900196c84aa4c89270 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 25 Nov 2023 11:52:20 +0100
Subject: [PATCH 02/13] test fixes

---
 test/resample.jl | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/test/resample.jl b/test/resample.jl
index a4364168..437a28a2 100644
--- a/test/resample.jl
+++ b/test/resample.jl
@@ -23,8 +23,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
 
     # Resample cea.tif using resample
     cea = Raster(raster_path; missingval=0x00)
-    raster_output = 
-    resample(cea; res=output_res, crs=output_crs, method)
+    raster_output = resample(cea; res=output_res, crs=output_crs, method)
     disk_output = resample(cea; res=output_res, crs=output_crs, method, filename="resample.tif")
 
     cea_permuted = permutedims(Raster(raster_path), (Y, X))
@@ -59,7 +58,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         @test cea == resampled
         # There is some floating point error here after Rasters -> GDAL -> Rasterss...
         # Should we correct it by detecting almost identical extent and using the original?
-        @test_broken extent(cea) = extent(resampled)
+        @test_broken extent(cea) == extent(resampled)
     end
 
     @testset "only `res` kw changes the array size predictably" begin
@@ -68,7 +67,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         @test crs(cea) == crs(resampled)
         @test size(dims(resampled, (X, Y))) == size(dims(cea, (X, Y))) .* 2
         # GDAL fp error see above
-        @test_broken extent(cea) = extent(resampled)
+        @test_broken extent(cea) == extent(resampled)
         resampled = resample(cea; res=(res, 2res))
         @test size(dims(resampled, (X, Y))) == (size(cea, X) .* 2, size(cea, Y))
         resampled = resample(cea; res=(X(2res), Y(res)))

From 69bce9e7ead5eaf1304b1d74d335e1ed6baf4baf Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 25 Nov 2023 12:35:47 +0100
Subject: [PATCH 03/13] fix output

---
 src/extensions.jl | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/extensions.jl b/src/extensions.jl
index c18601b2..94b87f81 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -171,11 +171,24 @@ Run `using ArchGDAL` to make this method available.
 
 ```jldoctest
 using Rasters, Rasters.LookupArrays
-dimz = X(Projected(90.0:0.1:120; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326))),
-       Y(Projected(0.0:0.1:50; sampling=Intervals(Start()), order = ForwardOrdered(), span = Regular(0.1), crs=EPSG(4326)))
+dimz = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326))),
+       Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326)))
 
 cs = cellsize(dimz)
 
+# output
+4×6 Raster{Float64,2} with dimensions: 
+  X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start} crs: EPSG,
+  Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start} crs: EPSG
+extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))
+missingval: missing
+crs: EPSG:4326
+parent:
+        0.0       10.0       20.0        30.0        40.0            50.0
+  90.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
+ 100.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
+ 110.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
+ 120.0  1.2332e6   1.1952e6   1.12048e6   1.01158e6   8.72085e5  706488.0
 ```
 $EXPERIMENTAL
 """

From 87c6d0a48dab5a981b2655d6df09579de45c3ce1 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 25 Nov 2023 12:38:57 +0100
Subject: [PATCH 04/13] fixes

---
 ext/RastersArchGDALExt/gdal_source.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index ee88ad84..8b56572f 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -77,7 +77,7 @@ function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
 )
     T = Missings.nonmissingtype(T)
     missingval = ismissing(missingval) ? RA._writeable_missing(T) : missingval
-    _create_with_driver(filename, dims, T, missingval; kw...) do _
+    _create_with_driver(filename, dims, T, missingval; kw...) do A
         _maybe_warn_south_up(A, verbose, "Creating South-up. Use `reverse(ydim)` first to write conventional North-up")
         nothing
     end

From 7cd298839389aa3ee38139416a908dc8d463993c Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 25 Nov 2023 12:53:45 +0100
Subject: [PATCH 05/13] f

---
 ext/RastersArchGDALExt/gdal_source.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index 8b56572f..3553efab 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -64,7 +64,7 @@ function Base.write(
     RA.check_can_write(filename, force)
     A1 = _maybe_correct_to_write(A)
     _create_with_driver(filename, dims(A1), eltype(A1), missingval(A1); _block_template=A1, kw...) do dataset
-        _maybe_warn_south_up(A, verbose, "Writing South-up. Use `reverse(raster; dims=Y)` first to write conventional North-up")
+        verbose && _maybe_warn_south_up(A, verbose, "Writing South-up. Use `reverse(x; dims=Y)` first to write conventional North-up")
         open(A1; write=true) do O
             AG.RasterDataset(dataset) .= parent(O)
         end
@@ -78,7 +78,7 @@ function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
     T = Missings.nonmissingtype(T)
     missingval = ismissing(missingval) ? RA._writeable_missing(T) : missingval
     _create_with_driver(filename, dims, T, missingval; kw...) do A
-        _maybe_warn_south_up(A, verbose, "Creating South-up. Use `reverse(ydim)` first to write conventional North-up")
+        verbose && _maybe_warn_south_up(dims, verbose, "Creating South-up. Use `reverse(x; dims=Y)` first to write conventional North-up")
         nothing
     end
 
@@ -86,7 +86,7 @@ function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
 end
 
 _maybe_warn_south_up(A, verbose, msg) = 
-    verbose && lookup(A, Y, msg) isa AbstractProjected && order(A, Y, msg) isa ForwardOrdered && @warn msg
+    verbose && lookup(A, Y) isa AbstractProjected && order(A, Y) isa ForwardOrdered && @warn msg
 
 function RA._open(f, ::Type{GDALsource}, filename::AbstractString; write=false, kw...)
     # Check the file actually exists because the GDAL error is unhelpful

From 6c6c2767349d80342f0fee958c513936bb892bcb Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 25 Nov 2023 13:03:21 +0100
Subject: [PATCH 06/13] tweaks

---
 ext/RastersArchGDALExt/gdal_source.jl | 9 +++++----
 test/methods.jl                       | 3 ++-
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index 3553efab..ca7a5929 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -77,16 +77,17 @@ function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
 )
     T = Missings.nonmissingtype(T)
     missingval = ismissing(missingval) ? RA._writeable_missing(T) : missingval
-    _create_with_driver(filename, dims, T, missingval; kw...) do A
-        verbose && _maybe_warn_south_up(dims, verbose, "Creating South-up. Use `reverse(x; dims=Y)` first to write conventional North-up")
+    _create_with_driver(filename, dims, T, missingval; kw...) do _
+        verbose && _maybe_warn_south_up(dims, verbose, "Creating a South-up raster. Use `reverse(x; dims=Y)` first to write conventional North-up")
         nothing
     end
 
     return Raster(filename; source=GDALsource, name, lazy, dropband=!hasdim(dims, Band))
 end
 
-_maybe_warn_south_up(A, verbose, msg) = 
-    verbose && lookup(A, Y) isa AbstractProjected && order(A, Y) isa ForwardOrdered && @warn msg
+function _maybe_warn_south_up(A, verbose, msg) 
+    verbose && lookup(A, Y) isa AbstractSampled && order(A, Y) isa ForwardOrdered && @warn msg
+end
 
 function RA._open(f, ::Type{GDALsource}, filename::AbstractString; write=false, kw...)
     # Check the file actually exists because the GDAL error is unhelpful
diff --git a/test/methods.jl b/test/methods.jl
index 52fdc6ac..eaaeee25 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -7,7 +7,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
 
 A = [missing 7.0f0; 2.0f0 missing]
 B = [1.0 0.4; 2.0 missing]
-ga = Raster(A, (X, Y); missingval=missing) 
+ga = Raster(A, (X(1.0:1:2.0), Y(1.0:1:2.0)); missingval=missing) 
 st = RasterStack((a=A, b=B), (X, Y); missingval=(a=missing,b=missing))
 
 pointvec = [(-20.0, 30.0),
@@ -52,6 +52,7 @@ gaMi = replace_missing(ga)
     @test all(map(values(replace_missing(st, NaN32)), (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
         all(x .=== y)
     end)
+    ga
     dNaN = replace_missing(ga, NaN32; filename="test.tif")
     @test all(isequal.(dNaN, [NaN32 7.0f0; 2.0f0 NaN32]))
     rm("test.tif")

From b258c8d9e07905b4cfecf03cf3e494bea73e31b9 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 25 Nov 2023 13:44:33 +0100
Subject: [PATCH 07/13] f

---
 test/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/methods.jl b/test/methods.jl
index eaaeee25..7db80598 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -69,7 +69,7 @@ end
     @test parent(boolmask(ga)) isa BitMatrix
     @test boolmask(ga99) == [false true; true false]
     @test boolmask(gaNaN) == [false true; true false]
-    @test dims(boolmask(ga)) == (X(NoLookup(Base.OneTo(2))), Y(NoLookup(Base.OneTo(2))))
+    @test dims(boolmask(ga)) === dims(ga)
     x = boolmask(polygon; res=1.0) 
     @test x == trues(X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing)))
     @test parent(x) isa BitMatrix

From ce9f8b9617ef72b3fa796464de5d81019a33b729 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 2 Dec 2023 10:51:00 +0100
Subject: [PATCH 08/13] fix test

---
 test/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/methods.jl b/test/methods.jl
index 7db80598..b070d02d 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -89,7 +89,7 @@ end
     @test all(missingmask(ga) .=== [missing true; true missing])
     @test all(missingmask(ga99) .=== [missing true; true missing])
     @test all(missingmask(gaNaN) .=== [missing true; true missing])
-    @test dims(missingmask(ga)) == (X(NoLookup(Base.OneTo(2))), Y(NoLookup(Base.OneTo(2))))
+    @test dims(missingmask(ga)) == dims(ga)
     @test missingmask(polygon; res=1.0) == fill!(Raster{Union{Missing,Bool}}(undef, X(Projected(-20:1.0:-1.0; crs=nothing)), Y(Projected(10.0:1.0:29.0; crs=nothing))), true)
     x = missingmask([polygon, polygon]; collapse=false, res=1.0)
     @test eltype(x) == Union{Bool,Missing}

From 564496180ca20ed7cb5597ec263bba5bfd07a652 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sat, 2 Dec 2023 18:57:45 +0100
Subject: [PATCH 09/13] doctest

---
 src/extensions.jl | 17 +++++++++--------
 1 file changed, 9 insertions(+), 8 deletions(-)

diff --git a/src/extensions.jl b/src/extensions.jl
index 94b87f81..4dff7411 100644
--- a/src/extensions.jl
+++ b/src/extensions.jl
@@ -25,7 +25,7 @@ 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 into a single array.
 
-If projections can be converted for each axis independently, it may 
+If projections can be converted for each axis independently, it may
 be faster and more accurate to use [`reproject`](@ref).
 
 Run `using ArchGDAL` to make this method available.
@@ -80,6 +80,7 @@ b = plot(resample(A; to=B))
 
 savefig(a, "docs/build/resample_example_before.png");
 savefig(b, "docs/build/resample_example_after.png"); nothing
+
 # output
 ```
 
@@ -98,12 +99,12 @@ resample(args...; kw...) = throw_extension_error(resample, "ArchGDAL", :RastersA
 """
     warp(A::AbstractRaster, flags::Dict; kw...)
 
-Gives access to the GDALs `gdalwarp` method given a `Dict` of 
+Gives access to the GDALs `gdalwarp` method given a `Dict` of
 `flag => value` arguments that can be converted to strings, or vectors
 where multiple space-separated arguments are required.
 
 Arrays with additional dimensions not handled by GDAL (other than `X`, `Y`, `Band`)
-are sliced, warped, and then combined to match the original array dimensions. 
+are sliced, warped, and then combined to match the original array dimensions.
 These slices will *not* be written to disk and loaded lazily at this stage -
 you will need to do that manually if required.
 
@@ -158,26 +159,26 @@ warp(args...; kw...) = throw_extension_error(warp, "ArchGDAL", :RastersArchGDALE
 """
     cellsize(x)
 
-Gives the approximate size of each cell in square km. 
+Gives the approximate size of each cell in square km.
 This function works for any projection, using an algorithm for polygons on a sphere. It approximates the true size to about 0.1%, depending on latitude.
 
 Run `using ArchGDAL` to make this method available.
 
 # Arguments
 
-- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions. 
+- `x`: A `Raster` or a `Tuple` of `X` and `Y` dimensions.
 
 ## Example
 
-```jldoctest
-using Rasters, Rasters.LookupArrays
+```julia
+using Rasters, ArchGDAL, Rasters.LookupArrays
 dimz = X(Projected(90.0:10.0:120; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326))),
        Y(Projected(0.0:10.0:50; sampling=Intervals(Start()), order=ForwardOrdered(), span=Regular(10.0), crs=EPSG(4326)))
 
 cs = cellsize(dimz)
 
 # output
-4×6 Raster{Float64,2} with dimensions: 
+4×6 Raster{Float64,2} with dimensions:
   X Projected{Float64} 90.0:10.0:120.0 ForwardOrdered Regular Intervals{Start} crs: EPSG,
   Y Projected{Float64} 0.0:10.0:50.0 ForwardOrdered Regular Intervals{Start} crs: EPSG
 extent: Extent(X = (90.0, 130.0), Y = (0.0, 60.0))

From c24280884f6c8c10e390fb16ce46c8800bbf1217 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sun, 3 Dec 2023 16:36:48 +0100
Subject: [PATCH 10/13] fp broken tests are not always broken

---
 test/resample.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/resample.jl b/test/resample.jl
index 437a28a2..c8d2e10b 100644
--- a/test/resample.jl
+++ b/test/resample.jl
@@ -58,7 +58,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         @test cea == resampled
         # There is some floating point error here after Rasters -> GDAL -> Rasterss...
         # Should we correct it by detecting almost identical extent and using the original?
-        @test_broken extent(cea) == extent(resampled)
+        # @test_broken extent(cea) == extent(resampled)
     end
 
     @testset "only `res` kw changes the array size predictably" begin
@@ -67,7 +67,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         @test crs(cea) == crs(resampled)
         @test size(dims(resampled, (X, Y))) == size(dims(cea, (X, Y))) .* 2
         # GDAL fp error see above
-        @test_broken extent(cea) == extent(resampled)
+        # @test_broken extent(cea) == extent(resampled)
         resampled = resample(cea; res=(res, 2res))
         @test size(dims(resampled, (X, Y))) == (size(cea, X) .* 2, size(cea, Y))
         resampled = resample(cea; res=(X(2res), Y(res)))

From 76d3f282479fdef60ca9fa898afa4368131b4343 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sun, 3 Dec 2023 21:30:30 +0100
Subject: [PATCH 11/13] only reorder sampled

---
 ext/RastersArchGDALExt/gdal_source.jl | 5 +++--
 test/runtests.jl                      | 1 +
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index ca7a5929..7ca89a86 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -549,8 +549,9 @@ _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = permutedims(A, di
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim}) = permutedims(A, dims)
 
 _maybe_restore_from_gdal(A, dims::Tuple) = reorder(permutedims(A, dims), dims) 
-_maybe_restore_from_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = reorder(A, dims)
-_maybe_restore_from_gdal(A, dims::Tuple{<:XDim,<:YDim}) = reorder(A, dims)
+function _maybe_restore_from_gdal(A, dims::Union{Tuple{<:XDim,<:YDim,<:Band},Tuple{<:XDim,<:YDim}}) 
+    all(map(l -> l isa AbstractSampled, lookup(A, (XDim, YDim)))) ? reorder(A, dims) : A
+end
 
 #= Geotranforms ########################################################################
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 4c801816..0bb713ba 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -32,6 +32,7 @@ if !haskey(ENV, "CI")
     @time @safetestset "rasterdatasources" begin include("sources/rasterdatasources.jl") end
     @time @safetestset "smap" begin include("sources/smap.jl") end
 end
+
 if !Sys.iswindows()
     # GDAL Environment vars need to be set manually for windows, so skip for now
     @time @safetestset "gdal" begin include("sources/gdal.jl") end

From 0b557bdaa125e8018fa28a06bf6db63e8f8c7da6 Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sun, 10 Dec 2023 12:57:56 +0100
Subject: [PATCH 12/13] bugfix _maybe_reorder

---
 ext/RastersArchGDALExt/gdal_source.jl | 12 ++++++++----
 1 file changed, 8 insertions(+), 4 deletions(-)

diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index 7ca89a86..7edfb93a 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -548,11 +548,15 @@ _maybe_permute_to_gdal(A, dims::Tuple) = A
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = permutedims(A, dims)
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim}) = permutedims(A, dims)
 
-_maybe_restore_from_gdal(A, dims::Tuple) = reorder(permutedims(A, dims), dims) 
-function _maybe_restore_from_gdal(A, dims::Union{Tuple{<:XDim,<:YDim,<:Band},Tuple{<:XDim,<:YDim}}) 
-    all(map(l -> l isa AbstractSampled, lookup(A, (XDim, YDim)))) ? reorder(A, dims) : A
-end
+_maybe_restore_from_gdal(A, dims::Tuple) = _maybe_reorder(permutedims(A, dims), dims) 
+_maybe_restore_from_gdal(A, dims::Union{Tuple{<:XDim,<:YDim,<:Band},Tuple{<:XDim,<:YDim}}) = 
+    _maybe_reorder(A, dims) 
 
+function _maybe_reorder(A, dims)
+    reduce(dims; init=A) do a, d
+        (lookup(d) isa AbstractSampled && lookup(A, d) isa AbstractSampled) ? reorder(A, d) : A
+    end
+end
 #= Geotranforms ########################################################################
 
 See https://lists.osgeo.org/pipermail/gdal-dev/2011-July/029449.html

From d0debfed27017c5ccef4a02f01a90aeef624ceca Mon Sep 17 00:00:00 2001
From: rafaqz <rafaelschouten@gmail.com>
Date: Sun, 10 Dec 2023 13:53:14 +0100
Subject: [PATCH 13/13] actually fix it

---
 ext/RastersArchGDALExt/gdal_source.jl | 45 ++++++++++++++-------------
 test/resample.jl                      |  1 +
 2 files changed, 25 insertions(+), 21 deletions(-)

diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl
index 7edfb93a..1efee6ab 100644
--- a/ext/RastersArchGDALExt/gdal_source.jl
+++ b/ext/RastersArchGDALExt/gdal_source.jl
@@ -3,7 +3,7 @@ const AG = ArchGDAL
 const GDAL_LOCUS = Start()
 
 # drivers supporting the gdal Create() method to directly write to disk
-const GDAL_DRIVERS_SUPPORTING_CREATE = ("GTiff", "HDF4", "KEA", "netCDF", "PCIDSK", "Zarr", "MEM"#=...=#) 
+const GDAL_DRIVERS_SUPPORTING_CREATE = ("GTiff", "HDF4", "KEA", "netCDF", "PCIDSK", "Zarr", "MEM"#=...=#)
 
 # order is equal to https://gdal.org/user/virtual_file_systems.html
 const GDAL_VIRTUAL_FILESYSTEMS = "/vsi" .* (
@@ -58,7 +58,7 @@ Write a `Raster` to file using GDAL.
 Returns `filename`.
 """
 function Base.write(
-    filename::AbstractString, ::Type{GDALsource}, A::AbstractRaster{T}; 
+    filename::AbstractString, ::Type{GDALsource}, A::AbstractRaster{T};
     force=false, verbose=true, kw...
 ) where T
     RA.check_can_write(filename, force)
@@ -85,7 +85,7 @@ function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
     return Raster(filename; source=GDALsource, name, lazy, dropband=!hasdim(dims, Band))
 end
 
-function _maybe_warn_south_up(A, verbose, msg) 
+function _maybe_warn_south_up(A, verbose, msg)
     verbose && lookup(A, Y) isa AbstractSampled && order(A, Y) isa ForwardOrdered && @warn msg
 end
 
@@ -105,7 +105,7 @@ function RA._open(f, ::Type{GDALsource}, filename::AbstractString; write=false,
             end
         end
     end
-    if write 
+    if write
         # Pass the OF_UPDATE flag to GDAL
         AG.readraster(RA.cleanreturn ∘ f, filename; flags=AG.OF_UPDATE)
     else
@@ -145,12 +145,12 @@ function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
     # otherwise use Transformed index, with an affine map.
     if _isalligned(gt)
         xstep = gt[GDAL_WE_RES]
-        if xstep > 0 
+        if xstep > 0
             xmin = gt[GDAL_TOPLEFT_X]
             xmax = gt[GDAL_TOPLEFT_X] + xstep * (xsize - 1)
             xorder = ForwardOrdered()
         else
-            xmin = gt[GDAL_TOPLEFT_X] + xstep 
+            xmin = gt[GDAL_TOPLEFT_X] + xstep
             xmax = gt[GDAL_TOPLEFT_X] + xstep * xsize
             xorder = ReverseOrdered()
         end
@@ -158,7 +158,7 @@ function DD.dims(raster::AG.RasterDataset, crs=nothing, mappedcrs=nothing)
         xorder = xstep > 0 ? ForwardOrdered() : ReverseOrdered()
 
         ystep = gt[GDAL_NS_RES] # Usually a negative number
-        if ystep > 0 
+        if ystep > 0
             ymax = gt[GDAL_TOPLEFT_Y]
             ymin = gt[GDAL_TOPLEFT_Y] + ystep * (ysize - 1)
             yorder = ForwardOrdered()
@@ -230,10 +230,10 @@ end
 # Create a Raster from a dataset
 RA.Raster(ds::AG.Dataset; kw...) = Raster(AG.RasterDataset(ds); kw...)
 function RA.Raster(ds::AG.RasterDataset;
-    crs=crs(ds), 
+    crs=crs(ds),
     mappedcrs=nothing,
     dims=dims(ds, crs, mappedcrs),
-    refdims=(), 
+    refdims=(),
     name=Symbol(""),
     metadata=metadata(ds),
     missingval=missingval(ds),
@@ -341,7 +341,7 @@ function _check_driver(filename::AbstractString, driver)
     return driver
 end
 
-# Handle creating a dataset with any driver, 
+# Handle creating a dataset with any driver,
 # applying the function `f` to the created dataset
 function _create_with_driver(f, filename, dims, T, missingval;
     options=Dict{String,String}(), driver="", _block_template=nothing, kw...
@@ -351,7 +351,7 @@ function _create_with_driver(f, filename, dims, T, missingval;
     x, y = map(DD.dims(dims, (XDim, YDim))) do d
         maybeshiftlocus(Start(), RA.nolookup_to_sampled(d))
     end
-    newdims = hasdim(dims, Band()) ? (x, y, DD.dims(dims, Band)) : (x, y) 
+    newdims = hasdim(dims, Band()) ? (x, y, DD.dims(dims, Band)) : (x, y)
     nbands = hasdim(dims, Band) ? length(DD.dims(dims, Band())) : 1
 
     driver = _check_driver(filename, driver)
@@ -465,7 +465,7 @@ end
 
 # Set the properties of an ArchGDAL Dataset to match
 # the dimensions and missingval of a Raster
-_set_dataset_properties!(ds::AG.Dataset, A) = 
+_set_dataset_properties!(ds::AG.Dataset, A) =
     _set_dataset_properties!(ds, dims(A), missingval(A))
 function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval)
     # We cant write mixed Points/Intervals, so default to Intervals if mixed
@@ -484,7 +484,7 @@ function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval)
     y = DD.maybeshiftlocus(GDAL_LOCUS, y)
 
     # Set GDAL AREA_OR_POINT metadata
-    area_or_point = sampling(x) isa Points ? "Point" : "Area" 
+    area_or_point = sampling(x) isa Points ? "Point" : "Area"
     AG.GDAL.gdalsetmetadataitem(dataset, "AREA_OR_POINT", area_or_point, "")
 
     # Set crs if it exists, converting crs to WKT
@@ -532,8 +532,8 @@ end
 _extensiondriver(filename::Nothing) = "MEM"
 function _extensiondriver(filename::AbstractString)
     # TODO move this check to ArchGDAL
-    if filename == "/vsimem/tmp" 
-        "MEM" 
+    if filename == "/vsimem/tmp"
+        "MEM"
     elseif splitext(filename)[2] == ".tif"
         # Force GTiff as the default for .tif because COG cannot do `create` yet
         "GTiff"
@@ -548,13 +548,16 @@ _maybe_permute_to_gdal(A, dims::Tuple) = A
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim,<:Band}) = permutedims(A, dims)
 _maybe_permute_to_gdal(A, dims::Tuple{<:XDim,<:YDim}) = permutedims(A, dims)
 
-_maybe_restore_from_gdal(A, dims::Tuple) = _maybe_reorder(permutedims(A, dims), dims) 
-_maybe_restore_from_gdal(A, dims::Union{Tuple{<:XDim,<:YDim,<:Band},Tuple{<:XDim,<:YDim}}) = 
-    _maybe_reorder(A, dims) 
+_maybe_restore_from_gdal(A, dims::Tuple) = _maybe_reorder(permutedims(A, dims), dims)
+_maybe_restore_from_gdal(A, dims::Union{Tuple{<:XDim,<:YDim,<:Band},Tuple{<:XDim,<:YDim}}) =
+    _maybe_reorder(A, dims)
 
 function _maybe_reorder(A, dims)
-    reduce(dims; init=A) do a, d
-        (lookup(d) isa AbstractSampled && lookup(A, d) isa AbstractSampled) ? reorder(A, d) : A
+    if all(map(l -> l isa AbstractSampled, lookup(dims, (XDim, YDim)))) &&
+        all(map(l -> l isa AbstractSampled, lookup(A, (XDim, YDim))))
+        reorder(A, dims)
+    else
+        A
     end
 end
 #= Geotranforms ########################################################################
@@ -572,7 +575,7 @@ adfGeoTransform[5] /* n-s pixel resolution (negative value) */
 =#
 
 # These function are defined in ext/RastersCoordinateTransformationsExt.jl
-const USING_COORDINATETRANSFORMATIONS_MESSAGE = 
+const USING_COORDINATETRANSFORMATIONS_MESSAGE =
     "Run `using CoordinateTransformations` to load affine transformed rasters"
 
 function RA.dims2geotransform(x::XDim, y::YDim)
diff --git a/test/resample.jl b/test/resample.jl
index c8d2e10b..04576505 100644
--- a/test/resample.jl
+++ b/test/resample.jl
@@ -137,6 +137,7 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
         rev_x_interval_dims = interval_dims[1], reverse(interval_dims[2])
         test_dims = (point_dims, interval_dims, rev_x_point_dims, rev_x_interval_dims, rev_y_point_dims, rev_y_interval_dims)
         ds_fwd = point_dims; f = identity
+        ds_fwd = point_dims; f = reverse
 
         for ds_fwd in test_dims, f in (identity, reverse)
             ds = f(ds_fwd)