# Interpolating along a single dimension of a multi dimensional array for particular points
**URL:** https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308
**Category:** General Usage
**Created:** [September 30, 2019, 4:16am UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308 "2019-09-30T04:16:41Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [September 30, 2019, 4:16am UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/1 "2019-09-30T04:16:41Z")
I have some simple linear interpolation code in a package. It would be nice to support spline etc and not write that myself - but it has to be extremely fast and non-allocating at the point of use. It can take as long as it likes to define an interpolator object as that happens outside the main loop.
The problem is to interpolate along a single dimension (time) of a 3 dimensional long/lat/time array. This happens one lat/long position at a time only a small (but changing) fraction of the lat/long grid is actually calculated, for performance reasons.
Can Interpolations.jl or another package do this without creating a view along the time axis?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 30, 2019, 4:25am UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/2 "2019-09-30T04:25:47Z")
You can broadcast over [GitHub - PumasAI/DataInterpolations.jl: A library of data interpolation and smoothing functions](https://github.com/PumasAI/DataInterpolations.jl) . More can be done but that should be a good start. I don’t think we implemented multiple 1D interpolations at one with this yet, but it would be a good idea.
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [September 30, 2019, 4:39am UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/3 "2019-09-30T04:39:00Z")
Thanks I wasn’t aware of that package, it looks great.
Unfortunately we can’t use `broadcast`, the key optimisation in [DynamicGrids.jl](https://github.com/cesaraustralia/DynamicGrids.jl) is to avoid unnecessary work by not broadcasting the whole array. It would be good to do something like this where the specified dimension(s) are interpolated:
```julia
u = rand(5, 5, 5)
t = 0:4
interp = LinearInterpolation(u, t; dims=3)
interp(1, 2, 3.5)
```
As a side point DataInterpolation.jl could play really well with DimensionalData.jl if it was generalised like this, you could specify the interpolation dimension instead of `t` and `dims`.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 30, 2019, 4:51am UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/4 "2019-09-30T04:51:58Z")
that looks like a fairly reasonable suggestion
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [September 30, 2019, 4:59am UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/5 "2019-09-30T04:59:54Z")
It would also be good to get an intermediate object that can apply an interpolation, like:
```julia
u = rand(5, 5, 5)
t = 0:4
interp = LinearInterpolation(u, t; dims=3)
partialinterp = interp(3.5)
partialinterp(1, 2)
```
Then the calculation of indices could be done before the actual interpolation, outside of main loops.
But I’m not sure how widely this kind of selective/non-broadcasted version would be required.
---
### Author:  [@tim.holy](https://discourse.julialang.org/u/tim.holy)
#### Post date: [September 30, 2019, 8:16pm UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/6 "2019-09-30T20:16:39Z")
Interpolations should be just fine doing this kind of thing. It’s not well documented, but Interpolations (now) does interpolation in the naive way, which IMO is a strength. For example, this is how the internals for linear interpolation in 1d actually work:
```julia
julia> using Interpolations
# Build a "weighted adjacent index" for indexing into any array: compute
# 0.2*a[j] + 0.8*a[j+1], with j=3
julia> widx = Interpolations.WeightedAdjIndex(3, (0.2, 0.8))
Interpolations.WeightedAdjIndex{2,Float64}(3, (0.2, 0.8))
julia> r = 1:10
1:10
julia> r[widx] # you can index arbitrary arrays with these weird indexes
3.8000000000000003
# Now do something kinda silly, but it shows you're not constrained in what you can do
julia> widx = Interpolations.WeightedAdjIndex(3, (0.0, 0.0, 0.0, 0.7))
Interpolations.WeightedAdjIndex{4,Float64}(3, (0.0, 0.0, 0.0, 0.7))
julia> r[widx]
4.199999999999999
julia> 0.7*r[6] # the 4 weights apply to positions 3, 4, 5, and 6, hence the r[6]
4.199999999999999
# Try a 2d example
julia> A = reshape(1:20, 4, 5)
4×5 reshape(::UnitRange{Int64}, 4, 5) with eltype Int64:
1 5 9 13 17
2 6 10 14 18
3 7 11 15 19
4 8 12 16 20
# Try indexing the first dimension with an integer and the 2nd with a weighted index
julia> i, j = 2, Interpolations.WeightedAdjIndex(2, (0.3, 0.7))
(2, Interpolations.WeightedAdjIndex{2,Float64}(2, (0.3, 0.7)))
julia> A[i,j]
8.8
julia> 0.3*A[i,2] + 0.7*A[i,3] # see, same as linear interpolation along dimension 2
8.8
# Interpolate in both dimensions, this time trying out the other index type, WeightedArbIndex.
# This type allows the index locations to be arbitrary, rather than incremental. It is used internally to
# implement things like periodic boundary conditions, reflecting boundary conditions, etc.
julia> i = Interpolations.WeightedArbIndex((1, 3), (-0.1, 100.0))
Interpolations.WeightedArbIndex{2,Float64}((1, 3), (-0.1, 100.0))
julia> A[i, j]
979.2199999999999
# Indeed, this is the full 2d interpolation
julia> A[1,2]*(-0.1*0.3) + A[1,3]*(-0.1*0.7) + A[3,2]*(100*0.3) + A[3,3]*(100*0.7)
979.22
```
As you can see, Interpolations puts all the magic into the index objects, not the arrays. This makes it very flexible. Moreover, if you put stuff inside `@inbounds` then the generated code should be very efficient.
Interpolations does include a `NoInterp` mode as an exported interface for doing these kinds of things, but the real magic is in the indexes.
---
### Author:  [@Raf](https://discourse.julialang.org/u/Raf)
#### Post date: [September 30, 2019, 11:56pm UTC](https://discourse.julialang.org/t/interpolating-along-a-single-dimension-of-a-multi-dimensional-array-for-particular-points/29308/7 "2019-09-30T23:56:34Z")
That’s perfect, thanks for the detailed example. I agree that kind naive composability is a real strength. Most of the interpolation I need day to day (in dynamic ecological models) are only for some subset of an array, optimizations on broadcast just cant compete with not doing 90% of the work.
It might be worth putting this in the docs.
Edit: I’m wondering if splines work in a similar way? It would be good to be able to structure a package so that users can easily switch between rounding, linear interpolation and splines - depending on their accuracy/performance tradeoff.