match to-dims with resample (#657)

* use size not step, maybe set dims

* error for irregular dimensions

* better check before rebuilding with `to`

* add a test with some outrageous dims

* remove an unnecessary `all`

* allow and test rasters with more than 2 dims

* remove some accidental code (oops)

* import `@dim` and `YDim`

* finally fix tests

ext/RastersArchGDALExt/resample.jl

@@ -17,6 +17,9 @@ function resample(A::RasterStackOrArray;
     # Method
     flags[:r] = method
 
+    # check if only to has been provided before overwriting arguments
+    onlyto = !isnothing(to) && !isnothing(dims(to)) && !(to isa Extents.Extent) && isnothing(res) && isnothing(size) && isnothing(crs)
+
     # Extent
     if to isa Extents.Extent || isnothing(to) || isnothing(dims(to))
         to = isnothing(to) || to isa Extents.Extent ? to : GeoInterface.extent(to)
@@ -26,12 +29,14 @@ function resample(A::RasterStackOrArray;
             flags[:te] = [xmin, ymin, xmax, ymax]
         end
     else
-        all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` must have both `XDim` and `YDim` dimensions to resize with GDAL"))
+        all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` must have both `XDim` and `YDim` dimensions to resample with GDAL"))
 
         # Set res from `to` if it was not already set
         if isnothing(res) && isnothing(size)
-            xres, yres = map(abs ∘ step, span(to, (XDim, YDim)))
-            flags[:tr] = [yres, xres]
+            todims = dims(to, (XDim, YDim))
+            isregular(todims) || throw(ArgumentError("`to` has irregular dimensions. Provide regular dimensions, or explicitly provide `res` or `size`."))
+            ysize, xsize = length.(todims)
+            flags[:ts] = [ysize, xsize]
         end
         (xmin, xmax), (ymin, ymax) = bounds(to, (XDim, YDim))
         flags[:te] = [xmin, ymin, xmax, ymax]
@@ -89,11 +94,17 @@ function resample(A::RasterStackOrArray;
     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(A))
-    else
-        return setcrs(resampled, crs)
+    if !isnothing(crs)
+        resampled = setcrs(resampled, crs)
     end
+
+    # if only to is provided and it has dims, make sure dims are the exact same 
+    if onlyto
+        newdims = (commondims(to, XDim, YDim)..., otherdims(A, (XDim, YDim))...)
+        resampled = rebuild(resampled; dims =newdims)
+    end
+
+    return resampled
 end
 
-_size_and_res_error() = throw(ArgumentError("Include only `size` or `res` keywords, not both"))
+_size_and_res_error() = throw(ArgumentError("Include only `size` or `res` keywords, not both"))

test/resample.jl

@@ -1,6 +1,8 @@
 using Rasters, ArchGDAL, GeoInterface, Extents
 using Test
 using Rasters.Lookups
+import DimensionalData: @dim, YDim
+@dim Lat YDim "Lat"
 
 include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
 
@@ -157,4 +159,21 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
             @test dims(resampled_res) == dims(resampled_size) == dims(raster)
         end
     end
+
+    @testset "dimensions matcha after resampling with only `to`" begin
+        # some weird dimensions
+        to = (
+            Lat(Projected(1:10, span = Regular(1 + eps()), crs = nothing, order = ForwardOrdered(), sampling = Intervals(Center()))),
+            X(Sampled([1,2,3], span = Regular(1), order = ForwardOrdered(), sampling = Intervals(Start())))
+        )
+
+        resampled = resample(cea; to)
+        @test dims(resampled) == to
+
+        # test with 3d
+        resampled_3D = resample(cat(cea, cea; dims = Z(1:2)); to)
+        @test length(dims(resampled_3D)) == 3
+        @test dims(resampled_3D, (1,2)) == to
+        @test dims(resampled_3D, Z) == Z(1:2)
+    end
 end