diff --git a/doc/api/index.rst b/doc/api/index.rst
index 591784e82c9..75dedbf737d 100644
--- a/doc/api/index.rst
+++ b/doc/api/index.rst
@@ -174,14 +174,6 @@ Crossover analysis with x2sys
     x2sys_init
     x2sys_cross
 
-Input/output
-------------
-
-.. autosummary::
-    :toctree: generated
-
-    load_dataarray
-
 GMT Defaults
 ------------
 
diff --git a/pygmt/__init__.py b/pygmt/__init__.py
index 8adecdf429a..67c9c583d5a 100644
--- a/pygmt/__init__.py
+++ b/pygmt/__init__.py
@@ -24,7 +24,6 @@
 from pygmt import datasets
 from pygmt._show_versions import __commit__, __version__, show_versions
 from pygmt.figure import Figure, set_display
-from pygmt.io import load_dataarray
 from pygmt.session_management import begin as _begin
 from pygmt.session_management import end as _end
 from pygmt.src import (
diff --git a/pygmt/io.py b/pygmt/io.py
deleted file mode 100644
index 0e7c560ae21..00000000000
--- a/pygmt/io.py
+++ /dev/null
@@ -1,66 +0,0 @@
-"""
-PyGMT input/output (I/O) utilities.
-"""
-
-import warnings
-
-import xarray as xr
-
-
-# TODO(PyGMT>=0.20.0): Remove pygmt.io.load_dataarray
-def load_dataarray(filename_or_obj, **kwargs):
-    """
-    Open, load into memory, and close a DataArray from a file or file-like object
-    containing a single data variable.
-
-    This is a thin wrapper around :py:func:`xarray.open_dataarray`. It differs
-    from :py:func:`xarray.open_dataarray` in that it loads the DataArray into
-    memory, gets GMT specific metadata about the grid via
-    :py:meth:`GMTDataArrayAccessor`, closes the file, and returns the
-    DataArray. In contrast, :py:func:`xarray.open_dataarray` keeps the file
-    handle open and lazy loads its contents. All parameters are passed directly
-    to :py:func:`xarray.open_dataarray`. See that documentation for further
-    details.
-
-    .. deprecated:: v0.16.0
-       The 'pygmt.io.load_dataarray' function will be removed in v0.20.0. Please use
-       `xarray.load_dataarray(..., engine='gmt', raster_kind='grid')` instead if you
-       were reading grids using the engine='netcdf'; otherwise use `raster_kind='image'`
-       if you were reading multi-band images using engine='rasterio'.
-
-    Parameters
-    ----------
-    filename_or_obj : str or pathlib.Path or file-like or DataStore
-        Strings and Path objects are interpreted as a path to a netCDF file
-        or an OpenDAP URL and opened with python-netCDF4, unless the filename
-        ends with .gz, in which case the file is gunzipped and opened with
-        scipy.io.netcdf (only netCDF3 supported). Byte-strings or file-like
-        objects are opened by scipy.io.netcdf (netCDF3) or h5py (netCDF4/HDF).
-
-    Returns
-    -------
-    datarray : xarray.DataArray
-        The newly created DataArray.
-
-    See Also
-    --------
-    xarray.open_dataarray
-    """
-    msg = (
-        "The 'pygmt.io.load_dataarray' function will be removed in v0.20.0. Please use "
-        "`xarray.load_dataarray(..., engine='gmt', raster_kind='grid')` instead if you "
-        "were reading grids using the engine='netcdf'; otherwise use "
-        "`raster_kind='image'` if you were reading multi-band images using "
-        "engine='rasterio'."
-    )
-    warnings.warn(message=msg, category=FutureWarning, stacklevel=1)
-
-    if "cache" in kwargs:
-        msg = "'cache' has no effect in this context."
-        raise TypeError(msg)
-
-    with xr.open_dataarray(filename_or_obj, **kwargs) as dataarray:
-        result = dataarray.load()
-        _ = result.gmt  # load GMTDataArray accessor information
-
-    return result
diff --git a/pygmt/tests/test_io.py b/pygmt/tests/test_io.py
deleted file mode 100644
index 24829ff2ae0..00000000000
--- a/pygmt/tests/test_io.py
+++ /dev/null
@@ -1,42 +0,0 @@
-"""
-Test input/output (I/O) utilities.
-"""
-
-import numpy as np
-import pytest
-import xarray as xr
-from pygmt.enums import GridRegistration, GridType
-from pygmt.helpers import GMTTempFile
-from pygmt.io import load_dataarray
-
-pytest.importorskip("netCDF4")
-
-
-# TODO(PyGMT>=0.20.0): Remove test_io_load_dataarray
-def test_io_load_dataarray():
-    """
-    Check that load_dataarray works to read a netCDF grid with GMTDataArrayAccessor
-    information loaded.
-    """
-    with GMTTempFile(suffix=".nc") as tmpfile:
-        rng = np.random.default_rng()
-        grid = xr.DataArray(
-            data=rng.random((2, 2)), coords=[[0.1, 0.2], [0.3, 0.4]], dims=("x", "y")
-        )
-        grid.to_netcdf(tmpfile.name)
-
-        with pytest.warns(FutureWarning):
-            dataarray = load_dataarray(tmpfile.name)
-
-        assert dataarray.gmt.gtype is GridType.CARTESIAN
-        assert dataarray.gmt.registration is GridRegistration.PIXEL
-        # this would fail if we used xr.open_dataarray instead of load_dataarray
-        dataarray.to_netcdf(tmpfile.name)
-
-
-def test_io_load_dataarray_cache():
-    """
-    Check that load_dataarray fails when the cache argument is used.
-    """
-    with pytest.raises(TypeError):
-        _ = load_dataarray("somefile.nc", cache=True)