How to grid bathymetry

Brian Sanderson

What sort of grid should we use?

Bathymetric data is available along irregular transects in some of New South Wales tidal estuaries and intermitently closing and opening lagoons (ICOL's). In addition, a digital elevation model (DEM) has been constructed to describe topography over much of New South Wales. The DEM has a spatial resolution of 25 m and is on a regular grid (uniform spacing). Coastline information is also available for many tidal estuaries and ICOL's. These data are typically expressed in terms of eastings and northings.

There are many potential uses for bathymetric data.
Logically, the bathymetric data might be used to fill in that part of the DEM that is below sealevel. The accompaning figure represents each bathymetric sounding in Hastings River-system as a blue dot. Coastlines are shown in cyan. The thick black line demarks an area within which bathymetry might be sensibly interpolated onto the DEM grid --- in order to essentially extend the topographic DEM below waterlevel. Even in this area, however, there are locations where:
  • measurements are sparse so interpolation onto a 25 m grid requires an algorithm that adds knowledge in addition to the information content of the bathymetric measurements themselves.
  • some channels have cross-sectional length scales that are not well-resolved on a uniform 25 m resolution grid.
My approach is, indeed, to represent the bathymetric information on a uniform grid in such areas. In this particular example, most of the important features can be represented on a uniform grid at a horizontal resolution of 12.5 m.

One can't, however, directly interpolate onto such a grid because for many locations the bathymetric data is very sparse when we consider such a fine-scale grid. Let us put the discussion of algorithms for gridding to one side while we consider other properties of the bathymetric data sets.

The figure to the left shows another portion of the system of rivers that make up Hastings River Estuary.
  • Mostly, bathymetry data consists of measurements across sections of rivers. Within a cross-section, depth might be measured every 1 to 5 metres. The distance between cross-sections can range from hundreds of metres to kilometres. Any algorithm to interpolate depths between measured cross-sections runs into several problems:
    1. The river may bend between measured cross-sections.
    2. Orientation of cross-sections varies between measured cross-sections and adjacent measured cross-sections can have totally different orientations.
    3. Channel width varies between measured cross-sections.
    4. Patterns of sediment deposition and erosion in a meandering river result in a loss of coherence between measured cross-sections.
  • Generally, bathymetry varies more rapidly in the cross-channel sense than in the along-channel sense. Bathymetry in channels is fundamentally inhomogeneous.
  • Creeks represent quasi one-dimensional geometric objects that occupy a small area within a large topographic domain.
In order to represent fine-scale bathymetric information on a regular grid, the grid would have to be of very high spatial resolution. At most locations, the high spatial resolution would be redundant. Yet the high resolution is an intrinsic requirement to represent the very features of the bathymetry that are of interest.

The regular grid

A DEM can be represented as a matrix with some uniform spacing between points. Thus, to specify the DEM we need only have a matrix of elevations, a corner coordinate, and the grid spacing. The positions of points on the regular grid can be derived from the corner coordinate and grid spacing.

The primary difficulty of a regular grid is that it can only record bathymetry (or topography) up to some resolution defined by grid spacing. Some areas may be adequately resolved with widely spaced grids ---- others may require closely spaced grids. In order to resolve high resolution features, the resolution must be high throughout the regular grid. The number of data points increases as the square of grid resolution which imposes practical limitations upon this approach.

The strained grid

One way of representing bathymetry of a narrow meandering creek is to define a coordinate system that follows the thalweg. Thus, depths become defined upon an x-y grid that is strained so that the x-axis follows the along-channel direction of the creek. To represent the bathymetry of a meandering creek on a strained grid we require three matrices: one for depths, one for x-coordinates, and one for y-coordinates. This scales as three times the data storage required for a regular grid of the same resolution --- but is far more efficient when high resolution information is only required on part of the domain of a DEM.

Thus, topography and areas with intrinsically low-resolution bathymetry may be represented using a DEM on a regular grid. High resolution information can be over-layed for meandering narrow creeks which would be represented on a strained coordinate system.

Existing bathymetric data are often inadequate at channel junctions (this should be very obvious from the plot above). At some later time I will write specifications for measuring channel junctions and methodologies for analysing such data. Channel junctions can be well represented on a strained-grid if appropriate data are available. They can also be well-represented on a regular grid that spans a local area with high resolution.

A polygon mesh

Coordinates can be defined as vertices on a mesh made up of polygons. In a regular grid the polygons are all squares of equal sides. In general, however, any sort of polygon can be used and the polygons do not need to have uniform size or shape. For example, it is quite possible to divide a complex area (such as a meandering river) into a mesh of triangles with different sizes and shapes. Small triangles could be used where the bathymetry has intrinsic small spatial scales and large triangles where the bathymetry varies slowly with respect to position.

The strained-grid, discussed above, can be considered as a special mesh of triangles that have some constraints imposed upon them. The constraints make the strained-grid coordinate system algorithmically easy to relate to the way bathymetric measurements have been made.

More general triangular (or polygon) grid can be valuable for some types of modelling. Once bathymetry represented on the strained grid above it is not so difficult to represent it on any other grid.

The multi-scale regular-grid

If we move our perspective from representing bathymetry/topography at a grid-point to representing bathymetry/topography averaged over a rectangular area then we obtian a much more powerful approach --- which I believe to be optimal. Thus, some domain can be divided into squares with sides of length L. A matrix can be used to represent average altitude/depth in each square. So far, we have something that is not that much different from a DEM.

If a bathymetry/topography varies little within a square then there is not much more to say about that square. On the other hand, bathymetry/topography might vary greatly within some square (or a bunch of squares). For those squares where bathymetry/topography varies greatly, we define higher-resolution sub-squares of length L/N where N is some small integer (eg. 2 or 3). The lower resolution square would contain a pointer to a matrix of sub-squares and each sub-square would have its own average altitude/depth.

The process can be made recursive so that those sub-squares where bathymetry is very complex point to a matrix of sub-sub-squares that locally represent bathymetry on a length scale L/N2.

If L=100 and N=3 then after four levels of recursion one achieves a horizontal resolution of 1.2 m. Those who have solved elliptic equations using multi-grid methods will immediately see that this representation of data is a slight generalization of data-storage techniques used by the code for the Full Multi-Grid Solver found in Numerical Recipes in Fortran 90: The art of parallel scientific computing, W.H. Press et al (Cambridge University Press 1996). Advanced Matlab users may see how "Structures" might be put to good use for such data storage. Generally, any respectable computing language that supports pointers, arrays, and derived data types will serve --- and many of us eat this sort of stuff for breakfast. C++ programmers might digress with an obscure object-oriented attack on the problem, be careful to keep such people on a very tight leash.

The multi-scale regular grid enables the best features of the regular grid commonly used for DEM --- but without its major liability.


The bathymetric data are often gappy. Generally, however, bathymetric measurements have been made in a manner that captures the most essential information at highest resolution (field-workers, take a bow). It is reasonable to assume that bathymetry varies more rapidly across rivers than along rivers --- and this is how measurements are made. (There are many other implicit assumptions that I won't get into, but which those who look at bathymetric data will quickly figure out for themselves.)
  1. A narrow meandering river occupies a small amount of area within a large spatial domain. In such circumstances, the strained-grid is a good way of representing bathymetry.
  2. Regions in which the water-body fills much of the area, should also be represented on a regular grid --- ie to fill in the water part of a DEM. Note, however, if a DEM has a spatial resolution of 50 m then it simply can't represent the bathymetry of a 10 m wide creek that meanders through it.
  3. Once both of the above are completed, one could use the resulting data base to calculate bathymetry on a multi-scale regular-grid. There is probably no need to do this immediately, but it is worth thinking about for the longer term.
  4. Ultimately, if high-resolution measurement methods are ever used to measure bathymetry, then it would be most effective to incorporate such measurements directly into the bathymetry/topography on a multi-scale regular-grid.

It is not difficult to use 1 and 2 above (or 3) to define topography/bathyemtry on any arbitrary grid (perhaps made from polygon vertices) as may be required for some sorts of modelling.

How to measure bathymetry in narrow estuaries

First and foremost, googleearth can often give you a rough picture of the channels and sand-bars in the typically shallow Australian estuaries. Use googleearth to qualitatively pre-plan your survey.
The figure to the right shows a portion of the Hastings River estuary. This portion is characterized as being long, narrow, and winding. Going across the estuary we would see bottom-depth varying rapidly. Going along the estuary thalweg we would see bottom-depth varying slowly. Geology has asymmetries that should be exploited when making bathymetric measurements.
  • Bathymetry should be measured with high resolution along a cross-section. Obviously, if the channel is only 10 metres wide the resolution should be 1 metre or better --- whereas a 500 metre wide channel might be well resolved by depth soundings at 5 metre intervals.
  • Distance between cross-sections can be much larger than than the spacing of bathymetric measurements along a cross-section. In general, distance between cross-sections should be adequate to resolve bends in the river, islands, changes in channel width, and branching channels.
Blue sections represent measurements made in the Hastings River. Clearly the river bend is poorly resolved and it would be advisable to measure depth along additional cross-sections marked in red. The additional green cross-sections are also advised in order to resolve the near-bend bathymetry at the scale of the bend.
The figure to the left shows a portion of the Hastings River estuary. The blue transects (actual measurements) resolve the major bend in the estuary but do not resolve features in the channel boundaries. Additional sections shown in red might be contemplated. Such detail may be important since changes in the coastline can indicate sand-bars and abrupt changes in the channel configuration. Googleearth can help you decide the importance of such coastline irregularities. Often they are of no major significance.
Branching channels can represent a major headache for the person charged with turning raw measurements into gridded bahtymetry. The blue transects (actual measurements) obviously do not resolve the bathymetry at the channel junction plotted to the right. Additional transects, plotted in red, are required. Different field-workers do all sorts of weird and wonderful things when confronted with a junction. The red transects at the junction should always be done. Googleearth may indicate channel structure requiring additional transects --- that's OK, but still do the basics. If the field worker feels compelled to do something weird and wonderful --- go ahead, but still do the basics.
Long-narrow islands can sometimes be well-resolved using across-channel transects. Ideally, these transects should be contiguous right across the island. At the very least, transects on one side of the channel should line up with those on the other side. Also, transects should be made immediately upstream and downstream of the island. Required measurements are shown on the figure to the left.
More chunky islands in lagoons, like that shown on the right can be well-resolved by measuring bathymetry along concentric paths that circle the island --- as shown to the right in red. Concentric circles are NOT a good way of doing things in narrow channels --- even if the island is chunky. They may, however, be useful in other areas (such as coastal lagoons) where the geometry of the bathymetry is more isotropic (ie, non-directional).

Algorithms for gridding bathymetric measurements

This section is under construction

Loading the data

Suppose we have measurements of bathymetry in a channel. Often the measurements will be made along a ships track --- which may zig-zag across the channel. Data are most naturally stored in a file with each line having an x,y,z coordinate. (The raw data would have had a time stamp on each line which we would have used to correct z for tides.) Often a digital elevation model of topography (DEM) and coastlines are available. The plot on the right shows positions of measurements:
  • Black dots show the DEM.
  • Cyan dots with connecting line show the coastline data.
  • Blue dots with connecting line show the bathymetric measurements.
The script ld_topog_bathy.m provides a template for loading up the data sets and plotting them in a fashion that is convenient to work with. The figure on the right was created by zooming in on the plot created by ld_topog_bathy.m

Typically bathymetric measurements have higher resolution along the ships track (the figure to the right is anomolous in this regard) and low resolution along tracks. We know that the bathymetry is less variable along-channel than across channel so a human can fill in the gaps with artificial data points as follows.

Selecting transects and filling in with artificial data --- simple case

The script line_sections.m can be used to select indices that define the end points of transects. The indices are saved in a matrix called isave. Many isave matrices are usually generated for a large and complex data set so they are saved in files ISAVE/isave#.mat where # is a number.

It it easy to get confused about what has and has not been done so we have created a script plt_isave.m which plots out the transects as shown in the figure. Note, the first point is plotted in magenta. First points must always be on the same river bank.

The Matlab script fill_in.m reads in the isave matrices and uses them define transects. It then interolated between neighbouring transects. A sneaky parametric interpolation is used which essentally folds neighbouring transects together. Thus we end up with a high density of points as shown in red on the figure to the right. This high-density pseudo-data is save in ISAVE/xyz_isave#.mat for later use.

Note, there is no "best" method for interpolation between transects. For the example shown, fill_in.m is appropriate. If we had transects in concentric circles about an island (or a spiral transect about an island) then fill_in.m would also work well providing values in isave were sensibly selected. To understand fill_in.m, you should read the commented code.

Later we shall that other interpolation methods are more appropriate for constructing fill-in points under other circumstances.

Selecting transects and filling in with artificial data --- curved channel

Consider a circumstance where stransects are widely separated on a narrow curving channel that has relatively uniform width. Separation is large compared to the channel width. Neighbouring transects could have quite different orientations and the straight line between transects might go over land. The simple parametric interpolation used above would clearly fail. Instead, we must define a thalweg (line along the centre of the channel) and interpolate along the thalweg.

The basic approach is to select transects and save the indices for transect end points in ISAVE/isave#.mat as above. Next one can use the function guide_thalweg.m to automatically calculate the vectors xt, yt containing thalweg coordinates for a given isave. The automated calculation is not always adequate --- see the red line in the plot opposite.

  • If only a few points need to be adjusted, then this can be done using the function correct_thalweg.m
  • If many points need to be adjusted it may be best to use mk_xy_thalweg.m to manually construct new thalweg vectors. In this case, the automatically generated thalweg can still act as a valuable guide.
Having constructed a satisfactory thalweg, save xt and yt in ISAVE/xy_thalweg#.mat --- where # corresponds to isave#.mat.
Next we use mk_cross_sections.m which loads topograhy and bathymetry and cycles through selected ISAVE/xy_thalweg#.mat and ISAVE/isave#.mat to:
  • Interpolates thalweg coordinates xt, yt to every metre along the thalweg. This creates a rather long vector which is used in a brute force manner by mk_strained_grid.m. (The next generation of code will be more elegant.)
  • Calls the function cross_section.m to create arrays xr, yr, and depth where each row gives coordinates and depths along a transect. Positions along a transect have uniform spacing. (The spacing varies from transect to transect because different transects have different lengths and each transect is divided into M=50 equal intervals.)
  • Saves xr, yr, depth, xt, yt into transects#.mat for later use by mk_strained_grid.m.
  • Uses the function z2WH.m to convert each transect into functions of channel depth H and channel width W as tabulated functions of water level zL=[2.5:-.1:-2.6]
  • Stores lookup tables of H and W at 50 metre intervals along the thalweg. These are used by the one-dimensional hydrodynamic estuary model.
The results saved in transects#.mat contain xr, yr positions of transects shown in magenta on the figure to the right. These transects are interpolated onto a strained grid by the script mk_strained_grid.m. The strained grid follows the thalweg and rotates the coordinate system as it does so using calls to the function rot_interp.m. In this manner we obtain bathymetry at the cyan transects (here every 25 m along the thalweg --- the thalweg is plotted with a dotted line). Along-channel resolution of 25 m resolves the channel bends reasonably well and the cross-channel direction is resolved at much higher resolution (50 points span the width of the creek). Obviously, the strained coordinate system only makes sense if the banks of the creek are smoothly varying between measured transects. The figure shows a 5 m contour (dark blue) which follows the bend. Clearly, we would want to interpolate to more closely spaced transects if we were to use psuedo-data on the strained grid as input for the calculation of depth at say 5 metre resolution on a uniform grid.

Interpolating onto a regular grid

After going through all of the above processes we will have: In principle, the Matlab function gridddata.m can interpolate irregular data onto a uniform grid. In practice, this is only possible if the data are mathematically well-behaved and not too numerous.

We have created a script call_griddata.m which selects data relevant to a point on the target grid and them makes a call to griddata.m in order to obtain an interpolated value. Note, call_griddata.m makes a call to remdup.m which removes duplicate points. Comments in call_griddata.m explain how you should save the gridded topography+bathymetry. Output from call_griddata.m is still likely to contain NaN values at points where there are no measurements to reasonably interpolate from.

The script removeNaN.m uses linear interpolations to remove the NaN values from the gridded topography+bathymetry. After you run this script, the final answer will be in the variable z1. It is recommended that z1 never be stored --- because it contains dishonest values (obtained by gross interpolation).

Finally, if new data come to hand then they could be patched into the existing gridded data using something like post_fix.m.

Using bathymetry