This is the command greensplinegmt that can be run in the OnWorks free hosting provider using one of our multiple free online workstations such as Ubuntu Online, Fedora Online, Windows online emulator or MAC OS online emulator

**PROGRAM:**

**NAME**

greenspline - Interpolate using Green's functions for splines in 1-3 dimensions

**SYNOPSIS**

**greenspline**[

__table__] [ [

**1**|

**2**|

**3**|

**4**|

**5**,]

__gradfile__] [ [

**n**|

**v**]

__cut__[/

__file__] ] [

__mode__] [

__grdfile__] [

__xinc__[/

__yinc__[/

__zinc__]] ] [ ] [

__nodefile__] [

__az__|

__x/y/z__] [

__west__/

__east__/

__south__/

__north__[/

__zmin__/

__zmax__][

**r**]

] [

**c|t|l|r|p|q**[

__pars__] ] [

__maskgrid__] [ [

__level__] ] [ ] [

**-b**<binary> ] [

**-d**<nodata> ] [

**-f**<flags> ] [

**-h**<headers> ] [

**-o**<flags> ] [

**-x**[[-]

__n__] ] [

**-:**[

**i**|

**o**] ]

**Note:**No space is allowed between the option flag and the associated arguments.

**DESCRIPTION**

**greenspline**uses the Green's function G(

**x**;

**x'**) for the chosen spline and geometry to

interpolate data at regular [or arbitrary] output locations. Mathematically, the solution

is composed as

__w__(

**x**) = sum {

__c__(

__i__) G(

**x'**;

**x**(

__i__))}, for

__i__= 1,

__n__, the number of data points

{

**x**(

__i__),

__w__(

__i__)}. Once the

__n__coefficients

__c__(

__i__) have been found the sum can be evaluated at any

output point

**x**. Choose between minimum curvature, regularized, or continuous curvature

splines in tension for either 1-D, 2-D, or 3-D Cartesian coordinates or spherical surface

coordinates. After first removing a linear or planar trend (Cartesian geometries) or mean

value (spherical surface) and normalizing these residuals, the least-squares matrix

solution for the spline coefficients

__c__(

__i__) is found by solving the

__n__by

__n__linear system

__w__(

__j__) = sum-over-

__i__{

__c__(

__i__) G(

**x**(

__j__);

**x**(

__i__))}, for

__j__= 1,

__n__; this solution yields an exact

interpolation of the supplied data points. Alternatively, you may choose to perform a

singular value decomposition (SVD) and eliminate the contribution from the smallest

eigenvalues; this approach yields an approximate solution. Trends and scales are restored

when evaluating the output.

**REQUIRED** **ARGUMENTS**

None.

**OPTIONAL** **ARGUMENTS**

__table__The name of one or more ASCII [or binary, see

**-bi**] files holding the

**x**,

__w__data

points. If no file is given then we read standard input instead.

**-A[1|2|3|4|5,]**

__gradfile__

The solution will partly be constrained by surface gradients

**v**=

__v__*

**n**, where

__v__is

the gradient magnitude and

**n**its unit vector direction. The gradient direction may

be specified either by Cartesian components (either unit vector

**n**and magnitude

__v__

separately or gradient components

**v**directly) or angles w.r.t. the coordinate axes.

Specify one of five input formats:

**0**: For 1-D data there is no direction, just

gradient magnitude (slope) so the input format is

__x__,

__gradient__. Options 1-2 are for

2-D data sets:

**1**: records contain

__x__,

__y__,

__azimuth__,

__gradient__(

__azimuth__in degrees is

measured clockwise from the vertical (north) [Default]).

**2**: records contain

__x__,

__y__,

__gradient__,

__azimuth__(

__azimuth__in degrees is measured clockwise from the vertical

(north)). Options 3-5 are for either 2-D or 3-D data:

**3**: records contain

**x**,

__direction(s)__,

