Better resample (#425)

* update resample

* add size kw to resample and test

* bugfix resample to geom

Author: rafaqz

src/array.jl

@@ -65,6 +65,7 @@ function crs(obj)
         nothing
     end
 end
+crs(::Nothing) = nothing
 crs(dim::Dimension) = crs(lookup(dim))
 
 """

src/methods/mask.jl

@@ -6,10 +6,11 @@ const TO_KEYWORD = """
 """
 const SIZE_KEYWORD = """
 - `size`: the size of the output array, as a `Tuple{Int,Int}` or single `Int` for a square.
-    Only required when `to` is not used or is an `Extents.Extent`, otherwise `size`.
+    Only required when `to` is not used or is an `Extents.Extent`, and `res` is not used.
 """
 const RES_KEYWORD = """
 - `res`: the resolution of the dimensions, a `Real` or `Tuple{<:Real,<:Real}`.
+    Only required when `to` is not used or is an `Extents.Extent`, and `size` is not used.
 """
 const CRS_KEYWORD = """
 - `crs`: a `crs` which will be attached to the resulting raster when `to` not passed

src/methods/resample.jl

@@ -1,24 +1,23 @@
 """
-	resample(x, resolution::Number; crs, method)
-	resample(x; to, method)
-    resample(xs...; to=first(xs), method)
+	resample(x; to, size, res, method)
+    resample(xs...; to=first(xs), size, res, method)
 
 `resample` uses `ArchGDAL.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`.
 
 # Arguments
 
-- `x`: the object to resample.
-- `resolution`: a `Number` specifying the resolution for the output.
-    If the keyword argument `crs` (described below) is specified, `resolution` must be in units of the `crs`.
+- `x`: the object/s to resample.
 
 # Keywords
 
-- `to`: an `AbstractRaster` whose resolution, crs and bounds will be snapped to.
-    For best results it should roughly cover the same extent, or be a subset of `A`.
-- `crs`: A `GeoFormatTypes.GeoFormat` such as `EPSG(x)` or `WellKnownText(string)` specifying an
-    output crs (`A` will be reprojected to `crs` in addition to being resampled). Defaults to `crs(A)`
+- `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`: A `GeoFormatTypes.GeoFormat` coordinate reference system for the output raster, 
+    such as `EPSG(x)` or `WellKnownText(string)`. Defaults to `crs(A)`.
 - `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).
@@ -38,12 +37,9 @@ or to snap to the bounds, resolution and crs of the object `to`.
 
     Where NODATA values are set to `missingval`.
 
-Notes:
-- `missingval` of `missing` does not work with GDAL. Use `replace_missing(A, newmissingval)` to
-  assign a missing value before using `resample` if the current value is `missing`. This will be
-  automated in future versions.
-- GDAL may cause some unexpected changes in the data, such as returning a reversed dimension or
-  changing the `crs` datatype from `EPSG` to `WellKnownText`.
+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).
 
 # Example
 
@@ -73,49 +69,107 @@ savefig(b, "build/resample_example_after.png"); nothing
 $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...)
     map(x -> resample(x, args...; kw...), ser)
 end
 function resample(xs::Union{Tuple,NamedTuple}; to=first(xs), kw...)
     map(x -> resample(x; to, kw...), xs)
 end
-function resample(A::RasterStackOrArray, resolution::Number; kw...)
-    dimres = map(d -> DD.basetypeof(d)(resolution), dims(A, (XDim, YDim)))
-    resample(A, dimres; kw...)
-end
-function resample(A::RasterStackOrArray, pairs::Pair...; kw...)
-    resample(A, map((D, x) -> D(x), pairs); kw...)
-end
-function resample(A::RasterStackOrArray, resdims::DimTuple;
-    crs::GeoFormat=crs(A), method=:near, kw...
+function resample(x::RasterStackOrArray; 
+    # We need to combine the `size` and `res` keywords with 
+    # the extent in extent2dims, even if we already have dims.
+    to=nothing, res=nothing, crs=nothing, size=nothing, method=:near, kw...
 )
-    wkt = convert(String, convert(WellKnownText, crs))
-    res = map(val, dims(resdims, (YDim, XDim)))
-    flags = Dict(
-        :t_srs => wkt,
-        :tr => res,
-        :r => method,
-    )
-    return warp(A, flags; kw...)
-end
-function resample(A::RasterStackOrArray; to, method=:near, kw...)
-    all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` mush have both XDim and YDim dimensions to resize with GDAL"))
-    if sampling(to, XDim) isa Points
-        to = set(to, dims(to, XDim) => Intervals(Start()))
+    (isnothing(size) || isnothing(res)) || _size_and_res_error()
+
+    # Flags to send to `warp`, then to GDAL
+    flags = Dict{Symbol,Any}()
+
+    # Method
+    flags[:r] = method
+
+    # Extent
+    if to isa Extents.Extent || isnothing(dims(to))
+        to = isnothing(to) || to isa Extents.Extent ? to : GeoInterface.extent(to)
+        if !isnothing(to)
+            # Get the extent of geometries
+            (xmin, xmax), (ymin, ymax) = to[(:X, :Y)]
+            flags[:te] = [xmin, ymin, xmax, ymax]
+        end
+    else
+        all(hasdim(to, (XDim, YDim))) || throw(ArgumentError("`to` mush 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)
+            xres, yres = map(abs ∘ step, span(to, (XDim, YDim)))
+            flags[:tr] = [yres, xres]
+        end
+        (xmin, xmax), (ymin, ymax) = bounds(to, (XDim, YDim))
+        flags[:te] = [xmin, ymin, xmax, ymax]
     end
-    if sampling(to, YDim) isa Points
-        to = set(to, dims(to, YDim) => Intervals(Start()))
+
+    # CRS
+    crs = if isnothing(crs) 
+        if to isa Extents.Extent
+            nothing
+        else
+            # get crs from `to` or `x` if none was passed in
+            isnothing(Rasters.crs(to)) ? Rasters.crs(x) : Rasters.crs(to)
+        end
+    else
+        crs
+    end
+    if !isnothing(crs)
+        wkt = convert(String, convert(WellKnownText, crs))
+        flags[:t_srs] = wkt
+        if isnothing(Rasters.crs(x))
+            @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
 
-    wkt = convert(String, convert(WellKnownText, crs(to)))
-    xres, yres = map(abs ∘ step, span(to, (XDim, YDim)))
-    (xmin, xmax), (ymin, ymax) = bounds(to, (XDim, YDim))
-    flags = Dict(
-        :t_srs => wkt,
-        :tr => [yres, xres],
-        :te => [xmin, ymin, xmax, ymax],
-        :r => method,
-    )
-    return warp(A, flags; kw...)
+    # Resolution
+    if !isnothing(res)
+        xres, yres = if res isa Real
+            res, res
+        elseif res isa Tuple{<:Dimension{<:Real},<:Dimension{<:Real}}
+            map(val, dims(res, (YDim, XDim)))
+        elseif res isa Tuple{<:Real,<:Real}
+            reverse(res)
+        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
+
+    # Size
+    if !isnothing(size)
+        xsize, ysize = if size isa Int
+            size, size
+        elseif size isa Tuple{<:Dimension{Int},<:Dimension{Int}}
+            map(val, dims(size, (YDim, XDim)))
+        elseif size isa Tuple{Int,Int}
+            reverse(size)
+        else
+            throw(ArgumentError("`size` must be a `Int`, or a 2 `Tuple` of `Int` or `Dimension`s wrapping `Int`. Got $size"))
+        end
+        flags[:ts] = [ysize, xsize]
+    end
+
+    # resample with `warp`
+    resampled = warp(x, flags; kw...)
+
+    # Return crs to the original type, from GDAL it will always be WellKnownText
+    if isnothing(crs)
+        return setcrs(resampled, Rasters.crs(x))
+    else
+        return setcrs(resampled, crs)
+    end
 end

src/utils.jl

@@ -137,7 +137,7 @@ function _extent2dims(to::Extents.Extent, size::Nothing, res::Nothing, crs)
     isnothing(res) && throw(ArgumentError("Pass either `size` or `res` keywords or a `Tuple` of `Dimension`s for `to`."))
 end
 function _extent2dims(to::Extents.Extent, size, res, crs)
-    isnothing(res) || throw(ArgumentError("Both `size` and `res` keywords are passed, but only one can be used"))
+    isnothing(res) || _size_and_res_error()
 end
 function _extent2dims(to::Extents.Extent{K}, size::Nothing, res::Real, crs) where K
     tuple_res = ntuple(_ -> res, length(K))
@@ -199,6 +199,8 @@ _progress(args...; kw...) = ProgressMeter.Progress(args...; color=:blue, barlen=
 # Function barrier for splatted vector broadcast
 @noinline _do_broadcast!(f, x, args...) = broadcast!(f, x, args...)
 
+_size_and_res_error() = throw(ArgumentError("Both `size` and `res` keywords are passed, but only one can be used"))
+
 _type_missingval(::Type{T}) where T = typemin(T)
 _type_missingval(::Type{T}) where T<:Unsigned = typemax(T) 
 

test/methods.jl

@@ -702,25 +702,25 @@ end
 
     output_res = 0.0027
     output_crs = EPSG(4326)
-    resample_method = "near"
+    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", "$(resample_method)"]) do warped
+                                "-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, output_res; crs=output_crs, method=resample_method)
-    disk_output = resample(cea, output_res; crs=output_crs, method=resample_method, filename="resample.tif")
+    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, Band))
-    permuted_output = resample(cea_permuted, output_res; crs=output_crs, method=resample_method)
+    permuted_output = resample(cea_permuted, output_res; crs=output_crs, method)
 
     # Compare ArchGDAL, resample and permuted resample 
     @test AG_output ==
@@ -744,4 +744,70 @@ end
         @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
 end