Expand a given bathymetry / DEM through interpolation towards 0 depth at coastline polygon

Question

I have a grid of a bathymetry of a fjord (grayscale in the image below). The grid doesn't reach all the way to the coastline (blue polygon) because it originates from sonar surveys where shallow regions could not be reached, obviously.

Now I would like to expand the grid outwards to the polygon border by interpolating between the outermost depths of of the raster layer and 0 depth at the coastline.

It is essentially the same as in this question from 2018 by @Jordi, which did not receive an answer unfortunately.

What is an efficient way to do this, preferably in QGIS (or any other open source tool).

Answer 1

One approach would be to create a raster grid that includes the zero values at pixels that fall on the polygon border (and a bit beyond), and then use one of the interpolation functions in QGIS to interpolate the null values (the blue areas).

There are a few ways that you could build out the raster. Here is one example, using one letter codes for the intermediate results (all non provided parameters left at default values):

0. starting with the incomplete DEM raster as A and the (multi)polygon vector layer surrounding the complete region as B

1. Vector > Geoprocessing Tools > Buffer --> output = C

   • Input layer: B
   • Distance: e.g. 10 x the grid resolution of A

   Dummy example:

   

2. from Processing Toolbox Create constant raster layer --> output = D

   • Extend: Calculate from layer: C
   • Pixel size: grid resolution of A
   • Target CRS: same as in A
   • Constant value: 0

3. Vector > Geoprocessing Tools > Difference --> output = E

   • Input layer: C
   • Overlay layer: B

4. Raster > Extraction > Clip by mask --> output = F

   • Input layer: D
   • Mask: E
   • NoData value: same as in A

5. Raster > Miscellaneous > Merge --> output = G

   • Input layers: A, F
   • Input NoData value: same as in A
   • Output NoData value: same as in A

At this stage you have the a raster layer, which is the original A raster plus a "ring" of 0-value pixels in the area of the buffer. In between the original DEM pixels and the 0-buffer there are NoData-pixels, as well as outside the buffer. Now interpolate the (inner) missing values (eg. GDAL > Fill nodata, or GRASS > r.fillnulls or WhiteBoxTools > FillMissingData (WBT needs to be installed separately). You may need to tweak some details depending on the chosen algorithm and on whether your coastline polygon surrounds the original grid or not.

6. "Interpolation method of choice" --> output = H
   • Input: G
   • Algorithm specific parameters influence the interpolated values

The various interpolation functions in QGIS don't generally perform that well when required to interpolate large areas. Below is the dummy example above, filled using r.fillnulls with bicubic interpolation. It's possible that playing with the parameters may yield better results.

7. Raster > Extraction > Clip by mask --> output = I

   • Input layer: H
   • Mask: B
   • NoData value: same as in A

8. right-click on I in Layers panel Export > Save As... to save to your preferred raster format.

Answer 2

One other suitable open source tool is GMT. The "surface" algorithm makes pretty nice results and it supports so called soft breaklines https://docs.generic-mapping-tools.org/latest/surface.html.

Use x, y, z data in the breakline file as a ‘soft breakline’. A ‘soft breakline’ is a line whose vertices will be used to constrain the nearest grid nodes without any further interpolation. A coastline or a lake shore are good examples of ‘soft breaklines’. Multi-segments files are accepted. If your lines do not have z-values or you wish to override those with a constant z-value, then append +zlevel to the filename. If no value is given then we default to 0.

As x, y, z data for example POLYGONZ type shapefiles which have 0 as z value can be used. The hardest thing with GMT is to learn the command line syntax.

The GRASS function v.surf.rst is quite similar to the GMT surface https://grass.osgeo.org/grass82/manuals/v.surf.rst.html. This function needs the shoreline as a raster mask:

User can either use r.mask to set a mask or specify a raster map in mask option, which will be used as a mask. The approximation is skipped for cells which have zero or NULL value in mask. NULL values will be assigned to these cells in all output raster maps. Data points are checked for identical points and points that are closer to each other than the given dmin are removed. If sparsely digitized contours or isolines are used as input, additional points are computed between each 2 points on a line if the distance between them is greater than specified dmax. Parameter zmult allows user to rescale the values used for approximation (useful e.g. for transformation of elevations given in feet to meters, so that the proper values of slopes and curvatures can be computed).

You can also try to extract the vertices of the shoreline as points with 0 depth and append them into your sonar data. Then you can use your current interpolation. The drawback compared to breakline or mask is that the interpolation algorithm does not know that the shoreline points are somehow special and therefore the algorithms tend to overshoots the interpolation over the dry land, and zero depths are not necessarily kept as zeroes.