Skip to content
1 change: 1 addition & 0 deletions doc/api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -90,6 +90,7 @@ Operations on grids:

grdclip
grdcut
grdfill
grdfilter
grdtrack

Expand Down
1 change: 1 addition & 0 deletions pygmt/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
grd2cpt,
grdclip,
grdcut,
grdfill,
grdfilter,
grdinfo,
grdtrack,
Expand Down
1 change: 1 addition & 0 deletions pygmt/src/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
from pygmt.src.grdclip import grdclip
from pygmt.src.grdcontour import grdcontour
from pygmt.src.grdcut import grdcut
from pygmt.src.grdfill import grdfill
from pygmt.src.grdfilter import grdfilter
from pygmt.src.grdimage import grdimage
from pygmt.src.grdinfo import grdinfo
Expand Down
77 changes: 77 additions & 0 deletions pygmt/src/grdfill.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,77 @@
"""
grdfill - Fill blank areas from a grid.
"""

import xarray as xr
from pygmt.clib import Session
from pygmt.helpers import (
GMTTempFile,
build_arg_string,
fmt_docstring,
kwargs_to_strings,
use_alias,
)


@fmt_docstring
@use_alias(
A="mode",
G="outgrid",
R="region",
)
@kwargs_to_strings(R="sequence")
def grdfill(grid, **kwargs):
r"""
Fill blank areas from a grid file.

Read a grid that presumably has unfilled holes that the user
wants to fill in some fashion. Holes are identified by NaN values but
this criteria can be changed. There are several different algorithms that
can be used to replace the hole values.

Full option list at :gmt-docs:`grdfill.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.
mode : str
Specify the hole-filling algorithm to use. Choose from **c** for
constant fill and append the constant value, **n** for nearest
neighbor (and optionally append a search radius in
pixels [default radius is :math:`r^2 = \sqrt{{ X^2 + Y^2 }}`,
where (*X,Y*) are the node dimensions of the grid]).
{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``)
"""
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("grdfill", 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
49 changes: 49 additions & 0 deletions pygmt/tests/test_grdfill.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,49 @@
"""
Tests for grdclip.
"""
import os

import numpy as np
import pytest
import xarray as xr
from pygmt import grdfill, grdinfo
from pygmt.datasets import load_earth_relief
from pygmt.helpers import GMTTempFile


@pytest.fixture(scope="module", name="grid")
def fixture_grid():
"""
Load the grid data from the sample earth_relief file and set value(s) to
NaN.
"""
grid = load_earth_relief(registration="pixel", region=[-5, 5, -5, 5])
grid[3:5, 3:5] = np.nan
grid[5:7, 5:7] = np.inf
return grid


def test_grdfill_dataarray_out(grid):
"""
grdfill with a DataArray output.
"""
result = grdfill(grid=grid, mode="c20")
# check information of the output grid
assert isinstance(result, xr.DataArray)
assert result[4, 4] == 20
assert result[5, 5] == np.inf
assert not result.isnull().all() # check that no NaN values exists
assert result.gmt.gtype == 1 # Geographic grid
assert result.gmt.registration == 1 # Pixel registration


def test_grdfill_file_out(grid):
"""
grdfill with an outgrid set.
"""
with GMTTempFile(suffix=".nc") as tmpfile:
result = grdfill(grid=grid, mode="c20", 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()
assert result == "-5 5 -5 5 -5130.5 inf 1 1 10 10 1 1"