pygmt.grdmask: Improve docs for "inside", "outside" and "edge" parameters

[yvonnefroehlich]

Description of the problem

pygmt.grdmask allows to cut out polygons from a grid. Using inside= "NaN" sets the nodes inside the grid to NaN. However the nodes outside of the polygon are changed to zero (right map). This is not the case for using outside="NaN": Here the nodes outside the polygon are set to NaN and the nodes inside the polygon are unchanged (middle map).

Minimal Complete Verifiable Example

import numpy as np
import pygmt

region = [125, 135, 25, 36]
polygon = np.array([[129, 31], [134, 31],  [134, 35],  [129, 35],  [129, 31],])

grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s")

mask_out = pygmt.grdmask(region=region, data=polygon, spacing="30s", outside="NaN")
mask_in = pygmt.grdmask(region=region, data=polygon, spacing="30s", inside="NaN")
grid_mask_out = grid * mask_out
grid_mask_in = grid * mask_in

args_basemap = {"region": region, "projection": "M12c", "frame": True}


fig = pygmt.Figure()
pygmt.makecpt(cmap="oleron", series=[-2000, 2000])

# Plot the elevation grid
fig.basemap(**args_basemap)
fig.grdimage(grid=grid, cmap=True)

fig.shift_origin(xshift="+w+2c")

# Plot the masked elevation grid outside
fig.basemap(**args_basemap)
fig.grdimage(grid=grid_mask_out, cmap=True)

fig.shift_origin(xshift="+w+2c")

# Plot the masked elevation grid inside
fig.basemap(**args_basemap)
fig.grdimage(grid=grid_mask_in, cmap=True)

fig.colorbar(frame=True)
fig.show()

Full error message

No error message occurrs. For the right map, the elevation grid should remain unchanged in the area outside the square.

System information

pygmt.show_versions()
PyGMT information:
  version: v0.19.0.dev151
System information:
  python: 3.14.2 | packaged by conda-forge | (main, Dec  6 2025, 11:17:27) [MSC v.1944 64 bit (AMD64)]
  executable: C:\ProgramData\Anaconda3\envs\pygmt_env_dev\python.exe
  machine: Windows-10-10.0.19045-SP0
Dependency information:
  numpy: 2.4.1
  pandas: 2.3.3
  xarray: 2025.12.0
  packaging: 25.0
  contextily: 1.7.0
  geopandas: 1.1.2
  IPython: 9.9.0
  pyarrow: 22.0.0
  rioxarray: 0.20.0
  gdal: 3.12.1
  ghostscript: 10.06.0
GMT library information:
  version: 6.6.0
  padding: 2
  share dir: C:/Program Files (x86)/gmt6/share
  plugin dir: C:/ProgramData/Anaconda3/envs/pygmt_env_dev/Library/bin/gmt_plugins
  library path: C:/ProgramData/Anaconda3/envs/pygmt_env_dev/Library/bin/gmt.dll
  cores: 4
  grid layout: rows
  image layout: 
  binary version: 6.6.0

[seisman]

I think this is the expected behavior:

import numpy as np
import pygmt

region = [125, 135, 25, 36]
polygon = np.array([[129, 31], [134, 31],  [134, 35],  [129, 35],  [129, 31],])

grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s")

mask_out = pygmt.grdmask(region=region, data=polygon, spacing="1d", outside="NaN")
mask_in = pygmt.grdmask(region=region, data=polygon, spacing="1d", inside="NaN")

print(mask_out)

print(mask_in)

The output is:

<xarray.DataArray 'z' (lat: 12, lon: 11)> Size: 528B
array([[nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan],
       [nan, nan, nan, nan,  0., nan, nan, nan, nan,  0., nan],
       [nan, nan, nan, nan,  0.,  1.,  1.,  1.,  1.,  0., nan],
       [nan, nan, nan, nan,  0.,  1.,  1.,  1.,  1.,  0., nan],
       [nan, nan, nan, nan,  0.,  1.,  1.,  1.,  1.,  0., nan],
       [nan, nan, nan, nan,  0.,  1.,  1.,  1.,  1.,  0., nan],
       [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]],
      dtype=float32)
