diff --git a/doc/api/index.rst b/doc/api/index.rst index c2974884a65..28683d3b75c 100644 --- a/doc/api/index.rst +++ b/doc/api/index.rst @@ -70,6 +70,7 @@ Operations on grids: .. autosummary:: :toctree: generated + grdcut grdinfo grdtrack diff --git a/pygmt/__init__.py b/pygmt/__init__.py index 1e4b4ba16af..963cbf39d48 100644 --- a/pygmt/__init__.py +++ b/pygmt/__init__.py @@ -19,6 +19,7 @@ from .sampling import grdtrack from .mathops import makecpt from .modules import config, info, grdinfo, which +from .gridops import grdcut from . import datasets diff --git a/pygmt/gridops.py b/pygmt/gridops.py new file mode 100644 index 00000000000..400a8d6f219 --- /dev/null +++ b/pygmt/gridops.py @@ -0,0 +1,119 @@ +""" +GMT modules for grid operations +""" + +import xarray as xr + + +from .clib import Session +from .helpers import ( + build_arg_string, + fmt_docstring, + kwargs_to_strings, + GMTTempFile, + use_alias, + data_kind, + dummy_context, +) +from .exceptions import GMTInvalidInput + + +@fmt_docstring +@use_alias( + G="outgrid", + R="region", + J="projection", + N="extend", + S="circ_subregion", + Z="z_subregion", +) +@kwargs_to_strings(R="sequence") +def grdcut(grid, **kwargs): + """ + Extract subregion from a grid. + + Produce a new *outgrid* file which is a subregion of *grid*. The + subregion is specified with *region*; the specified range must not exceed + the range of *grid* (but see *extend*). If in doubt, run + :meth:`pygmt.grdinfo` to check range. Alternatively, define the subregion + indirectly via a range check on the node values or via distances from a + given point. Finally, you can give *projection* for oblique projections to + determine the corresponding rectangular *region* setting that will give a + grid that fully covers the oblique domain. + + Full option list at :gmt-docs:`grdcut.html` + + {aliases} + + Parameters + ---------- + grid : str + The name of the input grid file. + outgrid : str or None + The name of the output netCDF file with extension .nc to store the grid + in. + {J} + {R} + extend : bool or int or float + Allow grid to be extended if new *region* exceeds existing boundaries. + Give a value to initialize nodes outside current region. + circ_subregion : str + ``'lon/lat/radius[unit][+n]'``. + Specify an origin (*lon* and *lat*) and *radius*; append a distance + *unit* and we determine the corresponding rectangular region so that + all grid nodes on or inside the circle are contained in the subset. + If **+n** is appended we set all nodes outside the circle to NaN. + z_subregion : str + ``'[min/max][+n|N|r]'``. + Determine a new rectangular region so that all nodes outside this + region are also outside the given z-range [-inf/+inf]. To indicate no + limit on *min* or *max* only, specify a hyphen (-). Normally, any NaNs + encountered are simply skipped and not considered in the + range-decision. Append **+n** to consider a NaN to be outside the given + z-range. This means the new subset will be NaN-free. Alternatively, + append **+r** to consider NaNs to be within the data range. In this + case we stop shrinking the boundaries once a NaN is found [Default + simply skips NaNs when making the range decision]. Finally, if your + core subset grid is surrounded by rows and/or columns that are all + NaNs, append **+N** to strip off such columns before (optionally) + considering the range of the core subset for further reduction of the + area. + + Returns + ------- + ret: xarray.DataArray or None + Return type depends on whether the *outgrid* parameter is set: + + - xarray.DataArray if *outgrid* is not set + - None if *outgrid* is set (grid output will be stored in *outgrid*) + """ + kind = data_kind(grid) + + with GMTTempFile(suffix=".nc") as tmpfile: + with Session() as lib: + if kind == "file": + file_context = dummy_context(grid) + elif kind == "grid": + raise NotImplementedError( + "xarray.DataArray is not supported as the input grid yet!" + ) + # file_context = lib.virtualfile_from_grid(grid) + # See https://github.com/GenericMappingTools/gmt/pull/3532 + # for a feature request. + else: + raise GMTInvalidInput("Unrecognized data type: {}".format(type(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("grdcut", arg_str) + + if outgrid == tmpfile.name: # if user did not set outgrid, return DataArray + with xr.open_dataarray(outgrid) as dataarray: + result = dataarray.load() + else: + result = None # if user sets an outgrid, return None + + return result diff --git a/pygmt/tests/test_grdcut.py b/pygmt/tests/test_grdcut.py new file mode 100644 index 00000000000..dfd06edc919 --- /dev/null +++ b/pygmt/tests/test_grdcut.py @@ -0,0 +1,64 @@ +""" +Tests for grdcut +""" +import os +import numpy as np +import pytest +import xarray as xr + +from .. import grdcut, grdinfo +from ..datasets import load_earth_relief +from ..exceptions import GMTInvalidInput +from ..helpers import GMTTempFile + + +def test_grdcut_file_in_file_out(): + "grduct an input grid file, and output to a grid file" + with GMTTempFile(suffix=".nc") as tmpfile: + result = grdcut("@earth_relief_01d", outgrid=tmpfile.name, region="0/180/0/90") + assert result is None # return value is None + assert os.path.exists(path=tmpfile.name) # check that outgrid exists + result = grdinfo(tmpfile.name, C=True) + assert result == "0 180 0 90 -8592.5 5559 1 1 181 91\n" + + +def test_grdcut_file_in_dataarray_out(): + "grdcut an input grid file, and output as DataArray" + outgrid = grdcut("@earth_relief_01d", region="0/180/0/90") + assert isinstance(outgrid, xr.DataArray) + # check information of the output grid + # the '@earth_relief_01d' is in pixel registration, so the grid range is + # not exactly 0/180/0/90 + assert outgrid.coords["lat"].data.min() == 0.0 + assert outgrid.coords["lat"].data.max() == 90.0 + assert outgrid.coords["lon"].data.min() == 0.0 + assert outgrid.coords["lon"].data.max() == 180.0 + assert outgrid.data.min() == -8592.5 + assert outgrid.data.max() == 5559 + assert outgrid.sizes["lat"] == 91 + assert outgrid.sizes["lon"] == 181 + + +def test_grdcut_dataarray_in_file_out(): + "grdcut an input DataArray, and output to a grid file" + # Not supported yet. + # See https://github.com/GenericMappingTools/gmt/pull/3532 + + +def test_grdcut_dataarray_in_dataarray_out(): + "grdcut an input DataArray, and output to a grid file" + # Not supported yet. + # See https://github.com/GenericMappingTools/gmt/pull/3532 + + +def test_grdcut_dataarray_in_fail(): + "Make sure that grdcut fails correctly if DataArray is the input grid" + with pytest.raises(NotImplementedError): + grid = load_earth_relief() + grdcut(grid, region="0/180/0/90") + + +def test_grdcut_fails(): + "Check that grdcut fails correctly" + with pytest.raises(GMTInvalidInput): + grdcut(np.arange(10).reshape((5, 2)))