pygmt.grdgradient: Add parameters 'normalize'/'norm_ambient'/'norm_amp'/'norm_sigma'/'norm_offset' for normalization

examples/gallery/images/grdgradient_shading.py

@@ -1,21 +1,23 @@
 """
-Calculating grid gradient with custom ``azimuth`` and ``normalize`` parameters
-==============================================================================
+Calculating grid gradient with custom azimuth and normalization parameters
+==========================================================================
 
-The :func:`pygmt.grdgradient` function calculates the gradient of a grid file.
-As input, :func:`pygmt.grdgradient` gets an :class:`xarray.DataArray` object
-or a path string to a grid file. It then calculates the respective gradient
-and returns an :class:`xarray.DataArray` object. The example below sets two
-main parameters:
+The :func:`pygmt.grdgradient` function calculates the gradient of a grid file. As input,
+:func:`pygmt.grdgradient` gets an :class:`xarray.DataArray` object or a path string to a
+grid file. It then calculates the respective gradient and returns an
+:class:`xarray.DataArray` object.
 
-- ``azimuth`` to set the illumination light source direction (0° is North,
-  90° is East, 180° is South, 270° is West).
-- ``normalize`` to enhance the three-dimensional sense of the topography.
+The ``normalize`` parameter calculates the azimuthal gradient of each point along a
+certain azimuth angle, then adjusts the brightness value of the color according to the
+positive/negative of the azimuthal gradient and the amplitude of each point. The example
+below shows how to customize the gradient by setting azimuth and normalization
+parameters.
 
-The ``normalize`` parameter calculates the azimuthal gradient of each point
-along a certain azimuth angle, then adjusts the brightness value of the color
-according to the positive/negative of the azimuthal gradient and the amplitude
-of each point.
+- ``azimuth`` sets the illumination light source direction (0° is North, 90° is East,
+  180° is South, 270° is West).
+- ``normalize`` sets the normalization method (e.g., Cauchy or Laplace distribution)
+- ``norm_amp`` sets the amplitude of the normalization
+- more parameters are available to further enhance the 3-D sense of the topography
 """
 
 # %%
@@ -28,7 +30,7 @@
 fig = pygmt.Figure()
 
 # Define a colormap to be used for topography
-pygmt.makecpt(cmap="gmt/terra", series=[-7000, 7000])
+pygmt.makecpt(cmap="gmt/terra", series=(-7000, 7000))
 
 # Define figure configuration
 pygmt.config(FONT_TITLE="10p,5", MAP_TITLE_OFFSET="1p", MAP_FRAME_TYPE="plain")
@@ -41,24 +43,28 @@
     sharex="b",
     sharey="l",
 ):
