Skip to content

Commit 24b8be1

Browse files
committed
address review
1 parent f81ba8b commit 24b8be1

File tree

1 file changed

+38
-4
lines changed

1 file changed

+38
-4
lines changed

content/tutorials/numpy_integration/grass_numpy_integration.qmd

Lines changed: 38 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,9 @@ description: >
99
on an example using Landlab modeling framework.
1010
image: thumbnail.webp
1111
other-links:
12+
- text: "Intro to GRASS Python API"
13+
href: "https://grass.osgeo.org/grass-devel/manuals/python_intro.html"
14+
icon: link
1215
- text: "Check out new xarray-grass interface"
1316
href: "https://github.com/lrntct/xarray-grass"
1417
icon: link
@@ -55,11 +58,12 @@ import os
5558
import sys
5659
import subprocess
5760
61+
# for some environments, this path setup can be skipped
5862
sys.path.append(
5963
subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
6064
)
6165
62-
# import GRASS python packages
66+
# import GRASS Python packages
6367
import grass.script as gs
6468
import grass.jupyter as gj
6569
from grass.tools import Tools
@@ -80,9 +84,10 @@ tools = Tools()
8084
```
8185

8286
Since we are generating artificial data, we need to specify the dimensions (number of rows and columns).
83-
We also need to let GRASS know the actual coordinates; we will do that by setting
87+
We also need to let GRASS know the actual coordinates we want to use; we will do all that by setting
8488
the [computational region](https://grass.osgeo.org/grass-stable/manuals/g.region.html). Lower-left (south-west) corner of the data will be at the coordinates (0, 0),
8589
the coordinates of the upper-right (nort-east) corner are number of rows times cell resolution and number of columns times cell resolution.
90+
We chose a cell resolution of 10 horizontal units, which would impact certain terrain analyses, such as slope computation, and certain hydrology processes.
8691

8792
```{python}
8893
rows = 200
@@ -108,7 +113,7 @@ fractal = tools.r_surf_fractal(output=np.array, seed=6)
108113

109114
Directly passing NumPy arrays to GRASS tools and receiving them back is a new feature in GRASS v8.5.
110115
If you work with older versions of GRASS, you can use
111-
[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):
116+
[grass.script.array](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array):
112117

113118
```{python}
114119
import grass.script as gs
@@ -168,7 +173,7 @@ grid = RasterModelGrid((rows, cols), xy_spacing=resolution)
168173
grid.add_field("topographic__elevation", fractal.ravel(), at="node")
169174
```
170175

171-
Run the erosion of the landscape, extract the resulting NumPy array and visualize it:
176+
Run the erosion of the landscape:
172177

173178
```{python}
174179
from landlab.components import LinearDiffuser
@@ -182,7 +187,11 @@ for i in range(100):
182187
fa.run_one_step()
183188
fsp.run_one_step(dt=1.0)
184189
ld.run_one_step(dt=1.0)
190+
```
191+
192+
Extract the resulting NumPy array and visualize it:
185193

194+
```{python}
186195
elevation = grid.at_node['topographic__elevation'].reshape(grid.shape)
187196
188197
plt.figure()
@@ -233,6 +242,8 @@ to create a GRASS array object with the dimensions given by the current computat
233242
Then we copy the NumPy array and write the raster map as GRASS native raster.
234243

235244
```{python}
245+
import grass.script.array as ga
246+
236247
# Create a new grasss.script.array object
237248
grass_elevation = ga.array()
238249
# Copy values from elevation array
@@ -241,6 +252,29 @@ grass_elevation[...] = elevation
241252
grass_elevation.write("elevation")
242253
```
243254

255+
Let's visualize the new GRASS elevation map with streams by setting its color
256+
scheme to `srtm_percent` and rendering it in 3D using the
257+
[grass.jupyter.Map3D](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.jupyter.html?highlight=map3d#module-grass.jupyter.map3d) class:
258+
259+
```{python}
260+
tools.r_colors(map="elevation", color="srtm_percent")
261+
elevation_3dmap = gj.Map3D()
262+
elevation_3dmap.render(
263+
elevation_map="elevation",
264+
resolution_fine=1,
265+
vline="streams",
266+
vline_width=3,
267+
vline_color="blue",
268+
vline_height=3,
269+
position=[0.50, 0.00],
270+
height=5600,
271+
perspective=10,
272+
)
273+
elevation_3dmap.show()
274+
```
275+
276+
![Resulting elevation with streams in 3D](thumbnail.webp)
277+
244278
::: {.callout-note title="What about N-dimensional arrays?"}
245279

246280
For 3D arrays, you can use [grass.script.array3d](https://grass.osgeo.org/grass-stable/manuals/libpython/grass.script.html#script.array.array3d) package.

0 commit comments

Comments
 (0)