more gdal cleanup

ext/RastersArchGDALExt/gdal_source.jl

@@ -85,16 +85,11 @@ function Base.write(
 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,
+    driver="", lazy=true, 
+    options=Dict{String,String}(), _block_template=nothing
 )
-    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
+    x, y = nolookup_to_sampled(DD.dims(dims, (XDim, YDim)))
     T = Missings.nonmissingtype(T)
 
     if ismissing(missingval)
@@ -110,27 +105,16 @@ function RA.create(filename, ::Type{GDALsource}, T::Type, dims::DD.DimTuple;
         newdims = (x, y)
     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
+    create_kw = (width=length(x), height=length(y), nbands=nbands, dtype=T)
+
+    _create_with_driver(filename, driver, create_kw; _block_template, options) do dataset
+        _set_dataset_properties!(dataset, A)
     end
+
     if hasdim(dims, Band)
-        return Raster(filename; source=GDALsource, lazy)
+        return Raster(filename; source=GDALsource, name, lazy)
     else
-        return view(Raster(filename; source=GDALsource, lazy), Band(1))
+        return view(Raster(filename; source=GDALsource, name, lazy), Band(1))
     end
 end
 
@@ -167,11 +151,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)
@@ -283,6 +269,7 @@ function RA.Raster(ds::AG.RasterDataset;
     lazy=false,
     dropband=false
 )
+    @show dims
     args = dims, refdims, name, metadata, missingval
     filelist = AG.filelist(ds)
     raster = if lazy && length(filelist) > 0
@@ -330,9 +317,8 @@ function AG.Dataset(f::Function, A::AbstractRaster; kw...)
     end
 end
 function AG.RasterDataset(f::Function, A::AbstractRaster; 
-    filename=nothing, driver="",
+    filename=nothing, driver="", options, _block_template=A
 )
-    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]
@@ -343,6 +329,7 @@ function AG.RasterDataset(f::Function, A::AbstractRaster;
 
     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)),
@@ -352,7 +339,8 @@ function AG.RasterDataset(f::Function, A::AbstractRaster;
     if filename == nothing 
         filename = ""
     end
-    _create_with_driver(filename, driver, kw; _block_template=A_p) do dataset
+
+    return _create_with_driver(filename, driver, kw; _block_template) do dataset
         _set_dataset_properties!(dataset, dims(A_p), missingval(A_p))
         rds = AG.RasterDataset(dataset)
         open(A_p) do a
@@ -389,8 +377,7 @@ _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)))
+        a -> RA.nolookup_to_sampled(a)
 end
 
 # Write a Raster to disk using GDAL
@@ -545,7 +532,9 @@ function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval)
     end
 
     # Set the geotransform from the updated lookups
-    AG.setgeotransform!(dataset, RA.dims2geotransform(x, y))
+    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
@@ -621,13 +610,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)

ext/RastersArchGDALExt/warp.jl

@@ -15,21 +15,27 @@ function warp(st::AbstractRasterStack, flags::Dict; filename=nothing, suffix=key
 end
 
 function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw...)
-    @show flags
-    @show sampling.(lookup(A, (X, Y))) bounds(A, (X, Y))
+
+            @show "start"
     A1 = _set_gdalwarp_sampling(A)
-    @show sampling.(lookup(A1, (X, Y))) bounds(A1, (X, Y)) 
+    # @show sampling.(lookup(A1, (X, Y))) bounds(A1, (X, Y)) 
     filename = RA._maybe_add_suffix(filename, suffix)
     flagvect = reduce([flags...]; init=[]) 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)
+            @show "here"
     out = AG.Dataset(A1; filename=tempfile, kw...) do dataset
+        A2 = Raster(dataset)
+        # @show sampling.(lookup(A2, (X, Y))) bounds(A2, (X, Y)) 
         AG.gdalwarp([dataset], flagvect; warp_kw...) do warped
+            # @show dims(AG.RasterDataset(warped))
+            @show "before"
             # Read the raster lazily, dropping Band if there is none in `A`
             raster = Raster(warped; lazy=true, dropband=!hasdim(A, Band()))
-            @show sampling.(lookup(raster, (X, Y))) bounds(raster, (X, Y)) bounds(raster, Y)[2]
+            @show "after"
+            # @show sampling.(lookup(raster, (X, Y))) bounds(raster, (X, Y)) bounds(raster, Y)[2]
 
             # Either read the MEM dataset, or get the filename as a FileArray
             # And permute the dimensions back to what they were in A
@@ -39,9 +45,9 @@ function _warp(A::AbstractRaster, flags::Dict; filename=nothing, suffix="", kw..
             return isnothing(filename) ? read(p_raster) : p_raster
         end
     end
-    @show sampling.(lookup(out, (X, Y))) bounds(out, (X, Y))
+    # @show sampling.(lookup(out, (X, Y))) bounds(out, (X, Y))
     x = _reset_gdalwarp_sampling(out, A)
-    @show sampling.(lookup(x, (X, Y))) bounds(x, (X, Y))
+    # @show sampling.(lookup(x, (X, Y))) bounds(x, (X, Y))
     println()
     return x
 end
@@ -54,12 +60,12 @@ _stringvect(x) = [string(x)]
 
 function _set_gdalwarp_sampling(A)
     x = if sampling(A, X) isa Points
-        convertlookup(Projected, dims(A, X))
+        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
-        convertlookup(Projected, dims(A, Y))
+        DD.maybeshiftlocus(Start(), set(convertlookup(Projected, dims(A, Y)), Intervals(Center())))
     else
         DD.maybeshiftlocus(Start(), convertlookup(Projected, dims(A, Y)))
     end
@@ -68,14 +74,14 @@ end
 
 function _reset_gdalwarp_sampling(A, template)
     x = if sampling(template, X) isa Points
-        lookup(A, X)
+        set(DD.maybeshiftlocus(Center(), lookup(A, X)), Points())
     else
-        DD.maybeshiftlocus(locus(template, X), lookup(A, Y))
+        DD.maybeshiftlocus(locus(template, X), lookup(A, X))
     end
     y = if sampling(template, Y) isa Points
-        lookup(A, Y)
+        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)
+    return set(A, X => x, Y => y)
 end

src/utils.jl

@@ -17,8 +17,8 @@ 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)
+noloopup_to_sampled(A) = rebuild(A; dims=nolookup_to_sampled(dims(A)))
+function nolookup_to_sampled(dims::DimTuple)
     map(dims) do d
         lookup(d) isa NoLookup ? set(d, Sampled) : d
     end

test/methods.jl

@@ -392,125 +392,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

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