POCS_SURFACE(9) POCS_SURFACE(9)1 July 1992

NAMEpocs_surface - gridding algorithm using POCS inversion method

SYNOPSISpocs_surface[xyzfile] [-Goutputfile.grd] [-Ix_inc[/y_inc]] [-Rxmin/xmax/ymin/ymax] [-V] [-D] [-Nmax_iterations] [-L] [-Cconvergence_limit] [-Ttolerance_factor] [-Emodel_spectrum_expression]

DESCRIPTIONpocs_surface reads randomly-spaced (x,y,z) triples from standard input [orxyzfile] and produces a binary netCDF .grdformat 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

ARGUMENTSxyzfileascii 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 defaultexpression is '1/(1+((0.2*kr/kxt)^2))', which is flatfor 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).

-DOmit initialization using surface & use only a random noise distribution

- 2 - Formatted: March 4, 1999

POCS_SURFACE(9) POCS_SURFACE(9)1 July 1992

-VRun in verbose mode. Prints out error for each iteration.

-LEnforce 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 calleddata4_grd.grd, try:

pocs_surface data4 -GPdata4.grd -I0.1 -R0/5/0/5 -V

BUGSIf we knew where the bugs were, don't you think we would fix them?

SEE ALSOgmtsystem, surface, pocs_helper

REFERENCESMenke, 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