Repository: GenericMappingTools/pygmt

Issue:
pygmt.grdmask: Create mask grid from polygons or point coverage
*This issue serves as the central place for discussing and tracking the implementation of the `grdmask` method in PyGMT. The issue will be closed when the initial implementation is complete. Progress is tracked at https://github.com/orgs/GenericMappingTools/projects/3.*

## Documentation

- GMT: https://docs.generic-mapping-tools.org/dev/grdmask.html
- GMT.jl: https://www.generic-mapping-tools.org/GMTjl_doc/documentation/modules/grdmask
- PyGMT: https://www.pygmt.org/dev/api/generated/pygmt.grdmask.html

## GMT Option Flags and Modifiers

☑️: *Implemented*; ⬜: *To be implemented/discussed*; ~~Strikethrough~~: *Won't implement*.

- [ ] `-A`: `steps` - Steps for polygon path construction
- [ ] `-C`: `clobber` - Clobber mode for selecting polygons (f|l|o|u)
- [x] `-G`: `outgrid` #4463
- [x] `-I`: `spacing` #4463
- [x] `-N`: `outside`/`edge`/`inside`/`id_start` #4463
- [x] `-R`: `region` #4463
- [ ] `-S`: `search_radius` - Search radius for point coverage mode
- [x] `-V`: `verbose` #4463
- [ ] `-r`: `registration` - Grid registration (gridline or pixel)
- [ ] `-x`: `cores` - Number of cores for multi-threading
- [ ] `-a`: `aspatial` - Aspatial column associations
- [ ] `-bi`: `binary_in` - Binary input format
- [ ] `-di`: `nodata_in` - Replace nodata values with NaN
- [ ] `-e`: `pattern` - Data record pattern matching
- [ ] `-f`: `coltypes` - Data types of input/output columns
- [ ] `-g`: `gaps` - Data gaps and line breaks
- [ ] `-h`: `header` - Header record handling
- [ ] `-i`: `incols` - Input column selection and transformations
- [ ] `-j`: `distcalc` - Distance calculation method
- [ ] `-n`: `interpolation` - Interpolation settings
- [ ] `-q`: `inrows` - Input row selection
- [ ] `-w`: `wrap` - Cyclic coordinate conversion
- [ ] ~~`--PAR=value`~~: Use `pygmt.config` instead

## Related GMT configurations

- `PROJ_ELLIPSOID`: Affects distance calculations for geographic coordinates
- `PROJ_GEODESIC`: Method used for geodesic calculations
- `PROJ_MEAN_RADIUS`: Used for spherical distance calculations

## Linked Pull Requests

