# Best way to interpolate bathymetry with a directional component?

I have a bathymetry dataset that has a small section missing. I need to fill in this section.

The bathymetry in this region has some strongly directional features (e.g. a canyon that runs through the missing region). Using a standard interpolation function would weaken these features. Is there an appropriate method to estimate the missing values that gives more weight to known values along a given axis?

At the moment the data is a GeoTIFF, and I'm using QGIS, but I can convert as necessary. I have straightforward access to MATLAB and R, as well as all the GRASS and SAGA tools included in QGIS.

• Asking this here because the application is firmly in earth science. I realise that there might be more answers over on gis.stackexchange; if folk can't answer here, I may request migration. Jan 5 '18 at 10:24
• Ultimately you have to use a weighted pseudo-distance in the covariance matrix. I have used some of the approaches in: sciencedirect.com/science/article/pii/S0278434301000073. You can add anisotropy in the covariance equations by weighting one direction more than the other. Jan 5 '18 at 16:30
• A package that does an okay job with this and it is relatively efficient is: dfo-mpo.gc.ca/Library/326156.pdf Jan 5 '18 at 16:36
• @aretxabaleta thanks for all the links. The last one doesn't work - I think it's a logged-in link. Can you try again, or give the title (or ideally, DOI) of the item? Jan 7 '18 at 10:24
• Sorry. Try onlinelibrary.wiley.com/doi/10.1029/97JC00697/pdf The title is: "Adaptation of a fast optimal interpolation algorithm to the mapping of oceanographic data" by Menemenlis et al. (JGR, 1997). Jan 7 '18 at 13:36

I had the same problem with surface DEM's. The best solution I found in Matlab was using the scatteredInterpolant class, it is inbuilt in Matlab. It allows Natural neighbour interpolation (that is a class of weighted distance interpolation as suggested in previous comments). The usage is like this:

% IF X and Y are the coords of each grid cell (you can use meshgrid to get
% those)
invalidValue=NaN; % or -99999 or whatever is in your case
gapIdxs= DEM==invalidValue; %Indexes of the void cells
interpolant = scatteredInterpolant(X(~gapIdxs),Y(~gapIdxs),DEM(~gapIdxs),'natural');
% interpolant is now a function, and we can evaluate it at any point you
% want, in this case we evaluate it in the void cells
interpolatedValues=interpolant(X(gapIdxs),Y(gapIdxs));
% Now we patch the DEM
patchedDEM=DEM;
patchedDEM(gapIdxs)=interpolatedValues;


That produce a result way better than cubic interpolation in my experience. If that is not enough, before you start doing manual picking of features you can explore Kriging interpolation. Finally if you what to force directionality, you can rotate your data (using matlab imrotate) so that your main direction is vertical and then interpolate column by column (using matlab interp1).

• I'm often a fan of natural neighbour interpolation, but in this case it'll still give equal weight to information from all directions, and thus would tend to "fill in" the canyon in the west of the image. You make a good point that with the small north-south size of the gap I have to fill, it may be good enough to simply to a 1-D interpolation across the gap in the north-south direction. If I can't get my head around the more sophisticated approach in the comments above, I'll probably do this! Jan 7 '18 at 10:41
• @SemidiurnalSimon I don't think it would fill much the canyon, because in that area it will give canyon pixels a much larger weight that that of pixels outside the canyon. If you upload the sample data I can give it a try. Jan 7 '18 at 15:22
• Hmm, you make a fair point - given the scales and resolution involved here, an omnidirectional approach with a sufficiently small search radius is probably OK. A criticism that I have of natural neighbour is that with dense data, it forces a very small "input" radius, but here that might not be a bad thing ;-) Jan 9 '18 at 7:21