diff --git a/examples/gallery/images/grdgradient_shading.py b/examples/gallery/images/grdgradient_shading.py index 14c2567d1a4..a2fc6978433 100644 --- a/examples/gallery/images/grdgradient_shading.py +++ b/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"), diff --git a/pygmt/src/grdgradient.py b/pygmt/src/grdgradient.py index 6ac656f181d..58174231739 100644 --- a/pygmt/src/grdgradient.py +++ b/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, diff --git a/pygmt/tests/test_grdgradient.py b/pygmt/tests/test_grdgradient.py index 40b27414c17..e5d9616a3fb 100644 --- a/pygmt/tests/test_grdgradient.py +++ b/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.