diff --git a/docs/source/reference/geotiff.rst b/docs/source/reference/geotiff.rst
index e3eab382f..23be3d85e 100644
--- a/docs/source/reference/geotiff.rst
+++ b/docs/source/reference/geotiff.rst
@@ -216,7 +216,9 @@ Reading
 =======
 ``open_geotiff`` is the single read entry point. The backend follows the
 parameters: ``gpu=True`` returns a CuPy-backed array, ``chunks=N`` returns a
-lazy dask array, and a ``.vrt`` source reads a mosaic.
+lazy dask array, and a ``.vrt`` source reads a mosaic. It is re-exported at
+the top level, so ``from xrspatial import open_geotiff`` and
+``from xrspatial.geotiff import open_geotiff`` both work.
 
 .. autosummary::
     :toctree: _autosummary
@@ -228,7 +230,8 @@ Writing
 ``to_geotiff`` is the single write entry point (``gpu=True`` or CuPy data
 selects the GPU path; a ``.vrt`` output path writes tiles plus an index).
 Writing to a ``.vrt`` path is how you produce a VRT mosaic; the underlying
-index emitter is internal.
+index emitter is internal. Like ``open_geotiff``, it is re-exported at the
+top level, so ``from xrspatial import to_geotiff`` also works.
 
 .. autosummary::
     :toctree: _autosummary
diff --git a/xrspatial/__init__.py b/xrspatial/__init__.py
index 106459f5f..506f0bb1b 100644
--- a/xrspatial/__init__.py
+++ b/xrspatial/__init__.py
@@ -141,6 +141,8 @@
 from xrspatial.zonal import stats as zonal_stats  # noqa
 from xrspatial.zonal import hypsometric_integral  # noqa
 from xrspatial.zonal import suggest_zonal_canvas as suggest_zonal_canvas  # noqa
+from xrspatial.geotiff import open_geotiff  # noqa
+from xrspatial.geotiff import to_geotiff  # noqa
 from xrspatial.reproject import merge  # noqa
 from xrspatial.reproject import reproject  # noqa
 from xrspatial.utils import rechunk_no_shuffle  # noqa
diff --git a/xrspatial/geotiff/tests/parity/test_api_consolidation.py b/xrspatial/geotiff/tests/parity/test_api_consolidation.py
index f296494d8..e1208d7c2 100644
--- a/xrspatial/geotiff/tests/parity/test_api_consolidation.py
+++ b/xrspatial/geotiff/tests/parity/test_api_consolidation.py
@@ -128,3 +128,16 @@ def test_to_geotiff_gpu_matches_direct_backend(tmp_path):
     _write_geotiff_gpu(arr, direct_path, compression="deflate")
     np.testing.assert_array_equal(
         open_geotiff(via_path).values, open_geotiff(direct_path).values)
+
+
+def test_top_level_reexports_are_the_subpackage_functions():
+    """``from xrspatial import open_geotiff, to_geotiff`` (issue #3005).
+
+    The two public entry points are re-exported from the top-level
+    package so they import the same way as every other public function.
+    Both spellings must resolve to the same object.
+    """
+    import xrspatial
+
+    assert xrspatial.open_geotiff is open_geotiff
+    assert xrspatial.to_geotiff is to_geotiff