This vignette provides a basic overview of the functions in R package
viewscape.
The basic viewshed analysis can be accessed through calling the
compute_viewshed. The two needed objects to compute the
viewshed are a digital surface model (DSM) and a viewpoint.
library(viewscape)#Load in DSM
test_dsm <- terra::rast(system.file("test_dsm.tif",
package ="viewscape"))
#Load in the viewpoints
# Load points (.shp file)
test_viewpoints <- sf::read_sf(system.file("test_viewpoints.shp",
package = "viewscape"))
test_viewpoint <- test_viewpoints[2,]
#Compute viewshed
output <- viewscape::compute_viewshed(dsm = test_dsm,
viewpoints = test_viewpoint,
offset_viewpoint = 3,
r = 800,
method = 'view_tree')# overlap viewshed on DSM
output_r <- viewscape::visualize_viewshed(output, outputtype = 'raster')
terra::plot(test_dsm, axes=FALSE, box=FALSE, legend = FALSE)
terra::plot(output_r, add=TRUE, col = "red", axes=FALSE, box=FALSE, legend = FALSE)
terra::plot(test_viewpoint, add = TRUE, col = "blue", axes=FALSE, box=FALSE, legend = FALSE)Angles are measured clockwise from north: 0 = north, 90 = east, −90 = west, ±180 = south.
# north-facing 120° arc (−60° to 60°)
sector <- viewscape::fov_mask(output, c(-60, 60))
terra::plot(test_dsm, axes=FALSE, box=FALSE, legend = FALSE)
terra::plot(viewscape::visualize_viewshed(sector, outputtype = 'raster'),
axes=FALSE, box=FALSE, legend = FALSE, add = TRUE, col = "red")
terra::plot(test_viewpoint, add = TRUE, col = "blue", axes=FALSE, box=FALSE, legend = FALSE)Visual Magnitude (VM) quantifies how much of the observer’s visual field each visible cell occupies. It combines three factors: the surface area of the cell (resolution²), the angle between the cell’s surface normal and the line of sight, and the squared distance from the observer. Cells that are large, face the observer directly, and are nearby have the highest visual magnitude. The method follows Chamberlain & Meitner (2013).
vm <- viewscape::visual_magnitude(output, test_dsm)Plot the raw VM raster. Higher values (warmer colours) indicate cells that dominate the visual field:
terra::plot(test_dsm, axes = FALSE, box = FALSE, legend = FALSE,
main = "Visual magnitude over DSM", col = gray.colors(256, start = 0.1, end = 0.9))
terra::plot(sqrt(vm), add=TRUE, alpha=0.8, axes=FALSE, box=FALSE)
terra::plot(test_viewpoint, add = TRUE, col = "red",
pch = 10, axes = FALSE, box = FALSE, legend = FALSE)Summarise the distribution of VM values across the viewshed:
vm_vals <- terra::values(vm, na.rm = TRUE)
summary(vm_vals)
# Total visual magnitude (sums to ~1 for a full hemisphere in theory)
sum(vm_vals)
# Identify the top 10% most visually dominant cells
threshold <- quantile(vm_vals, 0.90)
vm_top10 <- terra::classify(vm, cbind(0, threshold, NA)) # mask low values
terra::plot(test_dsm, axes = FALSE, box = FALSE, legend = FALSE,
main = "Top 10% most visually dominant cells",
col = gray.colors(256, start = 0.1, end = 0.9))
terra::plot(vm_top10, add = TRUE, col = "firebrick",
axes = FALSE, box = FALSE, legend = FALSE)
terra::plot(test_viewpoint, add = TRUE, col = "blue",
pch = 16, axes = FALSE, box = FALSE, legend = FALSE)VM can also be computed after subsetting the viewshed with
fov_mask() to restrict the analysis to a specific field of
view:
# Visual magnitude within a 120° north-east arc (0° to 120°)
sector <- viewscape::fov_mask(output, c(0, 120))
vm_sector <- viewscape::visual_magnitude(sector, test_dsm)
terra::plot(vm_sector, axes = FALSE, box = FALSE,
main = "Visual magnitude — 120° north-east arc")
terra::plot(test_viewpoint, add = TRUE, col = "blue",
pch = 16, axes = FALSE, box = FALSE, legend = FALSE)# Compute viewsheds
output <- viewscape::compute_viewshed(dsm = test_dsm,
viewpoints = test_viewpoints,
offset_viewpoint = 10,
parallel = TRUE,
workers = 1)# Use plot all viewsheds on DSM
par(mfrow=c(3,3))
for(i in 1:length(output)) {
each <- output[[i]]
raster_data <- viewscape::visualize_viewshed(each, outputtype="raster")
terra::plot(test_dsm, axes=FALSE, box=FALSE, legend = FALSE)
terra::plot(raster_data, add=TRUE, col = "red", axes=FALSE, box=FALSE, legend = FALSE)
}intervis_network() tests every pair of viewpoints for
mutual line of sight and returns the results either as an adjacency
matrix or as a spatial network of lines.
Each entry [i, j] is 1 if viewpoint
i can see viewpoint j, 0 if not,
and NA on the diagonal. Because observer heights can differ
between viewpoints, visibility is not always symmetric:
[i, j] may differ from [j, i].
# Binary N x N adjacency matrix
mat <- viewscape::intervis_network(
viewpoints = test_viewpoints,
dsm = test_dsm,
offset_viewpoint = 10,
r = 800
)
matSetting output = "lines" returns an sf data
frame of LINESTRING geometries — one per visible pair. The
mutual column is TRUE when both viewpoints can
see each other.
net <- viewscape::intervis_network(
viewpoints = test_viewpoints,
dsm = test_dsm,
offset_viewpoint = 10,
r = 800,
output = "lines"
)
netPlot the network over the DSM to map which viewpoints are intervisible:
terra::plot(test_dsm, axes = FALSE, box = FALSE, legend = FALSE,
main = "Intervisibility network")
# All visible links in grey
terra::plot(net["geometry"], add = TRUE, col = "grey60", lwd = 1)
# Mutually visible links in blue
terra::plot(net[net$mutual, "geometry"], add = TRUE, col = "#2171b5", lwd = 2)
# Viewpoints
terra::plot(test_viewpoints, add = TRUE, col = "red", pch = 16, cex = 1.2)Pass a vector of heights to offset_viewpoint to model
different observer elevations at each point — for example, observers on
rooftops of varying height:
n_vp <- nrow(sf::st_coordinates(test_viewpoints))
# cycle through a set of representative heights, one per viewpoint
heights <- rep(c(1.7, 6, 12), length.out = n_vp)
mat_heights <- viewscape::intervis_network(
viewpoints = test_viewpoints,
dsm = test_dsm,
offset_viewpoint = heights,
r = 1600
)
mat_heightsThe function of view depth analysis can calculate two different metrics: the furthest distance and standard deviation of distances. To calculate view depth, there are two needed objects: the DSM that was used to get viewshed and result from viewshed analysis.
The function of extent analysis can calculate the total area of viewshed and needs the DSM that was used to get viewshed and result from viewshed analysis.
The following function can calculate the area of ground surface and standard deviation of elevations within a viewshed. The function needs a DSM and a DEM/DTM to calculate the metrics.
#Load in DSM
test_dsm <- terra::rast(system.file("test_dsm.tif",
package ="viewscape"))
# Load DTM
test_dtm <- terra::rast(system.file("test_dtm.tif",
package ="viewscape"))
# Load canopy raster
test_canopy <- terra::rast(system.file("test_canopy.tif",
package ="viewscape"))
# Load building footprints raster
test_building <- terra::rast(system.file("test_building.tif",
package ="viewscape"))# calculate metrics given the viewshed, canopy, and building footprints
test_metrics <- viewscape::calculate_viewmetrics(output[[1]],
test_dsm,
test_dtm,
list(test_canopy, test_building))
test_metricscalculate_diversity() calculates the proportion of each type of land use/ cover within a viewshed to get the Shannon Diversity Index.
# load landuse raster
test_landuse <- terra::rast(system.file("test_landuse.tif",
package ="viewscape"))# the Shannon Diversity Index (SDI)
test_diversity <- viewscape::calculate_diversity(output[[1]],
test_landuse,
proportion = TRUE)
# SDI and The proportion of each type of land use
test_diversitycalculate_feature is to calculate the proportion of a feature (including trees, buildings, parking, and roads) within the viewshed. This function can be applied to
# calculate the percentage of canopy coverage
test_canopy_proportion <- viewscape::calculate_feature(viewshed = output[[1]],
feature = test_canopy,
type = 2,
exclude_value=0)
test_canopy_proportionpano_view() generates a 360° equirectangular (or
cylindrical) panorama from a viewpoint using raycasting through the DSM.
The output type depends on the inputs:
method |
semantic |
result |
|---|---|---|
"equirectangular" (default) |
NULL |
3-band RGB colour panorama (auto-downloads ESRI WorldImagery via greenSD) |
"equirectangular" |
land-cover raster | 2-layer: depth (m) + land-cover codes |
"cylindrical" |
NULL |
single-layer distance panorama |
"cylindrical" |
any raster | 2-layer: depth + semantic class |
When no semantic layer is supplied,
pano_view() automatically fetches ESRI WorldImagery tiles
for the DSM area via greenSD::get_tile_green() and projects
the satellite RGB colours onto the panoramic rays. Sky pixels are filled
with a configurable sky-blue colour.
# Requires the greenSD package:
# devtools::install_github("billbillbilly/greenSD")
# Default: equirectangular + ESRI WorldImagery → 3-band RGB panorama
pano_rgb <- viewscape::pano_view(
dsm = test_dsm,
vpt = test_viewpoint,
h = 3, # observer height above ground (m)
plot = TRUE
)
# pano_rgb is a 3-band SpatRaster (R, G, B, 0–255)
# Use plotRGB() to display it at any time:
terra::plotRGB(pano_rgb, r = 1, g = 2, b = 3)Change the sky colour or the panorama dimensions:
pano_custom <- viewscape::pano_view(
dsm = test_dsm,
vpt = test_viewpoint,
h = 6,
pano_dim = c(512, 1024), # taller panorama for more vertical detail
sky_color = c(200, 220, 255), # pale-blue sky
heading = 90, # face east
plot = TRUE
)When a semantic raster is supplied,
pano_view() encodes the land-cover or land-use code at each
ray-hit cell instead of the satellite colour. The result has two layers:
a depth layer (distance in metres to the first surface hit) and a
semantic layer (the code from the input raster).
test_landuse <- terra::rast(system.file("test_landuse.tif", package = "viewscape"))
pano_lc <- viewscape::pano_view(
dsm = test_dsm,
vpt = test_viewpoint,
h = 5,
semantic = test_landuse, # land-use raster with class codes
plot = TRUE # plots depth (left) and land-use codes (right)
)
# Access individual layers
depth_pano <- pano_lc[["equirectangular_depth"]]
lc_pano <- pano_lc[["equirectangular_semantic"]]
terra::plot(lc_pano, main = "Land use codes in panoramic view")Binary rasters (values 0/1 — e.g. building footprints, canopy masks) are detected automatically and treated as vertical extrusions: the semantic value fills the full column of voxels up to the DSM surface, matching the voxel-style projection of Lu et al. (2020).
test_building <- terra::rast(system.file("test_building.tif", package = "viewscape"))
pano_bld <- viewscape::pano_view(
dsm = test_dsm,
vpt = test_viewpoint,
h = 6,
semantic = test_building, # binary building-footprint mask
plot = TRUE
)test_building <- terra::rast(system.file("test_building.tif", package = "viewscape"))
pano_bld <- viewscape::pano_view(
dsm = test_dsm,
vpt = test_viewpoint,
h = 6,
semantic = test_canopy, # binary tree canopy mask
plot = TRUE
)