Add visibility analysis tutorial (r.viewshed) by Valyrian-Code · Pull Request #135 · OSGeo/grass-tutorials

Valyrian-Code · 2026-07-04T04:41:19Z

Adds a new intermediate tutorial on visibility analysis with
r.viewshed — a topic
not currently covered in the tutorial set.

What it covers

A single-observer viewshed, drawn over shaded relief, and measuring the visible area.
The effect of observer height — ground level (~6.4 km² visible) vs. a 40 m tower (~46.6 km²).
Limiting the search with max_distance (a 3 km radius).
Cumulative visibility from five sites summed with r.series, as a starting point for site optimization.

Details

Uses the standard North Carolina sample dataset (elevation DEM) and the grass.tools API (GRASS 8.5), consistent with the recent NumPy/Landlab tutorial.
All code was run against GRASS 8.5; the figures are real outputs of the workflow. The tutorial renders cleanly with Quarto (eval: false, following the other Python tutorials).
Categories: raster, terrain, visualization, intermediate, Python.
Cross-links to the Terrain/DEMs and Modeling Movement tutorials.

Feedback on scope, wording, or the choice of observer sites is very welcome.

New intermediate tutorial covering viewshed/visibility analysis: - a single-observer viewshed and measuring the visible area - the effect of observer height (ground vs 40 m tower) and max_distance - cumulative visibility from multiple sites with r.series Uses the North Carolina sample dataset and the grass.tools API (GRASS 8.5). Categories: raster, terrain, visualization, intermediate, Python.

Copilot

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

veroandreo · 2026-07-12T15:10:45Z

The example developed is pretty basic, and the images are not visually appealing. They are too big, the legend and scale bar overlap. I would like to read a more compelling story, i.e., the observer in the center of the map is easy, but why are they there? Perhaps you could pick a different dataset or a different area and still make the tutorial reproducible by using r.in.usgs or i.eodag, and place the observer and the towers in places that actually mean something, not just in the center of the region. Please have a look at https://github.com/ncsu-geoforall-lab/grass-gis-workshop-foss4gna-2023, esp. notebooks 2 and 3, which cover similar analysis.

Address review feedback: replace the synthetic center-of-map observer with a real story. The tutorial now uses three historic fire lookout towers in the Nantahala National Forest (Wayah Bald, Wesser Bald, Albert Mountain, verified coordinates on the Appalachian Trail) and downloads real USGS 1 arc-second elevation via the r.in.usgs addon into a new EPSG:32617 project, with a pinned tile for reproducibility. Content changes: - observer height: standing on the bald vs the 16 m cab (165 vs 277 km2) - Earth curvature with -c at a lookout's 32 km range, and why the effect is small in ridge-limited terrain - target_elevation as a rising smoke column (277 to 412 km2 detectable) - three-tower network summed with r.series, coverage reported as a table from r.report JSON output (44% of the region) - new worked section on siting a fourth tower: the best candidate recovers under 4% of the blind area because the gaps are valleys Figures regenerated at 600 px with labeled towers (d.text), shaded relief bases, non-overlapping legend and scale bar, and zoomed extents for the single-tower maps. Get-started link now points to the Python quick start; callouts are spread out; area estimation commands are shown in the text.

Valyrian-Code · 2026-07-12T16:56:12Z

Thanks @veroandreo, this was the right push. I reworked the tutorial from scratch around a real story: the historic fire lookout towers of the Nantahala National Forest (Wayah Bald, Wesser Bald, and Albert Mountain, all still standing along the Appalachian Trail near Franklin, NC). The terrain now comes from the USGS via the r.in.usgs addon into a fresh EPSG:32617 project, so it is fully reproducible: one pinned tile of about 58 MB, with a callout on previewing the download with -i first.

The observers now mean something. These towers were sited for visibility, so each section answers a question their builders actually faced:

