Combining SRTM height, water body and ETOPO bathymetry data

Data source

We have the three different sources of data :

An SRTM HGT height file covers 1 x 1 area with height values represented as regular grid of 1200x1200 intervals (1201x1201 points). Such file does not contain data about sea depth, as it was generated by the Shuttle radar measuring heights from land or sea surface. To be able to distinguish between land or water space, SWBD data is available - shore contours.

Before processing, SRTM and SWBD data are corrected :

Our purpose is to combine all the data into a regular HGT file grid..

Solution

Incorporation of depth is straightforward : once we get shore contours, we simply replace HGT height values by depth values everywhere inside water areas. At points where low-resolution ETOPO data is grossly inaccurate (e.g. provides positive land height where it should be negative), a default water depth is taken.

Original height data may substantially differ from much more accurate shore contours. Now we need to bring HGT data close to shore contour data.

Suppose we have a regular grid of 1200x1200 intervals (1201x1201 points) for geographic area of size 1 by 1. At each point we have a height value. So, a cell resolution for this file is, roughly speaking, 100 by 100 metres.

Let use represent height function by piecewise linear hat basis in the form

H(u,v) = Σi,j=0,1200Pi,j Ri,j(u,v)

where

Suppose we have a number of constraint points (they are water countour points) which should be moved to level H = 0. Call them Sf, f=1,...,h; they have parameter values (uf,vf), f=1,...h. These points will be moved to target points Tf, f=1,...,h with zero height.

In order to move Sf to Tf, f=1,...,h, we need to compute displacement δi,j of each control point Pi,j so that the following constraints are satisfied :

Tf = Σi,j=0,1200(Pi,j + δi,j) Ri,j(uf,vf) = Sf + Σi,j=0,1200δi,j Ri,j(uf,vf)

f = 1,...,h

Let us minimize displacements with constraints introduced with Lagrange multipliers

L = Σi,j=0,1200i,j|2 + Σf=1..hλf(Tf - Sf - Σi,j=0,1200δi,j Ri,j(uf,vf))

where λf, f = 1,...,h are Lagrange multipliers. After differentiating this function by displacements δ and Lagrange multipliers λ and equating to zero we get a system of two equations

Tf = Sf + Σi,j=0,1200δi,j Ri,j(uf,vf)  f = 1,...,h

δi,j = 0.5 Σf=1,h λf Ri,j(uf,vf)  0≤i≤1200, 0≤j≤1200.

In the matrix form

D=RTδ    (1)
δ=Rλ    (2)

By substituting equation (2) into equation (1)

D=RTRλ

from where we can get Lagrange multipliers by solving a system of linear equations. Further, we get the required grid node displacements from equation (2).

Nodes of HGT file grid are displaced by these values to make cell intersections with plane Z = 0 close to specified shore countour lines.

References

1. Direct manipulation of FFD : efficient explicit solutions and decomposible multiple point constraints. Shi-Min Hu, Hui Zhang, Chiew-Lan Tai, Jia-Guang Sun. Springer-Verlag 2001.