Skip to content

Commit

Permalink
Updated algorithm from https://github.com/STBrinkmann/protoVS
Browse files Browse the repository at this point in the history
  • Loading branch information
@STBrinkmann
STBrinkmann committed Sep 7, 2022
1 parent d98cce5 commit 3f780f8
Show file tree
Hide file tree
Showing 33 changed files with 886 additions and 1,122 deletions.
3 changes: 0 additions & 3 deletions NAMESPACE
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,16 +7,13 @@ export(vgvi_from_sf)
export(viewshed)
export(visibility_proportion)
export(visualizeWeights)
importFrom(doParallel,registerDoParallel)
importFrom(dplyr,as_tibble)
importFrom(dplyr,everything)
importFrom(dplyr,mutate)
importFrom(dplyr,n)
importFrom(dplyr,relocate)
importFrom(dplyr,rename)
importFrom(dplyr,select)
importFrom(foreach,"%dopar%")
importFrom(foreach,foreach)
importFrom(ggplot2,aes)
importFrom(ggplot2,geom_smooth)
importFrom(ggplot2,ggplot)
Expand Down
24 changes: 8 additions & 16 deletions R/RcppExports.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,22 +9,6 @@ LoS_reference <- function(x0_ref, y0_ref, r, nc_ref) {
.Call(`_GVI_LoS_reference`, x0_ref, y0_ref, r, nc_ref)
}

VGVI_cpp <- function(dsm, dsm_values, greenspace, greenspace_values, x0, y0, radius, h0, fun, m, b) {
.Call(`_GVI_VGVI_cpp`, dsm, dsm_values, greenspace, greenspace_values, x0, y0, radius, h0, fun, m, b)
}

viewshed_cpp_sum1 <- function(dsm, dsm_values, x0, y0, radius, h0) {
.Call(`_GVI_viewshed_cpp_sum1`, dsm, dsm_values, x0, y0, radius, h0)
}

viewshed_cpp_sum2 <- function(dsm, dsm_values, x0, y0, radius, h0) {
.Call(`_GVI_viewshed_cpp_sum2`, dsm, dsm_values, x0, y0, radius, h0)
}

viewshed_cpp <- function(dsm, dsm_values, x0, y0, radius, h0) {
.Call(`_GVI_viewshed_cpp`, dsm, dsm_values, x0, y0, radius, h0)
}

viewshed_distance_analysis_cpp <- function(dsm, dsm_values, x0, y0, radius, h0, ncores = 1L, display_progress = FALSE) {
.Call(`_GVI_viewshed_distance_analysis_cpp`, dsm, dsm_values, x0, y0, radius, h0, ncores, display_progress)
}
Expand All @@ -33,3 +17,11 @@ viewshed_and_greenness_distance_analysis_cpp <- function(dsm, dsm_values, greens
.Call(`_GVI_viewshed_and_greenness_distance_analysis_cpp`, dsm, dsm_values, greenspace, greenspace_values, x0, y0, radius, h0, ncores, display_progress)
}

VGVI_cpp <- function(dsm, dsm_values, greenspace, greenspace_values, x0, y0, h0, radius, fun, m, b, ncores = 1L, display_progress = FALSE) {
.Call(`_GVI_VGVI_cpp`, dsm, dsm_values, greenspace, greenspace_values, x0, y0, h0, radius, fun, m, b, ncores, display_progress)
}

viewshed_cpp <- function(dsm, dsm_values, x0, y0, h0, radius, display_progress = FALSE) {
.Call(`_GVI_viewshed_cpp`, dsm, dsm_values, x0, y0, h0, radius, display_progress)
}

