diff --git a/test/resample.jl b/test/resample.jl new file mode 100644 index 000000000..a43641686 --- /dev/null +++ b/test/resample.jl @@ -0,0 +1,151 @@ +using Rasters, ArchGDAL, GeoInterface, Extents +using Test +using Rasters.LookupArrays + +include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl")) + +@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 + + @testset "resample to the same size or resolution leaves raster unchanged" begin + res = 2 + ys = LinRange((90 - res / 2):-res:(-90 + res / 2)) + xs = LinRange((-180 + res / 2):res:(180 - res / 2)) + + point_dims = Y(ys), X(xs) + interval_dims = Y(ys; sampling=Intervals(Center())), X(xs; sampling=Intervals(Center())) + rev_y_point_dims = Y(reverse(ys)), X(xs) + rev_y_interval_dims = reverse(interval_dims[1]), interval_dims[2] + rev_x_point_dims = Y(ys), X(reverse(xs)) + 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 + + for ds_fwd in test_dims, f in (identity, reverse) + ds = f(ds_fwd) + raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN) + resampled_res = resample(raster; res) + resampled_size = resample(raster; size=size(raster), method=:near) + @test resampled_res == resampled_size == raster + @test dims(resampled_res) == dims(resampled_size) == dims(raster) + end + end +end