1.1 Compute single viewshed

#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)

1.2 Subset viewshed using by specifying the field of view

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)

1.3 Compute visual magnitude

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)

1.4 Compute the viewshed for multiple viewpoints

# 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)
}

1.5 Compute an intervisibility network

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.

# Binary N x N adjacency matrix
mat <- viewscape::intervis_network(
  viewpoints        = test_viewpoints,
  dsm               = test_dsm,
  offset_viewpoint  = 10,
  r                 = 800
)
mat

net <- viewscape::intervis_network(
  viewpoints        = test_viewpoints,
  dsm               = test_dsm,
  offset_viewpoint  = 10,
  r                 = 800,
  output            = "lines"
)
net

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)

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_heights

2.1 Calculate the metrics of viewshed

The 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_metrics

2.2 Calculate land use/cover diversity

calculate_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_diversity

2.3 calculate a single feature

calculate_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_proportion

3.1 Colour panorama (equirectangular, default)

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
)

3.2 Land-cover / land-use panorama (equirectangular + semantic)

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
)

3.3 Multi-layer semantic panorama

A single DSM often hides ground-level features: buildings, roads, sidewalks, or lawn under tree canopy are invisible to the raycaster because the canopy surface is what the DSM records. Passing a list of feature rasters solves this — each layer is extruded from the ground up to the DSM surface, so a building or sidewalk under canopy still occupies its column of space and is rendered in the panorama.

Earlier list elements take priority where layers overlap. Include whichever layers you have data for and omit the rest.

test_building <- terra::rast(system.file("test_building.tif", package = "viewscape"))
test_canopy   <- terra::rast(system.file("test_canopy.tif",   package = "viewscape"))

pano_multi <- viewscape::pano_view(
  dsm      = test_dsm,
  vpt      = test_viewpoint,
  h        = 6,
  semantic = list(
    buildings = test_building,   # highest priority
    # roads     = test_road,     # include if available
    # sidewalks = test_sidewalk, # include if available
    # lawn      = test_lawn,     # include if available
    canopy    = test_canopy      # lowest priority
  ),
  plot = TRUE
)
# semantic layer values: 0 = background, 1 = buildings, 2 = canopy, NA = sky
# with all five layers: 1 = buildings, 2 = roads, 3 = sidewalks, 4 = lawn, 5 = canopy

4.1 Preparing the layers

The understory functions need three aligned rasters: a DSM, a bare-earth DEM, and a binary canopy mask. Two things about the shipped layers are worth getting right, because both change the answer.

test_canopy.tif is fractional cover on a grid offset from the DSM, so it needs resampling before use. More subtly, a cell can have high fractional cover while the DSM there records the road surface rather than a crown — a crown edge overhanging a cell whose first return came from the ground. Treating such a cell as canopy puts a thin solid slab across the road, because the crown prism spans from the ground up to a canopy top that is barely above it. Require a real crown in the height model as well as cover.

This dataset’s CRS (EPSG:2253) is also in feet, while every height and distance argument is in metres and converted internally. The functions report the converted observer height in the CRS unit so the conversion can be checked.

test_dsm <- terra::rast(system.file("test_dsm.tif", package = "viewscape"))
test_dtm <- terra::rast(system.file("test_dtm.tif", package = "viewscape"))

canopy_frac   <- terra::rast(system.file("test_canopy.tif",   package = "viewscape"))
building_frac <- terra::rast(system.file("test_building.tif", package = "viewscape"))

chm <- test_dsm - test_dtm            # canopy height model, in feet here

# cover AND a crown that the height model actually sees
test_canopy_bin <- terra::ifel(
  terra::resample(canopy_frac, test_dsm, method = "bilinear") > 0.5 & chm > 10, 1, 0)
test_building_bin <- terra::ifel(
  terra::resample(building_frac, test_dsm, method = "bilinear") > 0.5 & chm > 10, 1, 0)

sf::st_crs(test_dsm)$units            # "ft" — arguments are still in metres

4.2 Finding a viewpoint on the road and under the canopy

The motivating case is a viewpoint snapped to a road that happens to sit beneath an overhanging crown. All three conditions are checkable from the shipped data: test_landuse.tif is parcel zoning in which class 0 is the road right-of-way, the canopy mask says there is a crown, and the CHM says how much of one.

landuse <- terra::resample(
  terra::rast(system.file("test_landuse.tif", package = "viewscape")),
  test_dsm, method = "near")

road <- terra::ifel(!is.na(landuse) & landuse == 0, 1, 0)

