diff --git a/NAMESPACE b/NAMESPACE index 6aae3de..fe1ce11 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,7 +7,6 @@ export(vgvi_from_sf) export(viewshed) export(visibility_proportion) export(visualizeWeights) -importFrom(doParallel,registerDoParallel) importFrom(dplyr,as_tibble) importFrom(dplyr,everything) importFrom(dplyr,mutate) @@ -15,8 +14,6 @@ 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) diff --git a/R/RcppExports.R b/R/RcppExports.R index fa5b311..6bf1bef 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -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) } @@ -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) +} + diff --git a/R/sf_to_rast.R b/R/sf_to_rast.R index 7cdb009..582aa2b 100644 --- a/R/sf_to_rast.R +++ b/R/sf_to_rast.R @@ -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)) } @@ -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) } @@ -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) @@ -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, diff --git a/R/vgvi_from_sf.R b/R/vgvi_from_sf.R index 3ab4dfa..20a5f78 100644 --- a/R/vgvi_from_sf.R +++ b/R/vgvi_from_sf.R @@ -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? #' @@ -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 @@ -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 @@ -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") } @@ -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) @@ -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") @@ -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) } \ No newline at end of file diff --git a/R/viewshed.R b/R/viewshed.R index 7923531..c308008 100644 --- a/R/viewshed.R +++ b/R/viewshed.R @@ -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), @@ -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 diff --git a/README.Rmd b/README.Rmd index 7393967..451f1ed 100644 --- a/README.Rmd +++ b/README.Rmd @@ -18,7 +18,7 @@ knitr::opts_chunk$set( [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.5068835.svg)](https://doi.org/10.5281/zenodo.5068835) [![Sample Data](https://badgen.net/badge/Sample%20Data/10.5281%252Fzenodo.5061257/blue?)](https://doi.org/10.5281/zenodo.5061257) [![jonnyhuck Python Library](https://badgen.net/badge/jonnyhuck/Python%20Library/blue?icon=github)](https://github.com/jonnyhuck/green-visibility-index) -The `GVI` R package helps researchers compute the Greenness Visibility Index (GVI) presented by [Labib, Huck and Lindley (2021)](https://doi.org/10.1016/j.scitotenv.2020.143050). The GVI is calculated using a Digital Surface Model (DSM), Digital Terrain Model (DTM) and Greenness Raster. `GVI` is written in C++ to provide fast and light weighted functionality. +The `GVI` R package helps researchers compute the Greenness Visibility Index (GVI) presented by [Labib, Huck and Lindley (2021)](https://doi.org/10.1016/j.scitotenv.2020.143050) and includes optimisation of the underlying algortihm proposed by [Brinkmann et al. (2022)](https://doi.org/10.5194/agile-giss-3-27-2022). The GVI is calculated using a Digital Surface Model (DSM), Digital Terrain Model (DTM) and Greenness Raster. `GVI` is written in C++ to provide fast and light weighted functionality. - [Functions](#functions) @@ -32,16 +32,16 @@ The `GVI` R package helps researchers compute the Greenness Visibility Index (GV GVI can also be computed in Python, by [Jonny Huck](https://github.com/jonnyhuck/green-visibility-index) -## Functions +## Functions {#functions} | Function | Description | -|---------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------| +|-------------------------|-----------------------------------------------| | [*viewshed*](#viewshed) | Computes the viewshed of a single point on a Digital Surface Model | | [*visualizeWeights*](#visualize-weights) | Helper function, to adjust spatial weight parameters in the *vgvi* and *vgvi_from_sf* functions | | [*vgvi_from_sf*](#viewshed-greenness-visibility-index-vgvi-from-sf) | Computes the *viewshed* and calculates the proportion of visible greenspace; Supports multiple points, lines or polygons | | sf_to_rast | Function for computing a continuous raster map from a sf object by interpolating the sf values using Inverse Distance Weighting (IDW). | -## Installation +## Installation {#installation} 1. Install [R](https://cran.r-project.org/) @@ -57,12 +57,12 @@ GVI can also be computed in Python, by [Jonny Huck](https://github.com/jonnyhuck - Install `GVI` - `GVI` is still in very active development. Therefore, the package is also not on CRAN yet. You can install the latest version of `GVI` from GitHub with: + `GVI` is still in very active development. Therefore, the package is also not on CRAN yet. You can install the latest version of `GVI` from GitHub with:\ `remotes::install_git("https://github.com/STBrinkmann/GVI")` -## Methods +## Methods {#methods} -### Viewshed +### Viewshed {#viewshed} The `viewshed` function computes a binary viewshed of a single point on a Digital Surface Model (DSM) raster. A radial buffer is applied on the observer position, and visibility is beeing calculated usig a C++ implementation of Bresenham's line algorithm [[Bresenham 1965](https://doi.org/10.1147/sj.41.0025), [Bresenham 1977](https://doi.org/10.1145/359423.359432)] and simple geometry. The result of the `viewshed` function is a radial raster where 0 = no-visible and 1 = visible area. @@ -74,7 +74,9 @@ Based on a viewshed and a binary greenspace raster, all visible points are class $$\begin{align*} f(x) = \cfrac{1}{1 + e^{b(x-m)}} && \text{(logistic)}\\ f(x) = \cfrac{1}{1 + (bx^{m})} && \text{(exponential)} \end{align*} $$ -### Visualize Weights +The full algorithm has been described in [Brinkmann et al. (2022)](https://doi.org/10.5194/agile-giss-3-27-2022). + +### Visualize Weights {#visualize-weights} The `visualizeWeights` function helps to adjust spatial weight parameters *m* and *b* used in the `vgvi` and `vgvi_from_sf` functions. @@ -82,11 +84,11 @@ The `visualizeWeights` function helps to adjust spatial weight parameters *m* an *Example output of the* `visualizeWeight` *function to compare and parameterize the decay weights of a logistic (left) and an exponential (right) function.* -### Viewshed Greenness Visibility Index (VGVI) from sf +### Viewshed Greenness Visibility Index (VGVI) from sf {#viewshed-greenness-visibility-index-vgvi-from-sf} The `vgvi_from_sf` function combines the `viewshed` function to compute a raster of visible cells, and calculates the VGVI using a binary GreenSpace map. The input of the observer location is an [sf](https://r-spatial.github.io/sf/) object (supports points, lines and polygons). If this object is a line or polygon feature, it will be converted to points. For every point, first the viewshed, than the VGVI is being calculated. The output of the `vgvi_from_sf` function is a single sf point object containing the VGVI for every point. This function supports parallel computation on Windows, Linux and MacOS, however it is highly recommended to use Linux or MacOS. -## Examples +## Examples {#examples} ### Data @@ -147,7 +149,7 @@ vgvi_from_sf(observer = observer, m = 0.5, b = 8, mode = "logit")$VGVI ``` -The output of \~0.52 indicates, that \~52% of the visible area, within a 200 meters radius, is green. +The output of \~0.53 indicates, that \~53% of the visible area, within a 200 meters radius, is green. ### 2. Road Network @@ -181,7 +183,7 @@ Often researchers need to compute the VGVI for a large study area, therefore we ## About -### Citation +### Citation {#citation} Run this command to get info on how to cite this package. @@ -197,8 +199,8 @@ Dr. S.M. Labib (Author) e-mail: [sml80\@medschl.cam.ac.uk](mailto:sml80@medschl. ### Paper authors -S.M. Labib (1, 2\*) -Jonny J. Huck (1) +S.M. Labib (1, 2\*)\ +Jonny J. Huck (1)\ Sarah Lindley (1) 1: Department of Geography, School of Environment, Education and Development (SEED), University of Manchester, Arthur Lewis building (1st Floor), Oxford Road, Manchester M13 9PL, United Kingdom. @@ -211,6 +213,8 @@ Sarah Lindley (1) Bresenham, J.E. (1965): Algorithm for computer control of a digital plotter. IBM Systems Journal, vol. 4, no. 1, pp. 25-30, doi: [10.1147/sj.41.0025](https://doi.org/10.1147/sj.41.0025). -Bresenham, Jack (1977): A linear algorithm for incremental digital display of circular arcs. Commun. ACM 20, 2 (Feb. 1977), 100–106. doi: [10.1145/359423.359432](https://doi.org/10.1145/359423.359432). +Bresenham, Jack (1977): A linear algorithm for incremental digital display of circular arcs. Commun. ACM 20, 2 (Feb. 1977), 100--106. doi: [10.1145/359423.359432](https://doi.org/10.1145/359423.359432). + +Brinkmann, S. T., Kremer, D., and Walker, B. B. (2022): Modelling eye-level visibility of urban green space: Optimising city-wide point-based viewshed computations through prototyping, AGILE GIScience Ser., 3, 27, [https://doi.org/10.5194/agile-giss-3-27-2022](https://doi.org/10.5194/agile-giss-3-27-2022,). 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. *Science of The Total Environment* 755 (February): 143050. doi: [j.scitotenv.2020.143050](https://doi.org/10.1016/j.scitotenv.2020.143050). diff --git a/README.md b/README.md index 9fd1efa..b80501b 100644 --- a/README.md +++ b/README.md @@ -9,7 +9,9 @@ Library](https://badgen.net/badge/jonnyhuck/Python%20Library/blue?icon=github)]( The `GVI` R package helps researchers compute the Greenness Visibility Index (GVI) presented by [Labib, Huck and Lindley -(2021)](https://doi.org/10.1016/j.scitotenv.2020.143050). The GVI is +(2021)](https://doi.org/10.1016/j.scitotenv.2020.143050) and includes +optimisation of the underlying algortihm proposed by [Brinkmann et +al. (2022)](https://doi.org/10.5194/agile-giss-3-27-2022). The GVI is calculated using a Digital Surface Model (DSM), Digital Terrain Model (DTM) and Greenness Raster. `GVI` is written in C++ to provide fast and light weighted functionality. @@ -29,12 +31,12 @@ Huck](https://github.com/jonnyhuck/green-visibility-index) ## Functions -| Function | Description | -|-----------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------| -| [*viewshed*](#viewshed) | Computes the viewshed of a single point on a Digital Surface Model | -| [*visualizeWeights*](#visualize-weights) | Helper function, to adjust spatial weight parameters in the *vgvi* and *vgvi\_from\_sf* functions | -| [*vgvi\_from\_sf*](#viewshed-greenness-visibility-index-vgvi-from-sf) | Computes the *viewshed* and calculates the proportion of visible greenspace; Supports multiple points, lines or polygons | -| sf\_to\_rast | Function for computing a continuous raster map from a sf object by interpolating the sf values using Inverse Distance Weighting (IDW). | +| Function | Description | +|---------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------------------| +| [*viewshed*](#viewshed) | Computes the viewshed of a single point on a Digital Surface Model | +| [*visualizeWeights*](#visualize-weights) | Helper function, to adjust spatial weight parameters in the *vgvi* and *vgvi_from_sf* functions | +| [*vgvi_from_sf*](#viewshed-greenness-visibility-index-vgvi-from-sf) | Computes the *viewshed* and calculates the proportion of visible greenspace; Supports multiple points, lines or polygons | +| sf_to_rast | Function for computing a continuous raster map from a sf object by interpolating the sf values using Inverse Distance Weighting (IDW). | ## Installation @@ -85,7 +87,10 @@ prominence of an object in space with increasing distance from the observer. Currently two options are supported, a logistic and an exponential function. -![\\begin{align\*} f(x) = \\cfrac{1}{1 + e^{b(x-m)}} && \\text{(logistic)}\\\\ f(x) = \\cfrac{1}{1 + (bx^{m})} && \\text{(exponential)} \\end{align\*} ](https://latex.codecogs.com/svg.latex?%5Cbegin%7Balign%2A%7D%20f%28x%29%20%3D%20%5Ccfrac%7B1%7D%7B1%20%2B%20e%5E%7Bb%28x-m%29%7D%7D%20%26%26%20%5Ctext%7B%28logistic%29%7D%5C%5C%20f%28x%29%20%3D%20%5Ccfrac%7B1%7D%7B1%20%2B%20%28bx%5E%7Bm%7D%29%7D%20%26%26%20%5Ctext%7B%28exponential%29%7D%20%5Cend%7Balign%2A%7D%20 "\begin{align*} f(x) = \cfrac{1}{1 + e^{b(x-m)}} && \text{(logistic)}\\ f(x) = \cfrac{1}{1 + (bx^{m})} && \text{(exponential)} \end{align*} ") +![\begin{align\*} f(x) = \cfrac{1}{1 + e^{b(x-m)}} && \text{(logistic)}\\\\ f(x) = \cfrac{1}{1 + (bx^{m})} && \text{(exponential)} \end{align\*}](https://latex.codecogs.com/png.image?%5Cdpi%7B110%7D&space;%5Cbg_white&space;%5Cbegin%7Balign%2A%7D%20f%28x%29%20%3D%20%5Ccfrac%7B1%7D%7B1%20%2B%20e%5E%7Bb%28x-m%29%7D%7D%20%26%26%20%5Ctext%7B%28logistic%29%7D%5C%5C%20f%28x%29%20%3D%20%5Ccfrac%7B1%7D%7B1%20%2B%20%28bx%5E%7Bm%7D%29%7D%20%26%26%20%5Ctext%7B%28exponential%29%7D%20%5Cend%7Balign%2A%7D "\begin{align*} f(x) = \cfrac{1}{1 + e^{b(x-m)}} && \text{(logistic)}\\ f(x) = \cfrac{1}{1 + (bx^{m})} && \text{(exponential)} \end{align*}") + +The full algorithm has been described in [Brinkmann et +al. (2022)](https://doi.org/10.5194/agile-giss-3-27-2022). ### Visualize Weights @@ -187,10 +192,10 @@ vgvi_from_sf(observer = observer, dsm_rast = DSM, dtm_rast = DEM, greenspace_rast = GreenSpace, max_distance = 200, observer_height = 1.7, m = 0.5, b = 8, mode = "logit")$VGVI -#> [1] 0.5199976 +#> [1] 0.5255536 ``` -The output of \~0.52 indicates, that \~52% of the visible area, within a +The output of \~0.53 indicates, that \~53% of the visible area, within a 200 meters radius, is green. ### 2. Road Network @@ -242,7 +247,7 @@ citation("GVI") #> Brinkmann, S.T. and Labib, S.M. (2021). GVI: R package for computing #> VGVI for Spatial Raster. doi: 10.5281/zenodo.5060325. #> -#> A BibTeX entry for LaTeX users is +#> Ein BibTeX-Eintrag für LaTeX-Benutzer ist #> #> @Manual{, #> title = {GVI: R package for computing VGVI for Spatial Raster}, @@ -285,6 +290,12 @@ Bresenham, Jack (1977): A linear algorithm for incremental digital display of circular arcs. Commun. ACM 20, 2 (Feb. 1977), 100–106. doi: [10.1145/359423.359432](https://doi.org/10.1145/359423.359432). +Brinkmann, S. T., Kremer, D., and Walker, B. B. (2022): Modelling +eye-level visibility of urban green space: Optimising city-wide +point-based viewshed computations through prototyping, AGILE GIScience +Ser., 3, 27, +[https://doi.org/10.5194/agile-giss-3-27-2022](https://doi.org/10.5194/agile-giss-3-27-2022,). + 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. *Science of The Total Environment* 755 diff --git a/man/distance_analysis.Rd b/man/distance_analysis.Rd index 6378748..c6302a8 100644 --- a/man/distance_analysis.Rd +++ b/man/distance_analysis.Rd @@ -9,7 +9,7 @@ distance_analysis( dsm_rast, dtm_rast, greenspace_rast = NULL, - max_distance = 1:800, + max_distance = 800, observer_height = 1.7, raster_res = NULL, summarise = FALSE, @@ -26,7 +26,7 @@ distance_analysis( \item{greenspace_rast}{optional; \code{\link[terra]{rast}} of the binary Greenspace mask. Values of the Greenspace mask must be 1 for Green and 0 for No-Green} -\item{max_distance}{numeric; Maximum buffer distance. If max_distance is a vector, the Viewshed Distance Analysis is conducted for multiple distance values.} +\item{max_distance}{numeric; Maximum buffer distance} \item{observer_height}{numeric > 0; Height of the observer (e.g. 1.7 meters)} diff --git a/man/figures/README-unnamed-chunk-4-1.png b/man/figures/README-unnamed-chunk-4-1.png index 44eeded..398e7b1 100644 Binary files a/man/figures/README-unnamed-chunk-4-1.png and b/man/figures/README-unnamed-chunk-4-1.png differ diff --git a/man/sf_to_rast.Rd b/man/sf_to_rast.Rd index ef9722f..d2c7c82 100644 --- a/man/sf_to_rast.Rd +++ b/man/sf_to_rast.Rd @@ -13,6 +13,7 @@ sf_to_rast( beta = 2, raster_res = NULL, spacing = raster_res, + na_only = FALSE, cores = 1, progress = FALSE ) diff --git a/man/vgvi_from_sf.Rd b/man/vgvi_from_sf.Rd index 1d9a45d..5215b15 100644 --- a/man/vgvi_from_sf.Rd +++ b/man/vgvi_from_sf.Rd @@ -17,7 +17,6 @@ vgvi_from_sf( b = 8, mode = c("logit", "exponential"), cores = 1, - chunk_size = 10000, folder_path = NULL, progress = FALSE ) @@ -48,8 +47,6 @@ The \code{spacing} parameter sets the distance in between the points on the line \item{cores}{numeric; The number of cores to use, i.e. at most how many child processes will be run simultaneously} -\item{chunk_size}{numeric; Chunk size for parallelization. See ‘Details’} - \item{folder_path}{optional; Folder path to where the output should be saved continuously. Must not inklude a filename extension (e.g. '.shp', '.gpkg').} \item{progress}{logical; Show progress bar and computation time?} @@ -71,7 +68,4 @@ The CRS (\code{\link[sf]{st_crs}}) needs to have a metric unit! 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. } diff --git a/src/GVI.dll b/src/GVI.dll index c94dd93..3d54947 100644 Binary files a/src/GVI.dll and b/src/GVI.dll differ diff --git a/src/GVI.so b/src/GVI.so deleted file mode 100644 index 264185f..0000000 Binary files a/src/GVI.so and /dev/null differ diff --git a/src/IDW_cpp.cpp b/src/IDW_cpp.cpp index 300f6b2..cdff781 100644 --- a/src/IDW_cpp.cpp +++ b/src/IDW_cpp.cpp @@ -52,16 +52,15 @@ NumericVector IDW_cpp(S4 &rast, const NumericVector &x, // Progress bar Progress pb(x.size(), display_progress); - + // Main loop: Loop over all values of the raster x #if defined(_OPENMP) omp_set_num_threads(ncores); #pragma omp parallel for shared(out) #endif for(int j = 0; j < x.size(); j++){ - if (!pb.is_aborted()) { - Progress::check_abort(); - + if ( pb.increment() ) { + if(na_only == false || (na_only && NumericVector::is_na(x[j]))){ // 1. Convert j to row/col and X/Y coordinates // row col from cell @@ -84,6 +83,7 @@ NumericVector IDW_cpp(S4 &rast, const NumericVector &x, // Distance const double dist = sqrt((x_j-sf_x[i])*(x_j-sf_x[i]) + (y_j-sf_y[i])*(y_j-sf_y[i])); + // If distance <= 0, use a small value (resolution / 4) instead if(dist <= 0){ d.push_back(rast_info.res/4); z.push_back(sf_z[i]); @@ -150,8 +150,6 @@ NumericVector IDW_cpp(S4 &rast, const NumericVector &x, } else { out[j] = x[j]; } - - pb.increment(); } } diff --git a/src/IDW_cpp.o b/src/IDW_cpp.o index bf069c5..533fe9c 100644 Binary files a/src/IDW_cpp.o and b/src/IDW_cpp.o differ diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3adc996..81b0b8c 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -46,30 +46,9 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// VGVI_cpp -Rcpp::NumericVector VGVI_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, Rcpp::S4& greenspace, const Rcpp::NumericVector& greenspace_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0, const int fun, const double m, const double b); -RcppExport SEXP _GVI_VGVI_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP greenspaceSEXP, SEXP greenspace_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP, SEXP funSEXP, SEXP mSEXP, SEXP bSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::S4& >::type dsm(dsmSEXP); - Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type dsm_values(dsm_valuesSEXP); - Rcpp::traits::input_parameter< Rcpp::S4& >::type greenspace(greenspaceSEXP); - Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type greenspace_values(greenspace_valuesSEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type x0(x0SEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type y0(y0SEXP); - Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); - Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type h0(h0SEXP); - Rcpp::traits::input_parameter< const int >::type fun(funSEXP); - Rcpp::traits::input_parameter< const double >::type m(mSEXP); - Rcpp::traits::input_parameter< const double >::type b(bSEXP); - rcpp_result_gen = Rcpp::wrap(VGVI_cpp(dsm, dsm_values, greenspace, greenspace_values, x0, y0, radius, h0, fun, m, b)); - return rcpp_result_gen; -END_RCPP -} -// viewshed_cpp_sum1 -Rcpp::IntegerVector viewshed_cpp_sum1(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0); -RcppExport SEXP _GVI_viewshed_cpp_sum1(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP) { +// viewshed_distance_analysis_cpp +Rcpp::IntegerMatrix viewshed_distance_analysis_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0, const int ncores, const bool display_progress); +RcppExport SEXP _GVI_viewshed_distance_analysis_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP, SEXP ncoresSEXP, SEXP display_progressSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -79,77 +58,69 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type y0(y0SEXP); Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type h0(h0SEXP); - rcpp_result_gen = Rcpp::wrap(viewshed_cpp_sum1(dsm, dsm_values, x0, y0, radius, h0)); + Rcpp::traits::input_parameter< const int >::type ncores(ncoresSEXP); + Rcpp::traits::input_parameter< const bool >::type display_progress(display_progressSEXP); + rcpp_result_gen = Rcpp::wrap(viewshed_distance_analysis_cpp(dsm, dsm_values, x0, y0, radius, h0, ncores, display_progress)); return rcpp_result_gen; END_RCPP } -// viewshed_cpp_sum2 -Rcpp::IntegerVector viewshed_cpp_sum2(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0); -RcppExport SEXP _GVI_viewshed_cpp_sum2(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP) { +// viewshed_and_greenness_distance_analysis_cpp +Rcpp::IntegerMatrix viewshed_and_greenness_distance_analysis_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, Rcpp::S4& greenspace, const Rcpp::NumericVector& greenspace_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0, const int ncores, const bool display_progress); +RcppExport SEXP _GVI_viewshed_and_greenness_distance_analysis_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP greenspaceSEXP, SEXP greenspace_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP, SEXP ncoresSEXP, SEXP display_progressSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::S4& >::type dsm(dsmSEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type dsm_values(dsm_valuesSEXP); + Rcpp::traits::input_parameter< Rcpp::S4& >::type greenspace(greenspaceSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type greenspace_values(greenspace_valuesSEXP); Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type x0(x0SEXP); Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type y0(y0SEXP); Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type h0(h0SEXP); - rcpp_result_gen = Rcpp::wrap(viewshed_cpp_sum2(dsm, dsm_values, x0, y0, radius, h0)); + Rcpp::traits::input_parameter< const int >::type ncores(ncoresSEXP); + Rcpp::traits::input_parameter< const bool >::type display_progress(display_progressSEXP); + rcpp_result_gen = Rcpp::wrap(viewshed_and_greenness_distance_analysis_cpp(dsm, dsm_values, greenspace, greenspace_values, x0, y0, radius, h0, ncores, display_progress)); return rcpp_result_gen; END_RCPP } -// viewshed_cpp -Rcpp::IntegerVector viewshed_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0); -RcppExport SEXP _GVI_viewshed_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP) { +// VGVI_cpp +std::vector VGVI_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, Rcpp::S4& greenspace, const Rcpp::NumericVector& greenspace_values, const Rcpp::IntegerVector& x0, const Rcpp::IntegerVector& y0, const Rcpp::NumericVector& h0, const int radius, const int fun, const double m, const double b, const int ncores, const bool display_progress); +RcppExport SEXP _GVI_VGVI_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP greenspaceSEXP, SEXP greenspace_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP h0SEXP, SEXP radiusSEXP, SEXP funSEXP, SEXP mSEXP, SEXP bSEXP, SEXP ncoresSEXP, SEXP display_progressSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::S4& >::type dsm(dsmSEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type dsm_values(dsm_valuesSEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type x0(x0SEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type y0(y0SEXP); - Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); + Rcpp::traits::input_parameter< Rcpp::S4& >::type greenspace(greenspaceSEXP); + Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type greenspace_values(greenspace_valuesSEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerVector& >::type x0(x0SEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerVector& >::type y0(y0SEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type h0(h0SEXP); - rcpp_result_gen = Rcpp::wrap(viewshed_cpp(dsm, dsm_values, x0, y0, radius, h0)); - return rcpp_result_gen; -END_RCPP -} -// viewshed_distance_analysis_cpp -Rcpp::IntegerMatrix viewshed_distance_analysis_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0, const int ncores, const bool display_progress); -RcppExport SEXP _GVI_viewshed_distance_analysis_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP, SEXP ncoresSEXP, SEXP display_progressSEXP) { -BEGIN_RCPP - Rcpp::RObject rcpp_result_gen; - Rcpp::RNGScope rcpp_rngScope_gen; - Rcpp::traits::input_parameter< Rcpp::S4& >::type dsm(dsmSEXP); - Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type dsm_values(dsm_valuesSEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type x0(x0SEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type y0(y0SEXP); Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); - Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type h0(h0SEXP); + Rcpp::traits::input_parameter< const int >::type fun(funSEXP); + Rcpp::traits::input_parameter< const double >::type m(mSEXP); + Rcpp::traits::input_parameter< const double >::type b(bSEXP); Rcpp::traits::input_parameter< const int >::type ncores(ncoresSEXP); Rcpp::traits::input_parameter< const bool >::type display_progress(display_progressSEXP); - rcpp_result_gen = Rcpp::wrap(viewshed_distance_analysis_cpp(dsm, dsm_values, x0, y0, radius, h0, ncores, display_progress)); + rcpp_result_gen = Rcpp::wrap(VGVI_cpp(dsm, dsm_values, greenspace, greenspace_values, x0, y0, h0, radius, fun, m, b, ncores, display_progress)); return rcpp_result_gen; END_RCPP } -// viewshed_and_greenness_distance_analysis_cpp -Rcpp::IntegerMatrix viewshed_and_greenness_distance_analysis_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, Rcpp::S4& greenspace, const Rcpp::NumericVector& greenspace_values, const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, const int radius, const Rcpp::NumericVector& h0, const int ncores, const bool display_progress); -RcppExport SEXP _GVI_viewshed_and_greenness_distance_analysis_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP greenspaceSEXP, SEXP greenspace_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP radiusSEXP, SEXP h0SEXP, SEXP ncoresSEXP, SEXP display_progressSEXP) { +// viewshed_cpp +std::vector viewshed_cpp(Rcpp::S4& dsm, const Rcpp::NumericVector& dsm_values, const Rcpp::IntegerVector& x0, const Rcpp::IntegerVector& y0, const Rcpp::NumericVector& h0, const int radius, const bool display_progress); +RcppExport SEXP _GVI_viewshed_cpp(SEXP dsmSEXP, SEXP dsm_valuesSEXP, SEXP x0SEXP, SEXP y0SEXP, SEXP h0SEXP, SEXP radiusSEXP, SEXP display_progressSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< Rcpp::S4& >::type dsm(dsmSEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type dsm_values(dsm_valuesSEXP); - Rcpp::traits::input_parameter< Rcpp::S4& >::type greenspace(greenspaceSEXP); - Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type greenspace_values(greenspace_valuesSEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type x0(x0SEXP); - Rcpp::traits::input_parameter< const Rcpp::IntegerVector >::type y0(y0SEXP); - Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerVector& >::type x0(x0SEXP); + Rcpp::traits::input_parameter< const Rcpp::IntegerVector& >::type y0(y0SEXP); Rcpp::traits::input_parameter< const Rcpp::NumericVector& >::type h0(h0SEXP); - Rcpp::traits::input_parameter< const int >::type ncores(ncoresSEXP); + Rcpp::traits::input_parameter< const int >::type radius(radiusSEXP); Rcpp::traits::input_parameter< const bool >::type display_progress(display_progressSEXP); - rcpp_result_gen = Rcpp::wrap(viewshed_and_greenness_distance_analysis_cpp(dsm, dsm_values, greenspace, greenspace_values, x0, y0, radius, h0, ncores, display_progress)); + rcpp_result_gen = Rcpp::wrap(viewshed_cpp(dsm, dsm_values, x0, y0, h0, radius, display_progress)); return rcpp_result_gen; END_RCPP } @@ -157,12 +128,10 @@ END_RCPP static const R_CallMethodDef CallEntries[] = { {"_GVI_IDW_cpp", (DL_FUNC) &_GVI_IDW_cpp, 12}, {"_GVI_LoS_reference", (DL_FUNC) &_GVI_LoS_reference, 4}, - {"_GVI_VGVI_cpp", (DL_FUNC) &_GVI_VGVI_cpp, 11}, - {"_GVI_viewshed_cpp_sum1", (DL_FUNC) &_GVI_viewshed_cpp_sum1, 6}, - {"_GVI_viewshed_cpp_sum2", (DL_FUNC) &_GVI_viewshed_cpp_sum2, 6}, - {"_GVI_viewshed_cpp", (DL_FUNC) &_GVI_viewshed_cpp, 6}, {"_GVI_viewshed_distance_analysis_cpp", (DL_FUNC) &_GVI_viewshed_distance_analysis_cpp, 8}, {"_GVI_viewshed_and_greenness_distance_analysis_cpp", (DL_FUNC) &_GVI_viewshed_and_greenness_distance_analysis_cpp, 10}, + {"_GVI_VGVI_cpp", (DL_FUNC) &_GVI_VGVI_cpp, 13}, + {"_GVI_viewshed_cpp", (DL_FUNC) &_GVI_viewshed_cpp, 7}, {NULL, NULL, 0} }; diff --git a/src/RcppExports.o b/src/RcppExports.o index 20fe148..4433c98 100644 Binary files a/src/RcppExports.o and b/src/RcppExports.o differ diff --git a/src/bresenham.o b/src/bresenham.o index aaa7dda..b83d33f 100644 Binary files a/src/bresenham.o and b/src/bresenham.o differ diff --git a/src/distance_analysis.cpp b/src/distance_analysis.cpp new file mode 100644 index 0000000..b5251e1 --- /dev/null +++ b/src/distance_analysis.cpp @@ -0,0 +1,233 @@ +#include +#include "rsinfo.h" +#include "rasterutils.h" +#include "bresenham.h" +#include "integrate.h" + +using namespace Rcpp; + +#ifdef _OPENMP +#include +#endif +// [[Rcpp::plugins(openmp)]] +// [[Rcpp::depends(RcppProgress)]] +#include + +// [[Rcpp::export]] +Rcpp::IntegerMatrix viewshed_distance_analysis_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, + const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, + const int radius, const Rcpp::NumericVector & h0, + const int ncores=1, const bool display_progress=false) +{ + // Cells from x0,y0 + Rcpp::IntegerVector x0_o = x0-1; + Rcpp::IntegerVector y0_o = y0-1; + const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); + + // Basic raster information + const RasterInfo ras(dsm); + + // Parameters + const int r = (int)(radius/ras.res); + const int nc_ref = (2*r)+1; + const int x0_ref = r, y0_ref = x0_ref; + const int c0_ref = y0_ref*nc_ref + x0_ref; + const int s = input_cells.size(); + + // Output + Rcpp::IntegerMatrix output( r*s,3 ); + for(int k = 0; k < s; k++) { + for(int d = 0; d dsm_values[input_cells[k]]){ + // Viewshed Analysis + const int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); + +#if defined(_OPENMP) + omp_set_num_threads(ncores); +#pragma omp parallel for shared(output) +#endif + for(int i = 0; i < (r*8); i++){ + if (!p.is_aborted()) { + Progress::check_abort(); + + double max_tan = -9999.0; + + for(int j = 0; j < r; j++){ + const int los_ref_cell = los_ref_vec[i*r + j]; + if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ + // Project reference cells to actual cell value + int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); + int row = trunc(cell/ras.ncol); + int col = cell - (row * ras.ncol); + int dcol = abs(col-x0_o[k]); + + if(!(cell < 0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ + // Compute tangents of LoS_vec and update output + double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); + double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); + + int d = round(distance_traveled) - 1; + d = (d < 0) ? 0 : d; + d = (d >= r) ? (r-1) : d; + + output( (d*s + k),1 ) += 1; + //output( d,1 ) += 1; + + if(this_tan > max_tan){ + max_tan = this_tan; + + output( (d*s + k),2 ) += 1; + //output( d,2 ) += 1; + } + } + } else { + break; + } + } + } + } + } + p.increment(); + } + + return output; +} + + +// [[Rcpp::export]] +Rcpp::IntegerMatrix viewshed_and_greenness_distance_analysis_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, + Rcpp::S4 &greenspace, const Rcpp::NumericVector &greenspace_values, + const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, + const int radius, const Rcpp::NumericVector & h0, + const int ncores=1, const bool display_progress=false) +{ + // Cells from x0,y0 + Rcpp::IntegerVector x0_o = x0-1; + Rcpp::IntegerVector y0_o = y0-1; + const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); + + // Basic raster information + const RasterInfo ras(dsm); + const RasterInfo gs_ras(greenspace); + + // Parameters + const int r = (int)(radius/ras.res); + const int nc_ref = (2*r)+1; + const int x0_ref = r, y0_ref = x0_ref; + const int c0_ref = y0_ref*nc_ref + x0_ref; + const int s = input_cells.size(); + + // Output + // v1: Distance + // V2: Total number of cells per distance + // V3: Number of all visible cells per distance + // V4: Total number of green cells per distance + // V5: Number of green visible cells per distance + Rcpp::IntegerMatrix output( r*s,5 ); + for(int k = 0; k < s; k++) { + for(int d = 0; d dsm_values[input_cells[k]]){ + // Viewshed Analysis + const int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); + + for(int i = 0; i < (r*8); i++){ + if (!p.is_aborted()) { + Progress::check_abort(); + + double max_tan = -9999.0; + + for(int j = 0; j < r; j++){ + const int los_ref_cell = los_ref_vec[i*r + j]; + if(los_ref_cell != NA_INTEGER){ //!Rcpp::IntegerVector::is_na(los_ref_cell) + // Project reference cells to actual cell value + int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); + + // row/col + int row = trunc(cell/ras.ncol); + int col = cell - (row * ras.ncol); + int dcol = abs(col-x0_o[k]); + + // X/Y coordinates + double x_j = ras.xmin + (col + 0.5) * ras.res; + double y_j = ras.ymax - (row + 0.5) * ras.res; + + // Cell of greenspace raster + // OpenMP can't call cellFromXY2 function without causing R to crash. + // Therefore, I'll calcualte the cell here manually + int cell_gs = cellFromXY2(gs_ras, x_j, y_j); + + + if(!(cell < 0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ + // Compute tangents of LoS_vec and update output + double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); + double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); + + int d = round(distance_traveled) - 1; + d = (d < 0) ? 0 : d; + d = (d >= r) ? (r-1) : d; + + output( (d*s + k),1 ) += 1; + if(cell_gs != NA_INTEGER && greenspace_values[cell_gs] == 1){ + output( (d*s + k),3 ) += 1; + } + + if(this_tan > max_tan){ + max_tan = this_tan; + + output( (d*s + k),2 ) += 1; + if(cell_gs != NA_INTEGER && greenspace_values[cell_gs] == 1){ + output( (d*s + k),4 ) += 1; + } + } + } + } else { + break; + } + } + } + } + } + p.increment(); + } + + return output; +} \ No newline at end of file diff --git a/src/distance_analysis.o b/src/distance_analysis.o new file mode 100644 index 0000000..dd0dbe2 Binary files /dev/null and b/src/distance_analysis.o differ diff --git a/src/eta_progress_bar.hpp b/src/eta_progress_bar.hpp new file mode 100644 index 0000000..16afdca --- /dev/null +++ b/src/eta_progress_bar.hpp @@ -0,0 +1,180 @@ +/* + * eta_progress_bar.hpp + * + * A custom ProgressBar class to display a progress bar with time estimation + * + * Author: clemens@nevrome.de + * + */ +#ifndef _RcppProgress_ETA_PROGRESS_BAR_HPP +#define _RcppProgress_ETA_PROGRESS_BAR_HPP + +#include +#include +#include +#include +#include + +#include "progress_bar.hpp" + +// for unices only +#if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__) +#include +#endif + +class ETAProgressBar: public ProgressBar{ + public: // ====== LIFECYCLE ===== + + /** + * Main constructor + */ + ETAProgressBar() { + _max_ticks = 50; + _finalized = false; + _timer_flag = true; + } + + ~ETAProgressBar() { + } + + public: // ===== main methods ===== + + void display() { + REprintf("0%% 10 20 30 40 50 60 70 80 90 100%%\n"); + REprintf("[----|----|----|----|----|----|----|----|----|----|\n"); + flush_console(); + } + + // update display + void update(float progress) { + + // stop if already finalized + if (_finalized) return; + + // start time measurement when update() is called the first time + if (_timer_flag) { + _timer_flag = false; + // measure start time + time(&start); + } else { + + // measure current time + time(&end); + + // calculate passed time and remaining time (in seconds) + double pas_time = std::difftime(end, start); + double rem_time = (pas_time / progress) * (1 - progress); + + // convert seconds to time string + std::string time_string = _time_to_string(rem_time); + + // create progress bar string + std::string progress_bar_string = _current_ticks_display(progress); + + // ensure overwriting of old time info + int empty_length = time_string.length(); + std::string empty_space = std::string(empty_length, ' '); + + // merge progress bar and time string + std::stringstream strs; + strs << "|" << progress_bar_string << "| " << time_string << empty_space; + std::string temp_str = strs.str(); + char const* char_type = temp_str.c_str(); + + // print: remove old and replace with new + REprintf("\r"); + REprintf("%s", char_type); + + // finalize display when ready + if(progress == 1) { + _finalize_display(); + } + } + } + + void end_display() { + update(1); + } + + protected: // ==== other instance methods ===== + + // convert double with seconds to time string + std::string _time_to_string(double seconds) { + + int time = (int) seconds; + + int hour = 0; + int min = 0; + int sec = 0; + + hour = time / 3600; + time = time % 3600; + min = time / 60; + time = time % 60; + sec = time; + + std::stringstream time_strs; + if (hour != 0) time_strs << hour << "h "; + if (min != 0) time_strs << min << "min "; + if (sec != 0) time_strs << sec << "s "; + std::string time_str = time_strs.str(); + + return time_str; + } + + // update the ticks display corresponding to progress + std::string _current_ticks_display(float progress) { + + int nb_ticks = _compute_nb_ticks(progress); + + std::string cur_display = _construct_ticks_display_string(nb_ticks); + + return cur_display; + } + + // construct progress bar display + std::string _construct_ticks_display_string(int nb) { + + std::stringstream ticks_strs; + for (int i = 0; i < (_max_ticks - 1); ++i) { + if (i < nb) { + ticks_strs << "*"; + } else { + ticks_strs << " "; + } + } + std::string tick_space_string = ticks_strs.str(); + + return tick_space_string; + } + + // finalize + void _finalize_display() { + if (_finalized) return; + + REprintf("\n"); + flush_console(); + _finalized = true; + } + + // compute number of ticks according to progress + int _compute_nb_ticks(float progress) { + return int(progress * _max_ticks); + } + + // N.B: does nothing on windows + void flush_console() { +#if !defined(WIN32) && !defined(__WIN32) && !defined(__WIN32__) + R_FlushConsole(); +#endif + } + + private: // ===== INSTANCE VARIABLES ==== + int _max_ticks; // the total number of ticks to print + bool _finalized; + bool _timer_flag; + time_t start,end; + +}; + +#endif diff --git a/src/integrate.cpp b/src/integrate.cpp index 7518847..33133fe 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -13,11 +13,11 @@ double f2(double x, double m, double b) { double integrate(double lower, double upper, int n, int fun, double m, double b){ double sum=0; - Rcpp::NumericVector x(n+1); - Rcpp::NumericVector y(n+1); + std::vector x(n+1); + std::vector y(n+1); double h = (upper-lower)/n; //get the width of the subintervals for (int i = 0; i <= n; i++) { //loop to evaluate x0,...xn and y0,...yn - x[i]=lower+i*h; //and store them in arrays + x[i]=lower+i*h; //and store them in arrays if ( fun == 1){ y[i]=f1(x[i], m, b); } else { diff --git a/src/integrate.o b/src/integrate.o index 58dbbe7..9a49d41 100644 Binary files a/src/integrate.o and b/src/integrate.o differ diff --git a/src/protoVS.dll b/src/protoVS.dll new file mode 100644 index 0000000..ca31fe4 Binary files /dev/null and b/src/protoVS.dll differ diff --git a/src/protoVS.so b/src/protoVS.so new file mode 100644 index 0000000..e9d675b Binary files /dev/null and b/src/protoVS.so differ diff --git a/src/rasterutils.cpp b/src/rasterutils.cpp index 5b7d4ba..b3467fd 100644 --- a/src/rasterutils.cpp +++ b/src/rasterutils.cpp @@ -67,48 +67,50 @@ Rcpp::IntegerMatrix colRowFromCell(const Rcpp::IntegerVector &cell, const int nc return result; } -Rcpp::NumericMatrix xyFromCell(Rcpp::S4 &raster, const Rcpp::IntegerVector &cell) { - RasterInfo ras(raster); +std::vector> xyFromCell(RasterInfo ras, const std::vector &cell) { + //RasterInfo ras(raster); int n = cell.size(); - Rcpp::NumericMatrix out(n,2); - std::fill( out.begin(),out.end(),NA_REAL ); + std::vector> out; + out.resize(n, std::vector(2, NA_REAL)); + for (int i = 0; i= ras.ncell)) continue; int row = cell[i] / ras.ncol; int col = cell[i] - (row * ras.ncol); - out(i,0) = ras.xmin + (col + 0.5) * ras.res; - out(i,1) = ras.ymax - (row + 0.5) * ras.res; + out[i][0] = ras.xmin + (col + 0.5) * ras.res; + out[i][1] = ras.ymax - (row + 0.5) * ras.res; } return out; } -Rcpp::NumericMatrix xyFromCell(Rcpp::S4 &raster, int cell) { - const Rcpp::IntegerVector vcell = {cell}; - return xyFromCell(raster, vcell); +std::vector> xyFromCell(RasterInfo ras, const int cell) { + std::vector vcell; + vcell.push_back(cell); + return xyFromCell(ras, vcell); } -Rcpp::IntegerVector cellFromXY(Rcpp::S4 &raster, Rcpp::NumericMatrix xy) { +std::vector cellFromXY(RasterInfo ras, const std::vector> &xy) { // size of x and y should be the same - RasterInfo ras(raster); + //RasterInfo ras(raster); - int size = xy.nrow(); - Rcpp::IntegerVector cells(size); + int s = xy.size(); + std::vector cells(s); double yr_inv = ras.nrow / (ras.ymax - ras.ymin); double xr_inv = ras.ncol / (ras.xmax - ras.xmin); - for (int i = 0; i < size; i++) { - int row = floor((ras.ymax - xy(i,1)) * yr_inv); + for (int i = 0; i < s; i++) { + int row = floor((ras.ymax - xy[i][1]) * yr_inv); // points in between rows go to the row below // except for the last row, when they must go up - if (xy(i,1) == ras.ymin) { + if (xy[i][1] == ras.ymin) { row = ras.nrow-1 ; } - int col = floor((xy(i,0) - ras.xmin) * xr_inv); + int col = floor((xy[i][0] - ras.xmin) * xr_inv); // as for rows above. Go right, except for last column - if (xy(i,0) == ras.xmax) { + if (xy[i][0] == ras.xmax) { col = ras.ncol - 1 ; } if (row < 0 || row >= ras.nrow || col < 0 || col >= ras.ncol) { @@ -122,7 +124,7 @@ Rcpp::IntegerVector cellFromXY(Rcpp::S4 &raster, Rcpp::NumericMatrix xy) { } int cellFromXY2(RasterInfo ras, double x, double y) { - + int cell; double yr_inv = ras.nrow / (ras.ymax - ras.ymin); diff --git a/src/rasterutils.h b/src/rasterutils.h index 6d2f923..bead471 100644 --- a/src/rasterutils.h +++ b/src/rasterutils.h @@ -2,14 +2,15 @@ #define RASTERUTILS #include +#include "rsinfo.h" extern void recycle(std::vector &x, std::vector &y); extern Rcpp::IntegerVector cellFromColRowSensitive(Rcpp::S4 &raster, Rcpp::IntegerVector &rcpp_col, Rcpp::IntegerVector &rcpp_row); extern Rcpp::IntegerVector cellFromColRow(const Rcpp::IntegerVector &x, const Rcpp::IntegerVector &y, const int ncol); extern Rcpp::IntegerMatrix colRowFromCell(const Rcpp::IntegerVector &cell, const int ncol); -extern Rcpp::NumericMatrix xyFromCell(Rcpp::S4 &raster, const Rcpp::IntegerVector &cell); -extern Rcpp::NumericMatrix xyFromCell(Rcpp::S4 &raster, int cell); -extern Rcpp::IntegerVector cellFromXY(Rcpp::S4 &raster, Rcpp::NumericMatrix xy); +extern std::vector> xyFromCell(RasterInfo ras, const std::vector &cell); +extern std::vector> xyFromCell(RasterInfo ras, const int cell); +extern std::vector cellFromXY(RasterInfo ras, const std::vector> &xy); extern int cellFromXY2(RasterInfo ras, double x, double y); #endif \ No newline at end of file diff --git a/src/rasterutils.o b/src/rasterutils.o index 7dd12df..247c39d 100644 Binary files a/src/rasterutils.o and b/src/rasterutils.o differ diff --git a/src/vgvi.cpp b/src/vgvi.cpp index cebe57f..c7e7281 100644 --- a/src/vgvi.cpp +++ b/src/vgvi.cpp @@ -1,418 +1,12 @@ #include +#include #include "rsinfo.h" #include "rasterutils.h" #include "bresenham.h" #include "integrate.h" -#include "vgvi.h" using namespace Rcpp; -// [[Rcpp::export]] -Rcpp::NumericVector VGVI_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - Rcpp::S4 &greenspace, const Rcpp::NumericVector &greenspace_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0, - const int fun, const double m, const double b) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)(radius/ras.res); - const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - - // Output vector - Rcpp::NumericVector output(input_cells.size(), NA_REAL); - - // Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells - // Will be used as a reference for all input points - const Rcpp::IntegerVector los_ref_vec = LoS_reference(x0_ref, y0_ref, r, nc_ref); - - // Main loop - for(int k=0; k dsm_values[input_cells[k]]) { - Rcpp::NumericVector max_tan_vec(r); - int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); - - for(int i = 0; i < (r*8); i++){ - // Re-use tangents calculated from prior LoS - int k_i = 0; - if(i > 0){ - for(int j = 0; j < r; j++){ - k_i = j; - if(los_ref_vec[i*r + j] != los_ref_vec[(i-1)*r + j]){ - break; - } - } - } - double max_tan = (k_i > 1) ? max_tan_vec[k_i-1] : -9999.0; - - - for(int j = k_i; j < r; j++){ - const int los_ref_cell = los_ref_vec[i*r + j]; - if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ - // Project reference cells to actual cell value - const int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); - const int row = trunc(cell/ras.ncol); - const int col = cell - (row * ras.ncol); - const int dcol = abs(col-x0_o[k]); - - if(!(cell<0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - const double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - const double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - if(this_tan > max_tan){ - max_tan = this_tan; - viewshed[los_ref_cell] = cell+1; - } - } - max_tan_vec[j] = max_tan; - } else { - break; - } - } - } - } - - Rcpp::IntegerVector viewshed_na = Rcpp::na_omit(viewshed); - - - // B: Greenspace Visibility Index - const int cs = viewshed_na.size(); - - // Get XY coordinates of visible cells - const Rcpp::NumericMatrix dsm_xy = xyFromCell(dsm, viewshed_na); - const Rcpp::NumericMatrix xy0 = xyFromCell(dsm, input_cells[k]); - - // Calculate distance - Rcpp::IntegerVector dxy(cs); - for(int i = 0; i < cs; i++){ - double d = sqrt( ((xy0(0,0) - dsm_xy(i,0))*(xy0(0,0) - dsm_xy(i,0))) + ((xy0(0,1) - dsm_xy(i,1))*(xy0(0,1) - dsm_xy(i,1))) ); - int di = round(d); - if(di < 1) { - dxy[i] = 1; - } else { - dxy[i] = di; - } - } - - - // Intersect XY with greenspace mask - const Rcpp::IntegerVector greenspace_cells = cellFromXY(greenspace, dsm_xy); - Rcpp::NumericVector greenspace_cell_values(cs, 0.0); - for(int i = 0; i < cs; i++){ - int gs_cell = greenspace_cells[i]; - if(!Rcpp::IntegerVector::is_na(gs_cell)){ - double gs = greenspace_values[gs_cell]; - greenspace_cell_values[i] = Rcpp::NumericVector::is_na(gs) ? 0.0 : gs; - } - } - - // Get number of green visible cells and total visible cells per distance - const double max_d = max(dxy); - const Rcpp::IntegerVector dxy_seq = seq(1, max_d); - const int n = dxy_seq.size(); - Rcpp::IntegerVector visibleTotal(n); - Rcpp::IntegerVector visibleGreen(n); - - for(int i = 0; i < cs; i++){ - int this_dxy = dxy[i]; - visibleTotal[this_dxy-1] += 1; - visibleGreen[this_dxy-1] += greenspace_cell_values[i]; - } - for(int i = 0; i < n; i++){ - if(visibleTotal[i] == 0){ - visibleTotal[i] = 1; - } - } - - - // Proportion of visible green cells - if(max_d == 1){ - output[k] = visibleGreen[0]/visibleTotal[0]; - } - Rcpp::NumericVector raw_GVI(n); - for(int i = 0; i < n; i++){ - raw_GVI[i] = (double)visibleGreen[i] / visibleTotal[i]; - } - - // Normalize distance - Rcpp::NumericVector nDxy(n); - for(int i = 0; i < n; i++){ - nDxy[i] = dxy_seq[i] / (double)radius; //max_d; // TAKE radius instead and save max_d / r for upper limit - } - - // Calculate weights by taking the proportion of the integral of each step from the integral of the whole area - //double big_integral = integrate(0, 1, 200, fun, m, b); - const double min_dist = min(nDxy); - Rcpp::NumericVector decayWeights(n); - for(int i = 0; i < n; i++){ - double d = nDxy[i]; - decayWeights[i] = integrate(d-min_dist, d, 200, fun, m, b); - } - const double big_integral = sum(decayWeights); - - // Proportion of visible green - double vgvi_sum = 0.0; - for(int i = 0; i < n; i++){ - vgvi_sum += raw_GVI[i] * (decayWeights[i]/big_integral); - } - - output[k] = vgvi_sum; - } - - return output; -} - - -// [[Rcpp::export]] -Rcpp::IntegerVector viewshed_cpp_sum1(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)(radius/ras.res); - const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - - // Output vector - Rcpp::IntegerVector output(input_cells.size(), NA_INTEGER); - - // Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells - // Will be used as a reference for all input points - const Rcpp::IntegerVector los_ref_vec = LoS_reference(x0_ref, y0_ref, r, nc_ref); - - // Main loop - for(int k=0; k dsm_values[input_cells[k]]){ - // 3. Viewshed Analysis - Rcpp::NumericVector max_tan_vec(r); - int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); - - for(int i = 0; i < (r*8); i++){ - double max_tan = -9999.0; - - for(int j = 0; j < r; j++){ - const int los_ref_cell = los_ref_vec[i*r + j]; - if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ - // Project reference cells to actual cell value - const int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); - const int row = trunc(cell/ras.ncol); - const int col = cell - (row * ras.ncol); - const int dcol = abs(col-x0_o[k]); - - if(!(cell<0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - const double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - const double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - if(this_tan > max_tan){ - max_tan = this_tan; - viewshed[los_ref_cell] = cell+1; - } - } - max_tan_vec[j] = max_tan; - } else { - break; - } - } - } - } - Rcpp::IntegerVector viewshed_na = Rcpp::na_omit(viewshed); - output[k] = viewshed_na.size(); - } - - - - return output; -} - - -// [[Rcpp::export]] -Rcpp::IntegerVector viewshed_cpp_sum2(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)(radius/ras.res); - const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - - // Output vector - Rcpp::IntegerVector output(input_cells.size(), NA_INTEGER); - - // Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells - // Will be used as a reference for all input points - const Rcpp::IntegerVector los_ref_vec = LoS_reference(x0_ref, y0_ref, r, nc_ref); - const Rcpp::IntegerVector los_start = shared_LoS(r, los_ref_vec); - - // Main loop - for(int k=0; k dsm_values[input_cells[k]]){ - // 3. Viewshed Analysis - Rcpp::NumericVector max_tan_vec(r); - int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); - - for(int i = 0; i < (r*8); i++){ - // Re-use tangents calculated for prior LoS - int k_i = los_start[i]; - double max_tan = (k_i > 1) ? max_tan_vec[k_i-1] : -9999.0; - - - for(int j = k_i; j < r; j++){ - const int los_ref_cell = los_ref_vec[i*r + j]; - if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ - // Project reference cells to actual cell value - const int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); - const int row = trunc(cell/ras.ncol); - const int col = cell - (row * ras.ncol); - const int dcol = abs(col-x0_o[k]); - - if(!(cell<0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - const double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - const double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - if(this_tan > max_tan){ - max_tan = this_tan; - viewshed[los_ref_cell] = cell+1; - } - } - max_tan_vec[j] = max_tan; - } else { - break; - } - } - } - } - Rcpp::IntegerVector viewshed_na = Rcpp::na_omit(viewshed); - output[k] = viewshed_na.size(); - } - - - - return output; -} - -// [[Rcpp::export]] -Rcpp::IntegerVector viewshed_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)round(radius/ras.res); - const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - - - // Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells - // Will be used as a reference for all input points - const Rcpp::IntegerVector los_ref_vec = LoS_reference(x0_ref, y0_ref, r, nc_ref); - - int k=0; - - // Outpur: x0/y0 is always visible - Rcpp::IntegerVector viewshed(nc_ref*nr_ref, NA_INTEGER); - viewshed[c0_ref] = input_cells[k]+1; - - if(h0[k] > dsm_values[input_cells[k]]){ - // 3. Viewshed Analysis - Rcpp::NumericVector max_tan_vec(r); - int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); - - for(int i = 0; i < (r*8); i++){ - // Re-use tangents calculated for prior LoS - int k_i = 0; - if(i > 0){ - for(int j = 0; j < r; j++){ - k_i = j; - if(los_ref_vec[i*r + j] != los_ref_vec[(i-1)*r + j]){ - break; - } - } - } - double max_tan = (k_i > 1) ? max_tan_vec[k_i-1] : -9999.0; - - - for(int j = k_i; j < r; j++){ - const int los_ref_cell = los_ref_vec[i*r + j]; - if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ - // Project reference cells to actual cell value - const int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); - const int row = trunc(cell/ras.ncol); - const int col = cell - (row * ras.ncol); - const int dcol = abs(col-x0_o[k]); - - if(!(cell<0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - const double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - const double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - if(this_tan > max_tan){ - max_tan = this_tan; - viewshed[los_ref_cell] = cell+1; - } - } - max_tan_vec[j] = max_tan; - } else { - break; - } - } - } - } - - return Rcpp::na_omit(viewshed); -} - #ifdef _OPENMP #include @@ -420,527 +14,215 @@ Rcpp::IntegerVector viewshed_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_v // [[Rcpp::plugins(openmp)]] // [[Rcpp::depends(RcppProgress)]] #include +#include "eta_progress_bar.hpp" // [[Rcpp::export]] -Rcpp::IntegerMatrix viewshed_distance_analysis_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0, - const int ncores=1, const bool display_progress=false) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)(radius/ras.res); - const int nc_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - const int s = input_cells.size(); - - // Output - Rcpp::IntegerMatrix output( r*s,3 ); - for(int k = 0; k < s; k++) { - for(int d = 0; d dsm_values[input_cells[k]]){ - // Viewshed Analysis - const int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); - -#if defined(_OPENMP) - omp_set_num_threads(ncores); -#pragma omp parallel for shared(output) -#endif - for(int i = 0; i < (r*8); i++){ - if (!p.is_aborted()) { - Progress::check_abort(); - - double max_tan = -9999.0; - - for(int j = 0; j < r; j++){ - const int los_ref_cell = los_ref_vec[i*r + j]; - if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ - // Project reference cells to actual cell value - int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); - int row = trunc(cell/ras.ncol); - int col = cell - (row * ras.ncol); - int dcol = abs(col-x0_o[k]); - - if(!(cell < 0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - int d = round(distance_traveled) - 1; - d = (d < 0) ? 0 : d; - d = (d >= r) ? (r-1) : d; - - output( (d*s + k),1 ) += 1; - //output( d,1 ) += 1; - - if(this_tan > max_tan){ - max_tan = this_tan; - - output( (d*s + k),2 ) += 1; - //output( d,2 ) += 1; - } - } - } else { - break; - } - } - } - } - } - p.increment(); - } - - return output; -} - -// [[Rcpp::export]] -Rcpp::IntegerMatrix viewshed_and_greenness_distance_analysis_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - Rcpp::S4 &greenspace, const Rcpp::NumericVector &greenspace_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0, - const int ncores=1, const bool display_progress=false) +std::vector VGVI_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, + Rcpp::S4 &greenspace, const Rcpp::NumericVector &greenspace_values, + const Rcpp::IntegerVector &x0, const Rcpp::IntegerVector &y0, + const Rcpp::NumericVector &h0, const int radius, + const int fun, const double m, const double b, + const int ncores=1, const bool display_progress=false) { - // Cells from x0,y0 + // Cells from x0, y0 Rcpp::IntegerVector x0_o = x0-1; Rcpp::IntegerVector y0_o = y0-1; const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); // Basic raster information - const RasterInfo ras(dsm); + const RasterInfo dsm_ras(dsm); const RasterInfo gs_ras(greenspace); // Parameters - const int r = (int)(radius/ras.res); - const int nc_ref = (2*r)+1; + const int r = (int)round(radius/dsm_ras.res); + const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; const int x0_ref = r, y0_ref = x0_ref; const int c0_ref = y0_ref*nc_ref + x0_ref; - const int s = input_cells.size(); - - // Output - // v1: Distance - // V2: Total number of cells per distance - // V3: Number of all visible cells per distance - // V4: Total number of green cells per distance - // V5: Number of green visible cells per distance - Rcpp::IntegerMatrix output( r*s,5 ); - for(int k = 0; k < s; k++) { - for(int d = 0; d output(input_cells.size(), NA_REAL); - // Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells + // Protoptype of Line of Sight (LoS) paths: // Will be used as a reference for all input points const Rcpp::IntegerVector los_ref_vec = LoS_reference(x0_ref, y0_ref, r, nc_ref); + const Rcpp::IntegerVector los_start = shared_LoS(r, los_ref_vec); // Progress bar - Progress p(s, display_progress); + ETAProgressBar pb_ETA; + Progress pb(input_cells.size(), display_progress, pb_ETA); // Main loop #if defined(_OPENMP) omp_set_num_threads(ncores); -#pragma omp parallel for shared(output) +#pragma omp parallel for schedule(dynamic) shared(output) #endif - for(int k = 0; k < s; k++) - { - if(h0[k] > dsm_values[input_cells[k]]){ - // Viewshed Analysis - const int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); + for(int k = 0; k < input_cells.size(); ++k){ + if ( !pb.is_aborted() ) { + Progress::check_abort(); - for(int i = 0; i < (r*8); i++){ - if (!p.is_aborted()) { - Progress::check_abort(); - - double max_tan = -9999.0; + const int this_input_cell = input_cells[k]; + + // Viewshed + std::vector viewshed(nc_ref*nr_ref, NA_INTEGER); + + // Viewshed at x0/y0 is always visible + viewshed[c0_ref] = input_cells[k]+1; + + // A: Viewshed Analysis + // Eye-level height at x0/y0 > DSM height? + if(h0[k] > dsm_values[this_input_cell]){ + + // Parameter for projecting reference cell to true cell values + const int x = this_input_cell - c0_ref - r*(dsm_ras.ncol-nc_ref); + + // Vector for storing max tangent from previous LoS + std::vector max_tan_vec(r, -9999.0); + + // Iterate over all LoS paths + for(int i = 0; i < (r*8); ++i){ + // Re-use tangents calculated from prior LoS or assign -9999.0 + const int k_i = los_start[i]; + double max_tan = (k_i > 1) ? max_tan_vec[k_i-1] : -9999.0; - for(int j = 0; j < r; j++){ + // Iterate over all cells of this LoS path starting at k_i + for(int j = k_i; j < r; ++j){ + + // This LoS path cell const int los_ref_cell = los_ref_vec[i*r + j]; - if(los_ref_cell != NA_INTEGER){ //!Rcpp::IntegerVector::is_na(los_ref_cell) - // Project reference cells to actual cell value - int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); - - // row/col - int row = trunc(cell/ras.ncol); - int col = cell - (row * ras.ncol); - int dcol = abs(col-x0_o[k]); - - // X/Y coordinates - double x_j = ras.xmin + (col + 0.5) * ras.res; - double y_j = ras.ymax - (row + 0.5) * ras.res; - - // Cell of greenspace raster - // OpenMP can't call cellFromXY2 function without causing R to crash. - // Therefore, I'll calcualte the cell here manually - int cell_gs = cellFromXY2(gs_ras, x_j, y_j); + + if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ + // Project reference cell to true cell value + const int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(dsm_ras.ncol-nc_ref); + // DSM height at this LoS path cell + const double h_cell = dsm_values[cell]; + // Test if this LoS path cell is within DSM raster + const int row = trunc(cell/dsm_ras.ncol); + const int col = cell - (row * dsm_ras.ncol); + const int dcol = abs(col-x0_o[k]); - if(!(cell < 0 || cell > ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - int d = round(distance_traveled) - 1; - d = (d < 0) ? 0 : d; - d = (d >= r) ? (r-1) : d; - - output( (d*s + k),1 ) += 1; - if(cell_gs != NA_INTEGER && greenspace_values[cell_gs] == 1){ - output( (d*s + k),3 ) += 1; - } + if(!(cell<0 || cell > dsm_ras.ncell || Rcpp::NumericVector::is_na(h_cell) || dcol>r)){ + // Compute tangent of x0/y0 (observer location) and this LoS path cell + const double distance_traveled = sqrt( + (x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row) + ); + const double this_tan = (h_cell - h0[k]) / (distance_traveled); + // Update viewshed and max tangent if(this_tan > max_tan){ max_tan = this_tan; - - output( (d*s + k),2 ) += 1; - if(cell_gs != NA_INTEGER && greenspace_values[cell_gs] == 1){ - output( (d*s + k),4 ) += 1; - } + viewshed[los_ref_cell] = cell+1; } } + + max_tan_vec[j] = max_tan; } else { break; } } } } - } - p.increment(); - } - - return output; -} - -// Alternative version of viewshed. Slightly faster for smaller radius, but more RAM intensive. -Rcpp::IntegerVector viewshed_cpp2(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)(radius/ras.res); - const int l = r+1; - const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - int k=0; - - - // 1. Reference Bresenham Lines for the first eigth of circle - const Rcpp::IntegerMatrix bh_mat = bresenham_map(x0_ref, y0_ref, r, nc_ref); - - // 2. Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells - Rcpp::IntegerMatrix los_ref_mat = Rcpp::IntegerMatrix(l*8 - 8, r); - - - for(int i=0; i> dsm_xy = xyFromCell(dsm_ras, viewshed); + const std::vector> xy0 = xyFromCell(dsm_ras, input_cells[k]); + + // Calculate distance + std::vector dxy(cs); + for(int i = 0; i < cs; i++){ + double d = sqrt( ((xy0[0][0] - dsm_xy[i][0])*(xy0[0][0] - dsm_xy[i][0])) + ((xy0[0][1] - dsm_xy[i][1])*(xy0[0][1] - dsm_xy[i][1])) ); + int di = round(d); + if(di < 1) { + dxy[i] = 1; + } else { + dxy[i] = di; } } - } - } - - - // A: Tangents Map - const Rcpp::NumericVector tan_vec = tangentsMap(x0_o[k], y0_o[k], h0[k], ras.nrow, ras.ncol, r, dsm_values); - - // B: Visibility - Rcpp::LogicalVector visibility_vec(nc_ref*nr_ref); - - for(int j = 0; j < (l*8 - 8); j++){ - - // Compute tangents of LoS_vec and update visibility_vec - LoS(los_ref_mat, tan_vec, visibility_vec, j, r); - } - - - // Project reference cells to actual cell value - Rcpp::IntegerVector viewshed(nc_ref*nr_ref, NA_INTEGER); - viewshed[c0_ref] = input_cells[k]+1; - int x = input_cells[k] - c0_ref - r*(ras.ncol-nc_ref); - for(int j = 0; j < visibility_vec.size(); j++){ - if(visibility_vec[j]){ - const int cell = x + j + trunc(j/nc_ref)*(ras.ncol-nc_ref); - const int row = cell - (trunc(cell/ras.ncol) * ras.ncol); - const int drow = abs(row-x0_o[k]); - viewshed[j] = (cell<0 || cell > ras.ncell || - Rcpp::NumericVector::is_na(dsm_values[cell]) || - drow>r) ? NA_INTEGER : (cell+1); - } - } - - return Rcpp::na_omit(viewshed); -} - - - -Rcpp::NumericVector VGVI_cpp2(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - Rcpp::S4 &greenspace, const Rcpp::NumericVector &greenspace_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0, - const int fun, const double m, const double b, - const int ncores=1) -{ - // Cells from x0,y0 - Rcpp::IntegerVector x0_o = x0-1; - Rcpp::IntegerVector y0_o = y0-1; - const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); - - // Basic raster information - const RasterInfo ras(dsm); - - // Parameters - const int r = (int)(radius/ras.res); - const int l = r+1; - const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; - const int x0_ref = r, y0_ref = x0_ref; - const int c0_ref = y0_ref*nc_ref + x0_ref; - - // Output vector - Rcpp::NumericVector output(input_cells.size(), NA_REAL); - - // 1. Reference Bresenham Lines for the first eigth of circle - const Rcpp::IntegerMatrix bh_mat = bresenham_map(x0_ref, y0_ref, r, nc_ref); - - // 2. Reference matrix (los_mat): Project reference Bresenham Lines to all perimeter cells - // Will be used as a reference for all input points - Rcpp::IntegerVector los_ref_vec = Rcpp::IntegerVector(r*8 * r); - - for(int i=0; i greenspace_cells = cellFromXY(gs_ras, dsm_xy); + std::vector greenspace_cell_values(cs, 0.0); + for(int i = 0; i < cs; i++){ + int gs_cell = greenspace_cells[i]; + if(!Rcpp::IntegerVector::is_na(gs_cell)){ + double gs = greenspace_values[gs_cell]; + greenspace_cell_values[i] = Rcpp::NumericVector::is_na(gs) ? 0.0 : gs; } - } - } - } - - // 3. Main loop - for(int k=0; k ras.ncell || Rcpp::NumericVector::is_na(dsm_values[cell]) || dcol>r)){ - // Compute tangents of LoS_vec and update output - const double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); - const double this_tan = (dsm_values[cell] - h0[k]) / (distance_traveled); - - if(this_tan > max_tan){ - max_tan = this_tan; - viewshed[los_ref_cell] = cell+1; - } - } - } else { - break; - } + // Get number of green visible cells and total visible cells per distance + const double max_d = *std::max_element(dxy.begin(), dxy.end()) + 0.0; + + std::vector dxy_seq(max_d); + std::iota(dxy_seq.begin(), dxy_seq.end(), 1); + + const int n = dxy_seq.size(); + + std::vector visibleTotal(n); + std::vector visibleGreen(n); + + for(int i = 0; i < cs; i++){ + int this_dxy = (int)dxy[i]; + visibleTotal[this_dxy-1] += 1; + visibleGreen[this_dxy-1] += greenspace_cell_values[i]; + } + for(int i = 0; i < n; i++){ + if(visibleTotal[i] == 0){ + visibleTotal[i] = 1; } } - } - Rcpp::IntegerVector viewshed_na = Rcpp::na_omit(viewshed); - - - // B: Greenspace Visibility Index - const int cs = viewshed_na.size(); - - // Get XY coordinates of visible cells - const Rcpp::NumericMatrix dsm_xy = xyFromCell(dsm, viewshed_na); - const Rcpp::NumericMatrix xy0 = xyFromCell(dsm, input_cells[k]); - - // Calculate distance - Rcpp::IntegerVector dxy(cs); - for(int i = 0; i < cs; i++){ - double d = sqrt( ((xy0(0,0) - dsm_xy(i,0))*(xy0(0,0) - dsm_xy(i,0))) + ((xy0(0,1) - dsm_xy(i,1))*(xy0(0,1) - dsm_xy(i,1))) ); - int di = round(d); - if(di < 1) { - dxy[i] = 1; - } else { - dxy[i] = di; + + // Proportion of visible green cells + if(max_d == 1.0){ + output[k] = visibleGreen[0]/visibleTotal[0]; + continue; } - } - - - // Intersect XY with greenspace mask - const Rcpp::IntegerVector greenspace_cells = cellFromXY(greenspace, dsm_xy); - Rcpp::NumericVector greenspace_cell_values(cs, 0.0); - for(int i = 0; i < cs; i++){ - int gs_cell = greenspace_cells[i]; - if(!Rcpp::IntegerVector::is_na(gs_cell)){ - double gs = greenspace_values[gs_cell]; - greenspace_cell_values[i] = Rcpp::NumericVector::is_na(gs) ? 0.0 : gs; + std::vector raw_GVI(n); + for(int i = 0; i < n; i++){ + raw_GVI[i] = (double)visibleGreen[i] / visibleTotal[i]; } - } - - // Get number of green visible cells and total visible cells per distance - const double max_d = max(dxy); - const Rcpp::IntegerVector dxy_seq = seq(1, max_d); - const int n = dxy_seq.size(); - Rcpp::IntegerVector visibleTotal(n); - Rcpp::IntegerVector visibleGreen(n); - - for(int i = 0; i < cs; i++){ - int this_dxy = dxy[i]; - visibleTotal[this_dxy-1] += 1; - visibleGreen[this_dxy-1] += greenspace_cell_values[i]; - } - for(int i = 0; i < n; i++){ - if(visibleTotal[i] == 0){ - visibleTotal[i] = 1; + + + // Normalize distance + std::vector nDxy(n); + for(int i = 0; i < n; i++){ + nDxy[i] = dxy_seq[i] / (double)radius; } + + + // Calculate weights by taking the proportion of the integral + // of each step from the integral of the whole area + const double min_dist = *std::min_element(nDxy.begin(), nDxy.end()); + std::vector decayWeights(n); + for(int i = 0; i < n; i++){ + double d = nDxy[i]; + decayWeights[i] = integrate(d-min_dist, d, 200, fun, m, b); + } + const double big_integral = std::accumulate(decayWeights.begin(), decayWeights.end(), 0.0); + + + // Proportion of visible green + double vgvi_sum = 0.0; + for(int i = 0; i < n; i++){ + vgvi_sum += raw_GVI[i] * (decayWeights[i]/big_integral); + } + + output[k] = vgvi_sum; + pb.increment(); } - - - // Proportion of visible green cells - if(max_d == 1){ - output[k] = visibleGreen[0]/visibleTotal[0]; - } - Rcpp::NumericVector raw_GVI(n); - for(int i = 0; i < n; i++){ - raw_GVI[i] = (double)visibleGreen[i] / visibleTotal[i]; - } - - // Normalize distance - Rcpp::NumericVector nDxy(n); - for(int i = 0; i < n; i++){ - nDxy[i] = dxy_seq[i] / (double)radius; //max_d; // TAKE radius instead and save max_d / r for upper limit - } - - // Calculate weights by taking the proportion of the integral of each step from the integral of the whole area - //double big_integral = integrate(0, 1, 200, fun, m, b); - const double min_dist = min(nDxy); - Rcpp::NumericVector decayWeights(n); - for(int i = 0; i < n; i++){ - double d = nDxy[i]; - decayWeights[i] = integrate(d-min_dist, d, 200, fun, m, b); - } - const double big_integral = sum(decayWeights); - - // Proportion of visible green - double vgvi_sum = 0.0; - for(int i = 0; i < n; i++){ - vgvi_sum += raw_GVI[i] * (decayWeights[i]/big_integral); - } - - output[k] = vgvi_sum; } + return output; -} - - - - +} \ No newline at end of file diff --git a/src/vgvi.h b/src/vgvi.h deleted file mode 100644 index e3911f3..0000000 --- a/src/vgvi.h +++ /dev/null @@ -1,12 +0,0 @@ -#ifndef VGVI -#define VGVI - -#include - -extern Rcpp::NumericVector VGVI_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, - Rcpp::S4 &greenspace, const Rcpp::NumericVector &greenspace_values, - const Rcpp::IntegerVector x0, const Rcpp::IntegerVector y0, - const int radius, const Rcpp::NumericVector & h0, - const int fun, const double m, const double b); - -#endif \ No newline at end of file diff --git a/src/vgvi.o b/src/vgvi.o index 1820107..57c8c98 100644 Binary files a/src/vgvi.o and b/src/vgvi.o differ diff --git a/src/viewshed.cpp b/src/viewshed.cpp new file mode 100644 index 0000000..9c82815 --- /dev/null +++ b/src/viewshed.cpp @@ -0,0 +1,159 @@ +#include +#include +#include "rsinfo.h" +#include "rasterutils.h" +#include "bresenham.h" + +using namespace Rcpp; + +// [[Rcpp::depends(RcppProgress)]] +#include + +std::vector LoS_tangent(const std::vector &shared_LoS, + const Rcpp::NumericVector &dsm_values, + const int p, const int nc_ref, const int nc_ras, const int ncell_ras, + const int x0, const int y0, const int h0, const int r){ + const int n = shared_LoS.size(); + std::vector out(n, -9999.0); + + for(int i = 0; i < n; ++i){ + // 1. Get LoS cell + const int los_cell = shared_LoS[i]; + + // 2. Project LoS cell + const int cell = p + los_cell + trunc(los_cell/nc_ref)*(nc_ras-nc_ref); + + // 3. DSM height at this LoS path cell + const double h_cell = dsm_values[cell]; + + // 4. Test if this LoS path cell is within DSM raster + const int row = trunc(cell/nc_ras); + const int col = cell - (row * nc_ras); + const int dcol = abs(col-x0); + + if(!(cell<0 || cell > ncell_ras || Rcpp::NumericVector::is_na(h_cell) || dcol>r)){ + // Compute tangent of x0/y0 (observer location) and this LoS path cell + const double distance_traveled = sqrt((x0 - col)*(x0 - col) + (y0 - row)*(y0 - row)); + const double this_tan = (h_cell - h0) / (distance_traveled); + + // Update viewshed and max tangent + out[i] = this_tan; + } + } + + return out; +} + + + +// [[Rcpp::export]] +std::vector viewshed_cpp(Rcpp::S4 &dsm, const Rcpp::NumericVector &dsm_values, + const Rcpp::IntegerVector &x0, const Rcpp::IntegerVector &y0, + const Rcpp::NumericVector &h0, const int radius, + const bool display_progress=false) +{ + // Cells from x0, y0 + Rcpp::IntegerVector x0_o = x0-1; + Rcpp::IntegerVector y0_o = y0-1; + const Rcpp::IntegerVector input_cells = Rcpp::na_omit(cellFromColRowSensitive(dsm, x0_o, y0_o)); + + // Basic raster information + const RasterInfo ras(dsm); + + // Parameters + const int r = (int)round(radius/ras.res); + const int nc_ref = (2*r)+1, nr_ref = (2*r)+1; + const int x0_ref = r, y0_ref = x0_ref; + const int c0_ref = y0_ref*nc_ref + x0_ref; + + // Output vector + std::vector output; + + // Protoptype of Line of Sight (LoS) paths: + // Will be used as a reference for all input points + const Rcpp::IntegerVector los_ref_vec = LoS_reference(x0_ref, y0_ref, r, nc_ref); + const Rcpp::IntegerVector los_start = shared_LoS(r, los_ref_vec); + + // Progress bar + Progress pb(input_cells.size(), display_progress); + + // Main loop + for(int k = 0; k < input_cells.size(); ++k){ + + const int this_input_cell = input_cells[k]; + + // Viewshed + std::vector viewshed(nc_ref*nr_ref, NA_INTEGER); + + // Viewshed at x0/y0 is always visible + viewshed[c0_ref] = input_cells[k]+1; + + if (!pb.is_aborted()) { + + Progress::check_abort(); + pb.increment(); + + // Eye-level height at x0/y0 > DSM height? + if(h0[k] > dsm_values[this_input_cell]){ + + // Parameter for projecting reference cell to true cell values + const int x = this_input_cell - c0_ref - r*(ras.ncol-nc_ref); + + // Vector for storing max tangent from previous LoS + std::vector max_tan_vec(r, -9999.0); + + // Iterate over all LoS paths + for(int i = 0; i < (r*8); ++i){ + // Re-use tangents calculated from prior LoS or assign -9999.0 + const int k_i = los_start[i]; + double max_tan = (k_i > 1) ? max_tan_vec[k_i-1] : -9999.0; + + // Iterate over all cells of this LoS path starting at k_i + for(int j = k_i; j < r; ++j){ + + // This LoS path cell + const int los_ref_cell = los_ref_vec[i*r + j]; + + if(!Rcpp::IntegerVector::is_na(los_ref_cell)){ + + // Project reference cell to true cell value + const int cell = x + los_ref_cell + trunc(los_ref_cell/nc_ref)*(ras.ncol-nc_ref); + + // DSM height at this LoS path cell + const double h_cell = dsm_values[cell]; + + // Test if this LoS path cell is within DSM raster + const int row = trunc(cell/ras.ncol); + const int col = cell - (row * ras.ncol); + const int dcol = abs(col-x0_o[k]); + + if(!(cell<0 || cell > ras.ncell || Rcpp::NumericVector::is_na(h_cell) || dcol>r)){ + // Compute tangent of x0/y0 (observer location) and this LoS path cell + const double distance_traveled = sqrt((x0_o[k] - col)*(x0_o[k] - col) + (y0_o[k] - row)*(y0_o[k] - row)); + const double this_tan = (h_cell - h0[k]) / (distance_traveled); + + // Update viewshed and max tangent + if(this_tan > max_tan){ + max_tan = this_tan; + viewshed[los_ref_cell] = cell+1; + } + } + + max_tan_vec[j] = max_tan; + } else { + break; + } + } + } + } + } + + viewshed.erase(std::remove_if(std::begin(viewshed), + std::end(viewshed), + [](const int& value) { return Rcpp::IntegerVector::is_na(value); }), + std::end(viewshed)); + output = viewshed; + } + + return output; +} \ No newline at end of file diff --git a/src/viewshed.o b/src/viewshed.o new file mode 100644 index 0000000..01c0176 Binary files /dev/null and b/src/viewshed.o differ