# Interpolate along single dimension of multidimensional array

**URL:** <https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290>\
**Category:** General Usage\
**Tags:** interpolations\
**Created:** [September 14, 2022, 10:06pm UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290 "2022-09-14T22:06:16Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![snowgum](https://avatars.discourse-cdn.com/v4/letter/s/ecc23a/32.png) [@snowgum](https://discourse.julialang.org/u/snowgum)\
**Post date:** [September 14, 2022, 10:06pm UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/1 "2022-09-14T22:06:16Z")

</div>

I’m trying to convert a Python code to Julia that makes use of `scipy.interpolate.interp1d`, and I’m having difficulty finding the best way to do this in Julia. The `interp1d` function does 1D interpolation along a given dimension of multidimensional array.

I’ve been using `Interpolations.jl` but my approach does not seem the most elegant (see below). There is a [very similar question](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308) asked three years ago, but the accepted answer does not provide a straightforward answer.

Here’s a simple example of what I’m looking for:

```julia
# Test data
x = [0., 2., 5., 10., 11.]
y = repeat(x .^ 2, outer=(1, 3, 3))

# We want to interpolate y (a 3D array) along the first dimension,
# so that the output is another 3D array
itp = interp1d(x, y; dims=1, k=1)

# Expected results
itp(2.0)
3×3 Matrix{Float64}:
 4.0 4.0 4.0
 4.0 4.0 4.0
 4.0 4.0 4.0

itp(3)
3×3 Matrix{Float64}:
 11.0 11.0 11.0
 11.0 11.0 11.0
 11.0 11.0 11.0

```

My data are all in irregular grids.

This seems to work, but I wonder if it is the best approach:

```julia
itp = interpolate((x, 1:3, 1:3), y, Gridded(Linear()))

itp(3, 1:3, 1:3)
3×3 Matrix{Float64}:
 11.0 11.0 11.0
 11.0 11.0 11.0
 11.0 11.0 11.0

```

And it still seems that with `Gridded` the highest order interpolation supported is linear.  
Anyone else has a better approach for this type of interpolation?

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [September 15, 2022, 7:33am UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/2 "2022-09-15T07:33:26Z")

</div>

Just in case, there is an [Interp1d.jl](https://github.com/AtsushiSakai/Interp1d.jl) package from one of the Python contributors, with an accompanying [notebook](https://github.com/AtsushiSakai/Interp1d.jl/blob/master/examples/interpolate.ipynb) with examples.

---

<div class="post-metadata">

**Author:** ![snowgum](https://avatars.discourse-cdn.com/v4/letter/s/ecc23a/32.png) [@snowgum](https://discourse.julialang.org/u/snowgum)\
**Post date:** [September 15, 2022, 8:29am UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/3 "2022-09-15T08:29:01Z")

</div>

Thanks for the suggestion. I had looked at `Interp1d.jl` before, but it seems very alpha. Not registered package, no documentation, does not work in 3D or higher (only 1D or 2D) and does not even support linear interpolation (only nearest neighbour or previous).

---

<div class="post-metadata">

**Author:** ![snowgum](https://avatars.discourse-cdn.com/v4/letter/s/ecc23a/32.png) [@snowgum](https://discourse.julialang.org/u/snowgum)\
**Post date:** [September 15, 2022, 8:36am UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/4 "2022-09-15T08:36:13Z")

</div>

For reference, the solution above with `Interpolations.jl` seems to be about 50% slower than the Python implementation:

```julia
import numpy as np
from scipy.interpolate import interp1d

tmp = np.random.rand(100, 100, 100).astype('f')
x = np.sort(np.random.rand(100).astype('f'))
new_x = np.linspace(x.min(), x.max(), 100).astype('f')
itp = interp1d(x, tmp, axis=-1, kind='linear')

%timeit itp(new_x)
3.08 ms ± 10.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

```

Julia:

```julia
using Interpolations
using BenchmarkTools

tmp = rand(Float32, (100, 100, 100))
x = sort(rand(Float32, 100))
new_x = LinRange(minimum(x), maximum(x), 100)
itp = interpolate((x, 1:100, 1:100), tmp, Gridded(Linear()))

@benchmark itp(new_x, 1:100, 1:100)
BenchmarkTools.Trial: 1059 samples with 1 evaluation.
 Range (min … max): 4.434 ms … 6.687 ms ┊ GC (min … max): 0.00% … 31.21%
 Time (median): 4.619 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 4.724 ms ± 343.316 μs ┊ GC (mean ± σ): 2.09% ± 5.69%

      ▃██▅▂▂▁
  █▅▄▁████████▅▅▁▄▁▄▅▄▁▁▁▁▁▁▁▁▄▁▁▁▁▁▅▅▅▅▇▇▇▇▅▅▆▇▆▅▅▅▅▄▅▅▄▄▄▅▄ █
  4.43 ms Histogram: log(frequency) by time 6.19 ms <

 Memory estimate: 3.83 MiB, allocs estimate: 7.

```

---

<div class="post-metadata">

**Author:** ![aplavin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aplavin/32/222056_2.png) [@aplavin](https://discourse.julialang.org/u/aplavin)\
**Post date:** [September 15, 2022, 9:56am UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/5 "2022-09-15T09:56:10Z")

</div>

It’s common in Julia to make your function work on 1-d arrays (vectors) and then just `map` it over slices of n-d arrays.  
Suppose you already have interpolation working exactly as you need for vectors, implemented in `func(vector)`. Then `mapslices(func, A; dims=1)` should work, or `map(func, eachslice(...))` - depending on how you prefer the result to be shaped.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [September 15, 2022, 4:45pm UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/6 "2022-09-15T16:45:57Z")

</div>

Telling interpolations to not perform interpolation in the second two axes produces better performance on my machine:

```julia
julia> @benchmark itp(new_x, 1:100, 1:100)
BenchmarkTools.Trial: 376 samples with 1 evaluation.
 Range (min … max): 12.599 ms … 19.921 ms ┊ GC (min … max): 0.00% … 0.00%
 Time (median): 12.970 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 13.298 ms ± 1.044 ms ┊ GC (mean ± σ): 0.31% ± 1.15%

  ▄▇█▇▆▄▃▁▁▃ ▁
  ██████████▇██▆▇▄▆▇▆▁▄▁▁▁▁▆▇▄▁▆▁▄▁▁▁▁▁▄▁▁▁▁▁▄▄▁▆▁▁▁▁▁▁▁▁▄▄▁▄ ▆
  12.6 ms Histogram: log(frequency) by time 18.6 ms <

 Memory estimate: 3.83 MiB, allocs estimate: 7.

julia> itp2 = interpolate((x, 1:100, 1:100), tmp, (Gridded(Linear()), NoInterp(), NoInterp()));

julia> @benchmark itp2(new_x, 1:100, 1:100)
BenchmarkTools.Trial: 2001 samples with 1 evaluation.
 Range (min … max): 2.210 ms … 6.351 ms ┊ GC (min … max): 0.00% … 18.00%
 Time (median): 2.311 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 2.485 ms ± 484.836 μs ┊ GC (mean ± σ): 1.90% ± 5.77%

  ██▇▆▅▃▃▂▁ ▂▂ ▁▁ ▁
  █████████████▇▇▇██████▆▇▆▇▄▅▆▁▆▅▅▅▅▆▅▆▅▅▅▄▄▄▅▅▅▅▆▅▄▅▅▅▁▄▁▁▅ █
  2.21 ms Histogram: log(frequency) by time 4.88 ms <

 Memory estimate: 3.82 MiB, allocs estimate: 7.

```

---

<div class="post-metadata">

**Author:** ![snowgum](https://avatars.discourse-cdn.com/v4/letter/s/ecc23a/32.png) [@snowgum](https://discourse.julialang.org/u/snowgum)\
**Post date:** [September 15, 2022, 6:21pm UTC](https://discourse.julialang.org/t/interpolate-along-single-dimension-of-multidimensional-array/87290/7 "2022-09-15T18:21:57Z")

</div>

Thank you! That is a great suggestion. Forgot that by default it was also interpolating in the two other axes.