From 2053fadf76710a9cd0177bf54e0d1f2a6434c357 Mon Sep 17 00:00:00 2001 From: Rafael Schouten Date: Mon, 11 Dec 2023 23:50:55 +0100 Subject: [PATCH] Extend fixes (#568) * bugfix extend * bugfix * delete display * dont change series --- ext/RastersArchGDALExt/gdal_source.jl | 3 +- ext/RastersNCDatasetsExt/ncdatasets_source.jl | 4 +- src/methods/crop_extend.jl | 45 ++++++++++++------- test/methods.jl | 8 ++-- 4 files changed, 39 insertions(+), 21 deletions(-) diff --git a/ext/RastersArchGDALExt/gdal_source.jl b/ext/RastersArchGDALExt/gdal_source.jl index 3c006b8ce..6b740f141 100644 --- a/ext/RastersArchGDALExt/gdal_source.jl +++ b/ext/RastersArchGDALExt/gdal_source.jl @@ -349,6 +349,7 @@ function _create_with_driver(f, filename, dims, T, missingval; x, y = map(DD.dims(dims, (XDim, YDim))) do d maybeshiftlocus(Start(), RA.nolookup_to_sampled(d)) end + newdims = hasdim(dims, Band()) ? (x, y, DD.dims(dims, Band)) : (x, y) nbands = hasdim(dims, Band) ? length(DD.dims(dims, Band())) : 1 @@ -468,7 +469,7 @@ _set_dataset_properties!(ds::AG.Dataset, A) = function _set_dataset_properties!(dataset::AG.Dataset, dims::Tuple, missingval) # We cant write mixed Points/Intervals, so default to Intervals if mixed xy = DD.dims(dims, (X, Y)) - if any(x -> x isa Intervals, sampling.(xy)) && any(x -> x isa Points, sampling.(xy)) + if any(x -> x isa Intervals, map(sampling, xy)) && any(x -> x isa Points, map(sampling, xy)) dims = set(dims, X => Intervals, Y => Intervals) end # Convert the dimensions to `Projected` if they are `Converted` diff --git a/ext/RastersNCDatasetsExt/ncdatasets_source.jl b/ext/RastersNCDatasetsExt/ncdatasets_source.jl index 8c4acfe0f..8b1029b12 100644 --- a/ext/RastersNCDatasetsExt/ncdatasets_source.jl +++ b/ext/RastersNCDatasetsExt/ncdatasets_source.jl @@ -493,8 +493,8 @@ function _def_dim_var!(ds::NCD.Dataset, dim::Dimension) attrib = _attribdict(metadata(dim)) _ncd_set_axis_attrib!(attrib, dim) # Bounds variables - if span(dim) isa Explicit - bounds = val(span(dim)) + if sampling(dim) isa Intervals + bounds = Dimensions.dim2boundsmatrix(dim) boundskey = get(metadata(dim), :bounds, string(dimkey, "_bnds")) push!(attrib, "bounds" => boundskey) NCD.defVar(ds, boundskey, bounds, ("bnds", dimkey)) diff --git a/src/methods/crop_extend.jl b/src/methods/crop_extend.jl index aaf46f0d6..0ea5a777c 100644 --- a/src/methods/crop_extend.jl +++ b/src/methods/crop_extend.jl @@ -93,6 +93,7 @@ function _crop_to(x, to; kw...) end end _crop_to(A, to::RasterStackOrArray; dims=DD.dims(to), kw...) = _crop_to(A, DD.dims(to, dims); kw...) +_crop_to(x, to::Dimension; kw...) = _crop_to(x, (to,); kw...) function _crop_to(x, to::DimTuple; kw...) # We can only crop to sampled dims (e.g. not categorical dims like Band) sampled = reduce(to; init=()) do acc, d @@ -170,38 +171,52 @@ function _extend_to(x::RasterStackOrArray, to; kw...) isnothing(ext) && throw(ArgumentError("No dims or extent available on `to` object of type $(typeof(to))")) return _extend_to(x, ext; kw...) end +_extend_to(x::RasterStackOrArray, to::Dimension; kw...) = _extend_to(x, (to,); kw...) function _extend_to(A::AbstractRaster, to::DimTuple; - filename=nothing, suffix=nothing, touches=false + filename=nothing, suffix=nothing, touches=false, missingval=missingval(A) ) + others = otherdims(to, A) + # Allow not specifying all dimensions + to = (set(dims(A), map(=>, dims(A, to), to)...)..., others...) # Calculate the range of the old array in the extended array - ranges = _without_mapped_crs(A) do A + rangedims = _without_mapped_crs(A) do A _without_mapped_crs(to) do to - map(dims(A), to) do d, t - range = if touches - DD.selectindices(t, LA.Touches(bounds(d))) - else - DD.selectindices(t, LA.ClosedInterval(bounds(d)...)) - end - rebuild(d, range) + map(dims(A, to), to) do d, t + # Values must match exactly, so use `At` + DD.selectindices(t, At(first(d))):DD.selectindices(t, At(last(d))) end end end + others1 = otherdims(to, A) + final_to = (set(dims(A), map(=>, dims(A, to), to)...)..., others1...) # Create a new extended array - newA = create(filename, eltype(A), to; - suffix, parent=parent(A), missingval=missingval(A), + newA = create(filename, eltype(A), final_to; + suffix, parent=parent(A), missingval, name=name(A), metadata=metadata(A) ) + # Input checks + map(dims(A, to), dims(newA, to)) do d1, d2 + if lookup(d1) isa Union{AbstractSampled,NoLookup} + b1, b2 = bounds(d1), bounds(d2) + b1[1] >= b2[1] || throw(ArgumentError("Lower bound of $(basetypeof(d1)) lookup of `$(b2[1])` are not larger than the original `$(b1[1])`")) + b1[2] <= b2[2] || throw(ArgumentError("Upper bound of $(basetypeof(d2)) lookup of `$(b2[2])` is not larger than the original `$(b1[2])`")) + elseif lookup(d1) isa Categorical + map(lookup(d1)) do x + x in d2 || throw(ArgumentError("category $x not in new dimension")) + end + end + end # The missingval may have changed for disk-based arrays - if !isequal(missingval(A), missingval(newA)) - A = replace_missing(A, missingval(newA)) + if !isequal(missingval, Rasters.missingval(newA)) + A = replace_missing(A, Rasters.missingval(newA)) end open(newA; write=true) do O # Fill it with missing/nodata values - O .= missingval(O) + O .= Rasters.missingval(O) # Copy the original data to the new array # Somehow this is slow from disk? - broadcast_dims!(identity, view(O, ranges...), A) + broadcast_dims!(identity, view(O, rangedims...), A) end return newA end diff --git a/test/methods.jl b/test/methods.jl index 9de5443e3..3b6aa6222 100644 --- a/test/methods.jl +++ b/test/methods.jl @@ -432,15 +432,17 @@ end @test all(trimmed_r .=== [missing 1.0; 0.5 2.0]) end cropped = crop(ga; to=trimmed) + cropped1 = crop(crop(ga; to=dims(trimmed, X)); to=dims(trimmed, Y)) _, cropped2 = crop(trimmed, ga) cropped_r = crop(ga_r; to=trimmed_r) - @test all(cropped .=== trimmed) + @test all(cropped .=== cropped1 .=== trimmed) @test all(cropped_r .=== trimmed_r) extended = extend(cropped, ga)[1] extended_r = extend(cropped_r; to=ga_r) - @test all(extended .=== ga) - @test all(extended_r .=== ga_r) + extended1 = extend(extend(cropped; to=dims(ga, X)); to=dims(ga, Y)) extended_d = extend(cropped; to=ga, filename="extended.tif") + @test all(extended .=== extended1 .=== replace_missing(extended_d) .=== ga) + @test all(extended_r .=== ga_r) @test all(map(==, lookup(extended_d), lookup(extended))) @testset "to polygons" begin