diff --git a/doc/python_intro.md b/doc/python_intro.md
index b3b35720d0d..c9ebbe654b3 100644
--- a/doc/python_intro.md
+++ b/doc/python_intro.md
@@ -1,6 +1,7 @@
 ---
 authors:
     - Corey T. White
+    - Vaclav Petras
     - GRASS Development Team
 ---
 
@@ -189,61 +190,94 @@ such as *text_split* function and *comma_items* attribute.
 
 ### NumPy interface
 
-The GRASS Python API includes a NumPy interface that allows you to read
-and write raster data as NumPy arrays.
-This makes it easy to integrate GRASS with the broader Python scientific
-stack for advanced analysis and custom modeling.
-Using *[grass.script.array](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array)*
-and *[grass.script.array3d](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array3d)*,
-you can switch between GRASS raster maps and NumPy arrays, run GRASS tools,
-and perform array-based operations as needed.
-It works for rasters as well as for 3D rasters.
+The *grass.tools* API allows tools to accept and output NumPy arrays instead of
+rasters. Although using GRASS native rasters is faster, NumPy arrays allow for
+easy integration with other tools in the Python scientific stack for advanced
+analysis and custom modeling. Additional APIs are available to switch between
+GRASS raster maps and NumPy arrays, run GRASS tools, and perform array-based
+operations.
 
-This example shows a workflow for writing a NumPy array
-to a GRASS raster, running a GRASS tool, and loading the result as a NumPy array:
+This example shows a workflow which starts with a NumPy array, runs a GRASS tool,
+and finishes again with a NumPy array:
 
 ```python
 import numpy as np
+import grass.script as gs
 from grass.tools import Tools
-from grass.script import array as garray
 
-# Create a 100x100 sinusoidal elevation surface
-xx, yy = np.meshgrid(np.linspace(0, 1, 100), np.linspace(0, 1, 100))
-elevation_array = np.sin(xx) + np.cos(yy)
+# Create a 100x100 elevation surface
+xx, yy = np.meshgrid(np.linspace(0, 1, 100), np.linspace(-1, 1, 100))
+elevation = xx * np.exp(-8 * np.abs(yy))
 
-# Set the region to match the array dimensions and resolution
-tools = Tools()
-tools.g_region(n=elevation_array.shape[0], s=0,
-               e=elevation_array.shape[1], w=0, res=1)
+gs.create_project("xy_project")
+session = gs.setup.init("xy_project")
 
-# Write the NumPy array to a new GRASS raster map
-map2d = garray.array()
-map2d[:] = elevation_array
-map2d.write("elevation", overwrite=True)
+# Set the region to match the array dimensions and resolution
+tools = Tools(session=session)
+rows = elevation.shape[0]
+cols = elevation.shape[1]
+tools.g_region(s=0, n=rows, w=0, e=cols, res=1)
 
-# Compute e.g., flow accumulation
-tools.r_watershed(elevation="elevation", accumulation="accumulation")
+# Use the NumPy array with a GRASS tool to compute, e.g., flow accumulation,
+# and get the result as a new NumPy array.
+accumulation = tools.r_watershed(elevation=elevation, accumulation=np.array, flags="a")
 
-# Load as numpy array
-accumulation_array = garray.array("accumulation")
+accumulation.max()
 ```
 
-This example demonstrates reading an existing GRASS raster into a NumPy array,
-modifying the array, and writing the modified array back to a GRASS raster:
+An additional APIs, *[grass.script.array](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array)*
+and *[grass.script.array3d](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array3d)*,
+provide you with full control over read and writing data
+between GRASS raster maps and NumPy arrays, including 3D raster maps.
+The following example demonstrates reading an existing GRASS raster into a NumPy
+array, plotting the data, modifying the array, and writing the modified array
+back to a GRASS raster:
 
 ```python
 import numpy as np
 import seaborn as sns
-from grass.script import array as garray
+import matplotlib.pyplot as plt
+import grass.script.array as ga
+from grass.tools import Tools
+
+# Set computational region to match the whole raster
+# This is assuming we are already using a GRASS project and that the project
+# contains a raster called 'elevation'
+tools = Tools()
+tools.g_region(raster="elevation")
 
-# Read elevation as numpy array
-elev = garray.array(mapname="elevation")
+# Read elevation as a NumPy array
+elevation = ga.array(mapname="elevation")
 # Plot raster histogram
-sns.histplot(data=elev.ravel(), kde=True)
+sns.histplot(data=elevation.ravel(), kde=True)
+plt.show()
 # Modify values
-elev_2 *= 2
+elevation *= 2
 # Write modified array into a GRASS raster
-elev_2.write(mapname="elevation_2")
+elevation.write(mapname="elevation_2")
+```
+
+Finally, this example shows a workflow for creating a NumPy array
+from scratch and writing it as a GRASS raster:
+
+```python
+import numpy as np
+from grass.tools import Tools
+import grass.script.array as ga
+
+# Create an array filled with values
+x = np.full((2, 3), 10)
+
+# Set the region to match the array dimensions and resolution
+tools = Tools()
+rows = elevation.shape[0]
+cols = elevation.shape[1]
+tools.g_region(s=0, n=rows, w=0, e=cols, res=1)
+
+# Write the NumPy array to a new GRASS raster map
+map2d = ga.array()
+map2d[:] = x
+map2d.write("x")  # add overwrite=True for repeated runs
 ```
 
 ## Fine-grained data handling
diff --git a/python/grass/tools/session_tools.py b/python/grass/tools/session_tools.py
index fc1f30e6d02..e1c6f161ad9 100644
--- a/python/grass/tools/session_tools.py
+++ b/python/grass/tools/session_tools.py
@@ -75,8 +75,10 @@ class Tools:
     Raster input and outputs can be NumPy arrays:
 
     >>> import numpy as np
-    >>> tools.g_region(rows=2, cols=3)
-    >>> slope = tools.r_slope_aspect(elevation=np.ones((2, 3)), slope=np.ndarray)
+    >>> rows = 2
+    >>> cols = 3
+    >>> tools.g_region(s=0, n=rows, w=0, e=cols, res=1)
+    >>> slope = tools.r_slope_aspect(elevation=np.ones((rows, cols)), slope=np.ndarray)
     >>> tools.r_grow(
     ...     input=np.array([[1, np.nan, np.nan], [np.nan, np.nan, np.nan]]),
     ...     radius=1.5,
@@ -113,6 +115,7 @@ class Tools:
     and text outputs from the tool as the result object has the same
     attributes and functionality as without arrays:
 
+    >>> # In this case, the tool did not produce any text, so the text is empty.
     >>> result.text
     ''
     """