# Interpolations of a function that has NaN
**URL:** https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704
**Category:** Numerics
**Tags:** question, interpolations
**Created:** [March 7, 2023, 11:26pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704 "2023-03-07T23:26:25Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [March 7, 2023, 11:26pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/1 "2023-03-07T23:26:25Z")
I am trying to interpolate a 2d function that is defined only on a partial set of my parameters.
Where it is not defined, my 2d array has `NaN`s. Note that in my example the range of where the function is defined is not a rectangle, but it is a function of x, y.
When I perform the interpolation, I get a function that has all `NaN`s. Is there a way to overcome this issue?
I am using [Interpolations.jl](https://juliamath.github.io/Interpolations.jl/dev/), and my code is like the [example](https://juliamath.github.io/Interpolations.jl/dev/control/#Control-of-interpolation-algorithm) in the documentation.
Here is a minimal example which includes `NaN`s:
```julia
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]
A[1:2, :] .= NaN
A[:, 1:2] .= NaN
itp = interpolate(A, BSpline(Cubic(Line(OnGrid()))))
```
The result of `itp` is all `NaN`.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 8, 2023, 12:12am UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/2 "2023-03-08T00:12:14Z")
Instead of using a matrix with NaNs, one possibility is to take the remaining valid entries as scattered points. See the example below using ScatteredInterpolation.jl.
> **ScatteredInterpolation.jl code**
>
> ```julia
> using ScatteredInterpolation, StaticArrays, Plots
>
> 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]
> A[50, :] .= NaN
> A[:, 20] .= NaN
>
> CI = findall(!isnan, A)
> points = SVector{2,Float64}[]
> samples = Float64[]
> for ci in CI
> push!(points, SA[A_x1[ci[1]], A_x2[ci[2]]])
> push!(samples, A[ci])
> end
> points = reduce(hcat, points)
>
> nx1, nx2 = length(A_x1), length(A_x2)
> X1, X2 = vec(repeat(A_x1, nx2)), vec(repeat(A_x2', nx1))
> gridPoints = [X1 X2]'
>
> itp = ScatteredInterpolation.interpolate(Multiquadratic(), points, samples)
> interpolated = evaluate(itp, gridPoints)
> gridded = reshape(interpolated, nx1, nx2)
>
> plot(heatmap(A), heatmap(gridded))
>
> ```
The packages GMT.jl and DIVAnd.jl are also very good for this task: [see this other post](https://discourse.julialang.org/t/plot-3d-data-in-a-2d-contour-plot/71251).
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [March 9, 2023, 3:34pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/3 "2023-03-09T15:34:16Z")
Thanks for your response!
It has been useful, but I am still not sure how to deal with my problem.
My function f(x, y) is defined on some non-rectangular grid. Namely, there is a line y\_{edge}(x) (known only numerically) that gives the boundary of where my function f(x, y) is defined. So when I calculate the values f(x,y) on some 2d grid I put `NaN` where the function is not defined.
The best scenario for me would be to simply extrapolate the function f(x,y) where it is not defined.
I tried to fill all the `NaN` values with flat values outside the region where the function is defined and to do the interpolations, but that does not work very well.
If you have other ideas on how to do that, I would love to hear!
---
### Author:  [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)
#### Post date: [March 9, 2023, 4:00pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/4 "2023-03-09T16:00:21Z")
With DIVAnd you can define a mask for your grid onto which you want to interpolate and the data points can be scattered. So basically you defined a 2D grid and mask it for the regions where you know your function is not defined. Only THEN you ask DIVAnd to interpolate and the interpolation will taken into account the mask. So you already created that mask in your tests and only need to use the (scattered) data points you actually have for the interpolation.
Also is your edge function known at the resolution of the interpolation grid ? Otherwise you probably need to interpolate it first to the right resolution too.
Hope it makes sense, otherwise I will provide a MWE.
---
### Author:  [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)
#### Post date: [March 9, 2023, 4:34pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/5 "2023-03-09T16:34:27Z")
```julia
using DIVAnd
using PyPlot
using Statistics
NX=300
NY=300
len=0.1
ND=NX
# Fine grid
xg=collect(range(0,stop=1,length=NX))
yg=collect(range(0,stop=1,length=NY))
# Boundary on coarse grid
yedge=[0.1,0.2,0.3,0.2,0.2]
xedge=[0.1,0.3,0.5,0.7,0.9]
# Interpolate fine boundary
yedgefine,s=DIVAndrun(trues(NX),(ones(size(xg)) / (xg[2]-xg[1]),),(xg,),(xedge,),yedge.-mean(yedge),len,0.01)
yedgefine=yedgefine.+mean(yedge)
figure()
plot(xg,yedgefine,"-",xedge,yedge,"o")
figure()
# Sample a function
fun(x,y) = 2*(sin.(6x) * cos.(6y)) .+ (0 .- 1) .* x .* y
x = 0.5.+0.25.*randn(ND);
y = 0.5 .+ 0.25 .* randn(ND);
f = fun.(x,y)+0.1*randn(ND)
# create the grid for the interpolation
xi,yi= DIVAnd.ndgrid(xg,yg)
# all points are valid points except those below the edge
mask = trues(size(xi));
for j=1:NY
for i=1:NX
if yi[i,j]
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 10, 2023, 7:14pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/6 "2023-03-10T19:14:04Z")
@JM_Beckers, could you please verify your code as it is giving errors?
After correcting the `espilon2` typo, the subsequent command:
```julia
fi, s = DIVAndrun(mask, (pm,pn), (xi,yi), (x,y), f, (len,len), epsilon2)
```
throws:
```julia
ERROR: CanonicalIndexError: setindex! not defined for LazyGrids.GridAV{Float64, 1, 2}
Stacktrace:
[1] error_if_canonical_setindex(::IndexCartesian, ::LazyGrids.GridAV{Float64, 1, 2}, ::Int64, ::Int64)
```
Using Julia 1.8.5 on Windows and DIVAnd v2.7.9.
---
### Author:  [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)
#### Post date: [March 10, 2023, 9:25pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/7 "2023-03-10T21:25:14Z")
I think there is a clash for function `ndgrid` which is also defined in [Methods · LazyGrids.jl](https://juliaarrays.github.io/LazyGrids.jl/stable/methods/#LazyGrids.ndgrid-Tuple%7BInt64,%20Vararg%7BInt64%7D).
Probably you loaded that one ? I changed to `DIVAnd.ndgrid` to make sure you use the right version. Let me know if it works.
(I also corrected the typo, thanks)
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [March 10, 2023, 10:19pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/8 "2023-03-10T22:19:35Z")
@JM_Beckers, thank you, you were spot on.
---
### Author:  [@roi.holtzman](https://discourse.julialang.org/u/roi.holtzman)
#### Post date: [March 10, 2023, 10:47pm UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/9 "2023-03-10T22:47:34Z")
Thanks for this wonderful example!
I tried to implement it in my case, but it did not work for me.
I get the following error:
```julia
MethodError: no method matching DIVAndrun(::BitMatrix, ::Tuple{Matrix{Float64}, Matrix{Float64}}, ::Tuple{Matrix{Float64}, Matrix{Float64}}, ::Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, ::typeof(Sλq_2d_grid_interpolate_function), ::Tuple{Float64, Float64}, ::Float64)
```
For me the function `f` in your code is a function that is defined using `Interpolations.jl`.
To have a function and not an interpolations type, I defined
`f(x, y) = interpolations_f(x,y)`.
Still, it does not work.
Do you have any idea of how to make it work?
---
### Author:  [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)
#### Post date: [March 11, 2023, 1:29am UTC](https://discourse.julialang.org/t/interpolations-of-a-function-that-has-nan/95704/10 "2023-03-11T01:29:11Z")
No, `f` is an array containing the scattered data values located in `(x,y)`, where `x` and `y` are the arrays containing the two coordinate values of each value of `f`.
In the MWE they are created by a random set of coordinates and a know function,
```julia
x = 0.5.+0.25.*randn(ND);
y = 0.5 .+ 0.25 .* randn(ND);
f = fun.(x,y)+0.1*randn(ND)
```
but you can provide the arrays you want. Basically you tell DIVAnd in which points you know the function value and which is the value.