<p>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 <a href="https://en.wikipedia.org/wiki/Natural_neighbor_interpolation" rel="nofollow noreferrer">Natural neighbour interpolation</a> (that is a class of weighted distance interpolation as suggested in previous comments). The usage is like this:</p>

<pre><code>% 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;
</code></pre>

<p>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 <a href="https://en.wikipedia.org/wiki/Kriging" rel="nofollow noreferrer">Kriging interpolation</a>.
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).</p>

CC BY-SA 3.0
2018-01-07T10:41:35.993Z
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!
Semidiurnal Simon
CC BY-SA 3.0
2018-01-07T15:22:05.640Z
@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.
Camilo Rada
CC BY-SA 3.0
2018-01-09T07:21:24.617Z
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 ;-)
Semidiurnal Simon
CC BY-SA 3.0
2018-01-05T17:04:07.020Z
2018-01-05T17:04:07.020Z
Camilo Rada
<p>I have a bathymetry dataset that has a small section missing. I need to fill in this section.</p>

<p>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.</p>

<p><a href="https://i.sstatic.net/AELcz.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/AELcz.png" alt="enter image description here"></a></p>

<p>Is there an appropriate method to estimate the missing values that gives more weight to known values along a given axis?</p>

<p>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.</p>

CC BY-SA 3.0
2018-01-05T10:24:21.700Z
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.
Semidiurnal Simon
CC BY-SA 3.0
2018-01-05T16:30:01.760Z
Ultimately you have to use a weighted pseudo-distance in the covariance matrix. I have used some of the approaches in: http://www.sciencedirect.com/science/article/pii/S0278434301000073. You can add anisotropy in the covariance equations by weighting one direction more than the other.
arkaia
CC BY-SA 3.0
2018-01-05T16:36:29.990Z
A package that does an okay job with this and it is relatively efficient is: http://www.dfo-mpo.gc.ca/Library/326156.pdf
arkaia
CC BY-SA 3.0
2018-01-07T10:24:22.220Z
@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?
Semidiurnal Simon
CC BY-SA 3.0
2018-01-07T13:36:13.980Z
Sorry. Try http://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).
arkaia
CC BY-SA 3.0
2018-01-18T15:57:54.447Z
2D objective mapping is a straightforward method of interpolation that allows you to specify two orthogonal directions and their respective decorrelation scales.  Usually these parameters are calculated from the data, but you can plug in whatever you like.  The typical application is more mapping scattered data onto a rectangular grid, but you could use it here.
Ben51
<data-analysis><bathymetry>
Best way to interpolate bathymetry with a directional component?
Semidiurnal Simon
# Best way to interpolate bathymetry with a directional component?

(Asked by: username_0 on 2018-01-05)

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.

[![enter image description here](https://i.sstatic.net/AELcz.png)](https://i.sstatic.net/AELcz.png)

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.

**Tags:** data-analysis, bathymetry

> **Comments:**
> **username_0**: 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.  
> **username_1**: Ultimately you have to use a weighted pseudo-distance in the covariance matrix. I have used some of the approaches in: http://www.sciencedirect.com/science/article/pii/S0278434301000073. You can add anisotropy in the covariance equations by weighting one direction more than the other.  
> **username_1**: A package that does an okay job with this and it is relatively efficient is: http://www.dfo-mpo.gc.ca/Library/326156.pdf  
> **username_0**: @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?  
> **username_1**: Sorry. Try http://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).  
> **username_2**: 2D objective mapping is a straightforward method of interpolation that allows you to specify two orthogonal directions and their respective decorrelation scales.  Usually these parameters are calculated from the data, but you can plug in whatever you like.  The typical application is more mapping scattered data onto a rectangular grid, but you could use it here.  

---

## 1 answers

### Answer by username_3 on 2018-01-05. Score: 2

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](https://en.wikipedia.org/wiki/Natural_neighbor_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](https://en.wikipedia.org/wiki/Kriging).
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).

> **Comments:**
> **username_0**: 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!  
> **username_3**: @username_0 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.  
> **username_0**: 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 ;-)  

---