- [x] Initial feature implementation(#4463)
- [x] Add unit tests(#4463)
- [x] Add documentation examples(#4463)

## Related Issues and Discussions

## Additional acceptance criteria
- `pygmt.grdmask` wraps GMT `grdmask` for polygon masking with `data`, `region`, `spacing`, `outgrid`, `outside`, `edge`, `inside`, `id_start`, and `verbose` arguments.
- If `outgrid` is provided, the grid is written to that file and the function returns `None`; otherwise it returns an `xarray.DataArray` grid.
- The returned or written grid preserves standard PyGMT grid metadata, dimensions, registration, and GMT mask semantics for outside, edge, inside, and polygon ID values.
- Calling `grdmask` without the required region/spacing information raises the standard PyGMT parameter error.


Discussion hints:

Solution patch:
diff --git a/doc/api/index.rst b/doc/api/index.rst
index 5824dc328e1..5f67d019399 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -126,6 +126,7 @@ Operations on tabular data
     blockmedian
     blockmode
     filter1d
+    grdmask
     nearneighbor
     project
     select
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index f10c4fad8a5..8adecdf429a 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -45,6 +45,7 @@
     grdhisteq,
     grdinfo,
     grdlandmask,
+    grdmask,
     grdpaste,
     grdproject,
     grdsample,
diff --git a/pygmt/src/__init__.py b/pygmt/src/__init__.py
index 3c3fb2c39e3..b3a41b3f637 100644
--- a/pygmt/src/__init__.py
+++ b/pygmt/src/__init__.py
@@ -25,6 +25,7 @@
 from pygmt.src.grdimage import grdimage
 from pygmt.src.grdinfo import grdinfo
 from pygmt.src.grdlandmask import grdlandmask
+from pygmt.src.grdmask import grdmask
 from pygmt.src.grdpaste import grdpaste
 from pygmt.src.grdproject import grdproject
 from pygmt.src.grdsample import grdsample
diff --git a/pygmt/src/grdmask.py b/pygmt/src/grdmask.py
new file mode 100644
index 00000000000..d1ad20e6bb8
--- /dev/null
+++ b/pygmt/src/grdmask.py
@@ -0,0 +1,229 @@
+"""
+grdmask - Create mask grid from polygons or point coverage.
+"""
+
+from collections.abc import Sequence
+from typing import Literal
+
+import xarray as xr
+from pygmt._typing import PathLike
+from pygmt.alias import Alias, AliasSystem
+from pygmt.clib import Session
+from pygmt.exceptions import GMTParameterError, GMTValueError
+from pygmt.helpers import build_arg_list, fmt_docstring
+
+__doctest_skip__ = ["grdmask"]
+
+
+def _alias_option_N(  # noqa: N802
+    outside: float | None = None,
+    edge: float | Literal["z", "id"] | None = None,
+    inside: float | Literal["z", "id"] | None = None,
+    id_start: float | None = None,
+) -> Alias:
+    """
+    Return an Alias object for the -N option.
+
+    Builds the -N parameter string for grdmask based on the inside, edge, and outside
+    values. Handles special modes "z" (use z-value from polygon data) and "id" (use
+    running polygon ID).
+
+    Examples
+    --------
+    >>> def parse(**kwargs):
+    ...     return AliasSystem(N=_alias_option_N(**kwargs)).get("N")
+    >>> parse()
+    >>> parse(outside=1, edge=2, inside=3)
+    '1/2/3'
+    >>> parse(outside=3)
+    '3/0/1'
+    >>> parse(inside="z")
+    'z'
+    >>> parse(outside=1, inside="z")
+    'z/1'
+    >>> parse(edge="z", inside="z")
+    'Z'
+    >>> parse(inside="id")
+    'p'
+    >>> parse(edge="id", inside="id")
+    'P'
+    >>> parse(inside="id", id_start=5)
+    'p5'
+    >>> parse(edge="id", inside="id", id_start=10)
+    'P10'
+    >>> parse(edge="id", inside="id", id_start=5, outside=3)
+    'P5/3'
+    >>> parse(edge="id", id_start=5, outside=3)
+    Traceback (most recent call last):
+        ...
+    pygmt.exceptions.GMTValueError: ...
+    >>> parse(edge="z")
+    Traceback (most recent call last):
+        ...
+    pygmt.exceptions.GMTValueError: ...
+    >>> parse(inside="z", edge="id")
+    Traceback (most recent call last):
+        ...
+    pygmt.exceptions.GMTValueError: ...
+    >>> parse(inside="z", id_start=5)
+    Traceback (most recent call last):
+        ...
+    pygmt.exceptions.GMTValueError: ...
+    """
+    _inside_modes = {"z": "z", "id": "p"}
+
+    if id_start is not None and inside != "id":
+        raise GMTValueError(
+            inside,
+            description="value for parameter 'inside'",
+            reason="Parameter 'id_start' requires inside='id'.",
+        )
+
+    # In the special modes, 'edge' must be None or the same as 'inside'
+    if (edge in _inside_modes or inside in _inside_modes) and edge not in {
+        None,
+        inside,
+    }:
+        raise GMTValueError(
+            edge,
+            description="edge",
+            reason=f"inside={inside!r} and edge={edge!r} must be the same.",
+        )
+
+    # outside/edge/inside are all omitted: keep GMT default 0/0/1
+    if all(v is None for v in (outside, inside, edge)):
+        return Alias(None, name="mask_values")
+    # Build -N argument
+    if inside in _inside_modes:  # Mode: -Nz, -NZ, -Np, or -NP
+        mode = "z" if inside == "z" else "p"
+        if edge == inside:
+            mode = mode.upper()
+        # Append id_start if specified (only valid for "id" mode)
+        if id_start:
+            mode = f"{mode}{id_start}"
+        mask_values = mode if outside is None else [mode, outside]
+    else:  # Build the full mask with defaults for any missing values.
+        mask_values = [
+            0 if outside is None else outside,
+            0 if edge is None else edge,
+            1 if inside is None else inside,
+        ]
+    return Alias(mask_values, name="mask_values", sep="/", size=(2, 3))
+
+
+@fmt_docstring
+def grdmask(
+    data,
+    outgrid: PathLike | None = None,
+    spacing: Sequence[float | str] | None = None,
+    region: Sequence[float | str] | str | None = None,
+    outside: float | None = None,
+    edge: float | Literal["z", "id"] | None = None,
+    inside: float | Literal["z", "id"] | None = None,
+    id_start: float | None = None,
+    verbose: Literal["quiet", "error", "warning", "timing", "info", "compat", "debug"]
+    | bool = False,
+    **kwargs,
+) -> xr.DataArray | None:
+    """
+    Create mask grid from polygons or point coverage.
+
+    Reads one or more files containing polygon or data point coordinates, and creates a
+    grid where nodes that fall inside, on the edge, or outside the polygons (or within
+    the search radius from data points) are assigned values based on the ``outside``,
+    ``edge``, and ``inside`` parameters.
+
+    The mask grid can be used to mask out specific regions in other grids using
+    :func:`pygmt.grdmath` or similar tools. For masking based on coastline features,
+    consider using :func:`pygmt.grdlandmask` instead.
+
+    Full GMT docs at :gmt-docs:`grdmask.html`.
+
+    **Aliases**
+
+    .. hlist::
+       :columns: 3
+
+       - G = outgrid
+       - I = spacing
+       - N = outside, edge, inside, id_start
+       - R = region
+       - V = verbose
+
+    Parameters
+    ----------
+    data
+        Pass in either a file name to an ASCII data table, a 2-D $table_classes
+        containing the polygon(s) or data points. Input can be:
+
+        - **Polygon mode**: One or more files containing closed polygon coordinates
+        - **Point coverage mode**: Data points (used with ``search_radius`` parameter)
+    $outgrid
+    $spacing
+    outside
+    edge
+    inside
+        Set the value assigned to nodes outside, on the edge, or inside the polygons.
+        Can be any number, or one of ``None``, ``"NaN"``, and ``np.nan`` for NaN.
+
+        ``inside`` and ``edge`` can also be set to one of the following values:
+
+        - ``"z"``: Use the z-value from polygon data (segment header ``-Zzval``,
+          ``-Lheader``, or via ``-aZ=name``).
+        - ``"id"``: Use a running polygon ID number.
+
+        To treat edges as inside, use the same value as ``inside``.
+    id_start
+        The starting number for polygon IDs when ``inside="id"`` [Default is ``0``].
+        Only valid when ``inside="id"``.
+    $region
+    $verbose
+
+    Returns
+    -------
+    ret
+        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 the file set by
+          ``outgrid``)
+
+    Example
+    -------
+    >>> import pygmt
+    >>> import numpy as np
+    >>> # Create a simple polygon as a triangle
+    >>> polygon = np.array([[125, 30], [130, 30], [130, 35], [125, 30]])
+    >>> # Create a mask grid with 1 arc-degree spacing
+    >>> mask = pygmt.grdmask(data=polygon, spacing=1, region=[125, 130, 30, 35])
+    >>> mask.values
+    array([[0., 0., 0., 0., 0., 0.],
+           [0., 0., 1., 1., 1., 0.],
+           [0., 0., 0., 1., 1., 0.],
+           [0., 0., 0., 0., 1., 0.],
+           [0., 0., 0., 0., 0., 0.],
+           [0., 0., 0., 0., 0., 0.]], dtype=float32)
+    """
+    if kwargs.get("I", spacing) is None or kwargs.get("R", region) is None:
+        raise GMTParameterError(required=["region", "spacing"])
+
+    aliasdict = AliasSystem(
+        I=Alias(spacing, name="spacing", sep="/", size=2),
+        N=_alias_option_N(outside=outside, edge=edge, inside=inside, id_start=id_start),
+    ).add_common(
+        R=region,
+        V=verbose,
+    )
+    aliasdict.merge(kwargs)
+
+    with Session() as lib:
+        with (
+            lib.virtualfile_in(check_kind="vector", data=data) as vintbl,
+            lib.virtualfile_out(kind="grid", fname=outgrid) as voutgrd,
+        ):
+            aliasdict["G"] = voutgrd
+            lib.call_module(
+                module="grdmask",
+                args=build_arg_list(aliasdict, infile=vintbl),
+            )
+            return lib.virtualfile_to_raster(vfname=voutgrd, outgrid=outgrid)


Test patch:
diff --git a/pygmt/tests/test_grdmask.py b/pygmt/tests/test_grdmask.py
new file mode 100644
index 00000000000..4a9cc7d3c70
--- /dev/null
+++ b/pygmt/tests/test_grdmask.py
@@ -0,0 +1,170 @@
+"""
+Test pygmt.grdmask.
+"""
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+import xarray as xr
+from pygmt import grdmask
+from pygmt.enums import GridRegistration, GridType
+from pygmt.exceptions import GMTParameterError
+from pygmt.helpers import GMTTempFile
+
+
+@pytest.fixture(scope="module", name="polygon_data")
+def fixture_polygon_data():
+    """
+    Create a simple polygon for testing.
+    """
+    return np.array([[125, 30], [130, 30], [130, 35], [125, 30]])
+
+
+@pytest.fixture(scope="module", name="expected_grid")
+def fixture_expected_grid():
+    """
+    Load the expected grdmask grid result.
+    """
+    return xr.DataArray(
+        data=[
+            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 1.0, 1.0, 1.0, 0.0],
+            [0.0, 0.0, 0.0, 1.0, 1.0, 0.0],
+            [0.0, 0.0, 0.0, 0.0, 1.0, 0.0],
+            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+        ],
+        coords={
+            "x": [125.0, 126.0, 127.0, 128.0, 129.0, 130.0],
+            "y": [30.0, 31.0, 32.0, 33.0, 34.0, 35.0],
+        },
+        dims=["y", "x"],
+    )
+
+
+@pytest.fixture(scope="module", name="expected_grid_outside_only")
+def fixture_expected_grid_outside_only():
+    """
+    Load the expected grdmask grid result when only outside is set.
+    """
+    return xr.DataArray(
+        data=[
+            [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+            [3.0, 0.0, 1.0, 1.0, 1.0, 0.0],
+            [3.0, 3.0, 0.0, 1.0, 1.0, 0.0],
+            [3.0, 3.0, 3.0, 0.0, 1.0, 0.0],
+            [3.0, 3.0, 3.0, 3.0, 0.0, 0.0],
+            [3.0, 3.0, 3.0, 3.0, 3.0, 0.0],
+        ],
+        coords={
+            "x": [125.0, 126.0, 127.0, 128.0, 129.0, 130.0],
+            "y": [30.0, 31.0, 32.0, 33.0, 34.0, 35.0],
+        },
+        dims=["y", "x"],
+    )
+
+
+def test_grdmask_outgrid(polygon_data, expected_grid):
+    """
+    Creates a mask grid with an outgrid argument.
+    """
+    with GMTTempFile(suffix=".nc") as tmpfile:
+        result = grdmask(
+            data=polygon_data,
+            outgrid=tmpfile.name,
+            spacing=1,
+            region=[125, 130, 30, 35],
+        )
+        assert result is None  # return value is None
+        assert Path(tmpfile.name).stat().st_size > 0  # check that outgrid exists
+        temp_grid = xr.load_dataarray(tmpfile.name, engine="gmt", raster_kind="grid")
+        # Check grid properties
+        assert temp_grid.dims == ("y", "x")
+        assert temp_grid.gmt.gtype is GridType.CARTESIAN
+        assert temp_grid.gmt.registration is GridRegistration.GRIDLINE
+        # Check grid values
+        xr.testing.assert_allclose(a=temp_grid, b=expected_grid)
+
+
+@pytest.mark.benchmark
+def test_grdmask_no_outgrid(polygon_data, expected_grid):
+    """
+    Test grdmask with no set outgrid.
+    """
+    result = grdmask(data=polygon_data, spacing=1, region=[125, 130, 30, 35])
+    # Check grid properties
+    assert isinstance(result, xr.DataArray)
+    assert result.dims == ("y", "x")
+    assert result.gmt.gtype is GridType.CARTESIAN
+    assert result.gmt.registration is GridRegistration.GRIDLINE
+    # Check grid values
+    xr.testing.assert_allclose(a=result, b=expected_grid)
+
+
+def test_grdmask_custom_mask_values(polygon_data):
+    """
+    Test grdmask with custom outside, edge, inside values.
+    """
+    result = grdmask(
+        data=polygon_data,
+        spacing=1,
+        region=[125, 130, 30, 35],
+        outside=10,
+        edge=20,
+        inside=30,
+    )
+    assert isinstance(result, xr.DataArray)
+    # Check that the grid has the right dimensions
+    assert result.shape == (6, 6)
+    # Check that we have expected values
+    assert result.sel(x=125, y=35) == 10.0  # outside
+    assert result.sel(x=128, y=30) == 20.0  # edge
+    assert result.sel(x=129, y=31) == 30.0  # inside
+
+
+def test_grdmask_outside_only(polygon_data, expected_grid_outside_only):
+    """
+    Test grdmask when only outside is set.
+    """
+    result = grdmask(
+        data=polygon_data,
+        spacing=1,
+        region=[125, 130, 30, 35],
+        outside=3,
+    )
+
+    assert isinstance(result, xr.DataArray)
+    assert result.dims == ("y", "x")
+    assert result.gmt.gtype is GridType.CARTESIAN
+    assert result.gmt.registration is GridRegistration.GRIDLINE
+    xr.testing.assert_allclose(a=result, b=expected_grid_outside_only)
+
+
+def test_grdmask_fails():
+    """
+    Check that grdmask fails correctly when region and spacing are not given.
+    """
+    with pytest.raises(GMTParameterError):
+        grdmask(data=np.array([[0, 0], [1, 1], [1, 0], [0, 0]]))
+
+
+@pytest.mark.parametrize(
+    ("edge", "id_start"),
+    [(None, 5), ("id", 10), (None, -1), (None, 1.5)],
+)
+def test_grdmask_id_start_valid(polygon_data, edge, id_start):
+    """
+    Check that id_start works when inside='id'.
+    """
+    result = grdmask(
+        data=polygon_data,
+        spacing=1,
+        region=[125, 130, 30, 35],
+        edge=edge,
+        inside="id",
+        id_start=id_start,
+    )
+    assert isinstance(result, xr.DataArray)
+    unique_values = np.unique(result.values)
+    assert id_start in unique_values