3
$\begingroup$

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

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Semidiurnal Simon
    Commented Jan 5, 2018 at 10:24
  • 1
    $\begingroup$ 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. $\endgroup$
    – arkaia
    Commented Jan 5, 2018 at 16:30
  • 1
    $\begingroup$ A package that does an okay job with this and it is relatively efficient is: dfo-mpo.gc.ca/Library/326156.pdf $\endgroup$
    – arkaia
    Commented Jan 5, 2018 at 16:36
  • 1
    $\begingroup$ @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? $\endgroup$
    – Semidiurnal Simon
    Commented Jan 7, 2018 at 10:24
  • 2
    $\begingroup$ 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). $\endgroup$
    – arkaia
    Commented Jan 7, 2018 at 13:36

1 Answer 1

Reset to default
2
$\begingroup$

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).

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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! $\endgroup$
    – Semidiurnal Simon
    Commented Jan 7, 2018 at 10:41
  • $\begingroup$ @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. $\endgroup$
    – Camilo Rada
    Commented Jan 7, 2018 at 15:22
  • $\begingroup$ 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 ;-) $\endgroup$
    – Semidiurnal Simon
    Commented Jan 9, 2018 at 7:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.