Skip to content
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
104 changes: 69 additions & 35 deletions doc/python_intro.md
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
---
authors:
- Corey T. White
- Vaclav Petras
- GRASS Development Team
---

Expand Down Expand Up @@ -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")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missing comment here.

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
Expand Down
7 changes: 5 additions & 2 deletions python/grass/tools/session_tools.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
''
"""
Expand Down
Loading