__v__(

__direction(s)__in degrees are measured counter-clockwise from the

horizontal (and for 3-D the vertical axis).

**4**: records contain

**x**,

**v**.

**5**: records

contain

**x**,

**n**,

__v__. Append name of ASCII file with the surface gradients (following a

comma if a format is specified).

**-C[n|v]**

__cut__

**[/**

__file__

**]**

Find an approximate surface fit: Solve the linear system for the spline

coefficients by SVD and eliminate the contribution from all eigenvalues whose ratio

to the largest eigenvalue is less than

__cut__[Default uses Gauss-Jordan elimination

to solve the linear system and fit the data exactly]. Optionally, append /

__file__to

save the eigenvalue ratios to the specified file for further analysis. Finally, if

a negative

__cut__is given then /

__file__is required and execution will stop after saving

the eigenvalues, i.e., no surface output is produced. Specify

**-Cv**to use the

largest eigenvalues needed to explain

__cut__% of the data variance. Alternatively,

use

**-Cn**to select the

__cut__largest eigenvalues. If a

__file__is given with

**-Cv**then we

save the eigenvalues instead of the ratios.

**-D**

__mode__Sets the distance flag that determines how we calculate distances between data

points. Select

__mode__0 for Cartesian 1-D spline interpolation:

**-D**0 means (

__x__) in user

units, Cartesian distances, Select

__mode__1-3 for Cartesian 2-D surface spline

interpolation:

**-D**1 means (

__x__,

__y__) in user units, Cartesian distances,

**-D**2 for (

__x__,

__y__) in

degrees, Flat Earth distances, and

**-D**3 for (

__x__,

__y__) in degrees, Spherical distances in

km. Then, if PROJ_ELLIPSOID is spherical, we compute great circle arcs, otherwise

geodesics. Option

__mode__= 4 applies to spherical surface spline interpolation only:

**-D**4 for (

__x__,

__y__) in degrees, use cosine of great circle (or geodesic) arcs. Select

__mode__5 for Cartesian 3-D surface spline interpolation:

**-D**5 means (

__x__,

__y__,

__z__) in user

units, Cartesian distances.

**-G**

__grdfile__

Name of resulting output file. (1) If options

**-R**,

**-I**, and possibly

**-r**are set we

produce an equidistant output table. This will be written to stdout unless

**-G**is

specified. Note: for 2-D grids the

**-G**option is required. (2) If option

**-T**is

selected then

**-G**is required and the output file is a 2-D binary grid file. Applies

to 2-D interpolation only. (3) If

**-N**is selected then the output is an ASCII (or

binary; see

**-bo**) table; if

**-G**is not given then this table is written to standard

output. Ignored if

**-C**or

**-C**0 is given.

**-I**

__xinc__

**[/**

__yinc__

**[/**

__zinc__

**]]**

Specify equidistant sampling intervals, on for each dimension, separated by

slashes.

**-L**Do

__not__remove a linear (1-D) or planer (2-D) trend when

**-D**selects mode 0-3 [For

those Cartesian cases a least-squares line or plane is modeled and removed, then

restored after fitting a spline to the residuals]. However, in mixed cases with

both data values and gradients, or for spherical surface data, only the mean data

value is removed (and later and restored).

**-N**

__nodefile__

ASCII file with coordinates of desired output locations

**x**in the first column(s).

The resulting

__w__values are appended to each record and written to the file given in

**-G**[or stdout if not specified]; see

**-bo**for binary output instead. This option

eliminates the need to specify options

**-R**,

**-I**, and

**-r**.

**-Q**

__az__

**|**

__x/y/z__

Rather than evaluate the surface, take the directional derivative in the

__az__azimuth

and return the magnitude of this derivative instead. For 3-D interpolation, specify

the three components of the desired vector direction (the vector will be normalized

before use).

**-R**

__xmin__

**/**

__xmax__

**[/**

__ymin__

**/**

__ymax__

**[/**

__zminzmax__

**]]**

Specify the domain for an equidistant lattice where output predictions are

required. Requires

**-I**and optionally

**-r**.

__1-D:__Give

__xmin/xmax__, the minimum and maximum

__x__coordinates.

__2-D:__Give

__xmin/xmax/ymin/ymax__, the minimum and maximum

__x__and

__y__coordinates. These

may be Cartesian or geographical. If geographical, then

__west__,

__east__,

__south__, and

__north__specify the Region of interest, and you may specify them in decimal degrees

or in [+-]dd:mm[:ss.xxx][W|E|S|N] format. The two shorthands

**-Rg**and

**-Rd**stand for

global domain (0/360 and -180/+180 in longitude respectively, with -90/+90 in

latitude).

__3-D:__Give

__xmin/xmax/ymin/ymax/zmin/zmax__, the minimum and maximum

__x__,

__y__and

__z__

coordinates. See the 2-D section if your horizontal coordinates are geographical;

note the shorthands

**-Rg**and

**-Rd**cannot be used if a 3-D domain is specified.

**-Sc|t|l|r|p|q[**

__pars__

**]**

Select one of six different splines. The first two are used for 1-D, 2-D, or 3-D

Cartesian splines (see

**-D**for discussion). Note that all tension values are

expected to be normalized tension in the range 0 <

__t__< 1: (

**c**) Minimum curvature

spline [

__Sandwell__, 1987], (

**t**) Continuous curvature spline in tension [

__Wessel__

__and__

__Bercovici__, 1998]; append

__tension__[/

__scale__] with

__tension__in the 0-1 range and

optionally supply a length scale [Default is the average grid spacing]. The next is

a 1-D or 2-D spline: (

**l**) Linear (1-D) or Bilinear (2-D) spline; these produce

output that do not exceed the range of the given data. The next is a 2-D or 3-D

spline: (

**r**) Regularized spline in tension [

__Mitasova__

__and__

__Mitas__, 1993]; again, append

__tension__and optional

__scale__. The last two are spherical surface splines and both

imply

**-D**4: (

**p**) Minimum curvature spline [

__Parker__, 1994], (

**q**) Continuous curvature

spline in tension [

__Wessel__

__and__

__Becker__, 2008]; append

__tension__. The G(

**x'**;

**x'**) for the

last method is slower to compute (a series solution) so we pre-calculate values and

use cubic spline interpolation lookup instead. Optionally append

**+n**

__N__(an odd

integer) to change how many points to use in the spline setup [10001]. The finite

Legendre sum has a truncation error [1e-6]; you can lower that by appending

**+e**

__limit__

at the expense of longer run-time.

**-T**

__maskgrid__

For 2-D interpolation only. Only evaluate the solution at the nodes in the

__maskgrid__

that are not equal to NaN. This option eliminates the need to specify options

**-R**,

**-I**, and

**-r**.

**-V[**

__level__

**]**

**(more**

**...)**

Select verbosity level [c].

**-W**Expect data weights in the final input column, typically given as weight = 1 /

sigma, the data uncertainty. This results in a weighted least squares fit. Note

that this only has an effect if

**-CC**is used.

**-bi[**

__ncols__

**][t]**

**(more**

**...)**

Select native binary input. [Default is 2-4 input columns (

**x**,

__w__); the number depends

on the chosen dimension].

**-bo[**

__ncols__

**][**

__type__

**]**

**(more**

**...)**

Select native binary output.

**-d[i|o]**

__nodata__

**(more**

**...)**

Replace input columns that equal

__nodata__with NaN and do the reverse on output.

**-f[i|o]**

__colinfo__

**(more**

**...)**

Specify data types of input and/or output columns.

**-h[i|o][**

__n__

**][+c][+d][+r**

__remark__

**][+r**

__title__

**]**

**(more**

**...)**

Skip or produce header record(s).

**-i**

__cols__

**[l][s**

__scale__

**][o**

__offset__

**][,**

__...__

**]**

**(more**

**...)**

Select input columns (0 is first column).

**-o**

__cols__

**[,...]**

**(more**

**...)**

Select output columns (0 is first column).

**-r**

**(more**

**...)**

Set pixel node registration [gridline].

**-x[[-]**

__n__

**]**

**(more**

**...)**

Limit number of cores used in multi-threaded algorithms (OpenMP required).

**-^**

**or**

**just**

**-**

Print a short message about the syntax of the command, then exits (NOTE: on Windows

use just

**-**).

**-+**

**or**

**just**

**+**

Print an extensive usage (help) message, including the explanation of any

module-specific option (but not the GMT common options), then exits.

**-?**

**or**

**no**

**arguments**

Print a complete usage (help) message, including the explanation of options, then

exits.

**--version**

Print GMT version and exit.

**--show-datadir**

Print full path to GMT share directory and exit.

**1-D** **EXAMPLES**

To resample the

__x__,

__y__Gaussian random data created by

**gmtmath**and stored in 1D.txt,

requesting output every 0.1 step from 0 to 10, and using a minimum cubic spline, try

gmt math -T0/10/1 0 1 NRAND = 1D.txt

gmt psxy -R0/10/-5/5 -JX6i/3i -B2f1/1 -Sc0.1 -Gblack 1D.txt -K > 1D.ps

gmt greenspline 1D.txt -R0/10 -I0.1 -Sc -V | psxy -R -J -O -Wthin >> 1D.ps

To apply a spline in tension instead, using a tension of 0.7, try

gmt psxy -R0/10/-5/5 -JX6i/3i -B2f1/1 -Sc0.1 -Gblack 1D.txt -K > 1Dt.ps

gmt greenspline 1D.txt -R0/10 -I0.1 -St0.7 -V | psxy -R -J -O -Wthin >> 1Dt.ps

**2-D** **EXAMPLES**

To make a uniform grid using the minimum curvature spline for the same Cartesian data set

from Davis (1986) that is used in the GMT Technical Reference and Cookbook example 16, try

gmt greenspline table_5.11 -R0/6.5/-0.2/6.5 -I0.1 -Sc -V -D1 -GS1987.nc

gmt psxy -R0/6.5/-0.2/6.5 -JX6i -B2f1 -Sc0.1 -Gblack table_5.11 -K > 2D.ps

gmt grdcontour -JX6i -B2f1 -O -C25 -A50 S1987.nc >> 2D.ps

To use Cartesian splines in tension but only evaluate the solution where the input mask

grid is not NaN, try

gmt greenspline table_5.11 -Tmask.nc -St0.5 -V -D1 -GWB1998.nc

To use Cartesian generalized splines in tension and return the magnitude of the surface

slope in the NW direction, try

gmt greenspline table_5.11 -R0/6.5/-0.2/6.5 -I0.1 -Sr0.95 -V -D1 -Q-45 -Gslopes.nc

Finally, to use Cartesian minimum curvature splines in recovering a surface where the

input data is a single surface value (pt.d) and the remaining constraints specify only the

surface slope and direction (slopes.d), use

gmt greenspline pt.d -R-3.2/3.2/-3.2/3.2 -I0.1 -Sc -V -D1 -A1,slopes.d -Gslopes.nc

**3-D** **EXAMPLES**

To create a uniform 3-D Cartesian grid table based on the data in table_5.23 in Davis

(1986) that contains

__x__,

__y__,

__z__locations and a measure of uranium oxide concentrations (in

percent), try

gmt greenspline table_5.23 -R5/40/-5/10/5/16 -I0.25 -Sr0.85 -V -D5 -G3D_UO2.txt

**2-D** **SPHERICAL** **SURFACE** **EXAMPLES**

To recreate Parker's [1994] example on a global 1x1 degree grid, assuming the data are in

file mag_obs_1990.d, try

greenspline -V -Rg -Sp -D3 -I1 -GP1994.nc mag_obs_1990.d

To do the same problem but applying tension of 0.85, use

greenspline -V -Rg -Sq0.85 -D3 -I1 -GWB2008.nc mag_obs_1990.d

**CONSIDERATIONS**

1. For the Cartesian cases we use the free-space Green functions, hence no boundary

conditions are applied at the edges of the specified domain. For most applications

this is fine as the region typically is arbitrarily set to reflect the extent of your

data. However, if your application requires particular boundary conditions then you may

consider using

**surface**instead.

2. In all cases, the solution is obtained by inverting a

__n__x

__n__double precision matrix for

the Green function coefficients, where

__n__is the number of data constraints. Hence, your

computer's memory may place restrictions on how large data sets you can process with

__greenspline__

__<greenspline.html>__. Pre-processing your data with

__blockmean__

__<blockmean.html>__,

__blockmedian__

__<blockmedian.html>__, or

__blockmode__

__<blockmode.html>__is

recommended to avoid aliasing and may also control the size of

__n__. For information, if

__n__

= 1024 then only 8 Mb memory is needed, but for

__n__= 10240 we need 800 Mb. Note that

__greenspline__

__<greenspline.html>__is fully 64-bit compliant if compiled as such. For

spherical data you may consider decimating using

__gmtspatial__

__<gmtspatial.html>__nearest

neighbor reduction.

3. The inversion for coefficients can become numerically unstable when data neighbors are

very close compared to the overall span of the data. You can remedy this by

pre-processing the data, e.g., by averaging closely spaced neighbors. Alternatively,

you can improve stability by using the SVD solution and discard information associated

with the smallest eigenvalues (see

**-C**).

4. The series solution implemented for

**-Sq**was developed by Robert L. Parker, Scripps

Institution of Oceanography, which we gratefully acknowledge.

5. If you need to fit a certain 1-D spline through your data points you may wish to

consider

__sample1d__

__<sample1d.html>__instead. It will offer traditional splines with

standard boundary conditions (such as the natural cubic spline, which sets the

curvatures at the ends to zero). In contrast, greenspline's 1-D spline, as is

explained in note 1, does

__not__specify boundary conditions at the end of the data

domain.

**TENSION**

Tension is generally used to suppress spurious oscillations caused by the minimum

curvature requirement, in particular when rapid gradient changes are present in the data.

The proper amount of tension can only be determined by experimentation. Generally, very

smooth data (such as potential fields) do not require much, if any tension, while rougher

data (such as topography) will typically interpolate better with moderate tension. Make

sure you try a range of values before choosing your final result. Note: the regularized

spline in tension is only stable for a finite range of

__scale__values; you must experiment

to find the valid range and a useful setting. For more information on tension see the

references below.

**REFERENCES**

Davis, J. C., 1986,

__Statistics__

__and__

__Data__

__Analysis__

__in__

__Geology__, 2nd Edition, 646 pp., Wiley,

New York,

Mitasova, H., and L. Mitas, 1993, Interpolation by regularized spline with tension: I.

Theory and implementation,

__Math.__

__Geol.__,

**25**, 641-655.

Parker, R. L., 1994,

__Geophysical__

__Inverse__

__Theory__, 386 pp., Princeton Univ. Press,

Princeton, N.J.

Sandwell, D. T., 1987, Biharmonic spline interpolation of Geos-3 and Seasat altimeter

data,

__Geophys.__

__Res.__

__Lett.__,

**14**, 139-142.

Wessel, P., and D. Bercovici, 1998, Interpolation with splines in tension: a Green's

function approach,

__Math.__

__Geol.__,

**30**, 77-93.

Wessel, P., and J. M. Becker, 2008, Interpolation using a generalized Green's function for

a spherical surface spline in tension,

__Geophys.__

__J.__

__Int__,

**174**, 21-28.

Wessel, P., 2009, A general-purpose Green's function interpolator,

__Computers__

__&__

__Geosciences__,

**35**, 1247-1254, doi:10.1016/j.cageo.2008.08.012.

Use greensplinegmt online using onworks.net services