From 3183553890a5389f6f705663b7695d2dacfeb0ca Mon Sep 17 00:00:00 2001 From: rafaqz Date: Sun, 22 Oct 2023 14:23:31 +0200 Subject: [PATCH] more gdal cleanup --- ext/RastersArchGDALExt/gdal_source.jl | 72 ++++++++------- ext/RastersArchGDALExt/warp.jl | 30 ++++--- src/utils.jl | 4 +- test/data/grid_mapping_test.nc | Bin 612 -> 0 bytes test/methods.jl | 122 -------------------------- test/runtests.jl | 1 + 6 files changed, 56 insertions(+), 173 deletions(-) delete mode 100644 test/data/grid_mapping_test.nc diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl index f91190774..fde200607 100644 --- a/ext/RastersArchGDALExt/gdal_source.jl +++ b/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) diff --git a/ext/RastersArchGDALExt/warp.jl b/ext/RastersArchGDALExt/warp.jl index b17e09260..caf98a6aa 100644 --- a/ext/RastersArchGDALExt/warp.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index f180331e0..2de678321 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/data/grid_mapping_test.nc b/test/data/grid_mapping_test.nc deleted file mode 100644 index 5cdf46730ae29523e9225ac67d12545a4ef07696..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 612 zcmZ9J%}&EG41m+|&jtu~=E!>tfw&_f!~rQ-Yno-GNpZ4i+~~XTAUpC>NZ83DrbTLP z$6x$=cmGg{d`?LW=&L@n&G*&nU(cElLsFPKYN6|l)8&3d$}xSH8ci)XPU)akncz3$ z&RakkqC(|0LN_y0_Fh|vCN63veIq;cb2W{|JL%DA 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 2da7b3a66..1adc542e8 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