# road cells that sit under a crown
road_under_canopy <- terra::ifel(road > 0 & test_canopy_bin > 0, 1, 0)
terra::global(road_under_canopy, "sum", na.rm = TRUE)   # 5562 cells

# a point well inside the right-of-way with a crown overhead
vpt <- c(13293423, 282322.9)
terra::extract(chm, matrix(vpt, ncol = 2))              # 19 ft, about 5.8 m

4.3 The observer fix on its own

compute_viewshed() accepts canopy and dem. Canopy cells drop to bare earth so crowns stop obstructing, and — the important part — the observer is placed above bare earth wherever the viewpoint falls on canopy, while roofs and bridge decks keep their own elevation.

It is worth computing the bare-earth viewshed first as an upper bound, so the remaining reduction can be attributed to surface objects rather than to terrain.

# upper bound: what the terrain alone permits
terrain_only <- viewscape::compute_viewshed(
  dsm = test_dtm, viewpoints = vpt, offset_viewpoint = 1.7, r = 150)
sum(terrain_only@visible == 1)   # 154171 cells — the terrain is wide open

# the default: observer read off the DSM, i.e. standing on top of the crown
wrong <- viewscape::compute_viewshed(
  dsm = test_dsm, viewpoints = vpt, offset_viewpoint = 1.7, r = 150)
sum(wrong@visible == 1)          # 48281 cells, seen from above the treetops

# canopy cleared, observer on bare earth
right <- viewscape::compute_viewshed(
  dsm = test_dsm, viewpoints = vpt, offset_viewpoint = 1.7, r = 150,
  canopy = test_canopy_bin, dem = test_dtm)
sum(right@visible == 1)          # 600 cells, seen from 1.7 m under the canopy

4.4 The crown band: compute_viewshed_understory()

Clearing canopy to bare earth is all-or-nothing: a crown either blocks everything or nothing. In a 2.5D surface it has to be, because a cell obstructs everything below its elevation — writing dem + cbh into a canopy cell would wall off exactly the gap the observer looks through.

compute_viewshed_understory() drops that restriction by modelling each crown as a prism between its base and its top and testing the height of the sight line where it crosses the cell. A line that stays below the crown base passes underneath; one that has climbed above the crown top passes over. Two consequences follow that no single pitch threshold can reproduce: the slope at which a given crown starts to block scales as 1/distance, and visibility is not monotone in pitch.

Setting cbh = 0 makes the prism span the whole column, which reproduces an opaque canopy exactly — a useful reference point, and here a stark one.

opaque <- viewscape::compute_viewshed_understory(
  dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
  viewpoints = vpt, cbh = 0, r = 150)
sum(opaque$viewshed@visible == 1)     # 9 cells: walled in by the crowns

band <- viewscape::compute_viewshed_understory(
  dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
  viewpoints = vpt, cbh = 3, r = 150)
sum(band$viewshed@visible == 1)       # 612 cells: the street opens up

# band$viewshed is a Viewshed object, so the usual helpers still work
viewscape::visualize_viewshed(band$viewshed, plottype = "raster")

# band$class distinguishes what is seen:
#   0 invisible, 1 surface, 2 canopy only, 3 both surface and crown
terra::plot(band$class)

Nine cells against 612 is the whole argument for the method in one comparison. Both are far below the 154171 that the terrain alone allows, and that gap is real rather than an artefact: this DSM records a great deal of low vegetation and other surface clutter along the street, none of it tall enough to be classed canopy by the mask above, all of it genuinely opaque to an eye at 1.7 m.

Where the assumed crown base would reach the canopy top, the whole column is treated as foliage and the cell becomes opaque from the ground up, so hedges and shrubs need no special handling. Section 4.6 covers when that rule turns into a trap.

4.5 The hemispherical view: pano_view_understory()

The viewshed answers “is this cell visible”. Canopy exposure, sky view and any comparison against a segmented street-level photograph need the other question: “what does the observer see in this direction”. pano_view_understory() is the first-hit counterpart, and its classes are deliberately the ones a semantic segmentation gives you — sky, ground, built, canopy.

Fractions are weighted by solid angle, not counted in pixels: an equirectangular row near the zenith subtends cos(latitude) of a row at the horizon, so counting pixels badly overweights the poles. On open flat ground the function returns a sky view factor of exactly 1.

pv <- viewscape::pano_view_understory(
  dsm    = test_dsm,
  dem    = test_dtm,
  canopy = test_canopy_bin,
  built  = test_building_bin,
  vpt    = vpt,
  h      = 1.7,      # eye height in metres above bare earth
  cbh    = 3,
  max_dist = 120)

