Source code for pygmt.src.surface

"""
surface - Grid table data using adjustable tension continuous curvature
splines.
"""
from pygmt.clib import Session
from pygmt.helpers import (
    GMTTempFile,
    build_arg_string,
    fmt_docstring,
    kwargs_to_strings,
    use_alias,
)
from pygmt.io import load_dataarray

__doctest_skip__ = ["surface"]


[docs]@fmt_docstring @use_alias( C="convergence", G="outgrid", I="spacing", Ll="lower", Lu="upper", M="maxradius", R="region", T="tension", V="verbose", a="aspatial", b="binary", d="nodata", e="find", f="coltypes", h="header", i="incols", r="registration", w="wrap", ) @kwargs_to_strings(I="sequence", R="sequence") def surface(data=None, x=None, y=None, z=None, **kwargs): r""" Grid table data using adjustable tension continuous curvature splines. Surface reads randomly-spaced (x, y, z) triplets and produces gridded values z(x,y) by solving: .. math:: (1 - t)\nabla^2(z)+t\nabla(z) = 0 where :math:`t` is a tension factor between 0 and 1, and :math:`\nabla` indicates the Laplacian operator. Here, :math:`t = 0` gives the "minimum curvature" solution. Minimum curvature can cause undesired oscillations and false local maxima or minima (see Smith and Wessel, 1990), and you may wish to use :math:`t > 0` to suppress these effects. Experience suggests :math:`t \sim 0.25` usually looks good for potential field data and t should be larger (:math:`t \sim 0.35`) for steep topography data. :math:`t = 1` gives a harmonic surface (no maxima or minima are possible except at control data points). It is recommended that the user preprocess the data with :func:`pygmt.blockmean`, :func:`pygmt.blockmedian`, or :func:`pygmt.blockmode` to avoid spatial aliasing and eliminate redundant data. You may impose lower and/or upper bounds on the solution. These may be entered in the form of a fixed value, a grid with values, or simply be the minimum/maximum input data values. Natural boundary conditions are applied at the edges, except for geographic data with 360-degree range where we apply periodic boundary conditions in the longitude direction. Takes a matrix, (x, y, z) triplets, or a file name as input. Must provide either ``data`` or ``x``, ``y``, and ``z``. Full option list at :gmt-docs:`surface.html` {aliases} Parameters ---------- data : str or {table-like} Pass in (x, y, z) or (longitude, latitude, elevation) values by providing a file name to an ASCII data table, a 2-D {table-classes}. x/y/z : 1-D arrays Arrays of x and y coordinates and values z of the data points. {spacing} {region} outgrid : str Optional. The file name for the output netCDF file with extension .nc to store the grid in. convergence : float Optional. Convergence limit. Iteration is assumed to have converged when the maximum absolute change in any grid value is less than ``convergence``. (Units same as data z units). Alternatively, give limit in percentage of root-mean-square (rms) deviation by appending %. [Default is scaled to :math:`10^{{-4}}` of the rms deviation of the data from a best-fit (least-squares) plane.] This is the final convergence limit at the desired grid spacing; for intermediate (coarser) grids the effective convergence limit is divided by the grid spacing multiplier. maxradius : int or str Optional. After solving for the surface, apply a mask so that nodes farther than ``maxradius`` away from a data constraint are set to NaN [Default is no masking]. Append a distance unit (see :gmt-docs:`Units <surface.html#units>`) if needed. One can also select the nodes to mask by using the *n_cells*\ **c** form. Here *n_cells* means the number of cells around the node is controlled by a data point. As an example ``"0c"`` means that only the cell where the point lies is filled, ``"1c"`` keeps one cell beyond that (i.e. makes a 3x3 square neighborhood), and so on. lower : float or str Optional. Impose limits on the output solution. Parameter ``lower`` sets the lower bound. ``lower`` can be the name of a grid file with lower bound values, a fixed value, **d** to set to minimum input value, or **u** for unconstrained [Default]. Grid files used to set the limits may contain NaNs. In the presence of NaNs, the limit of a node masked with NaN is unconstrained. upper : float or str Optional. Impose limits on the output solution. Parameter ``upper`` sets the upper bound and can be the name of a grid file with upper bound values, a fixed value, **d** to set to maximum input value, or **u** for unconstrained [Default]. Grid files used to set the limits may contain NaNs. In the presence of NaNs, the limit of a node masked with NaN is unconstrained. tension : float or str [**b**\|\ **i**]. Optional. Tension factor[s]. These must be between 0 and 1. Tension may be used in the interior solution (above equation, where it suppresses spurious oscillations) and in the boundary conditions (where it tends to flatten the solution approaching the edges). Add **i**\ *tension* to set interior tension, and **b**\ *tension* to set boundary tension. If you do not prepend **i** or **b**, both will be set to the same value. [Default is 0 for both and gives minimum curvature solution.] {verbose} {aspatial} {binary} {nodata} {find} {coltypes} {header} {incols} {registration} {wrap} Returns ------- ret: xarray.DataArray or None Return type depends on whether the ``outgrid`` parameter is set: - :class:`xarray.DataArray`: if ``outgrid`` is not set - None if ``outgrid`` is set (grid output will be stored in file set by ``outgrid``) Example ------- >>> import pygmt >>> # Load a sample table of topography >>> topography = pygmt.datasets.load_sample_data( ... name="notre_dame_topography" ... ) >>> # Perform gridding of topography data >>> grid = pygmt.surface(data=topography, spacing=1, region=[0, 4, 0, 8]) """ with GMTTempFile(suffix=".nc") as tmpfile: with Session() as lib: # Choose how data will be passed into the module file_context = lib.virtualfile_from_data( check_kind="vector", data=data, x=x, y=y, z=z, required_z=True ) with file_context as infile: if (outgrid := kwargs.get("G")) is None: kwargs["G"] = outgrid = tmpfile.name # output to tmpfile lib.call_module( module="surface", args=build_arg_string(kwargs, infile=infile) ) return load_dataarray(outgrid) if outgrid == tmpfile.name else None