I am working with a dataset which is essentially a depth sample map: x and y coordinates on an irregular grid and the corresponding depth values. I need to calculate tons of depth profiles for different lines defined by two points `A`

and `B`

.

I thought interpolating would be an easy task, as I have done it so many times but I realised that the dataset is too large (millions of points) and I get an `OutOfMemoryError`

from Interpolations.jl.

Since this is a fairly common task to create these “depth profiles” (but not my expertise at all) I thought I quickly hack together a MWE and ask around if there is a sophisticated solution out there, before I naively hack together something

Below is my MWE. The issue is that `n_datapoints`

grows (at least) with O(n^2) in memory. My dataset has millions of points and the example below already starts to hang with a few thousands of entry (which is no surprise, given how the interpolator works).

**Is there any state of the art solution for this kind of problem which is maybe already implemented somewhere? Like “lazy interpolations” or so?** At least this is implemented in basically every map software which can calculate z profiles for a drawn line in real-time

Another approach would be: based on `A`

and `B`

, only select those datapoints which are in the road width of the line between them, however this requires some sorting and can become a bit of a hassle when dealing with lots of `A`

and `B`

pairs.

Here is the Pluto.jl notebook in case you prefer that:

2D Interpolation of Depth Profiles.jl (49.3 KB)

```
using LinearAlgebra
using CairoMakie
using Interpolations
n_datapoints = 100 # this is the sensitive parameter
n_profile_samplepoints = 50
x = range(-2, 2; length=n_datapoints)
y = range(-2, 2; length=n_datapoints)
f(z) = (z^3 - 3) / z
z = [angle(f(u + 1im * v)) for u in x, v in y]
# two points to define a z-cut
A = Point2f(1.0, 1.5)
B = Point2f(-1.5, -1.0)
# the distance so that we can scale the profile accordingly
distance = norm(A - B)
# these will be our sample points for the interpolator
coords = [(A[1] + t * (B[1] - A[1]), A[2] + t * (B[2] - A[1])) for t in range(0, 1; length=n_profile_samplepoints)]
itp = linear_interpolation((x, y), z)
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), [itp(p...) for p in coords])
fig
```