round(pv$metrics, 3)
#>              sky           ground            built           canopy
#>            0.043            0.600            0.003            0.354
#>              svf  canopy_overhead canopy_enclosure       mean_depth
#>            0.085            0.709            0.211           11.234

# six layers: depth, class, and the per-direction weights
names(pv$pano)
terra::plot(pv$pano[["class"]])

The ground fills 0.600 of the sphere — rather more than the 0.5 of an open hemisphere, because the observer is looking slightly down at a surface that curves away under the crowns — while canopy takes 0.354 and sky only 0.043. Note that canopy_enclosure (0.211) is much lower than canopy_overhead (0.709): enclosure is restricted to the band within 30° of the horizon, because that is where foliage acts as a wall rather than as a ceiling, and the two are different experiences of the same trees.

h defaults to 1.7 m here rather than the 6 m of pano_view(), and sky_threshold defaults to 0 rather than 3. That second default matters: the surface test counts a hit as soon as the surface comes within sky_threshold of the ray, so a 3 m tolerance with an eye at 1.7 m registers the ground at the first step in every direction and returns a panorama with no sky in it at all. The function warns if sky_threshold is not smaller than h.

4.6 Road corridors, and the cbh trap

An observer beneath a street tree is standing in an opening, and the crown over a roadway is pruned for vehicle clearance rather than hanging to the ground. Here the right-of-way is already available as land-use class 0, so it can go straight in as a clear_mask.

pv_road <- viewscape::pano_view_understory(
  dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
  built = test_building_bin, vpt = vpt, cbh = 3,
  clear_mask = road, max_dist = 120)

round(pv_road$metrics, 3)
#>              sky           ground            built           canopy
#>            0.305            0.604            0.003            0.088
#>              svf  canopy_overhead canopy_enclosure       mean_depth
#>            0.610            0.177            0.104           13.142

band_road <- viewscape::compute_viewshed_understory(
  dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
  viewpoints = vpt, cbh = 3, clear_mask = road, r = 150)
sum(band_road$viewshed@visible == 1)   # 617 cells, against 612 without

That single change is instructive. Ground visibility moves by 0.8%, from 612 cells to 617 — the sight lines that reach distant ground were already passing below the crown base. Sky view factor goes from 0.085 to 0.610, a sevenfold change. The corridor is a canopy-exposure decision, not a ground-visibility one, and it is worth reporting as such whenever both sets of metrics appear together.

It also over-corrects in the vertical. clear_mask deletes the crown prism outright over those cells, which is right at eye level and wrong at 6 m: it removes the arching branches, so straight up becomes sky. Where the crowns are tall enough, raising the crown base over the corridor with a cbh raster is the better instrument, since it keeps the canopy overhead while letting the street view through.

That is not available on this dataset, and the reason is a trap worth understanding. Crowns here top out around 19 ft, so an assumed crown base above that reaches the canopy top, the collapse rule fires, and those cells become opaque from the ground up. Visibility therefore rises with cbh and then falls:

for (m in c(2, 3, 4, 5)) {
  v <- viewscape::compute_viewshed_understory(
    dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
    viewpoints = vpt, cbh = m, r = 150)
  cat(sprintf("cbh = %d m: %d cells\n", m, sum(v$viewshed@visible == 1)))
}
#> cbh = 2 m: 470 cells
#> cbh = 3 m: 612 cells
#> cbh = 4 m: 172 cells     <- past the canopy top, cells turn opaque
#> cbh = 5 m: 57 cells

So cbh is only meaningful below the local canopy height. Check it against dsm - dem before choosing a value; the function warns when the collapse affects more than a tenth of the canopy cells. canopy_removal_r is the cruder version of the corridor idea — a radius around the viewpoint within which crowns contribute nothing — and being isotropic it also clears the flanking canopy where nothing should be cleared.

Most study areas have no right-of-way polygon to hand, but they all have OpenStreetMap centrelines. road_corridor_mask() buffers those into a corridor, and the query is driven straight off the DSM’s own footprint, so nothing has to be digitised. This needs the osmdata package (install.packages("osmdata")); OSM works in lon/lat, so the bounding box is reprojected to EPSG:4326 for the query and road_corridor_mask() transforms the result back.

# the DSM footprint, in lon/lat, as the Overpass query area
bb <- sf::st_bbox(sf::st_transform(sf::st_as_sfc(sf::st_bbox(test_dsm)), 4326))

