How to grid bathymetry
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
Coastline information is also available for many tidal estuaries and
These data are typically expressed in terms of eastings and northings.
There are many potential uses for bathymetric data.
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:
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.
- 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.
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.
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.
- 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
- The river may bend between measured cross-sections.
- Orientation of cross-sections varies between measured cross-sections
and adjacent measured cross-sections can have
totally different orientations.
- Channel width varies between measured cross-sections.
- 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.
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.
- The grid spacing is no longer regular,
but the local y-direction is defined as being orthogonal to the
local along-channel direction.
- Spacing of points along the thalweg
is at an interval fixed to capture the sparse bathymetry and
also the bends in the river. (Information about bends in the
river often is not present in the bathymetry --- rather it
is obtained from coastlines or topography.)
- Spacing of points along channel cross-sections is determined
by the local channel width and the number of points used
to span that width.
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
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
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.)
- 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.
- 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.
- 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.
- Ultimately, if high-resolution measurement methods are ever used
to measure bathymetry, then it would be most effective to
measurements directly into the bathymetry/topography on a
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
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.
- 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.
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:
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
- Black dots show the DEM.
- Cyan dots with connecting line show the coastline data.
- Blue dots with connecting line show the bathymetric measurements.
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
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.
Having constructed a satisfactory thalweg,
save xt and yt in ISAVE/xy_thalweg#.mat --- where # corresponds to isave#.mat.
- 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.
Next we use
loads topograhy and bathymetry and cycles through selected
ISAVE/xy_thalweg#.mat and ISAVE/isave#.mat to:
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
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.
- 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
- 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.
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.
- Topography on a regular grid at 25 m resolution (locations
in water being undefined).
- Bathymetric data as x,y,z information at sparse and irregular
- Coastline contours x,y,z=0
- Filled-in points which are high-density pseudo-bathymetric data
at irregular locations.
- Optionally, we may also have x,y,z information on a strained
grid that follows the thalweg of a creek.
We have created a script
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
Output from call_griddata.m is still likely to contain NaN
values at points where there are no measurements to reasonably
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
- hypsometry.m calculates hypsometry from
gridded bathymetry (either on a uniform or strained grid).
- z2WH.m to convert transects
into channel depth H and channel width W as
tabulated functions of water level zL for use by a
one-dimensional hydrodynamic estuary model.
- HW2hypsometry.m converts H, W lookup
tables to hypsometry --- this conversion can never be exact because
information is lost going from depth-transects to H, W tables.
Nevertheless, it's handy --- and uncertainty in the bathymetry
is often the bigger issue.