Observer height: standing on the bald vs the 16 m cab (165 vs 277 km² visible), which is why the towers were built.
Earth curvature with -c at a lookout's 32 km range, including why the effect is small in ridge-limited terrain.
target_elevation as a rising smoke column: detectable area grows from 277 to 412 km², since a lookout watches for smoke, not ground.
The three-tower network summed with r.series, with per-class coverage reported as a table from r.report JSON output (44% of the region).
A worked "where should a fourth tower go?" section: the best candidate summit recovers under 4% of the blind area, because the remaining gaps are valley bottoms that no ridgetop tower can see. That diminishing-returns result is, I think, the most interesting number in the tutorial.

On the images: regenerated at 600 px with labeled towers, shaded relief bases, non-overlapping legend and scale bar, and zoomed extents for the single-tower maps. The get-started link now points to the Python quick start, the callouts are spread out, and the area estimation commands are shown in the text.

Happy to adjust anything else.

veroandreo

This version reads much better and it is more interesting too. I leave some comments and suggestions here and there before merging.

veroandreo · 2026-07-15T12:55:24Z

Legend is not readable, please adjust the size

veroandreo · 2026-07-15T13:02:09Z

In all cases, labels overlap the tower symbol, try moving them a little bit to the right so it looks better.

veroandreo · 2026-07-15T13:02:51Z

+  Use r.viewshed and real USGS terrain to explore why fire lookout towers stand
+  where they do, how far they could see, and how well a network of them covered


Suggested change

Use r.viewshed and real USGS terrain to explore why fire lookout towers stand

where they do, how far they could see, and how well a network of them covered

Use r.viewshed and USGS terrain data to explore why fire lookout towers stand

where they do, how far they can see, and how well a network of them cover

veroandreo · 2026-07-15T13:07:23Z

+```bash
+# One-time: create the project and install the add-on
+grass -c EPSG:32617 ~/grassdata/nantahala -e
+grass ~/grassdata/nantahala/PERMANENT --exec g.extension extension=r.in.usgs
+```


You can use python for this too; see the create_project() function and how it is used in other tutorials

veroandreo · 2026-07-15T13:16:47Z

