Sebastian T. Brinkmann (31.3.2022)
For this analysis we used our protoVS R package, which has been built for this paper specifically. The functions used in this analysis are the same as in the GVI R package (Brinkmann and Labib 2021), that is being used for calculating VGVI.
library(protoVS)
library(terra)
library(sf)
library(tidyverse)For the viewshed experiment we used LiDAR derived a Digital Surface Model (DSM) and Digital Terrain Model (DTM) with 1m spatial resolution of the Vancouver metropolitan area (Natural Resources Canada 2019).
DTM <- rast("data/DTM_Vancouver_1m.tif")
DSM <- rast("data/DSM_Vancouver_1m.tif")Publicly available Land Use and Land Cove (LULC) data has been acquired by Metro Vancouver (31.11.2019) at 2m resolution and reclassified to a binary greenness raster (0=no-green; 1=green).
LULC <- rast("data/LULC_2014.tif")
rcl_mat <- matrix(c(1, 6, 0, # no vegetation
6, 11, 1, # vegetation
11, 14, 0), # no vegetation
ncol = 3, byrow = TRUE)
GS <- terra::classify(LULC, rcl = rcl_mat, include.lowest = TRUE, right = FALSE)To demonstrate runtime for a large-scale greenness visibility study we estimated city-wide VGVI at 5 m intervals, resulting in 17,329,345 observer locations. We will remove observer locations that are located inside buildings or on water based on the LULC raster.
# 5 m intervals
DSM_5 <- aggregate(DSM, 5)
observer <- xyFromCell(DSM_5, which(!is.na(DSM_5[]))) %>%
sfheaders::sf_point() %>%
st_sf(crs = crs(DSM_5)) %>%
dplyr::rename(geom = geometry)
# Buildings and waters
lulc_build_water <- (LULC %in% c(1, 12)) %>%
terra::mask(LULC)
# intersect with LULC
observer_lulc_vals <- unlist(terra::extract(lulc_build_water, y = st_coordinates(observer)),
use.names = FALSE)
valid_obs <- which(observer_lulc_vals == 0)
observer <- observer[valid_obs,]As suggested by Labib, Huck, and Lindley (2021) we assumed eye-level observer height of 1.7 m and used an exponential distance decay weight function with a viewing distance threshold of 800 m. Total runtime time using 20 CPU threads was 20.3 hours, computation time per observer was 84 milliseconds.The result of the city-wide VGVI calculation is mapped in the figure below.
vgvi_1 <- vgvi_from_sf(observer = observer,
dsm_rast = DSM, dtm_rast = DTM, greenspace_rast = GS,
max_distance = 800, observer_height = 1.7,
m = 1, b = 3, mode = "exponential", cores = 20, progress = TRUE)
Click the figure for an interactive version
To compare our algorithm against other results by Labib, Huck, and Lindley (2021) and Cimburova and Blumentrath (2022), respectively, we also computed VGVI with 5 m resolution elevation models and 100 m viewing distance, and 1 m resolution elevation models and 100 m viewing distance. As the number of observer locations differ, total runtime is not comparable. We therefore calculated computation time per point.
DSM/DTM resolution: 5 m
Viewing distance: 800 m
| Labib et al. 2021 | Brinkmann et al. 2022 | |
|---|---|---|
| Total runtime | 11.5 days | 73 minutes |
| Computation time per point | 0.8 seconds | 5 milliseconds |
DSM_5 <- aggregate(DSM, 5)
DTM_5 <- aggregate(DTM, 5)
vgvi_2 <- vgvi_from_sf(observer = observer,
dsm_rast = DSM_5, dtm_rast = DTM_5, greenspace_rast = GS,
max_distance = 800, observer_height = 1.7,
m = 1, b = 3, mode = "exponential", cores = 20, progress = TRUE)DSM/DTM resolution: 1 m
Viewing distance: 100 m
| Cimburova & Blumentrath (2022) | Brinkmann et al. 2022 | |
|---|---|---|
| Total runtime | 5.6 days | 33 minutes |
| Computation time per point | 80 milliseconds | 2 milliseconds |
vgvi_3 <- vgvi_from_sf(observer = observer,
dsm_rast = DSM, dtm_rast = DTM, greenspace_rast = GS,
max_distance = 100, observer_height = 1.7,
m = 1, b = 3, mode = "exponential", cores = 20, progress = TRUE)library(GVI)
# Convert points to raster using a parallel IDW algorithm
# We used IDW interpolation for visualisation purposes only
vgvi_rast <- sf_to_rast(observer = vgvi_1, v = "VGVI", aoi = aoi,
max_distance = 400, n = 10, raster_res = 5,
cores = 22, progress = TRUE)
# Smoothing
vgvi_rast_5 <- focal(vgvi_rast, 3, fun = "median", na.rm = TRUE)
# Water mask from LULC
water <- as.polygons(LULC == 12) %>%
st_as_sf() %>%
filter(LULC_2014 == 1) %>%
mutate(a = st_area(.) %>% as.numeric()) %>%
filter(a > 1000) %>%
vect()
# Remove water and mask by AOI
aoi <- vect("data/Vancouver.gpkg")
vgvi_rast_5 <- vgvi_rast_5 %>%
crop(aoi) %>%
mask(water, inverse = TRUE) %>%
mask(aoi)
# Using Jenks algorithm to reclassified raster with 9 classes
set.seed(1234)
this_jenks <- classInt::classIntervals(var = vgvi_rast_5 %>%
terra::values(mat = FALSE) %>%
na.omit() %>%
sample(100000),
n = 9, style = "fisher", warnLargeN = FALSE)
this_jenks$brks[c(1, 10)] <- c(vgvi_rast_5@ptr$range_min, vgvi_rast_5@ptr$range_max)
rcl_mat <- matrix(c(this_jenks$brks[1:9],
this_jenks$brks[2:10],
1:9),
ncol = 3, byrow = F)
vgvi_rast_5_jenks <- classify(vgvi_rast_5, rcl_mat, include.lowest=TRUE)
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}
vgvi_rast_5_jenks[] <- range01(vgvi_rast_5_jenks[])Brinkmann, Sebastian T., and S. M. Labib. 2021. “GVI: Greenness Visibility Index r Package.” Zenodo. https://doi.org/10.5281/zenodo.5068835.
Cimburova, Zofie, and Stefan Blumentrath. 2022. “Viewshed-Based Modelling of Visual Exposure to Urban Greenery – an Efficient GIS Tool for Practical Planning Applications.” Landscape and Urban Planning 222: 104395. https://doi.org/10.1016/j.landurbplan.2022.104395.
Labib, S. M., Jonny J. Huck, and Sarah Lindley. 2021. “Modelling and Mapping Eye-Level Greenness Visibility Exposure Using Multi-Source Data at High Spatial Resolutions.” The Science of the Total Environment 755 (Pt 1): 143050. https://doi.org/10.1016/j.scitotenv.2020.143050.
Metro Vancouver. 31.11.2019. “Land Cover Classification 2014 - 2m LiDAR (Raster).” http://www.metrovancouver.org/data.
Natural Resources Canada. 2019. “High Resolution Digital Elevation Model (HRDEM) - CanElevation Serie.” https://open.canada.ca/data/en/dataset/957782bf-847c-4644-a757-e383c0057995.