osm <- osmdata::opq(bbox = as.numeric(bb)) |>
  osmdata::add_osm_feature(key = "highway") |>
  osmdata::osmdata_sf()

roads <- osm$osm_lines

# half-widths per OSM highway class, in metres, with a fallback so an unlisted
# class does not arrive as NA
hw <- c(motorway = 12, trunk = 11, primary = 9, secondary = 8, tertiary = 7,
        residential = 6, unclassified = 5, service = 4,
        footway = 2, path = 2, cycleway = 2)
half_width <- unname(hw[roads$highway])
half_width[is.na(half_width)] <- 5

corridor_osm <- viewscape::road_corridor_mask(roads, test_dsm,
                                              half_width = half_width)
terra::plot(corridor_osm)

# and straight into either understory function
pv_osm <- viewscape::pano_view_understory(
  dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
  vpt = vpt, cbh = 3, clear_mask = corridor_osm, max_dist = 120)

Two practical notes. Drop the classes you do not want before buffering — highway = "footway" also returns pedestrian paths through parks, which are open space but not street corridors, and motorway in a rural extract can pull in a carriageway that dominates the mask. And reprojecting lon/lat to a State Plane CRS such as this dataset’s EPSG:2253 needs a NADCON5 grid shift file; if PROJ cannot find it locally and has no network access it degrades silently, dropping vertices rather than failing. road_corridor_mask() will warn that it discarded degenerate geometries, which is the symptom to watch for — the fix is a local proj-data installation, not a change to the call.

4.7 Foliage transmittance

Left opaque, a crown stops the first ray that touches it. That makes canopy an upper bound and sky a lower one. Setting foliage_k switches crowns from walls to attenuators: optical depth accumulates along the path travelled inside the prism, T = exp(-k * L), and each direction is partitioned between sky, foliage and surface rather than labelled with one of them.

pv_t <- viewscape::pano_view_understory(
  dsm = test_dsm, dem = test_dtm, canopy = test_canopy_bin,
  built = test_building_bin, vpt = vpt, h = 1.7, cbh = 3,
  max_dist = 120, foliage_k = 0.1)     # extinction coefficient, per metre

round(pv_t$metrics, 3)
#>              sky           ground            built           canopy
#>            0.229            0.603            0.005            0.163
#>              svf  canopy_overhead canopy_enclosure       mean_depth
#>            0.459            0.327            0.089           13.033

# the sky layer is now continuous, and directly comparable with the soft sky
# mask of a segmented street-level photograph
terra::plot(pv_t$pano[["sky"]])

Treat foliage_k as an effective parameter to be calibrated, not as a leaf area measurement. The turbid-medium relation k = G * LAD would put broadleaf crowns near 0.25 to 1.0 per metre, but the prism fills the whole column from crown base to canopy top, whereas a real crown is largely a shell and a DSM-derived canopy top over a patch is the upper envelope of several crowns. The model therefore overstates the path length through leaves, and the coefficient that reproduces observed sky has to absorb that: on a synthetic 32%-cover scene, k = 0.5 left a sky view factor of 0.006 and only k near 0.1 produced the 0.2 to 0.3 range a tree-lined street plausibly has. Fit it against segmented sky pixels and report the fitted value.

4.8 Which choices matter for which metric

The decisions above do not carry equal weight, and it is worth knowing which of your numbers depend on which assumption. All figures are for the road viewpoint above, with a 150 m analysis radius and a 120 m ray limit:

The observer fix is worth two orders of magnitude and is not optional. The crown band is worth another two. Beyond that, the corridor and the transmittance move sky view a great deal and ground visibility almost not at all — so structural visibility conclusions are robust to them, and canopy exposure conclusions are not.

4.9 Known limitations

• A crown is a vertical prism over its cell footprint, so rounded crown margins are wrong at roughly cell scale.
• The prism is fully opaque unless foliage_k is set, and foliage_k needs calibrating.
• cbh is not monotone in its effect once it exceeds the local canopy height, as section 4.6 shows. Check it against the canopy height model.
• Trunks are absent, since the column below the crown base is empty. Where stem positions are known they can be burned into canopy as non-canopy cells.
• No curvature or refraction correction is applied.
• compute_viewshed_understory() evaluates point-to-point visibility for every target independently, which is not what compute_viewshed() computes: the view-tree and reference-plane algorithms propagate a horizon between cells, so an isolated obstruction casts a shadow that persists and widens behind it. On the same surface with no canopy the understory function is the more conservative of the two at grazing incidence.