Coordinates:
  * lat      (lat) float64 96B 25.0 26.0 27.0 28.0 29.0 ... 33.0 34.0 35.0 36.0
  * lon      (lon) float64 88B 125.0 126.0 127.0 128.0 ... 133.0 134.0 135.0
Attributes:
    Conventions:   CF-1.7
    title:         
    history:       gmt grdmask @GMTAPI@-S-I-D-M-T-N-000000 -G@GMTAPI@-S-O-G-G...
    description:   
    actual_range:  [nan nan]
    long_name:     z
<xarray.DataArray 'z' (lat: 12, lon: 11)> Size: 528B
array([[ 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.,  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.,  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.,  0., nan, nan, nan, nan,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., nan, nan, nan, nan,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., nan, nan, nan, nan,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., nan, nan, nan, nan,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]],
      dtype=float32)
Coordinates:
  * lat      (lat) float64 96B 25.0 26.0 27.0 28.0 29.0 ... 33.0 34.0 35.0 36.0
  * lon      (lon) float64 88B 125.0 126.0 127.0 128.0 ... 133.0 134.0 135.0
Attributes:
    Conventions:   CF-1.7
    title:         
    history:       gmt grdmask @GMTAPI@-S-I-D-M-T-N-000000 -G@GMTAPI@-S-O-G-G...
    description:   
    actual_range:  [nan nan]
    long_name:     z

The key point is that, for the -N option, the argument is in the format of outside/edge/inside, which defaults to 0/0/1.

inside="NaN" is translated to -N0/0/NaN, and outside="NaN" is translated to -NNaN/0/1.

To achieve what you want, you need to set inside="NaN", outside=1, which gives:

[seisman]

I think we need to improve the documentation to clarify the default values.

Edit: Should be done in #4524.

[yvonnefroehlich]

Just extended my example exploring the outside parameter and found it has to be larger 0, then the outside nodes are not changed:

import numpy as np
import pygmt

region = [125, 135, 25, 36]
polygon = np.array([[129, 31], [134, 31],  [134, 35],  [129, 35],  [129, 31],])

grid = pygmt.datasets.load_earth_relief(region=region, resolution="30s")

mask_out = pygmt.grdmask(region=region, data=polygon, spacing="30s", outside="NaN")
mask_in = pygmt.grdmask(region=region, data=polygon, spacing="30s", inside="NaN")
mask_inout0 = pygmt.grdmask(region=region, data=polygon, spacing="30s", inside="NaN", outside=0)
mask_inout1 = pygmt.grdmask(region=region, data=polygon, spacing="30s", inside="NaN", outside=1)
grid_mask_out = grid * mask_out
grid_mask_in = grid * mask_in
grid_mask_inout0 = grid * mask_inout0
grid_mask_inout1 = grid * mask_inout1

args_basemap = {"region": region, "projection": "M12c", "frame": True}


fig = pygmt.Figure()
pygmt.makecpt(cmap="oleron", series=[-2000, 2000])

# Plot the elevation grid
fig.basemap(**args_basemap)
fig.grdimage(grid=grid, cmap=True)

fig.shift_origin(xshift="+w+2c")
0
# Plot the masked elevation grid outside=NaN
fig.basemap(**args_basemap)
fig.grdimage(grid=grid_mask_out, cmap=True)

fig.shift_origin(xshift="+w+2c")

# Plot the masked elevation grid inside=NaN
fig.basemap(**args_basemap)
fig.grdimage(grid=grid_mask_in, cmap=True)

fig.shift_origin(xshift="+w+2c")

# Plot the masked elevation grid inside=NaN, outside=0
fig.basemap(**args_basemap)
fig.grdimage(grid=grid_mask_inout0, cmap=True)

fig.shift_origin(xshift="+w+2c")

# Plot the masked elevation grid inside=NaN, outside=1
fig.basemap(**args_basemap)
fig.grdimage(grid=grid_mask_inout1, cmap=True)

fig.colorbar(frame=True)
fig.show()

[yvonnefroehlich]

The key point is that, for the -N option, the argument is in the format of outside/edge/inside, which defaults to 0/0/1.

Ahhh, I see 🙂.

I think we need to improve the documentation to clarify the default values.

Yes, we definitely need to specify this behaviour in more detail in the docs.