18 changes: 11 additions & 7 deletions R/sf_to_rast.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,9 +132,7 @@ sf_to_rast <- function(observer, v, aoi = NULL, max_distance = Inf, n = Inf, bet

#### 3. Prepare data for interpolation analysis ####
# observer
if(!is.null(aoi)){
observer <- observer[aoi,]
} else {
if(is.null(aoi)){
aoi <- st_as_sfc(sf::st_bbox(observer))
}

Expand All @@ -146,7 +144,7 @@ sf_to_rast <- function(observer, v, aoi = NULL, max_distance = Inf, n = Inf, bet
# mode
if(mode == 1) {
box_size <- (as.integer(max_distance/raster_res)) * (as.integer(max_distance/raster_res))
mode <- ifelse(box_size > n, 0, 1)
mode <- ifelse(box_size > n, 1, 0)
rm(box_size)
}

Expand All @@ -156,8 +154,14 @@ sf_to_rast <- function(observer, v, aoi = NULL, max_distance = Inf, n = Inf, bet
ymin = sf::st_bbox(aoi)[2],
ymax = sf::st_bbox(aoi)[4],
crs = sf::st_crs(aoi)$proj4string,
resolution = raster_res) %>%
terra::rasterize(terra::vect(observer), ., v, background = NA)
resolution = raster_res)
obs_cells <- terra::cellFromXY(iwd_rast, sf::st_coordinates(observer))

# Remove observer outside aoi
observer <- observer[which(!is.na(obs_cells)), ]
obs_cells <- na.omit(obs_cells)

iwd_rast[obs_cells] <- dplyr::pull(observer, v)

if (progress) setTxtProgressBar(pb, 2)

Expand All @@ -168,7 +172,7 @@ sf_to_rast <- function(observer, v, aoi = NULL, max_distance = Inf, n = Inf, bet

#### 4. IWD ####
if (progress) {
message("Computing IWD interpolation:")
message("\nComputing IWD interpolation:")
}

iwd_vals <- IDW_cpp(rast = iwd_cpp_rast, x = iwd_raster_vec,
Expand Down
96 changes: 22 additions & 74 deletions R/vgvi_from_sf.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#' @param b numeric; See ‘Details’
#' @param mode character; 'logit' or 'exponential'. See ‘Details’
#' @param cores numeric; The number of cores to use, i.e. at most how many child processes will be run simultaneously
#' @param chunk_size numeric; Chunk size for parallelization. See ‘Details’
#' @param folder_path optional; Folder path to where the output should be saved continuously. Must not inklude a filename extension (e.g. '.shp', '.gpkg').
#' @param progress logical; Show progress bar and computation time?
#'
Expand All @@ -29,9 +28,6 @@
#' The type of function, used for calculating the distance decay weights, can be defined with the \code{mode} parameter.
#' The argument 'logit' uses the logistic function, d = 1 / (1 + e^(b * (x - m))) and 'exponential' the exponential function d = 1 / (1 + (b * x^m)).
#' The decay function can be visualized using the \code{\link[GVI]{visualizeWeights}} function.
#'
#' Higher values of chunk_size increase computation time, but may also be more RAM intensive.
#' It is highly recommended to use a Linux or Mac system, for better parallel performance.
#'
#' @return sf_object containing the weighted VGVI values as POINT features, where 0 = no green cells are visible, and 1 = all of the visible cells are green.
#' @export
Expand Down Expand Up @@ -72,17 +68,13 @@
#' @importFrom methods is
#' @importFrom utils txtProgressBar
#' @importFrom utils setTxtProgressBar
#' @importFrom doParallel registerDoParallel
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#' @useDynLib GVI, .registration = TRUE

vgvi_from_sf <- function(observer, dsm_rast, dtm_rast, greenspace_rast,
max_distance = 800, observer_height = 1.7,
raster_res = NULL, spacing = raster_res,
m = 0.5, b = 8, mode = c("logit", "exponential"),
cores = 1, chunk_size = 10000,
folder_path = NULL, progress = FALSE) {
cores = 1, folder_path = NULL, progress = FALSE) {

#### 1. Check input ####
# observer
Expand Down Expand Up @@ -143,12 +135,12 @@ vgvi_from_sf <- function(observer, dsm_rast, dtm_rast, greenspace_rast,
}

# mode
if (is.character(mode) && (mode == c("logit", "exponential"))){
mode = 2
if (is.character(mode) && (mode == "logit")){
mode = 1
} else if (is.character(mode) && (mode == "exponential")){
mode = 2
} else if (is.character(mode) && (mode == "logit")){
mode = 1
} else if (is.character(mode) && (mode == c("logit", "exponential"))){
mode = 2
} else {
stop("mode must be logit or exponential")
}
Expand Down Expand Up @@ -250,15 +242,10 @@ vgvi_from_sf <- function(observer, dsm_rast, dtm_rast, greenspace_rast,
id = 1:dplyr::n()) %>%
dplyr::select(id, VGVI, dplyr::everything())

# Convert to list
observer_list <- suppressWarnings(
1:nrow(observer) %>%
split(seq(1, length(.), chunk_size))
)

# convert x0/y0 to col/row
x0 <- terra::colFromX(dsm_rast, x0)
y0 <- terra::rowFromY(dsm_rast, y0)
c0 <- terra::colFromX(dsm_rast, x0)
r0 <- terra::rowFromY(dsm_rast, y0)

if (progress) setTxtProgressBar(pb, 4)

Expand All @@ -274,71 +261,28 @@ vgvi_from_sf <- function(observer, dsm_rast, dtm_rast, greenspace_rast,
} else if (length(invalid_points) > 1) {
message(paste(length(invalid_points), "points have been removed, because they were outside of the DSM or DTM"))
}

invisible(gc())

#### 7. Calculate viewsheds and VGVI ####
if (progress) {
message(paste0("Computing VGVI for ", nrow(observer), ifelse(nrow(observer)>1, " points:", " point:")))
pb = txtProgressBar(min = 0, max = length(observer_list), initial = 0, style = 3)
#pb = txtProgressBar(min = 0, max = length(observer_list), initial = 0, style = 3)
start_time <- Sys.time()
}

if (cores > 1 && Sys.info()[["sysname"]] == "Windows") {
cl <- parallel::makeCluster(cores)
doParallel::registerDoParallel(cl)
}
vgvi_values <- VGVI_cpp(dsm = dsm_cpp_rast, dsm_values = dsm_vec,
greenspace = greenspace_cpp_rast, greenspace_values = greenspace_vec,
x0 = c0, y0 = r0, radius = max_distance, h0 = height_0_vec,
fun = mode, m = m, b = b, ncores = cores, display_progress = progress)

valid_values <- unlist(lapply(vgvi_values, is.numeric), use.names = FALSE)
observer[valid_values,2] <- vgvi_values[valid_values]

for (j in seq_along(observer_list)){
this_aoi <- observer[observer_list[[j]], ]
this_ids <- this_aoi$id

#### Calculate VGVI
if (cores > 1) {
if (Sys.info()[["sysname"]] == "Windows") {
# Calculate VGVI in parallel
par_fun <- function(i){
VGVI_cpp(dsm = dsm_cpp_rast, dsm_values = dsm_vec,
greenspace = greenspace_cpp_rast, greenspace_values = greenspace_vec,
x0 = x0[i], y0 = y0[i], radius = max_distance, h0 = height_0_vec[i],
fun = mode, m = m, b = b)
}

vgvi_list <- foreach::foreach(i=this_ids) %dopar% par_fun(i)
}
else {
vgvi_list <- parallel::mclapply(this_ids, function(i){
VGVI_cpp(dsm = dsm_cpp_rast, dsm_values = dsm_vec,
greenspace = greenspace_cpp_rast, greenspace_values = greenspace_vec,
x0 = x0[i], y0 = y0[i], radius = max_distance, h0 = height_0_vec[i],
fun = mode, m = m, b = b)},
mc.cores = cores, mc.preschedule = TRUE)
}
} else {
vgvi_list <- lapply(this_ids, function(i){
VGVI_cpp(dsm = dsm_cpp_rast, dsm_values = dsm_vec,
greenspace = greenspace_cpp_rast, greenspace_values = greenspace_vec,
x0 = x0[i], y0 = y0[i], radius = max_distance, h0 = height_0_vec[i],
fun = mode, m = m, b = b)})
}

vgvi_values <- unlist(vgvi_list, use.names = FALSE)

# Update observer
valid_values <- unlist(lapply(vgvi_values, is.numeric), use.names = FALSE)
observer[this_ids[valid_values],2] <- vgvi_values[valid_values]

if (!is.null(folder_path)) {
sf::st_write(observer[this_ids, ], folder_path, append = TRUE, quiet = T)
}

# Update ProgressBar
if (progress) setTxtProgressBar(pb, j)
}
if (cores > 1 && Sys.info()[["sysname"]] == "Windows") {
parallel::stopCluster(cl)
if (!is.null(folder_path)) {
sf::st_write(observer, folder_path, append = TRUE, quiet = T)
}


if (progress) {
time_dif <- round(cores * ((as.numeric(difftime(Sys.time(), start_time, units = "s"))*1000) / nrow(observer)), 2)
cat("\n")
Expand All @@ -359,5 +303,9 @@ vgvi_from_sf <- function(observer, dsm_rast, dtm_rast, greenspace_rast,

message(paste("Average time for a single point:", time_dif, "milliseconds"))
}

rm(dsm_cpp_rast, dsm_vec, greenspace_cpp_rast, greenspace_vec, c0, r0, height_0_vec)
invisible(gc())

return(observer)
}
13 changes: 7 additions & 6 deletions R/viewshed.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,15 +83,14 @@ viewshed <- function(observer, dsm_rast, dtm_rast,
rm(dsm_res)

# observer inside DSM
observer_cell <- terra::cellFromXY(object = dsm_rast, xy = sf::st_coordinates(observer))
if(is.na(observer_cell)) {
if(is.na(terra::cellFromXY(object = dsm_rast, xy = sf::st_coordinates(observer)))) {
stop("observer outside dsm_rast")
}

#### 2. Prepare Data for viewshed analysis ####
# Coordinates of start point
x0 <- terra::xFromCell(dsm_rast, observer_cell)
y0 <- terra::yFromCell(dsm_rast, observer_cell)
x0 <- sf::st_coordinates(observer)[1]
y0 <- sf::st_coordinates(observer)[2]

# AOI
output <- terra::rast(crs = terra::crs(dsm_rast),
Expand Down Expand Up @@ -126,8 +125,10 @@ viewshed <- function(observer, dsm_rast, dtm_rast,
dsm_cpp_rast <- terra::rast(dsm_rast_masked) %>% raster::raster()

# Apply viewshed (C++) function
viewshed <- viewshed_cpp(dsm_cpp_rast, dsm_vec, c0, r0, max_distance, height0)

viewshed <- viewshed_cpp(dsm = dsm_cpp_rast, dsm_values = dsm_vec,
x0 = c0, y0 = r0, h0 = height0,
radius = max_distance)

# Copy result of lineOfSight to the output raster
output[viewshed] <- 1

Expand Down
Loading

0 comments on commit 3f780f8

Please sign in to comment.