diff --git a/doc/api/index.rst b/doc/api/index.rst index fc69a95d1e2..1e74629c5af 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -92,6 +92,7 @@ Operations on grids: grdcut grdfill grdfilter + grdgradient grdtrack Crossover analysis with x2sys: diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 5cbea0fc853..57e223d5c14 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -37,6 +37,7 @@ grdcut, grdfill, grdfilter, + grdgradient, grdinfo, grdtrack, info, diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py index b4855a37737..4e79d0d96f7 100644 --- a/pygmt/src/__init__.py +++ b/pygmt/src/__init__.py @@ -15,6 +15,7 @@ from pygmt.src.grdcut import grdcut from pygmt.src.grdfill import grdfill from pygmt.src.grdfilter import grdfilter +from pygmt.src.grdgradient import grdgradient from pygmt.src.grdimage import grdimage from pygmt.src.grdinfo import grdinfo from pygmt.src.grdtrack import grdtrack diff --git a/pygmt/src/grdgradient.py b/pygmt/src/grdgradient.py new file mode 100644 index 00000000000..8dca2018c10 --- /dev/null +++ b/pygmt/src/grdgradient.py @@ -0,0 +1,118 @@ +""" +grdgradient - Compute directional gradients from a grid. +""" + +import xarray as xr +from pygmt.clib import Session +from pygmt.exceptions import GMTInvalidInput +from pygmt.helpers import ( + GMTTempFile, + args_in_kwargs, + build_arg_string, + fmt_docstring, + kwargs_to_strings, + use_alias, +) + + +@fmt_docstring +@use_alias( + A="azimuth", + D="direction", + E="radiance", + G="outgrid", + R="region", + V="verbose", + n="interpolation", +) +@kwargs_to_strings(A="sequence", E="sequence", R="sequence") +def grdgradient(grid, **kwargs): + r""" + Compute the directional derivative of the vector gradient of the data. + + Can accept ``azimuth``, ``direction``, and ``radiance`` input to create + the resulting gradient. + + Full option list at :gmt-docs:`grdgradient.html` + + {aliases} + + Parameters + ---------- + grid : str or xarray.DataArray + The file name of the input grid or the grid loaded as a DataArray. + outgrid : str or None + The name of the output netCDF file with extension .nc to store the grid + in. + azimuth : int or float or str or list + *azim*\ [/*azim2*]. + Azimuthal direction for a directional derivative; *azim* is the + angle in the x,y plane measured in degrees positive clockwise from + north (the +y direction) toward east (the +x direction). The + negative of the directional derivative, -[dz/dx\*sin(*azim*) + + dz/dy\*cos(\ *azim*)], is found; negation yields positive values + when the slope of z(x,y) is downhill in the *azim* direction, the + correct sense for shading the illumination of an image by a light + source above the x,y plane shining from the *azim* direction. + Optionally, supply two azimuths, *azim*/*azim2*, in which case the + gradients in each of these directions are calculated and the one + larger in magnitude is retained; this is useful for illuminating data + with two directions of lineated structures, e.g., *0*/*270* + illuminates from the north (top) and west (left). Finally, if *azim* + is a file it must be a grid of the same domain, spacing and + registration as *grid* that will update the azimuth at each output + node when computing the directional derivatives. + direction : str + [**a**][**c**][**o**][**n**]. + Find the direction of the positive (up-slope) gradient of the data. + To instead find the aspect (the down-slope direction), use **a**. + By default, directions are measured clockwise from north, as *azim* + in ``azimuth``. Append **c** to use conventional Cartesian angles + measured counterclockwise from the positive x (east) direction. + Append **o** to report orientations (0-180) rather than + directions (0-360). Append **n** to add 90 degrees to all angles + (e.g., to give local strikes of the surface). + radiance : str or list + [**m**\|\ **s**\|\ **p**]\ *azim/elev*\ [**+a**\ *ambient*][**+d**\ + *diffuse*][**+p**\ *specular*][**+s**\ *shine*]. + Compute Lambertian radiance appropriate to use with ``grdimage`` + and ``grdview``. The Lambertian Reflection assumes an ideal surface + that reflects all the light that strikes it and the surface appears + equally bright from all viewing directions. Here, *azim* and *elev* are + the azimuth and elevation of the light vector. Optionally, supply + *ambient* [0.55], *diffuse* [0.6], *specular* [0.4], or *shine* [10], + which are parameters that control the reflectance properties of the + surface. Default values are given in the brackets. Use **s** for a + simpler Lambertian algorithm. Note that with this form you only have + to provide azimuth and elevation. Alternatively, use **p** for + the Peucker piecewise linear approximation (simpler but faster + algorithm; in this case the *azim* and *elev* are hardwired to 315 + and 45 degrees. This means that even if you provide other values + they will be ignored.) + {R} + {V} + {n} + """ + with GMTTempFile(suffix=".nc") as tmpfile: + if not args_in_kwargs(args=["A", "D", "E"], kwargs=kwargs): + raise GMTInvalidInput( + """At least one of the following parameters must be specified: + azimuth, direction, or radiance""" + ) + with Session() as lib: + file_context = lib.virtualfile_from_data(check_kind="raster", data=grid) + with file_context as infile: + if "G" not in kwargs.keys(): # if outgrid is unset, output to tempfile + kwargs.update({"G": tmpfile.name}) + outgrid = kwargs["G"] + arg_str = " ".join([infile, build_arg_string(kwargs)]) + lib.call_module("grdgradient", arg_str) + + if outgrid == tmpfile.name: # if user did not set outgrid, return DataArray + with xr.open_dataarray(outgrid) as dataarray: + result = dataarray.load() + _ = result.gmt # load GMTDataArray accessor information + else: + result = None # if user sets an outgrid, return None + + return result diff --git a/pygmt/tests/test_grdgradient.py b/pygmt/tests/test_grdgradient.py new file mode 100644 index 00000000000..4b01038ab6b --- /dev/null +++ b/pygmt/tests/test_grdgradient.py @@ -0,0 +1,62 @@ +""" +Tests for grdgradient. +""" +import os + +import numpy.testing as npt +import pytest +from pygmt import grdgradient, grdinfo +from pygmt.datasets import load_earth_relief +from pygmt.exceptions import GMTInvalidInput +from pygmt.helpers import GMTTempFile + + +@pytest.fixture(scope="module", name="grid") +def fixture_grid(): + """ + Load the grid data from the sample earth_relief file. + """ + return load_earth_relief(resolution="01d", region=[-5, 5, -5, 5]) + + +def test_grdgradient_outgrid(grid): + """ + Test the azimuth and direction parameters for grdgradient with a set + outgrid. + """ + with GMTTempFile(suffix=".nc") as tmpfile: + result = grdgradient(grid=grid, outgrid=tmpfile.name, azimuth=10, direction="c") + assert result is None # return value is None + assert os.path.exists(path=tmpfile.name) # check that outgrid exists + result = ( + grdinfo(grid=tmpfile.name, force_scan="a", per_column="n").strip().split() + ) + npt.assert_allclose(float(result[4]), -0.0045060496) # min + npt.assert_allclose(float(result[5]), 0.0575332976) # max + # Check spherically weighted statistics below + npt.assert_allclose(float(result[10]), 0.000384754501283) # median + npt.assert_allclose(float(result[12]), 0.00285958005568) # mean + + +def test_grdgradient_no_outgrid(grid): + """ + Test the azimuth and direction parameters for grdgradient with no set + outgrid. + """ + temp_grid = grdgradient(grid=grid, azimuth=10, direction="c") + assert temp_grid.dims == ("lat", "lon") + assert temp_grid.gmt.gtype == 1 # Geographic grid + assert temp_grid.gmt.registration == 1 # Pixel registration + npt.assert_allclose(temp_grid.min(), -0.0045060496) + npt.assert_allclose(temp_grid.max(), 0.0575332976) + npt.assert_allclose(temp_grid.median(), 0.0004889865522272885) + npt.assert_allclose(temp_grid.mean(), 0.0028633063193410635) + + +def test_grdgradient_fails(grid): + """ + Check that grdgradient fails correctly when neither of azimuth, direction + or radiance is given. + """ + with pytest.raises(GMTInvalidInput): + grdgradient(grid=grid)