Nearest Neighbor Interpolation on Matrix to Match Another's Non-Zero Structure

Hello, Julia Community!

I am working with two matrices: C_values, which contains parameter values obtained from an optimization process, and V, which represents observed ice surface velocity (please see the attached image for visualization). C_values has dimensions that match V, but with many NaN entries that I would like to fill using nearest neighbor interpolation based on the non-zero structure of V.


C_values is a sparse representation where only a subset of cells have been calculated (left figure), whereas V is a dense matrix with varying velocities (right figure).

I am looking for an efficient way to perform nearest neighbor interpolation on C_values to ensure that it contains interpolated values at all non-zero points in V. Specifically, I need to interpolate NaN values in C_values only where there are non-zero values in V.

My question are:

  1. How can I implement nearest neighbor interpolation in this context using Julia?
  2. Is there an existing package or function that can handle this selective interpolation directly?

Any pointers or examples would be greatly appreciated!

Thank you in advance!

Do you mean it is a sparse matrix? Or just a regular matrix with NaN values? Either way to me it sounds like you could just loop over all indices and insert elements of V into C_values whenever they are NaN in the latter and nonzero on V? afterwards you could probably just interpolate the result using Interpolations.

for i in eachindex(V, C_values)
   if V[i] != 0 && isnan(C_values[i])
      C_values[i] = V[i]

and then just interpolation the resulting array?

This is a typical problem that is trivial to solve with the GeoStats.jl framework. Take a look at our InterpolateNeighbors transform with the NN model. You can define the CartesianGrid of interpolation, and run the interpolation in one line.

It is all explained in the GDSJL book:

Many other packages provide alternative interfaces to interpolation. You will benefit from our approach if you need to replace grids by more complicated domains in the future.

1 Like

Sorry for my late response, I had some other obligations coming up. I have read through the documents, but I struggle to get the interpolation working. I feel I just do not understand some things good enough.

So I have this matrix C from which we would like to sample the data, I should then define grid = CartesianGrid(size(C)) and somehow geo reference this right and then proceed with the interpolation?

I find it hard to see how this interpolation can be run in one line.
Could you maybe provide a minimal working example?

You have to understand two things:

  1. The domain of interpolation, which is the CartesianGrid with appropriate size, origin, spacing, etc. Check the docstring for details.
  2. The geospatial data, which is the GeoTable with known values of the variable at specific locations. Check the docstring of georef for details.

After you construct these two objects, say grid and geotable, then you can:

geotable |> InterpolateNeighbors(grid) # using nearest neighbors


geotable |> Interpolate(grid) # using all neighbors

You can visualize all these objects with viz and viewer as explained in the book, to make sure that they are where you expect them to be in the coordinate reference system.

1 Like

Okay thanks for the quick response and clarification!

Okay what I do not understand here, is how can I convert the 2D matrix to a suitable DataFrame for georef?

I am running the following, where C is a 95x137 matrix.

using GeoStats, DataFrames
df = DataFrame(C, :auto)
grid = CartesianGrid(size(C))
gt = georef(df, grid)
InterpolateNeighbors(gt, grid)

The above code does not work, because the dataframe df is not setup right.
I am to be honest not so familiar with the use of DataFrame structures within Julia…

Please share a MWE with the data. Have you seen the examples in the docs?

Yes I have seen this example, but I struggle on how to assign the coordinates to the points in the matrix in an automated way.

It is a bit hard to share the matrix itself, as there is a whole inversion process to obtain it. But this is the shape, where I would like to interpolate the white (NaN) values.

One possible solution is:

gtable = georef((; C)) # georeference array C into grid

inds = findall(!isnan, gtable.C) # find indices with data

data = gtable[inds, :] # view of geotable

data |> InterpolateNeighbors(gtable.geometry) # interpolate over grid

Okay I think this is exactly what I was looking for, except for the fact that it is mirrored.


My last question is if there is a straightforward way to extract the interpolated matrix as a matrix{Float64}?

It is mirrored because the convention is that the first dimension of the array is attached to the X axis, the second to the Y axis, etc. This also holds in 3D. You can easily result |> Rotate(pi/2) if you need a rotated visualization.

You can extract columns from geotables as arrays with asarray(result, "C")

Perfect, you really helped me a lot! Thanks for all the effort and the quick responses, it is really appreciated :smile:

1 Like

We introduced a new transform in GeoStats.jl v0.56 called InterpolateNaN to hide these low-level details from end-users. Here is an updated solution:

georef((; C)) |> InterpolateNaN()