-    # E.g. "0/90" illuminates light source from the North (top) and East
-    # (right), and so on
-    for azi in ["0/90", "0/300", "180/225"]:
-        # "e" and "t" are cumulative Laplace distribution and cumulative
-        # Cauchy distribution, respectively
-        # "amp" (e.g. 1 or 10) controls the brightness value of the color
-        for nor in ["t1", "e1", "t10", "e10"]:
-            # Making an intensity DataArray using azimuth and normalize
-            # parameters
-            shade = pygmt.grdgradient(grid=grid, azimuth=azi, normalize=nor)
-            fig.grdimage(
-                grid=grid,
-                shading=shade,
-                projection="M?",
-                frame=["a4f2", f"+tazimuth={azi}, normalize={nor}"],
-                cmap=True,
-                panel=True,
-            )
+    # Setting azimuth angles, e.g. (0, 90) illuminates light source from the North (top)
+    # and East (right).
+    for azi in [(0, 90), (0, 300), (180, 225)]:
+        # "cauchy"/"laplace" sets cumulative Cauchy/Laplace distribution, respectively.
+        for normalize in ("cauchy", "laplace"):
+            # amp (e.g., 1 or 10) controls the brightness value of the color.
+            for amp in (1, 10):
+                # Making an intensity DataArray using azimuth and normalize parameters
+                shade = pygmt.grdgradient(
+                    grid=grid, azimuth=azi, normalize=normalize, norm_amp=amp
+                )
+                fig.grdimage(
+                    grid=grid,
+                    shading=shade,
+                    projection="M?",
+                    frame=[
+                        "a4f2",
+                        f"+tazimuth={azi}, normalize={normalize}, amp={amp}",
+                    ],
+                    cmap=True,
+                    panel=True,
+                )
 
 fig.colorbar(
     position=Position("BC", cstype="outside"),

pygmt/src/grdgradient.py

@@ -9,26 +9,80 @@
 from pygmt._typing import PathLike
 from pygmt.alias import Alias, AliasSystem
 from pygmt.clib import Session
-from pygmt.exceptions import GMTParameterError
+from pygmt.exceptions import GMTInvalidInput, GMTParameterError
 from pygmt.helpers import build_arg_list, fmt_docstring, use_alias
 
 __doctest_skip__ = ["grdgradient"]
 
 
+def _alias_option_N(  # noqa: N802
+    normalize=False,
+    norm_amp=None,
+    norm_ambient=None,
+    norm_sigma=None,
+    norm_offset=None,
+):
+    """
+    Helper function to create the alias list for the -N option.
+
+    Examples
+    --------
+    >>> def parse(**kwargs):
+    ...     return AliasSystem(N=_alias_option_N(**kwargs)).get("N")
+    >>> parse(normalize=True)
+    ''
+    >>> parse(normalize="laplace")
+    'e'
+    >>> parse(normalize="cauchy")
+    't'
+    >>> parse(
+    ...     normalize="laplace",
+    ...     norm_amp=2,
+    ...     norm_offset=10,
+    ...     norm_sigma=0.5,
+    ...     norm_ambient=0.1,
+    ... )
+    'e2+a0.1+s0.5+o10'
+    >>> # Check for backward compatibility with old syntax
+    >>> parse(normalize="e2+a0.2+s0.5+o10")
+    'e2+a0.2+s0.5+o10'
+    """
+    _normalize_mapping = {"laplace": "e", "cauchy": "t"}
+    # Check for old syntax for normalize
+    if isinstance(normalize, str) and normalize not in _normalize_mapping:
+        if any(
+            v is not None and v is not False
+            for v in [norm_amp, norm_ambient, norm_sigma, norm_offset]
+        ):
+            msg = (
+                "Parameter 'normalize' using old syntax with raw GMT CLI string "
+                "conflicts with parameters 'norm_amp', 'norm_ambient', 'norm_sigma', "
+                "or 'norm_offset'."
+            )
+            raise GMTInvalidInput(msg)
+        _normalize_mapping = None
+
+    return [
+        Alias(normalize, name="normalize", mapping=_normalize_mapping),
+        Alias(norm_amp, name="norm_amp"),
+        Alias(norm_ambient, name="norm_ambient", prefix="+a"),
+        Alias(norm_sigma, name="norm_sigma", prefix="+s"),
+        Alias(norm_offset, name="norm_offset", prefix="+o"),
+    ]
+
+
 @fmt_docstring
-@use_alias(
-    D="direction",
-    N="normalize",
-    Q="tiles",
-    S="slope_file",
-    f="coltypes",
-    n="interpolation",
-)
-def grdgradient(
+@use_alias(D="direction", Q="tiles", S="slope_file", f="coltypes", n="interpolation")
+def grdgradient(  # noqa: PLR0913
     grid: PathLike | xr.DataArray,
     outgrid: PathLike | None = None,
     azimuth: float | Sequence[float] | None = None,
     radiance: Sequence[float] | str | None = None,
+    normalize: Literal["laplace", "cauchy"] | bool = False,
+    norm_amp: float | None = None,
+    norm_ambient: float | None = None,
+    norm_sigma: float | None = None,
+    norm_offset: float | None = None,
     region: Sequence[float | str] | str | None = None,
     verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"]
     | bool = False,
@@ -49,6 +103,12 @@ def grdgradient(
        - R = region
        - V = verbose
 
+    .. hlist::
+       :columns: 1
+
+       - N = normalize, norm_amp, **+a**: norm_ambient, **+s**: norm_sigma,
+         **+o**: norm_offset
+
     Parameters
     ----------
     $grid
@@ -101,31 +161,37 @@ def grdgradient(
         algorithm; in this case *azim* and *elev* are hardwired to 315
         and 45 degrees. This means that even if you provide other values
         they will be ignored.).
-    normalize : str or bool
-        [**e**\|\ **t**][*amp*][**+a**\ *ambient*][**+s**\ *sigma*]\
-        [**+o**\ *offset*].
-        The actual gradients :math:`g` are offset and scaled to produce
-        normalized gradients :math:`g_n` with a maximum output magnitude of
-        *amp*. If *amp* is not given, default *amp* = 1. If *offset* is not
-        given, it is set to the average of :math:`g`. The following forms are
-        supported:
-
-        - **True**: Normalize using math:`g_n = \mbox{amp}\
-          (\frac{g - \mbox{offset}}{max(|g - \mbox{offset}|)})`
-        - **e**: Normalize using a cumulative Laplace distribution yielding:
-          :math:`g_n = \mbox{amp}(1 - \
-          \exp{(\sqrt{2}\frac{g - \mbox{offset}}{\sigma}))}`, where
-          :math:`\sigma` is estimated using the L1 norm of
-          :math:`(g - \mbox{offset})` if it is not given.
-        - **t**: Normalize using a cumulative Cauchy distribution yielding:
-          :math:`g_n = \
-          \frac{2(\mbox{amp})}{\pi}(\tan^{-1}(\frac{g - \
-          \mbox{offset}}{\sigma}))` where :math:`\sigma` is estimated
-          using the L2 norm of :math:`(g - \mbox{offset})` if it is not
-          given.
-
-        As a final option, you may add **+a**\ *ambient* to add *ambient* to
-        all nodes after gradient calculations are completed.
+    normalize
+        Normalize the output gradients. Valid values are:
+
+        - ``False``: No normalization is done [Default].
+        - ``True``: Normalize using max absolute value.
+        - ``"laplace"``: Normalize using cumulative Laplace distribution.
+        - ``"cauchy"``: Normalize using cumulative Cauchy distribution.
+
+        The normalization process is controlled via the additional parameters
+        ``norm_amp``, ``norm_ambient``, ``norm_sigma``, and ``norm_offset``.
+
+        Let :math:`g` denote the actual gradients, :math:`g_n` the normalized gradients,
+        :math:`a` the maximum output magnitude (``norm_amp``), :math:`o` the offset
+        value (``norm_offset``), and :math:`\sigma` the sigma value (``norm_sigma``).
+        The normalization is computed as follows:
+
+        - ``True``: :math:`g_n = a (\frac{g - o}{\max(|g - o|)})`
+        - ``"laplace"``: :math:`g_n = a(1 - \exp(\sqrt{2}\frac{g - o}{\sigma}))`
+        - ``"cauchy"``: :math:`g_n = \frac{2a}{\pi}\arctan(\frac{g - o}{\sigma})`
+    norm_amp
+        Set the maximum output magnitude [Default is 1].
+    norm_ambient
+        The ambient value to add to all nodes after gradient calculations are completed
+        [Default is 0].
+    norm_offset
+        The offset value used in the normalization. If not given, it is set to the
+        average of :math:`g`.
+    norm_sigma
+        The *sigma* value used in the Laplace or Cauchy normalization. If not given,
+        it is estimated from the L1 norm of :math:`g-o` for Laplace or the L2 norm of
+        :math:`g-o` for Cauchy.
     tiles : str
         **c**\|\ **r**\|\ **R**.
         Control how normalization via ``normalize`` is carried out. When
@@ -137,10 +203,10 @@ def grdgradient(
         grid output is not needed for this run then do not specify
         ``outgrid``. For subsequent runs, just use **r** to read these
         values. Using **R** will read then delete the statistics file.
-    $region
     slope_file : str
         Name of output grid file with scalar magnitudes of gradient vectors.
         Requires ``direction`` but makes ``outgrid`` optional.
+    $region
     $verbose
     $coltypes
     $interpolation
@@ -179,6 +245,13 @@ def grdgradient(
     aliasdict = AliasSystem(
         A=Alias(azimuth, name="azimuth", sep="/", size=2),
         E=Alias(radiance, name="radiance", sep="/", size=2),
+        N=_alias_option_N(
+            normalize=normalize,
+            norm_amp=norm_amp,
+            norm_ambient=norm_ambient,
+            norm_sigma=norm_sigma,
+            norm_offset=norm_offset,
+        ),
     ).add_common(
         R=region,
         V=verbose,

pygmt/tests/test_grdgradient.py

@@ -8,7 +8,7 @@
 import xarray as xr
 from pygmt import grdgradient
 from pygmt.enums import GridRegistration, GridType
-from pygmt.exceptions import GMTParameterError
+from pygmt.exceptions import GMTInvalidInput, GMTParameterError
 from pygmt.helpers import GMTTempFile
 from pygmt.helpers.testing import load_static_earth_relief
 
@@ -73,6 +73,22 @@ def test_grdgradient_no_outgrid(grid, expected_grid):
     xr.testing.assert_allclose(a=result, b=expected_grid)
 
 
+def test_grdgradient_normalize_mixed_syntax(grid):
+    """
+    Test that mixed syntax for normalize in grdgradient raises GMTInvalidInput.
+    """
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, azimuth=10, normalize="t", norm_amp=2)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, azimuth=10, normalize="t1", norm_amp=2)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, azimuth=10, normalize="t1", norm_sigma=2)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, azimuth=10, normalize="e", norm_offset=5)
+    with pytest.raises(GMTInvalidInput):
+        grdgradient(grid=grid, azimuth=10, normalize="e", norm_ambient=0.1)
+
+
 def test_grdgradient_fails(grid):
     """
     Check that grdgradient fails correctly.