Combining SRTM height, water body and ETOPO bathymetry data
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..
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 multipliersL = Σi,j=0,1200|δi,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
By substituting equation (2) into equation (1)
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.
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.