From d29a863fb049fcff4aff7e085f7ef0585e51a891 Mon Sep 17 00:00:00 2001 From: rafaqz Date: Wed, 27 Dec 2023 01:36:20 +0100 Subject: [PATCH 1/4] add Interpolations.jl extension --- Project.toml | 5 ++- ext/DimensionalDataInterpolations.jl | 47 ++++++++++++++++++++++++++++ test/interpolations.jl | 16 ++++++++++ 3 files changed, 67 insertions(+), 1 deletion(-) create mode 100644 ext/DimensionalDataInterpolations.jl create mode 100644 test/interpolations.jl diff --git a/Project.toml b/Project.toml index f9486f601..b5830b094 100644 --- a/Project.toml +++ b/Project.toml @@ -22,10 +22,12 @@ TableTraits = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" [weakdeps] +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" [extensions] DimensionalDataMakie = "Makie" +DimensionalDataInterpolations = "Interpolations" [compat] Adapt = "2, 3.0, 4" @@ -77,6 +79,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -87,4 +90,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" [targets] -test = ["Aqua", "ArrayInterface", "BenchmarkTools", "ColorTypes", "Combinatorics", "CoordinateTransformations", "DataFrames", "Distributions", "Documenter", "ImageFiltering", "ImageTransformations", "CairoMakie", "OffsetArrays", "Plots", "Random", "SafeTestsets", "StatsPlots", "Test", "Unitful"] +test = ["Aqua", "ArrayInterface", "BenchmarkTools", "ColorTypes", "Combinatorics", "CoordinateTransformations", "DataFrames", "Distributions", "Documenter", "ImageFiltering", "ImageTransformations", "Interpolations", "CairoMakie", "OffsetArrays", "Plots", "Random", "SafeTestsets", "StatsPlots", "Test", "Unitful"] diff --git a/ext/DimensionalDataInterpolations.jl b/ext/DimensionalDataInterpolations.jl new file mode 100644 index 000000000..cdf4ce451 --- /dev/null +++ b/ext/DimensionalDataInterpolations.jl @@ -0,0 +1,47 @@ +module DimensionalDataMakie + +using DimensionalData, Interpolations +using DimensionalData.Dimensions, DimensionalData.LookupArrays +using Rasters + +const DD = DimensionalData + +function interp(A::AbstractDimArray; + to=nothing, res=nothing, size=nothing, degree=nothing +) + isempty(otherdims(to, dims(A))) || throw(DimensionMismatch("Cannot interpolate over dimensions not in the source array")) + # TODO handle permutations + shared_dims = commondims(A, to) + dest_dims = dims(Rasters._extent2dims(commondims(to, shared_dims); size, res, crs=crs(A)), shared_dims) + other_dims = otherdims(A, to) + degrees = if isnothing(degree) + map(_ -> Linear(), dims(A)) + elseif !(degree isa Tuple) + map(_ -> degree, dims(A)) + else + degree + end + + data = if isregular(A, (XDim, YDim)) # Do linear interpolation + interpmode = map(dims(A), degrees) do d, deg + BSpline(deg) + end + if isempty(other_dims) + itp = interpolate(A, BSpline(Cubic(Line(OnGrid())))) + sitp = scale(itp, map(parent, lookup(shared_dims))...) + sitp(map(parent, lookup(dest_dims))...) + else + for D in DimIndices(other_dims) + itp = interpolate(view(A, D...), BSpline(Cubic(Line(OnGrid())))) + sitp = scale(itp, map(parent, lookup(shared_dims))...) + sitp(map(parent, lookup(dest_dims))...) + end + end + else + error("interpolate is only implemented for regular grids") + end + + return Raster(data, dest_dims) +end + +end diff --git a/test/interpolations.jl b/test/interpolations.jl new file mode 100644 index 000000000..539d4568e --- /dev/null +++ b/test/interpolations.jl @@ -0,0 +1,16 @@ +using DimensionalData +using Interpolations + +f((x1, x2)) = log(x1+x2) +A = f.(DimPoints((X(1:.1:10), Y(1:.5:20)))) +to = rand(X(2:.3:7), Y(2:.3:17)) +out = DimensionalData.interp(A; to) + +A_x1 = 1:.1:10 +A_x2 = 1:.5:20 +f(x1, x2) = log(x1+x2) +A = [f(x1,x2) for x1 in A_x1, x2 in A_x2] +itp = interpolate(A, BSpline(Cubic(Line(OnGrid())))) +sitp = scale(itp, A_x1, A_x2) +sitp(5., 10.) # exactly log(5 + 10) +sitp([5.6, 5.2], [7.1, 7.1]) # approximately log(5.6 + 7.1) From 138b7777c43d5606720c33dd6e461e67f6ef115d Mon Sep 17 00:00:00 2001 From: rafaqz Date: Wed, 27 Dec 2023 01:45:47 +0100 Subject: [PATCH 2/4] da --- ext/DimensionalDataInterpolations.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/DimensionalDataInterpolations.jl b/ext/DimensionalDataInterpolations.jl index cdf4ce451..9285b5d03 100644 --- a/ext/DimensionalDataInterpolations.jl +++ b/ext/DimensionalDataInterpolations.jl @@ -12,7 +12,8 @@ function interp(A::AbstractDimArray; isempty(otherdims(to, dims(A))) || throw(DimensionMismatch("Cannot interpolate over dimensions not in the source array")) # TODO handle permutations shared_dims = commondims(A, to) - dest_dims = dims(Rasters._extent2dims(commondims(to, shared_dims); size, res, crs=crs(A)), shared_dims) + # TODO move `_extent2dims` to DimensionalData` + dest_dims = dims(Rasters._extent2dims(commondims(to, shared_dims); size, res), shared_dims) other_dims = otherdims(A, to) degrees = if isnothing(degree) map(_ -> Linear(), dims(A)) @@ -41,7 +42,7 @@ function interp(A::AbstractDimArray; error("interpolate is only implemented for regular grids") end - return Raster(data, dest_dims) + return DimArray(data, dest_dims) end end From 06b5cd7b433ca8c74ef839f8858a384faf014c4e Mon Sep 17 00:00:00 2001 From: rafaqz Date: Wed, 27 Dec 2023 01:46:45 +0100 Subject: [PATCH 3/4] clean up test --- test/interpolations.jl | 9 --------- 1 file changed, 9 deletions(-) diff --git a/test/interpolations.jl b/test/interpolations.jl index 539d4568e..f6a4cef49 100644 --- a/test/interpolations.jl +++ b/test/interpolations.jl @@ -5,12 +5,3 @@ f((x1, x2)) = log(x1+x2) A = f.(DimPoints((X(1:.1:10), Y(1:.5:20)))) to = rand(X(2:.3:7), Y(2:.3:17)) out = DimensionalData.interp(A; to) - -A_x1 = 1:.1:10 -A_x2 = 1:.5:20 -f(x1, x2) = log(x1+x2) -A = [f(x1,x2) for x1 in A_x1, x2 in A_x2] -itp = interpolate(A, BSpline(Cubic(Line(OnGrid())))) -sitp = scale(itp, A_x1, A_x2) -sitp(5., 10.) # exactly log(5 + 10) -sitp([5.6, 5.2], [7.1, 7.1]) # approximately log(5.6 + 7.1) From 1f8767bec3e2456aeaf0bbf9d12ac20b3d986e1a Mon Sep 17 00:00:00 2001 From: Daniel Boland Date: Mon, 24 Jun 2024 12:17:50 -0400 Subject: [PATCH 4/4] linear interpolation (#609) * linear interpolation * Apply suggestions from code review Co-authored-by: Rafael Schouten * Simpler implementation * compat * fixed typo missed = * Update ext/DimensionalDataInterpolations.jl * Update ext/DimensionalDataInterpolations.jl --------- Co-authored-by: Rafael Schouten --- Project.toml | 1 + ext/DimensionalDataInterpolations.jl | 49 ++++------------------------ test/interpolations.jl | 20 ++++++++++-- 3 files changed, 26 insertions(+), 44 deletions(-) diff --git a/Project.toml b/Project.toml index b5830b094..940e6322c 100644 --- a/Project.toml +++ b/Project.toml @@ -46,6 +46,7 @@ Documenter = "1" Extents = "0.1" ImageFiltering = "0.7" ImageTransformations = "0.10" +Interpolations = "0.15" IntervalSets = "0.5, 0.6, 0.7" InvertedIndices = "1" IteratorInterfaceExtensions = "1" diff --git a/ext/DimensionalDataInterpolations.jl b/ext/DimensionalDataInterpolations.jl index 9285b5d03..29f5997ae 100644 --- a/ext/DimensionalDataInterpolations.jl +++ b/ext/DimensionalDataInterpolations.jl @@ -1,48 +1,13 @@ -module DimensionalDataMakie +module DimensionalDataInterpolations using DimensionalData, Interpolations -using DimensionalData.Dimensions, DimensionalData.LookupArrays -using Rasters -const DD = DimensionalData - -function interp(A::AbstractDimArray; - to=nothing, res=nothing, size=nothing, degree=nothing -) - isempty(otherdims(to, dims(A))) || throw(DimensionMismatch("Cannot interpolate over dimensions not in the source array")) - # TODO handle permutations - shared_dims = commondims(A, to) - # TODO move `_extent2dims` to DimensionalData` - dest_dims = dims(Rasters._extent2dims(commondims(to, shared_dims); size, res), shared_dims) - other_dims = otherdims(A, to) - degrees = if isnothing(degree) - map(_ -> Linear(), dims(A)) - elseif !(degree isa Tuple) - map(_ -> degree, dims(A)) - else - degree - end - - data = if isregular(A, (XDim, YDim)) # Do linear interpolation - interpmode = map(dims(A), degrees) do d, deg - BSpline(deg) - end - if isempty(other_dims) - itp = interpolate(A, BSpline(Cubic(Line(OnGrid())))) - sitp = scale(itp, map(parent, lookup(shared_dims))...) - sitp(map(parent, lookup(dest_dims))...) - else - for D in DimIndices(other_dims) - itp = interpolate(view(A, D...), BSpline(Cubic(Line(OnGrid())))) - sitp = scale(itp, map(parent, lookup(shared_dims))...) - sitp(map(parent, lookup(dest_dims))...) - end - end - else - error("interpolate is only implemented for regular grids") - end - - return DimArray(data, dest_dims) +function Interpolations.linear_interpolation(A::AbstractDimArray; kw...) + linear_interpolation(DimensionalData.index(dims(A)), A; kw...) end +function Interpolations.cubic_spline_interpolation(A::AbstractDimArray; kw...) + cubic_spline_interpolation(DimensionalData.index(dims(A)), A; kw...) end + +end \ No newline at end of file diff --git a/test/interpolations.jl b/test/interpolations.jl index f6a4cef49..d9c523da7 100644 --- a/test/interpolations.jl +++ b/test/interpolations.jl @@ -2,6 +2,22 @@ using DimensionalData using Interpolations f((x1, x2)) = log(x1+x2) -A = f.(DimPoints((X(1:.1:10), Y(1:.5:20)))) +a = f.(DimPoints((X(1:1:10), Y(1:1.5:20)))) +b = f.(DimPoints((X(1:1:10), Y([1.0,3,7.5,10,15])))) to = rand(X(2:.3:7), Y(2:.3:17)) -out = DimensionalData.interp(A; to) +# out = DimensionalData.interp(A; to) + +itp_a = linear_interpolation(a) +itp_b = linear_interpolation(b) +itp_ex = linear_interpolation(a, extrapolation_bc=(Flat(), Linear())) + +itp_a(7.5,7.5) +itp_b(7.5,7.5) +itp_ex(10,13) +itp_ex(11,13) +itp_ex(10,14) + +itp_ca = cubic_spline_interpolation(a) +# itp_cb = cubic_spline_interpolation(b) #fails, as expected + +itp_ca(10,10) \ No newline at end of file