From 51898a50c1b512129a49c6291a0ff736035128cd Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 5 Sep 2025 09:14:27 -0400 Subject: [PATCH 1/3] docs: Show grass.tools with arrays before grass.script.array The main interface doc for Python now shows how GRASS tools are used with arrays through grass.tools before showing direct conversions between GRASS rasters and NumPy arrays with grass.script.array. --- doc/python_intro.md | 78 ++++++++++++++++++++++++++------------------- 1 file changed, 46 insertions(+), 32 deletions(-) diff --git a/doc/python_intro.md b/doc/python_intro.md index b3b35720d0d..9d70cbfc76d 100644 --- a/doc/python_intro.md +++ b/doc/python_intro.md @@ -189,61 +189,75 @@ 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 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) +elevation = np.sin(xx) + np.cos(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) - -# Write the NumPy array to a new GRASS raster map -map2d = garray.array() -map2d[:] = elevation_array -map2d.write("elevation", overwrite=True) +tools.g_region(n=elevation.shape[0], s=0, + e=elevation.shape[1], w=0, res=1) -# Compute e.g., flow accumulation -tools.r_watershed(elevation="elevation", accumulation="accumulation") - -# Load as numpy array -accumulation_array = garray.array("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) ``` -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, 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 grass.script.array as ga -# Read elevation as numpy array -elev = garray.array(mapname="elevation") +# Read elevation as a NumPy array +elev = ga.array(mapname="elevation") # Plot raster histogram sns.histplot(data=elev.ravel(), kde=True) # Modify values -elev_2 *= 2 +elev *= 2 # Write modified array into a GRASS raster -elev_2.write(mapname="elevation_2") +elev.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 a 100x100 sinusoidal elevation surface +x = np.full((2, 3), 10) + +# Set the region to match the array dimensions and resolution +tools = Tools() +tools.g_region(n=x.shape[0], s=0, e=x.shape[1], w=0, res=1) + +# Write the NumPy array to a new GRASS raster map +map2d = gs.array() +map2d[:] = x +map2d.write("x") # add overwrite=True for repeated runs ``` ## Fine-grained data handling From 82178412e9e9c8014ad79f9d6de02f62c772f1c8 Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 5 Sep 2025 09:54:47 -0400 Subject: [PATCH 2/3] Make the examples consistent and make them run --- doc/python_intro.md | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/doc/python_intro.md b/doc/python_intro.md index 9d70cbfc76d..d8b28460602 100644 --- a/doc/python_intro.md +++ b/doc/python_intro.md @@ -209,8 +209,7 @@ elevation = np.sin(xx) + np.cos(yy) # Set the region to match the array dimensions and resolution tools = Tools() -tools.g_region(n=elevation.shape[0], s=0, - e=elevation.shape[1], w=0, res=1) +tools.g_region(n=elevation.shape[0], s=0, e=elevation.shape[1], w=0, res=1) # Use the NumPy array with a GRASS tool to compute, e.g., flow accumulation, # and get the result as a new NumPy array. @@ -222,21 +221,29 @@ and *[grass.script.array3d](https://grass.osgeo.org/grass-stable/manuals/libpyth 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, modifying the array, and writing the modified array back to a GRASS raster: +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 +import matplotlib.pyplot as plt import grass.script.array as ga +from grass.tools import Tools + +# Set computational region to match the whole raster +tools = Tools() +tools.g_region(raster="elevation") # Read elevation as a NumPy array -elev = ga.array(mapname="elevation") +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 +elevation *= 2 # Write modified array into a GRASS raster -elev.write(mapname="elevation_2") +elevation.write(mapname="elevation_2") ``` Finally, this example shows a workflow for creating a NumPy array @@ -255,7 +262,7 @@ tools = Tools() tools.g_region(n=x.shape[0], s=0, e=x.shape[1], w=0, res=1) # Write the NumPy array to a new GRASS raster map -map2d = gs.array() +map2d = ga.array() map2d[:] = x map2d.write("x") # add overwrite=True for repeated runs ``` From 138534a47ac34a05d507096420c006a83333f72f Mon Sep 17 00:00:00 2001 From: Vaclav Petras Date: Fri, 5 Sep 2025 15:15:19 -0400 Subject: [PATCH 3/3] docs: Better use of computational region in NumPy examples --- doc/python_intro.md | 29 +++++++++++++++++++++-------- python/grass/tools/session_tools.py | 7 +++++-- 2 files changed, 26 insertions(+), 10 deletions(-) diff --git a/doc/python_intro.md b/doc/python_intro.md index d8b28460602..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 --- @@ -201,19 +202,27 @@ and finishes again with a NumPy array: ```python import numpy as np +import grass.script as gs from grass.tools import Tools -# Create a 100x100 sinusoidal elevation surface -xx, yy = np.meshgrid(np.linspace(0, 1, 100), np.linspace(0, 1, 100)) -elevation = 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)) + +gs.create_project("xy_project") +session = gs.setup.init("xy_project") # Set the region to match the array dimensions and resolution -tools = Tools() -tools.g_region(n=elevation.shape[0], s=0, e=elevation.shape[1], w=0, res=1) +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) # 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) +accumulation = tools.r_watershed(elevation=elevation, accumulation=np.array, flags="a") + +accumulation.max() ``` An additional APIs, *[grass.script.array](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array)* @@ -232,6 +241,8 @@ 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") @@ -254,12 +265,14 @@ import numpy as np from grass.tools import Tools import grass.script.array as ga -# Create a 100x100 sinusoidal elevation surface +# Create an array filled with values x = np.full((2, 3), 10) # Set the region to match the array dimensions and resolution tools = Tools() -tools.g_region(n=x.shape[0], s=0, e=x.shape[1], w=0, res=1) +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() 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 '' """