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:
How can I implement nearest neighbor interpolation in this context using Julia?
Is there an existing package or function that can handle this selective interpolation directly?
Any pointers or examples would be greatly appreciated!
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]
end
end
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.
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.
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?
The domain of interpolation, which is the CartesianGrid with appropriate size, origin, spacing, etc. Check the docstring for details.
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
or
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.
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…
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.
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
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")
