# 2D interpolation of large data to calculate depth profiles
**URL:** https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190
**Category:** Visualization
**Tags:** plotting, interpolations, makie
**Created:** [June 25, 2024, 7:59am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190 "2024-06-25T07:59:09Z")
**Posts on this page:** 13
**Page:** 2
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [June 25, 2024, 10:46am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/21 "2024-06-25T10:46:41Z")
Im surprised that you didnt find GeoStats.jl in your searches. We easily interpolate more than 200Million points in (3D) industrial applications.
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [June 25, 2024, 11:10am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/22 "2024-06-25T11:10:19Z")
In general, interpolation only depends on the local values. So in any methods used, subsampling to the local area (chosen using some nearest neighbor method - see NearestNeighbors.jl) and then interpolating using the methods previously used will work.
For the example in the OP:
```julia
# provide local interpolation
function local_itp(p)
pidx = searchsortedfirst.((x,y),p)
xrng, yrng = (:).(clamp.(pidx .- 5, 1, 100),
clamp.(pidx .+5, 1, 100))
itp = linear_interpolation((x[xrng],y[yrng]),z[xrng,yrng])
itp(p...)
end
```
```julia
# use the local interpolation
fig = Figure(size=(900, 400), fontsize=20)
axs = [Axis(fig[1, j], aspect=1) for j in 1:2]
cmap = :roma
contour!(axs[1], x, y, z, levels=30, colormap=cmap)
linesegments!(axs[1], [A, B], color=:red, linewidth=5)
lines!(axs[2],
range(0, distance; length=n_profile_samplepoints),
[local_itp(p) for p in coords])
fig
```
---
### Author:  [@tamasgal](https://discourse.julialang.org/u/tamasgal)
#### Post date: [June 25, 2024, 11:11am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/23 "2024-06-25T11:11:35Z")
Oh OK, sounds like something I’d need. I will have a look. I searched for 2D interpolations and thought Interpolations.jl would cut it…
I got it working with DIVAnd.jl but it seems that the interpolation is too coarse. Somehow it caps at 100x100 in the interpolation backend?
---
### Author:  [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)
#### Post date: [June 25, 2024, 11:13am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/24 "2024-06-25T11:13:14Z")
That is because it uses default settings. If you want higher resolutions, you can specify the underlying output grid using something as `xi=(collect(range(minimum(x),maximum(x),1000)),collect(range(minimum(y),maximum(y),1000))) ` as keyword and look at the general DIVAnd documentation
---
### Author:  [@tamasgal](https://discourse.julialang.org/u/tamasgal)
#### Post date: [June 25, 2024, 11:21am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/25 "2024-06-25T11:21:45Z")
> [@Dan](#):
>
> For the example in the OP:
>
> ```julia
> # provide local interpolation
> function local_itp(p)
> pidx = searchsortedfirst.((x,y),p)
> xrng, yrng = (:).(clamp.(pidx .- 5, 1, 100),
> clamp.(pidx .+5, 1, 100))
> itp = linear_interpolation((x[xrng],y[yrng]),z[xrng,yrng])
> itp(p...)
> end
>
> ```
Well I see my mistake that I forgot to emphasise that I am not able to create `z` (the 2D depth map) due to memory restrictions. This is the reason why things get tricky very quickly.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [June 25, 2024, 11:25am UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/26 "2024-06-25T11:25:38Z")
@tamasgal check the InterpolateNeighbors transform in the GeoStats.jl docs or in the book. It will do the job using the NearestNeighbors.jl approach.
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [June 25, 2024, 3:01pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/27 "2024-06-25T15:01:36Z")
> [@tamasgal](#):
>
> This assumes that I have `z`, a 2D matrix of my z-data, right?
From your responses your input consists of 14347538 scattered points represented by 3 arrays `x, y, z` of length 14347538 each. This does not look too big and no need for a matrix. Have I got it correct?
---
### Author:  [@tamasgal](https://discourse.julialang.org/u/tamasgal)
#### Post date: [June 25, 2024, 3:08pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/28 "2024-06-25T15:08:04Z")
Yep, that’s correct 🙂
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [June 25, 2024, 3:27pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/29 "2024-06-25T15:27:20Z")
But if you have it interpolated into a grid, you may have a use for the goodies in GMT’s [grdtrack](https://www.generic-mapping-tools.org/GMTjl_doc/documentation/modules/grdtrack/index.html#grdtrack) (stacked profiles, crossprofiles profiles, etc)
EDIT: Is this multi-beam data?
---
### Author:  [@tamasgal](https://discourse.julialang.org/u/tamasgal)
#### Post date: [June 25, 2024, 6:04pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/30 "2024-06-25T18:04:06Z")
Thanks, I will try that too!
No multi-beam, neutrinos through Earth from far outside our galaxy 😉
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [June 25, 2024, 8:47pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/31 "2024-06-25T20:47:19Z")
FWIW, ScatteredInterpolation.jl takes about 1 minute to run on my Windows laptop for 14347538 points. However, the results interpolated along the cyan line don’t look great in the example below:
> **ScatteredInterpolation code**
>
> ```julia
> # 1 - Generate input data
> peaks(x, y) = 3(1-x)^2*exp(-x^2-(y+1)^2) -10(x/5-x^3-y^5)*exp(-x^2-y^2) - exp(-(x+1)^2-y^2)/3
> N = 14347538
> x, y = 6*(rand(N) .- 0.5), 6*(rand(N) .- 0.5)
> z = peaks.(x, y)
>
> # 2 - Interpolate scattered data
> using ScatteredInterpolation
> # interpolate z function of x, y
> itp_z = interpolate(Shepard(), permutedims([x y]), z)
>
> coords = [[-3,-3], [3,3]] # end points of line along which interpolate data
> Ni = 200 # number of equidistant points to interpolate
> xi = LinRange(coords[1][1], coords[2][1], Ni)
> yi = LinRange(coords[1][2], coords[2][2], Ni)
> zi = evaluate(itp_z, permutedims([xi yi]))
>
> # 3 - Plot results
> using Plots; gr(dpi=600)
> p1 = contourf(sort(x[1:5000:end]), sort(y[1:5000:end]), peaks, ratio=1, lims=(-3,3))
> plot!(first.(coords), last.(coords), c=:cyan)
> scatter!(x[1:5000:end], y[1:5000:end], c=:red, msw=0, ms=1)
> p2 = plot(xi, peaks.(xi, yi), lw=5, c=:blues, label="Exact")
> plot!(xi, zi, lw=2, ls=:dash, c=:black, label="Shepard")
> plot(p1, p2, size=(1000,600))
>
> ```
_PS:
there are 5000x more points than displayed, for clarity purposes._
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [June 25, 2024, 9:57pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/32 "2024-06-25T21:57:15Z")
This thread lists a number of packages you could try (some of which have already been mentioned above): [2D Interpolation on an irregular grid](https://discourse.julialang.org/t/2d-interpolation-on-an-irregular-grid/105247)
---
### Author:  [@joa-quim](https://discourse.julialang.org/u/joa-quim)
#### Post date: [June 25, 2024, 10:14pm UTC](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190/33 "2024-06-25T22:14:11Z")
This takes ~5 sec on my computer (x,y,z from Rafael’s post above)
```julia
using GMT
@time G = gridit([x y z], limis=(-3,3,-3,3), method="nearneighbor", inc=0.1, search_radius=0.1)
4.803249 seconds (453 allocations: 328.435 MiB, 0.79% gc time)
```
but, if I understood it, the OP’s problem was memory available.
[Previous page](https://discourse.julialang.org/t/2d-interpolation-of-large-data-to-calculate-depth-profiles/116190.md?page=1)