Let surface use optimal region under the hood

doc/examples/images.dvc

@@ -1,5 +1,5 @@
 outs:
-- md5: f06948ccba06dea0f9a18b63a4d8739b.dir
-  size: 35022240
+- md5: df85591c00343ed16310813c991db900.dir
+  size: 35020779
   nfiles: 53
   path: images

doc/rst/source/surface.rst

@@ -22,7 +22,7 @@ Synopsis
 [ |-L|\ **l**\ *lower* ] [ |-L|\ **u**\ *upper* ]
 [ |-M|\ *max_radius* ]
 [ |-N|\ *max_iterations* ]
-[ |-Q| ]
+[ |-Q|\ [**r**] ]
 [ |-S|\ *search_radius*\ [**m**\|\ **s**] ]
 [ |-T|\ [**b**\|\ **i**]\ *tension_factor* ]
 [ |SYN_OPT-V| ]
@@ -47,22 +47,21 @@ Description
 -----------
 
 **surface** reads randomly-spaced (x,y,z) triples from standard input
-[or *table*] and produces a binary grid file of gridded values z(x,y) by
+[or *table*] and produces a binary file of gridded values z(x,y) by
 solving the differential equation (away from data points)
 
 .. math::
 
     (1 - t) \nabla ^2(z) + t \nabla (z) = 0,
 
 where *t* is a tension factor between 0 and 1, and :math:`\nabla` indicates the
-Laplacian operator. Here, *t* = 0 gives the "minimum curvature" solution which
-is equivalent to SuperMISP and the ISM packages. Minimum curvature can
-cause undesired oscillations and false local maxima or minima (See Smith
-and Wessel, 1990), and you may wish to use *t* > 0 to suppress these
+Laplacian operator. Here, *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 *t* > 0 to suppress these
 effects. Experience suggests *t* ~ 0.25 usually looks good for potential
 field data and *t* should be larger (*t* ~ 0.35) for steep topography data.
 *t* = 1 gives a harmonic surface (no maxima or minima are possible except
-at control data points). It is recommended that the user pre-process the
+at control data points). It is recommended that the user preprocess the
 data with :doc:`blockmean`, :doc:`blockmedian`, or :doc:`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
@@ -129,7 +128,7 @@ Optional Arguments
 .. _-D:
 
 **-D**\ *breakline*\ [**+z**\ [*level*]]
-    Use xyz data in the *breakline* file as a 'soft breakline'. A 'soft breakline'
+    Use *x, y, z* data in the *breakline* file as a 'soft breakline'. A 'soft breakline'
     is a line whose vertices will be used to constrain the nearest grid nodes without
     any further interpolation. A coastline or a lake shore are good examples of
     'soft breaklines'. Multi-segments files are accepted.  If your lines do not have
@@ -139,26 +138,27 @@ Optional Arguments
 .. _-L:
 
 **-Ll**\ *lower* and **-Lu**\ *upper*
-    Impose limits on the output solution. **l**\ *lower* sets the lower
+    Impose limits on the output solution. Directive **l**\ *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]. **u**\ *upper* sets the upper bound and
+    for unconstrained [Default]. Directive **u**\ *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.
+    with NaN is unconstrained. **Note**: Grids given via **-L** must be compatible
+    with the desired output domain and increments.
 
 .. _-M:
 
 **-M**\ *max_radius*
     After solving for the surface, apply a mask so that nodes farther
-    than *max_radius* away from a data constraint is set to NaN [no masking].
+    than *max_radius* away from a data constraint are set to NaN [no masking].
     Append a distance unit (see `Units`_) if needed.
     One can also select the nodes to mask by using the **-M**\ *n_cells*\ **c** form.
     Here *n_cells* means the number of cells around the node controlled by a data point. As an example
-    **-M0c** means that only the cell where point lies is filled, **-M1c** keeps one cell
-    beyond that (i.e. makes a 3x3 neighborhood), and so on.
+    **-M0c** means that only the cell where the point lies is filled, **-M1c** keeps one cell
+    beyond that (i.e. makes a 3x3 square neighborhood), and so on.
 
 .. _-N:
 
@@ -171,13 +171,18 @@ Optional Arguments
 
 .. _-Q:
 
-**-Q**
+**-Q**\ [**r**]
     Suggest grid dimensions which have a highly composite greatest
     common factor. This allows surface to use several intermediate steps
     in the solution, yielding faster run times and better results. The
     sizes suggested by |-Q| can be achieved by altering |-R| and/or
     |-I|. You can recover the |-R| and |-I| you want later by
     using :doc:`grdsample` or :doc:`grdcut` on the output of **surface**.
+    Alternatively, append **r** to have **surface** use the specified
+    **-R** setting exactly as given in the calculations [Default automatically
+    seeks a slightly larger region if that allows for more intermediate
+    steps to ensure the best possible convergence; the region is then trimmed
+    back to what was requested in |-R| upon output].
 
 .. _-S:
 
@@ -268,13 +273,21 @@ Examples
 
 To grid 5 by 5 minute gravity block means from the ASCII data in
 hawaii_5x5.xyg, using a *tension_factor* = 0.25, a
-*convergence_limit* = 0.1 milligal, writing the result to a file called
+*convergence_limit* = 0.1 mGal, writing the result to a file called
 hawaii_grd.nc, and monitoring each iteration, try:
 
    ::
 
     gmt surface hawaii_5x5.xyg -R198/208/18/25 -I5m -Ghawaii_grd.nc -T0.25 -C0.1 -Vi
 
+Notes
+-----
+
+While the region specified by **-R** determines your final output grid, internally
+we may use a slightly larger region that will allow for more intermediate grids
+(i.e., more common factors between *n_columns - 1* and *n_rows - 1*). This
+should allow for better convergence in the final solution.
+
 Gridding Geographic Data: Boundary Conditions
 ---------------------------------------------
 
@@ -289,7 +302,8 @@ will continue.  Because the result is a geographic grid, the GMT i/o machinery w
 interfere and detect inconsistencies at the pole points and replace all values along
 a pole with their mean value.  This will introduce further distortion into the
 grid near the poles.  We recommend you instead consider spherical gridding for global
-data sets; see :doc:`greenspline` (for modest data sets) or :doc:`sphinterpolate`.
+data sets; see :doc:`greenspline` (for modest data sets) or :doc:`sphinterpolate`, or
+project your data using a stereographic projection and grid the projected Cartesian data.
 
 Gridding Geographic Data: Setting Increments
 --------------------------------------------
@@ -327,10 +341,6 @@ specify more decimal places by editing the parameter
 the decimators or choose binary input and/or output using single or
 double precision storage.
 
-Note that only gridline registration is possible with **surface**. If
-you need a pixel-registered grid you can resample a gridline registered
-grid using :doc:`grdsample` |-T|.
-
 See Also
 --------