+
+# Download and import the National Elevation Dataset for this region
+tools.r_in_usgs(product="ned", ned_dataset="ned1sec", ned_release="all",
+                title_filter="20220725", output_name="elevation",


Suggested change

title_filter="20220725", output_name="elevation",

output_name="elevation",

Better remove this so the default dataset is downloaded, even if in the future the data itself is different, it will always download something.

veroandreo · 2026-07-15T13:34:35Z

+    detection_maps.append(out)
+
+tools.r_series(input=detection_maps, output="coverage_sum", method="sum")
+tools.r_mapcalc(expression="coverage = int(coverage_sum)")


why do you need this?

veroandreo · 2026-07-15T13:35:34Z

+| 2 | 214.8 | 13.7% |
+| 3 | 61.1 | 3.9% |
+
+The three towers between them watch about 692 km², or 44% of this 1,560 km² region. Only


Suggested change

The three towers between them watch about 692 km², or 44% of this 1,560 km² region. Only

The three towers together allow to watch about 692 km², or 44% of this 1,560 km² region. Only

veroandreo · 2026-07-15T13:37:16Z

+test a tower there.
+
+```{python}
+# cells no tower currently sees, and their elevation


Suggested change

# cells no tower currently sees, and their elevation

# cells no tower currently see, and their elevation

veroandreo · 2026-07-15T13:39:00Z

+gs.write_command("r.colors", map="covered_now", rules="-", stdin="1 150:180:200")
+gs.write_command("r.colors", map="recovered", rules="-", stdin="1 215:25:28")
+gs.write_command("v.in.ascii", input="-", output="candidate", separator="pipe",
+                 columns="x double precision, y double precision", x=1, y=2,
+                 stdin=f"{cx:.0f}|{cy:.0f}")


I think this can also be translated into tools API

veroandreo · 2026-07-15T13:39:51Z

+
+# Summary
+
+You have used


Suggested change

You have used

We have used

- create the project in Python with create_project() and install the add-on with g.extension, instead of a shell block - read the tower coordinates (and the candidate site) with StringIO so v.in.ascii runs through the Tools API - report visible areas with r.report instead of a custom helper - redraw the figures: tower labels shifted clear of the markers, larger legends, smaller file sizes - wording fixes from the review; comment the integer coverage step The NED tile stays pinned with title_filter for reproducible numbers; the download note explains how to drop it to fetch the current data.

Valyrian-Code · 2026-07-15T19:59:29Z

Thanks @veroandreo, glad the story landed! Pushed a round addressing the review:

Project setup in Python: create_project() plus g.extension to install the add-on, instead of the shell block.
Tools API for input: the tower coordinates and the candidate site now come from a StringIO buffer, so v.in.ascii runs through the Tools API.
r.report for areas: dropped the custom function and report visible areas with r.report (and the JSON output for the coverage table).
Figures: tower labels shifted clear of the markers, legends enlarged for readability, and smaller file sizes.
Applied the wording suggestions and commented the integer-coverage step.

One point I did not apply, and wanted to flag with data: removing title_filter. When I drop it, r.in.usgs with ned_release=all patches four NED vintages for this footprint, which shifts the terrain by up to ±45 m locally and, more importantly, changes the headline result. The ground-vs-tower comparison drops from ~70% more visible area to ~15%, and because the multi-vintage patch order is not deterministic the numbers stop being reproducible. So I kept the single pinned vintage for stable, reproducible figures, and the download note explains that a reader can drop title_filter to always fetch whatever NED is current.

Happy to switch to the unpinned default if you'd still prefer that; just wanted you to see the tradeoff first.

petrasovaa · 2026-07-15T20:38:41Z

One point I did not apply, and wanted to flag with data: removing title_filter. When I drop it, r.in.usgs with ned_release=all patches four NED vintages for this footprint, which shifts the terrain by up to ±45 m locally and, more importantly, changes the headline result. The ground-vs-tower comparison drops from ~70% more visible area to ~15%, and because the multi-vintage patch order is not deterministic the numbers stop being reproducible. So I kept the single pinned vintage for stable, reproducible figures, and the download note explains that a reader can drop title_filter to always fetch whatever NED is current.

The suggestion was to use the default current, not all. Reproducibility is anyway a problem because USGS does not guarantee the specific historical dataset will be available. I've seen them changing the data quite a bit. 45 m difference is strange.

Per review: drop ned_release=all and title_filter and let r.in.usgs fetch the current 1 arc-second tile (the default). For this footprint the current tile is identical to the previously pinned one, so all figures and reported areas are unchanged; the download note now says the numbers may shift slightly if USGS republishes the data.

Valyrian-Code · 2026-07-15T20:53:05Z

You're right, and thanks for the nudge, I had it wrong. I switched to the default current release (dropped both ned_release=all and title_filter):

tools.r_in_usgs(product="ned", ned_dataset="ned1sec", output_name="elevation",
                memory=1000, nprocs=2)

The ±45 m difference was an artifact of ned_release=all patching several vintages together, not something current does. I checked: for this footprint the current tile is byte-identical to the one I had pinned (elevation difference is 0 everywhere), so a single clean tile downloads and all the figures and reported areas stay exactly as they are now, just without the fragile pin. The download note now mentions that the numbers may shift a little if USGS republishes the data, which matches what you have seen with their updates.

Copilot AI review requested due to automatic review settings July 4, 2026 04:41

Copilot AI reviewed Jul 4, 2026

veroandreo reviewed Jul 12, 2026

View reviewed changes

veroandreo requested changes Jul 15, 2026

View reviewed changes

		Use r.viewshed and real USGS terrain to explore why fire lookout towers stand
		where they do, how far they could see, and how well a network of them covered

	title_filter="20220725", output_name="elevation",
	output_name="elevation",

	The three towers between them watch about 692 km², or 44% of this 1,560 km² region. Only
	The three towers together allow to watch about 692 km², or 44% of this 1,560 km² region. Only

	# cells no tower currently sees, and their elevation
	# cells no tower currently see, and their elevation

Uh oh!

Conversation

Valyrian-Code commented Jul 4, 2026

What it covers

Details

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Uh oh!

veroandreo commented Jul 12, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Valyrian-Code commented Jul 12, 2026

Uh oh!

veroandreo left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Valyrian-Code commented Jul 15, 2026

Uh oh!

petrasovaa commented Jul 15, 2026

Uh oh!

Valyrian-Code commented Jul 15, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants