POCS_SURFACE(9) POCS_SURFACE(9)
1 July 1992
NAME
pocs_surface - gridding algorithm using POCS inversion method
SYNOPSIS
pocs_surface [xyzfile] [-Goutputfile.grd] [-Ix_inc[/y_inc]]
[-Rxmin/xmax/ymin/ymax] [-V] [-D] [-Nmax_iterations] [-L]
[-Cconvergence_limit] [-Ttolerance_factor]
[-Emodel_spectrum_expression]
DESCRIPTION
pocs_surface reads randomly-spaced (x,y,z) triples from standard input
[or xyzfile] and produces a binary netCDF .grd format file of gridded
values z(x,y) by applying the method of Projection Onto Convex Sets
(POCS). Smoothness of the solution is determined by a model amplitude
spectrum, s(kx,ky), where (kx, ky) are wavenumbers. The model spectrum
is specified by an expression [-Emodel_spectrum_expression] or by the
default expression. It is recommended that the user check the radially
averaged spectrum of the xyzfile with pocs_helper in order to choose
an expression that is compatible with the data. Number of iterations
needed for convergence will vary according to amount and distribution
of data points and desired smoothness. The output is to some degree
dependent on the initial guess used to start the POCS iteration.
pocs_surface provides two possible initial guesses: (by default) A
smooth interpolation provided by the GMT 'surface' gridding program,
perturbed with 10% random noise; or (using the -D option) 100% random
noise. The output satisfies the data and has a spectrum that is less
than (1+tolerance)*s(kx,ky). With the [-L] option, the user can force
the spectrum to be greater than (1-tolerance)*s(kx,ky). However, this
constraint is not convex, so that convergence is not guaranteed.
pocs_surface is an enhancement of the LDEO-HIG-SIO GMTSYSTEM software
package, and uses that package's 'surface' gridding routine. The
Environmental variable GMT_HOME and ARCH must be set so that
$GMT_HOME/$ARCH/surface is the path of the 'surface' gridding program.
pocs_surface uses the Lamont Expression Parser 'parse' libaray and
UCAR's Network Common Data Format
ARGUMENTS
xyzfile
ascii input data file of 'x y z' values, white-space separated,
one set per line. If no file is specified, surface will read from
standard input.
-G -Ggrdfile output grd file name
-I -Idx/dy the grid increment. Default is dx=dy=0.1. If only dx is
specified, the dy=dx. Warning: dx and dy are decreased, if
necessary, to make the number of rows and columns in the image
integral powers of 2.
- 1 - Formatted: March 4, 1999
POCS_SURFACE(9) POCS_SURFACE(9)
1 July 1992
-R -Rxmin/xmax/ymin/ymax specifies the min/max coordinates of the
grid in user units. Default is 0/1/0/1.
-C -Cconvergence_limit. POCS iteration stops when the error is less
than the convergence limit. The error is defined as the mean
square (data_predicted-data) divided by the mean square data.
Default is 0.0001.
-T -Ttolerance_factor. The spectrum of the image is forced below
(1+tolerance)* model_spectrum, where the tolerance decreases with
the radial wavenumber, kr, by the formula tolerance =
tolerance_factor / (1+(0.2*kr/kxt)^2). Here kxt is an estimate of
the maximum wavenumber controlled by the data, that is, 1/(2*dx),
where dx is a typical spacing between points. Default is
tolerance_factor = 1.0.
-N -Niterations is the maximum number of iterations. Default = 20.
-E -Emodel_spectrum_expression determines value of model spectrum.
It is specified in the C-like syntax of the Lamont Expression
Parser, and may contain the following symbols:
kx: horizontal wavenumber
kxt: typical horizontal wavenumber of data, as above
kxny: nyquist horizontal wavenumber
ky: vertical wavenumber
kyt: typical vertical wavenumber of data
kyny: nyquist horizontal wavenumber
kr: radial wavenumber, i.e. sqrt(kx^2+ky^2)
theta: angular wavenumber, radians, i.e. atan2(ky,kx) The default
expression is '1/(1+((0.2*kr/kxt)^2))', which is flat for krkx_t, leading to an
image rich in small scale features.
Note 1: Quote the expression, so the shell doesn't try to parse
it!
Note 2: The overall amplitude of the spectrum is irrelevant. It
is adjusted to match the amplitude of the input data (using
Parseval's equality).
-D Omit initialization using surface & use only a random noise
distribution
- 2 - Formatted: March 4, 1999
POCS_SURFACE(9) POCS_SURFACE(9)
1 July 1992
-V Run in verbose mode. Prints out error for each iteration.
-L Enforce lower bound of model spectrum.
To make a 5 by 5 grid from the data in data4 with options
-Niterations=100 and -Vverbose, writing the result to a file called
data4_grd.grd, try:
pocs_surface data4 -GPdata4.grd -I0.1 -R0/5/0/5 -V
BUGS
If we knew where the bugs were, don't you think we would fix them?
SEE ALSO
gmtsystem, surface, pocs_helper
REFERENCES
Menke, W., Applications of the POCS inversion method to interpolating
topography and other geophysical fields, Geophys. Res. Lett. 18, 435-
438, 1991.
Gubin, L.G., B.T. Polyakm and E. Raik, The method of projections for
the common point of convex sets, USSR Comput. Math. and Phys. 7, 1-24,
1967.
Wessel, P., and W. H. F. Smith, 1992, The GMT-SYSTEM v. 2.1 Technical
Reference & Cookbook, SOEST/SIO.
Wessel, P., and W. H. F. Smith, 1991, Free software helps map and
display data, EOS Trans. AGU, 72, 441.
Smith, W. H. F, and P. Wessel, 1990, Gridding with continuous
curvature splines in tension, Geophysics, 55, 293-305.
- 3 - Formatted: March 4, 1999