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 Riversystem 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 crosssectional length scales that are not wellresolved
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 finescale 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 crosssection, depth might be measured
every 1 to 5 metres. The distance between crosssections can
range from hundreds of metres to kilometres. Any algorithm
to interpolate depths between measured crosssections runs into
several problems:
 The river may bend between measured crosssections.
 Orientation of crosssections varies between measured crosssections
and adjacent measured crosssections can have
totally different orientations.
 Channel width varies between measured crosssections.
 Patterns of sediment deposition and erosion in a meandering
river result in a loss of coherence between measured crosssections.
 Generally, bathymetry varies more rapidly in the crosschannel
sense than in the alongchannel sense. Bathymetry in channels
is fundamentally inhomogeneous.
 Creeks represent quasi onedimensional geometric objects that
occupy a small area within a large topographic domain.
In order to represent finescale 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 xy grid
that is strained so that the xaxis follows the alongchannel
direction of the creek.
 The grid spacing is no longer regular,
but the local ydirection is defined as being orthogonal to the
local alongchannel 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 crosssections is determined
by the local channel width and the number of points used
to span that width.
To represent the bathymetry of a meandering creek on a strained
grid we require three matrices: one for depths, one for xcoordinates,
and one for ycoordinates. 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
lowresolution bathymetry may be represented using a DEM on a regular grid.
High resolution information can be overlayed 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 strainedgrid
if appropriate data are available. They can also be wellrepresented
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 strainedgrid, discussed above, can be considered as a
special mesh of triangles that have some constraints imposed
upon them. The constraints make the strainedgrid 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 multiscale regulargrid
If we move our perspective from representing bathymetry/topography
at a gridpoint 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 higherresolution subsquares 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 subsquares and each subsquare would have its own average
altitude/depth.
The process can be made recursive so that those subsquares where
bathymetry is very complex point to a matrix of subsubsquares
that locally represent bathymetry on a length scale L/N^{2}.
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 multigrid methods will immediately see that
this representation of data is a slight generalization
of datastorage techniques used by the code for the Full MultiGrid 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 objectoriented attack on the problem, be careful to
keep such people on a very tight leash.
The multiscale regular grid enables the best features of the
regular grid commonly used for DEM  but without its major
liability.
Recommendations
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 (fieldworkers, 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 strainedgrid is
a good way of representing bathymetry.
 Regions in which the waterbody 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 multiscale regulargrid.
There is probably no need to do this immediately,
but it is worth thinking about for the longer term.
 Ultimately, if highresolution 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
multiscale regulargrid.
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 sandbars in the typically shallow Australian estuaries.
Use googleearth to qualitatively preplan 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 bottomdepth varying rapidly.
Going along the estuary thalweg we would see bottomdepth varying slowly.
Geology has asymmetries that should be exploited when making bathymetric
measurements.
 Bathymetry should be measured with high resolution along a crosssection.
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 crosssections can be much larger than than
the spacing of bathymetric measurements along a crosssection.
In general, distance between crosssections 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 crosssections
marked in red. The additional green crosssections are also advised
in order to resolve the nearbend 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 sandbars 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 fieldworkers 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.



Longnarrow islands can sometimes be wellresolved using
acrosschannel 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 wellresolved 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, nondirectional).


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
zigzag 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 alongchannel 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 highdensity pseudodata 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 fillin 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
onedimensional 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).
Alongchannel resolution of 25 m resolves the channel bends reasonably well
and the crosschannel 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 psuedodata 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:
 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
locations.
 Coastline contours x,y,z=0
 Filledin points which are highdensity pseudobathymetric data
at irregular locations.
 Optionally, we may also have x,y,z information on a strained
grid that follows the thalweg of a creek.
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 wellbehaved 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
 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
onedimensional hydrodynamic estuary model.
 HW2hypsometry.m converts H, W lookup
tables to hypsometry  this conversion can never be exact because
information is lost going from depthtransects to H, W tables.
Nevertheless, it's handy  and uncertainty in the bathymetry
is often the bigger issue.