diff --git a/doc/api/index.rst b/doc/api/index.rst index 74bd2f182c4..8e888419b57 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -94,6 +94,7 @@ Operations on grids: grdfilter grdgradient grdlandmask + grdproject grdsample grdtrack xyz2grd diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 25ad96a2127..96fe5cf908d 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -40,6 +40,7 @@ grdgradient, grdinfo, grdlandmask, + grdproject, grdsample, grdtrack, info, diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py index 2e1f8b56186..c7a15ce8715 100644 --- a/pygmt/src/__init__.py +++ b/pygmt/src/__init__.py @@ -19,6 +19,7 @@ from pygmt.src.grdimage import grdimage from pygmt.src.grdinfo import grdinfo from pygmt.src.grdlandmask import grdlandmask +from pygmt.src.grdproject import grdproject from pygmt.src.grdsample import grdsample from pygmt.src.grdtrack import grdtrack from pygmt.src.grdview import grdview diff --git a/pygmt/src/grdproject.py b/pygmt/src/grdproject.py new file mode 100644 index 00000000000..c8a5d539421 --- /dev/null +++ b/pygmt/src/grdproject.py @@ -0,0 +1,93 @@ +""" +grdproject - Forward and inverse map transformation of grids. +""" + +import xarray as xr +from pygmt.clib import Session +from pygmt.exceptions import GMTInvalidInput +from pygmt.helpers import ( + GMTTempFile, + build_arg_string, + fmt_docstring, + kwargs_to_strings, + use_alias, +) + + +@fmt_docstring +@use_alias( + G="outgrid", + J="projection", + I="inverse", + R="region", + V="verbose", + n="interpolation", + r="registration", +) +@kwargs_to_strings(R="sequence") +def grdproject(grid, **kwargs): + r""" + Change projection of gridded data between geographical and rectangular. + + This module will project a geographical gridded data set onto a + rectangular grid. If ``inverse`` is ``True``, it will project a + rectangular coordinate system to a geographic system. To obtain the value + at each new node, its location is inversely projected back onto the input + grid after which a value is interpolated between the surrounding input + grid values. By default bi-cubic interpolation is used. Aliasing is + avoided by also forward projecting the input grid nodes. If two or more + nodes are projected onto the same new node, their average will dominate in + the calculation of the new node value. Interpolation and aliasing is + controlled with the ``interpolation`` option. The new node spacing may be + determined in one of several ways by specifying the grid spacing, number + of nodes, or resolution. Nodes not constrained by input data are set to + NaN. The ``region`` parameter can be used to select a map region larger or + smaller than that implied by the extent of the grid file. + + {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. + inverse : bool + When set to ``True`` transforms grid from rectangular to + geographical [Default is False]. + {J} + {R} + {V} + {n} + {r} + + 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``) + """ + if "J" not in kwargs.keys(): + raise GMTInvalidInput("The projection must be specified.") + with GMTTempFile(suffix=".nc") as tmpfile: + 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("grdproject", 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_grdproject.py b/pygmt/tests/test_grdproject.py new file mode 100644 index 00000000000..851272a6115 --- /dev/null +++ b/pygmt/tests/test_grdproject.py @@ -0,0 +1,59 @@ +""" +Tests for grdproject. +""" +import os + +import numpy.testing as npt +import pytest +from pygmt import grdinfo, grdproject +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_grdproject_file_out(grid): + """ + grdproject with an outgrid set. + """ + with GMTTempFile(suffix=".nc") as tmpfile: + result = grdproject(grid=grid, projection="M10c", outgrid=tmpfile.name) + assert result is None # return value is None + assert os.path.exists(path=tmpfile.name) # check that outgrid exists + result = grdinfo(tmpfile.name, per_column=True).strip().split() + npt.assert_allclose(float(result[0]), 0) # x min + npt.assert_allclose(float(result[1]), 10) # x max + npt.assert_allclose(float(result[2]), 0, atol=1.0e-10) # y min + npt.assert_allclose(float(result[3]), 9.94585661273) # y max + npt.assert_allclose(float(result[4]), -5130.48193359) # min + npt.assert_allclose(float(result[5]), -152.585281372) # max + + +def test_grdproject_no_outgrid(grid): + """ + Test grdproject with no set outgrid. + """ + assert grid.gmt.gtype == 1 # Geographic grid + temp_grid = grdproject(grid=grid, projection="M10c") + assert temp_grid.dims == ("y", "x") + assert temp_grid.gmt.gtype == 0 # Rectangular grid + assert temp_grid.gmt.registration == 1 # Pixel registration + npt.assert_allclose(temp_grid.min(), -5130.482) + npt.assert_allclose(temp_grid.max(), -152.58528) + npt.assert_allclose(temp_grid.median(), -4578.288) + npt.assert_allclose(temp_grid.mean(), -4354.3296) + + +def test_grdproject_fails(grid): + """ + Check that grdproject fails correctly. + """ + with pytest.raises(GMTInvalidInput): + grdproject(grid=grid)