# Interpolation on non-rectangular domain

**URL:** <https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568>\
**Category:** Numerics\
**Created:** [December 9, 2020, 10:11pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568 "2020-12-09T22:11:03Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [December 9, 2020, 10:11pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/1 "2020-12-09T22:11:03Z")

</div>

Dear all,

I would like to ask for advice, how to interpolate function on a non-rectangular domain. I am working on a problem, where the domain of my function is 2 or 3 dimensional, where all variables are restricted to be in [0,1] and their sum has to be in [0,1], hence this problem produces a non-rectangular domain. Is there some convenient package for this type of problem? As far as I know (maybe I am wrong), Interpolations.jl works with rectangular domains.

Best,  
Honza

 ![Výstřižek](https://global.discourse-cdn.com/julialang/original/3X/6/d/6de5654fba6ef6902710c0a9503dde13eda9f365.png)

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [December 9, 2020, 10:18pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/2 "2020-12-09T22:18:43Z")

</div>

Sounds like interpolation on triangles and tetrahedral domains. Do you want piecewise linear or high order polynomial interpolation?

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 9, 2020, 10:19pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/3 "2020-12-09T22:19:45Z")

</div>

You can perform such interpolations with [GeoStats.jl](https://github.com/JuliaEarth/GeoStats.jl). If you are looking for a lightweight dependency, other members can suggest alternatives.

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [December 9, 2020, 10:21pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/4 "2020-12-09T22:21:55Z")

</div>

Higher-order would be nice, but linear is also perfectly ok.

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [December 9, 2020, 10:26pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/5 "2020-12-09T22:26:43Z")

</div>

Is there some tutorial about interpolations in your package? What I need is to interpolate values of my function on a non-rectangular grid, and then use this interpolation as a representation of this function in further computations (I am using it in an iterative algorithm that solves for the unknown function). Compile time is ok for me, I am more concerned about runtime.

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 9, 2020, 10:59pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/6 "2020-12-09T22:59:00Z")

</div>

Say you have arbitrary locations with measurements of a variable `z`:

```julia
julia> using GeoStats
        
julia> data = georef((z=rand(100),), rand(2,100))
100 PointSet{Float64,2}
  variables
    └─z (Float64)

```

And another set of arbitrary locations where you want to perform the interpolation:

```julia
julia> locs = PointSet(rand(2,100))
100 PointSet{Float64,2}
 0.8877317040268236 0.5228607963917429 … 0.7112810163212397
 0.5372754529043251 0.7050593626184458 0.6308336193731663

```

You can define the “interpolation” problem:

```julia
julia> problem = EstimationProblem(data, locs, :z)
2D EstimationProblem
  data: 100 SpatialData{Float64,2}
  domain: 100 PointSet{Float64,2}
  variables: z (Float64)

```

and pick one of the solvers like for example inverse distance weighting (IDW):

```julia
julia> solution = solve(problem, IDW())
2D EstimationSolution
  domain: 100 PointSet{Float64,2}
  variables: z

```

These implementations should have a decent performance as they are using KDTrees from NearestNeighbors.jl and other tricks. But I never used them in the context of a tight loop like it is usually done in function approximation inside finite element codes for example.

You may also find [ApproxFun.jl](https://github.com/JuliaApproximation/ApproxFun.jl) useful depending on what you want to do.

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [December 10, 2020, 12:04am UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/7 "2020-12-10T00:04:44Z")

</div>

@juliohm Thank you very much! Would this be efficient for scalar operations? Systems that I am solving are typically implicit systems (macroeconomic discrete time functional equations), and hence I need to use at each grid point nonlinear solver…

If I would wrap EstimationProblem and solve as a scalar function, would it be efficient, or would the overhead be too high?

Best,  
Honza

---

<div class="post-metadata">

**Author:** ![juliohm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juliohm/32/215266_2.png) [@juliohm](https://discourse.julialang.org/u/juliohm)\
**Post date:** [December 10, 2020, 1:08am UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/8 "2020-12-10T01:08:36Z")

</div>

You are welcome @Honza9723. I’d benchmark the code with some typical setup for your use case. You can easily quantify possible bottlenecks if they exist. Take a look at BenchmarkTools.jl. The problem setup is cost-less (it is just a wrapper pointing to the other variables). The solvers have a compilation time, but runtime should be fine. The major downside I see is that GeoStats.jl is a big project to have as a dependency, and we are constantly evolving the API.

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [December 10, 2020, 1:14am UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/9 "2020-12-10T01:14:02Z")

</div>

That sounds nice. Yeah, the compilation time of GeoStats.jl is quite large, but I can live with it. Thanks again!

---

<div class="post-metadata">

**Author:** ![jlchan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jlchan/32/10958_2.png) [@jlchan](https://discourse.julialang.org/u/jlchan)\
**Post date:** [December 10, 2020, 4:48am UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/10 "2020-12-10T04:48:30Z")

</div>

Sorry for the later reply, was actually working on the package I’m using for this answer. I think @juliohm’s approach is probably better suited (esp if you can’t choose the points at which you interpolate), but I figure I’d include this in case it’s helpful to anyone.

I have a package NodesAndModes.jl for high order polynomials and interpolation points on several shapes. A warning: I found a bug in the main branch, so I’m using a branch where I fixed this by making it more “Julian”.

You can get high order interpolation nodes for a Tri() or Tet() domain via

```julia
pkg> add NodesAndModes#dev_elem_types
julia> using NodesAndModes
julia> using Plots
julia> x,y = nodes(Tri(),25) # degree 25 interpolation nodes

```

This gives some “decent” interpolation nodes.

 ![Screen Shot 2020-12-09 at 10.37.35 PM](https://global.discourse-cdn.com/julialang/original/3X/4/5/451106b2239dc5930bc6d4954ca5be5e6b1adb6d.jpeg)

I do the actual interpolation using linear algebra (implicitly defined Lagrange basis functions). If I want to interpolate, say, a version of Runge’s function in 2D

 ![Screen Shot 2020-12-09 at 10.38.57 PM](https://global.discourse-cdn.com/julialang/original/3X/d/7/d7717eeb033b4e491433e6cadf364f61a43dcb9d.jpeg)  
with a degree N=20 interpolant, I can do this using

```julia
using NodesAndModes
N = 20 # polynomial degree
x,y = nodes(Tri(),N)

xp,yp = equi_nodes(Tri(),100) # plotting nodes
f(x,y) = 1/(1 + 5*((x+1/3)^2+(y+1/3)^2)) # Runge-ish function

function interp(elem,N,x_eval,y_eval,x,y,f_vals)
    coeffs = vandermonde(elem,N,x,y)\f_vals
    return vandermonde(elem,N,x_eval,y_eval)*coeffs
end

f_interp = interp(Tri(),N,xp,yp,x,y,f.(x,y))

err = maximum(abs.(f.(xp,yp)-f_interp))
@show err

```

which gives an error of 0.015041654247017144.

 ![Screen Shot 2020-12-09 at 10.38.38 PM](https://global.discourse-cdn.com/julialang/original/3X/e/2/e2f0d6cc3d6b0c517c87026723b60051edb26a35.jpeg)

Runge’s function is sort of the worst-case scenario. If you interpolate a smoother function like `f(x,y) = sin(pi*x)*sin(pi*y)`, then the error is 1.167277052793736e-9.

A few notes

- if the points at which you want to interpolate are fixed, you can just replace the interpolation `x,y` in the above with your node locations. However, due to the nature of high order interpolation, this can lead to pretty poor results.
- degree N \> 25 starts to get iffy since some of the polynomials get really large and roundoff effects come into play. For many “nice” functions, N between 15-20 is more than enough though.
- if your evaluation and interpolation points don’t change between iterations, you can reduce the interpolation procedure to just a matrix-vector product, which should be pretty fast.

Finally, for completeness, the plots were generated using  
`scatter(xp,yp,f_interp,zcolor=f_interp,msw=0,leg=false,title="Interpolant")` and `scatter(xp,yp,f.(xp,yp),zcolor=f.(xp,yp),msw=0,leg=false,title="Exact")`

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [December 10, 2020, 9:14am UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/11 "2020-12-10T09:14:42Z")

</div>

@jlchan Thank you! I will definitely try both approaches!

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [December 10, 2020, 9:46am UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/12 "2020-12-10T09:46:05Z")

</div>

FastTransforms.jl has a fast transform for triangle OPs:

[https://github.com/JuliaApproximation/FastTransforms.jl/blob/master/examples/triangle.jl](https://github.com/JuliaApproximation/FastTransforms.jl/blob/master/examples/triangle.jl)

Note it’s only technically interpolation if your function is a polynomial (the number of sample points is double the degrees of freedom), but for function _approximation_ it works great even with millions of unknowns

---

<div class="post-metadata">

**Author:** ![MikaelSlevinsky](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikaelslevinsky/32/38545_2.png) [@MikaelSlevinsky](https://discourse.julialang.org/u/MikaelSlevinsky)\
**Post date:** [December 10, 2020, 6:12pm UTC](https://discourse.julialang.org/t/interpolation-on-non-rectangular-domain/51568/13 "2020-12-10T18:12:08Z")

</div>

Note that docs are generated for examples via Literate.jl (for prettier output) [Calculus on the reference triangle · FastTransforms.jl](https://juliaapproximation.github.io/FastTransforms.jl/dev/generated/triangle/)