From bdc26667bdd5d02314b697423a5c101b6c9e6c4c Mon Sep 17 00:00:00 2001 From: Jan Caha Date: Mon, 16 Dec 2019 10:40:06 +0000 Subject: [PATCH] version 0.5.0 --- DESCRIPTION | 31 ++ LICENSE | 2 + MD5 | 26 ++ NAMESPACE | 38 +++ NEWS.md | 13 + R/RcppExports.R | 7 + R/create_grid.R | 90 +++++ R/create_raster.R | 42 +++ R/kde.R | 182 +++++++++++ R/utils-pipe.R | 11 + R/utils.R | 85 +++++ R/zzz.R | 3 + README.md | 40 +++ build/vignette.rds | Bin 0 -> 210 bytes inst/doc/SpatialKDE.R | 51 +++ inst/doc/SpatialKDE.Rmd | 100 ++++++ inst/doc/SpatialKDE.html | 381 ++++++++++++++++++++++ man/create_grid.Rd | 57 ++++ man/create_raster.Rd | 30 ++ man/kde.Rd | 65 ++++ man/pipe.Rd | 12 + src/RcppExports.cpp | 33 ++ src/kde.cpp | 136 ++++++++ tests/testthat.R | 4 + tests/testthat/test-create_grids_raster.R | 50 +++ tests/testthat/test-kde.R | 62 ++++ vignettes/SpatialKDE.Rmd | 100 ++++++ 27 files changed, 1651 insertions(+) create mode 100644 DESCRIPTION create mode 100644 LICENSE create mode 100644 MD5 create mode 100644 NAMESPACE create mode 100644 NEWS.md create mode 100644 R/RcppExports.R create mode 100644 R/create_grid.R create mode 100644 R/create_raster.R create mode 100644 R/kde.R create mode 100644 R/utils-pipe.R create mode 100644 R/utils.R create mode 100644 R/zzz.R create mode 100644 README.md create mode 100644 build/vignette.rds create mode 100644 inst/doc/SpatialKDE.R create mode 100644 inst/doc/SpatialKDE.Rmd create mode 100644 inst/doc/SpatialKDE.html create mode 100644 man/create_grid.Rd create mode 100644 man/create_raster.Rd create mode 100644 man/kde.Rd create mode 100644 man/pipe.Rd create mode 100644 src/RcppExports.cpp create mode 100644 src/kde.cpp create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test-create_grids_raster.R create mode 100644 tests/testthat/test-kde.R create mode 100644 vignettes/SpatialKDE.Rmd diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..685fc41 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,31 @@ +Package: SpatialKDE +Type: Package +Title: Kernel Density Estimation for Spatial Data +Version: 0.5.0 +URL: https://jancaha.github.io/SpatialKDE/index.html, + https://github.com/JanCaha/SpatialKDE +Authors@R: + person(given = "Jan", + family = "Caha", + role = c("aut", "cre"), + email = "jan.caha@outlook.com", + comment = c(ORCID = "0000-0003-0165-0606")) +Description: Calculate Kernel Density Estimation (KDE) for spatial data. + The algorithm is inspired by the tool 'Heatmap' from 'QGIS'. The method is described by: + Hart, T., Zandbergen, P. (2014) , + Nelson, T. A., Boots, B. (2008) , + Chainey, S., Tompson, L., Uhlig, S.(2008) . +License: MIT + file LICENSE +Encoding: UTF-8 +LazyData: true +RoxygenNote: 7.0.2 +VignetteBuilder: knitr +LinkingTo: Rcpp +Imports: Rcpp, sf, dplyr, glue, magrittr, rlang, methods, raster +Suggests: tmap, sp, knitr, testthat (>= 2.1.0) +NeedsCompilation: yes +Packaged: 2019-12-16 08:14:56 UTC; Caha +Author: Jan Caha [aut, cre] () +Maintainer: Jan Caha +Repository: CRAN +Date/Publication: 2019-12-16 11:40:06 UTC diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..8f9738f --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +YEAR: 2019 +COPYRIGHT HOLDER: Jan Caha diff --git a/MD5 b/MD5 new file mode 100644 index 0000000..aa1a0f0 --- /dev/null +++ b/MD5 @@ -0,0 +1,26 @@ +45dac6d9b8643dcc67e022fdc27f2026 *DESCRIPTION +82179e94a9c57510eeb0dfa196cd97cb *LICENSE +f9384b67978f362e89b618ef0b1b5a48 *NAMESPACE +156fd5b034e4172f87effa1677d5cc1b *NEWS.md +46ff90eda42db6d37a5a20a97cecc8d1 *R/RcppExports.R +e11ac6cfc676e4a3d0ebc35c4492e8de *R/create_grid.R +bb39c117653c3bf19fcde4d7aea1dd6c *R/create_raster.R +eefb7801ac5594fea1f9bad03ccfc4a8 *R/kde.R +ad3f3c64fcdb96e66b1859bd35c48478 *R/utils-pipe.R +ea034c24e4e96a741e65721835b733cf *R/utils.R +6e0dc28feda82e7c39bb6b1638ad036b *R/zzz.R +fc09e35d33437d635f6361e42c35696e *README.md +0e95615d823df27d4b2a94e7be3a749a *build/vignette.rds +512913a6d2935ccb9dc5d5dc73640404 *inst/doc/SpatialKDE.R +abf6cca3ed36a52d85d0c50a4c33a8d2 *inst/doc/SpatialKDE.Rmd +f70b055123458caa7ef0afa06fc47d4b *inst/doc/SpatialKDE.html +fc333f44de04b8cb7aec0ce9c2d3ff15 *man/create_grid.Rd +89e5d758105098f7ca26068280fd22ad *man/create_raster.Rd +00277ed05f610057238f0007553ded4e *man/kde.Rd +1f7896a1b866ff9ae89ba35be7c7b6f1 *man/pipe.Rd +c8e651984bf03903e92a85a54f4e79b4 *src/RcppExports.cpp +746686d03189a416b9c16fb0536c1489 *src/kde.cpp +70b197edb62a1a7250b6a1aa009a5ecd *tests/testthat.R +78b8bef43e80f100e302e047fabaa9af *tests/testthat/test-create_grids_raster.R +9ba68f18046da3b5bd6f4474671a677b *tests/testthat/test-kde.R +abf6cca3ed36a52d85d0c50a4c33a8d2 *vignettes/SpatialKDE.Rmd diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..bb1348c --- /dev/null +++ b/NAMESPACE @@ -0,0 +1,38 @@ +# Generated by roxygen2: do not edit by hand + +export("%>%") +export(create_grid_hexagonal) +export(create_grid_rectangular) +export(create_raster) +export(kde) +importFrom(Rcpp,evalCpp) +importFrom(dplyr,filter) +importFrom(dplyr,mutate) +importFrom(dplyr,select) +importFrom(glue,glue) +importFrom(glue,glue_collapse) +importFrom(magrittr,"%>%") +importFrom(methods,setMethod) +importFrom(raster,crs) +importFrom(raster,isLonLat) +importFrom(raster,raster) +importFrom(raster,values) +importFrom(raster,xyFromCell) +importFrom(rlang,.data) +importFrom(rlang,enquo) +importFrom(rlang,quo_text) +importFrom(sf,st_as_sf) +importFrom(sf,st_bbox) +importFrom(sf,st_buffer) +importFrom(sf,st_centroid) +importFrom(sf,st_convex_hull) +importFrom(sf,st_coordinates) +importFrom(sf,st_covered_by) +importFrom(sf,st_geometry) +importFrom(sf,st_geometry_type) +importFrom(sf,st_intersects) +importFrom(sf,st_is_longlat) +importFrom(sf,st_make_grid) +importFrom(sf,st_sf) +importFrom(sf,st_union) +useDynLib(SpatialKDE) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..8118129 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,13 @@ +# SpatialKDE 0.4.0 + +* rename of `create_raster_rectangular()` to `create_grid_rectangular()` and `create_raster_hexagonal()` to `create_grid_hexagonal()` to avoid confusion about the type of outcome. + +* fix of inner workings of `create_grid_rectangular()` and `create_grid_hexagonal()` + +* fix usage of kernell in function `kde()` + +* Added a `NEWS.md` file to track changes to the package. + +# SpatialKDE 0.3.2 + +* Website via pkgdown diff --git a/R/RcppExports.R b/R/RcppExports.R new file mode 100644 index 0000000..fb55986 --- /dev/null +++ b/R/RcppExports.R @@ -0,0 +1,7 @@ +# Generated by using Rcpp::compileAttributes() -> do not edit by hand +# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +kde_estimate <- function(fishnet, points, bw, kernel, scaled = FALSE, decay = 1) { + .Call('_SpatialKDE_kde_estimate', PACKAGE = 'SpatialKDE', fishnet, points, bw, kernel, scaled, decay) +} + diff --git a/R/create_grid.R b/R/create_grid.R new file mode 100644 index 0000000..4f18239 --- /dev/null +++ b/R/create_grid.R @@ -0,0 +1,90 @@ +#' Create grid +#' +#' Create grid of equaly spaced rectangles or hexagons. The distance between centre points +#' in both x and y dimension is equal to \code{cell_size}. The function is effectively a wrapper around +#' \code{\link[sf]{st_make_grid}} with a little bit of preprocessing including generation of grid only inside +#' \code{\link[sf]{st_convex_hull}}. +#' +#' @param geometry \code{\link[sf]{sf}} \code{data.frame} containing geometry which should be cover by +#' the grid. +#' @param cell_size \code{numeric} specifing the distance for equally spaced centers of polygons +#' (rectangular or hexagonal). +#' @param side_offset \code{numeric} specifing the side offset, distance added to the convex hull +#' of input geometry to generate grid for KDE. Good estimate is usually the same value as band width of KDE. +#' @param only_inside \code{logical} specifing if the grid cells should be generated only inside of the +#' geometry. Default value is \code{FALSE}. +#' +#' @return \code{\link[sf]{sf}} \code{data.frame}. +#' @export +#' +#' @describeIn create_grid Create rectangular grid +#' +#' @examples +#' library(sf) +#' nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% st_transform(32031) +#' grid <- create_grid_hexagonal(nc, cell_size = 100000) +#' grid <- create_grid_rectangular(nc, cell_size = 100000, only_inside = TRUE) +#' +create_grid_rectangular <- function(geometry, cell_size, side_offset = 0, only_inside = FALSE){ + + .create_grid(geometry, cell_size, side_offset, only_inside, square = TRUE) +} + +#' @export +#' +#' @describeIn create_grid Create hexagonal grid +create_grid_hexagonal <- function(geometry, cell_size, side_offset = 0, only_inside = FALSE){ + + .create_grid(geometry, cell_size, side_offset, only_inside, square = FALSE) +} + +#' @importFrom sf st_geometry st_union st_convex_hull st_buffer st_make_grid +#' @importFrom sf st_sf st_intersects st_covered_by +#' @importFrom rlang .data +#' @importFrom dplyr mutate filter select +.create_grid <- function(geometry, + cell_size, + side_offset = 0, + only_inside = FALSE, + square = TRUE) { + + .validate_sf(geometry) + + .validate_sideoffset(side_offset) + + .validate_cellsize(cell_size) + + if (!(typeof(only_inside) == "logical")) { + stop(glue::glue( + "Parameter `only_inside` must be \"logical\". Currently it is of type: `{typeof(only_inside)}`." + )) + } + + buff_convex_hull <- geometry %>% + sf::st_geometry() %>% + sf::st_union() %>% + sf::st_convex_hull() %>% + sf::st_buffer(side_offset) + + grid <- buff_convex_hull %>% + sf::st_make_grid(cellsize = cell_size, + what = "polygons", + square = square) %>% + sf::st_sf() + + if (only_inside) { + + grid <- grid %>% + dplyr::mutate(covered = as.numeric(sf::st_covered_by(grid, buff_convex_hull))) %>% + dplyr::filter(!is.na(.data$covered)) + + } else { + + grid <- grid %>% + dplyr::mutate(intersect = as.numeric(sf::st_intersects(grid, buff_convex_hull))) %>% + dplyr::filter(!is.na(.data$intersect)) + } + + grid %>% + dplyr::select() +} diff --git a/R/create_raster.R b/R/create_raster.R new file mode 100644 index 0000000..4ccb130 --- /dev/null +++ b/R/create_raster.R @@ -0,0 +1,42 @@ +#' Create raster +#' +#' Create raster of equaly spaced cells. The distance between centre of cells +#' in both x and y dimension is equal to \code{cell_size}. +#' +#' @param geometry \code{\link[sf]{sf}} \code{data.frame} containing geometry which should be cover by +#' the raster. +#' @param cell_size \code{numeric} specifing the distance for equally spaced cells. +#' @param side_offset \code{numeric} specifing the side offset, distance added to the convex hull +#' of input geometry to generate raster for KDE. Good estimate is usually the same value as band width of KDE. +#' +#' @return \code{\link[raster]{Raster-class}} +#' @export +#' +#' @importFrom sf st_convex_hull st_buffer st_geometry st_union st_as_sf +#' @importFrom raster raster +#' +#' @examples +#' library(sf) +#' nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% st_transform(32031) +#' raster <- create_raster(nc, cell_size = 100000) +#' +create_raster <- function(geometry, cell_size, side_offset = 0){ + + .validate_sf(geometry) + + .validate_sideoffset(side_offset) + + .validate_cellsize(cell_size) + + buffered_geometry <- geometry %>% + sf::st_geometry() %>% + sf::st_union() %>% + sf::st_convex_hull() %>% + sf::st_buffer(side_offset) %>% + sf::st_as_sf() + + raster <- raster::raster(buffered_geometry, + resolution = c(cell_size, cell_size)) + + raster +} diff --git a/R/kde.R b/R/kde.R new file mode 100644 index 0000000..16c1df5 --- /dev/null +++ b/R/kde.R @@ -0,0 +1,182 @@ +#' Kernel Density Estimation +#' +#' KDE for spatial data. The algorithm is heavily inspired by +#' \href{https://github.com/qgis/QGIS/blob/b3d2619976a69d7fb67b884492da491dfaba287c/src/analysis/raster/qgskde.cpp}{Heatmap tool} +#' in QGIS. The help for QGIS tools is provided \href{https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/qgis/interpolation.html#heatmap-kernel-density-estimation}{at the QGIS website}. +#' The a tutorial is provided \href{https://grindgis.com/software/heat-map-using-qgis}{here}. +#' +#' @details +#' \code{grid} parameter specifies output of the function. KDE is calculated on the specified \code{grid}. +#' If grid is \code{\link[raster]{Raster-class}} then outcome is also \code{\link[raster]{Raster-class}}. +#' If grid is \code{\link[sf]{sf}} \code{data.frame} then outcome is also \code{\link[sf]{sf}} \code{data.frame}. +#' +#' +#' @param points \code{\link[sf]{sf}} \code{data.frame} containing only POINTS. +#' @param band_width \code{numeric} specifing the band width for KDE. +#' @param decay \code{numeric} specifing the decay parameter for \code{"triangular"} kernel. For +#' other kernels besides \code{"triangular"} the parameter is not used. +#' @param kernel \code{character} specifing type of kernel to use. Available implemented kernels are +#' \code{"uniform", "quartic", "triweight", "epanechnikov", "triangular"}. Default is \code{"quartic"} and if +#' uknown kernel name is used it falls back to the default value. +#' @param scaled \code{logical} specifing if the output values should be scaled. Default value is +#' \code{FALSE}. +#' @param grid either \code{\link[sf]{sf}} \code{data.frame} (outcome of function +#' \code{\link{create_grid_rectangular}} or \code{\link{create_grid_hexagonal}}) or +#' \code{\link[raster]{Raster-class}} (outcome of function \code{\link{create_raster}}). +#' Does not have to be specified if \code{cell_size} is set. +#' @param cell_size \code{numeric} specifing the distance for equal spaced points. Must be +#' higher than 0. Can be left out if \code{grid} is provided as \code{grid} is used instead. +#' The code used to generate grid is \code{\link{create_grid_rectangular}(points, cell_size, band_width)}. +#' +#' @return either \code{\link[sf]{sf}} \code{data.frame} or \code{\link[raster]{Raster-class}} +#' depending on class of \code{grid} parameter. +#' @export +#' +#' @importFrom sf st_is_longlat st_bbox st_geometry st_union st_convex_hull st_buffer st_make_grid +#' @importFrom sf st_sf st_coordinates st_geometry_type st_centroid +#' @importFrom dplyr mutate +#' @importFrom glue glue glue_collapse +#' @importFrom Rcpp evalCpp +#' @importFrom rlang .data +#' +#' @useDynLib SpatialKDE +#' +#' @examples +#' library(sf) +#' nc <- st_read(system.file("shape/nc.shp", package="sf")) %>% st_transform(32031) +#' grid <- create_grid_hexagonal(nc, cell_size = 100000) +#' points <- st_sample(nc, 500) %>% st_as_sf() +#' kde_estimate_grid <- kde(points, band_width = 150000, grid = grid) +#' raster <- create_raster(nc, cell_size = 100000) +#' kde_estimate_raster <- kde(points, band_width = 150000, grid = raster) +#' +kde <- function(points, + band_width, + decay = 1, + kernel = "quartic", + scaled = FALSE, + grid, + cell_size){ + + available_kernels = c("uniform", "quartic", "triweight", "epanechnikov", "triangular") + + .validate_sf(points) + + .validate_projected(points) + + .validate_points(points) + + if (missing(grid) & missing(cell_size)) { + stop("Both variables `grid` and `cellsize` are not specified. Don't know how to create grid for KDE estimation.") + } + + if (!missing(cell_size)) { + .validate_cellsize(cell_size) + } + + if (missing(grid)) { + + grid <- create_grid_rectangular(points, cell_size, band_width) + + } + + .validate_bandwidth(band_width) + + + if (!(kernel %in% available_kernels)) { + warning(glue::glue("Unknown `kernel` used. The implemented kernels are: ", + glue::glue_collapse(available_kernels, sep = ", "), ". ", + "Using `quartic` kernel as default.")) + kernel = available_kernels[2] + } + + kde_calculated <- .kde(grid = grid, + points = points, + band_width = band_width, + decay = decay, + kernel = kernel, + scaled = scaled) + + return(kde_calculated) +} + + +.kde <- function(grid, + points, + band_width, + decay, + kernel, + scaled) { + UseMethod(".kde") +} + +#' @importFrom sf st_centroid st_coordinates +#' @importFrom dplyr mutate +.kde.sf <- function(grid, + points, + band_width, + decay, + kernel, + scaled) { + + .validate_sf(grid) + + .validate_projected(grid) + + cells_number <- nrow(grid) + + .warn_long_calculation(cells_number) + + if (all(unique(st_geometry_type(grid)) == "POINT")) { + grid_points <- grid + } else { + message("Using centroids instead of provided `grid` geometries to calculate KDE estimates.") + suppressWarnings( + grid_points <- grid %>% + sf::st_centroid(of_largest_polygon = TRUE) + ) + } + + kde_values <- kde_estimate(sf::st_coordinates(grid_points), + sf::st_coordinates(points), + bw = band_width, + kernel = kernel, + scaled = scaled, + decay = decay) + + grid <- grid %>% + dplyr::mutate(kde_value = kde_values) + + grid +} + +#' @importFrom raster xyFromCell values +#' @importFrom sf st_coordinates +#' @importFrom methods setMethod +setMethod(".kde", + "RasterLayer", + function(grid, + points, + band_width, + decay, + kernel, + scaled) { + + .validate_raster_projected(grid) + + cells_number <- length(grid) + + .warn_long_calculation(cells_number) + + kde_values <- kde_estimate(raster::xyFromCell(grid, 1:cells_number), + sf::st_coordinates(points), + bw = band_width, + kernel = kernel, + scaled = scaled, + decay = decay) + + raster::values(grid) <- kde_values + + grid + } +) diff --git a/R/utils-pipe.R b/R/utils-pipe.R new file mode 100644 index 0000000..333ce08 --- /dev/null +++ b/R/utils-pipe.R @@ -0,0 +1,11 @@ +#' Pipe operator +#' +#' See \code{magrittr::\link[magrittr]{\%>\%}} for details. +#' +#' @name %>% +#' @rdname pipe +#' @keywords internal +#' @export +#' @importFrom magrittr %>% +#' @usage lhs \%>\% rhs +NULL diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..6300a02 --- /dev/null +++ b/R/utils.R @@ -0,0 +1,85 @@ +#' @importFrom glue glue +.warn_long_calculation <- function(cells_number){ + if (cells_number > 50000) { + message(glue::glue("The number of cells (points) in the outcomes is large (`{cells_number}`). ", + "The calculation may take a while. To speed it up you can use larger `cell_size`.")) + } +} + +#' @importFrom rlang quo_text enquo +#' @importFrom glue glue +.validate_sf <- function(x){ + var_name <- rlang::quo_text(rlang::enquo(x)) + + if (!("sf" %in% class(x))) { + stop(glue::glue("Variable `{var_name}` must be of class `sf`, currently is has classes: ", + glue::glue_collapse(class(x), sep = ", "), ".")) + } +} + +#' @importFrom sf st_geometry_type +#' @importFrom glue glue +.validate_points <- function(points){ + + if (!(all(unique(sf::st_geometry_type(points)) == "POINT"))) { + stop(glue::glue("`Points` must be only of sfc `POINT`. Other types were found: ", + glue::glue_collapse(unique(sf::st_geometry_type(points)), sep = ", "), ".")) + } +} + +#' @importFrom glue glue +#' @importFrom rlang quo_text enquo +#' @importFrom sf st_is_longlat +.validate_projected <- function(x){ + var_name <- rlang::quo_text(rlang::enquo(x)) + + if (sf::st_is_longlat(x)) { + stop(glue::glue("Variable `{var_name}` layer must be projected. Cannot calculate KDE on geographical coordinates.")) + } +} + +#' @importFrom glue glue +#' @importFrom rlang quo_text enquo +#' @importFrom raster isLonLat crs +.validate_raster_projected <- function(x){ + var_name <- rlang::quo_text(rlang::enquo(x)) + + if (raster::isLonLat(raster::crs(x))) { + stop(glue::glue("Raster layer `{var_name}` must be projected. Cannot calculate KDE on geographical coordinates.")) + } +} + +#' @importFrom glue glue +.validate_bandwidth <- function(band_width){ + + if (!is.numeric(band_width)) { + if (band_width <= 0) { + stop(glue::glue("Band_width parameter must be numerical and higher than zero. ", + "Currently it is `{class(band_width)}` with value `{band_width}`.")) + } + } +} + +#' @importFrom glue glue +.validate_cellsize <- function(cell_size){ + + if (!is.numeric(cell_size) | cell_size <= 0) { + stop(glue::glue("Cell_size parameter must be numerical and higher than zero. ", + "Currently it is `{class(cell_size)}` with value `{cell_size}`.")) + } +} + +#' @importFrom sf st_geometry_type +.is_polygon <- function(geometry){ + + all(unique(sf::st_geometry_type(geometry)) %in% c("POLYGON", "MULTIPOLYGON")) +} + +#' @importFrom glue glue +.validate_sideoffset <- function(side_offset){ + + if (!is.numeric(side_offset) | side_offset < 0) { + stop(glue::glue("Side_offset parameter must be numerical and higher than zero. ", + "Currently it is `{class(side_offset)}` with value `{side_offset}`.")) + } +} diff --git a/R/zzz.R b/R/zzz.R new file mode 100644 index 0000000..db7a839 --- /dev/null +++ b/R/zzz.R @@ -0,0 +1,3 @@ +.onUnload <- function (libpath) { + library.dynam.unload("SpatialKDE", libpath) +} diff --git a/README.md b/README.md new file mode 100644 index 0000000..b59bafc --- /dev/null +++ b/README.md @@ -0,0 +1,40 @@ +# SpatialKDE + + +Github actions: +![R-CMD-check](https://github.com/JanCaha/SpatialKDE/workflows/R-CMD-check/badge.svg) +![Windows Release](https://github.com/JanCaha/SpatialKDE/workflows/Build%20Windows%20Binaries%20and%20Create%20Release/badge.svg) +![Web deploy](https://github.com/JanCaha/SpatialKDE/workflows/Pkgdown%20-%20build%20and%20deploy%20website/badge.svg) +![Lifecycle: experimental](https://img.shields.io/badge/Lifecycle-stable-green.svg) + +Travis build - R-devel: +[![Build Status](https://travis-ci.org/JanCaha/SpatialKDE.svg?branch=master)](https://travis-ci.org/JanCaha/SpatialKDE) + +Cran Status: +![Cran: version](https://www.r-pkg.org/badges/version/SpatialKDE) +![Cran: downloads](https://cranlogs.r-pkg.org/badges/grand-total/SpatialKDE) + + +R package to calculate spatial KDE. Inspired by the tool Heatmap tool from QGIS. Help for Heatmap tool can be found [here](https://docs.qgis.org/2.18/en/docs/user_manual/plugins/plugins_heatmap.html), the help is for older version of the tool, but the window of the tool looks relatively the same. + +## Documentation + +[Available here.](https://jancaha.github.io/SpatialKDE/) + +## Instalation + +If you have RTools (Windows), or any C++ compiler (Linux, macOS) installed then you can use: + +``` r +remotes::install_github("JanCaha/SpatialKDE") +``` + +### Compiled Window binaries + +For Windows the binaries are available from [release page](https://github.com/JanCaha/SpatialKDE/releases/). + +The instalation can be done from R using command, where you just replace __x.x.x__ with the current version (e.g. __0.1.0__): + +``` r +install.packages('SpatialKDE_x.x.x.zip', repos = NULL, type = "win.binary") +``` diff --git a/build/vignette.rds b/build/vignette.rds new file mode 100644 index 0000000000000000000000000000000000000000..c592467ba7bfb8c0ad82878e195219b553c97c32 GIT binary patch literal 210 zcmV;@04@I?iwFP!0000028-ZgU|?WkU}j@vU}6R`nT3G_8xRWsF(U&D11FH?3ob}3 z$xO`gc5&4U%1uF&6-CG@6qaTtXBU?w7L}kW;zu$fqa-&6O_~SU2B3q4!Ez8aK)}L; zq?R==F*mgs%>Z_n)PmH!6p-*g?D`r1qiNyrPOU7@FM{boQNspR!{U~ilM1spB(nt0 zVGeL{gEAR|(A|JxK8pJ}i}K6$V1D2Q1%CsO_zwi=9!<_kEJpVzS4v_@qF!1NP(R4s M03@ykbHD)r06z~@LI3~& literal 0 HcmV?d00001 diff --git a/inst/doc/SpatialKDE.R b/inst/doc/SpatialKDE.R new file mode 100644 index 0000000..c51accf --- /dev/null +++ b/inst/doc/SpatialKDE.R @@ -0,0 +1,51 @@ +## ----setup, include=FALSE----------------------------------------------------- +knitr::opts_chunk$set(echo = TRUE) + +## ---- message=FALSE----------------------------------------------------------- +library(SpatialKDE) +library(sp) +library(sf) +library(dplyr) +library(tmap) + +## ----------------------------------------------------------------------------- +data(meuse) + +meuse <- meuse %>% + st_as_sf(coords = c("x", "y"), dim = "XY") %>% + st_set_crs(28992) %>% + select() + +## ----------------------------------------------------------------------------- +cell_size <- 100 +band_width <- 150 + +## ----------------------------------------------------------------------------- +grid_meuse <- meuse %>% + create_grid_rectangular(cell_size = cell_size, side_offset = band_width) + +## ----------------------------------------------------------------------------- +kde <- meuse %>% + kde(band_width = band_width, kernel = "quartic", grid = grid_meuse) + +## ----------------------------------------------------------------------------- +tm_shape(kde) + + tm_polygons(col = "kde_value", palette = "viridis", title = "KDE Estimate") + +tm_shape(meuse) + + tm_bubbles(size = 0.1, col = "red") + +## ----------------------------------------------------------------------------- +raster_meuse <- meuse %>% + create_raster(cell_size = cell_size, side_offset = band_width) + +## ----------------------------------------------------------------------------- +kde <- meuse %>% + kde(band_width = band_width, kernel = "triweight", grid = raster_meuse) + +## ----------------------------------------------------------------------------- +tm_shape(kde) + + tm_raster(palette = "viridis", title = "KDE Estimate") + +tm_shape(meuse) + + tm_bubbles(size = 0.1, col = "red") + +tm_layout(legend.outside = TRUE) + diff --git a/inst/doc/SpatialKDE.Rmd b/inst/doc/SpatialKDE.Rmd new file mode 100644 index 0000000..1a2c1c2 --- /dev/null +++ b/inst/doc/SpatialKDE.Rmd @@ -0,0 +1,100 @@ +--- +title: "SpatialKDE quickstart" +author: "Jan Caha" +date: "`r Sys.time()`" +output: + rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::rmarkdown} + %\VignetteIndexEntry{SpatialKDE quickstart} + %\usepackage[utf8]{inputenc} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +## Inspiration + +SpatialKDE implements kernel density estimation for spatial data with all the necessary settings, including selection of bandwidth, kernel type and underlying grid (cell size and shape). The algorithm is based on [Heatmap tool](https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/qgis/interpolation.html#heatmap-kernel-density-estimation) from [QGIS](https://qgis.org/en/site/). The tool is adapted directly from [source code](https://github.com/qgis/QGIS/blob/b3d2619976a69d7fb67b884492da491dfaba287c/src/analysis/raster/qgskde.cpp). Example tutorial about the QGIS tool is [available here](https://grindgis.com/software/heat-map-using-qgis). + +## Example + +First we load all necessary packages. + +```{r, message=FALSE} +library(SpatialKDE) +library(sp) +library(sf) +library(dplyr) +library(tmap) +``` + +Then we load the example dataset and prepare it into expected format of `sf` `data.frame`. + +```{r} +data(meuse) + +meuse <- meuse %>% + st_as_sf(coords = c("x", "y"), dim = "XY") %>% + st_set_crs(28992) %>% + select() +``` + +Let's define variables necessary for KDE estimation, cell size of the resulting grid and band width of points. + +```{r} +cell_size <- 100 +band_width <- 150 +``` + +### Vector grid + +Now we can prepare grid for KDE estimation. We prepare rectangular grid (hexagonal is the second option) with given cell size which is slightly bigger than convex hull of the data. + +```{r} +grid_meuse <- meuse %>% + create_grid_rectangular(cell_size = cell_size, side_offset = band_width) +``` + +At this moment it is possible to calculate KDE using `kde()` function with specified settings. + +```{r} +kde <- meuse %>% + kde(band_width = band_width, kernel = "quartic", grid = grid_meuse) +``` + +The result can be visualized using `tmap` package. + +```{r} +tm_shape(kde) + + tm_polygons(col = "kde_value", palette = "viridis", title = "KDE Estimate") + +tm_shape(meuse) + + tm_bubbles(size = 0.1, col = "red") +``` + +### Raster + +Now we can prepare raster for KDE estimation. We prepare raster with given cell size which is slightly bigger than convex hull of the data. + +```{r} +raster_meuse <- meuse %>% + create_raster(cell_size = cell_size, side_offset = band_width) +``` + +At this moment it is possible to calculate KDE using `kde()` function with specified settings. + +```{r} +kde <- meuse %>% + kde(band_width = band_width, kernel = "triweight", grid = raster_meuse) +``` + +The result can be visualized using `tmap` package. + +```{r} +tm_shape(kde) + + tm_raster(palette = "viridis", title = "KDE Estimate") + +tm_shape(meuse) + + tm_bubbles(size = 0.1, col = "red") + +tm_layout(legend.outside = TRUE) +``` diff --git a/inst/doc/SpatialKDE.html b/inst/doc/SpatialKDE.html new file mode 100644 index 0000000..f33a03e --- /dev/null +++ b/inst/doc/SpatialKDE.html @@ -0,0 +1,381 @@ + + + + + + + + + + + + + + + +SpatialKDE quickstart + + + + + + + + + + + + + + + + + + + + + +

SpatialKDE quickstart

+

Jan Caha

+

2019-12-16 09:14:47

+ + + +
+

Inspiration

+

SpatialKDE implements kernel density estimation for spatial data with all the necessary settings, including selection of bandwidth, kernel type and underlying grid (cell size and shape). The algorithm is based on Heatmap tool from QGIS. The tool is adapted directly from source code. Example tutorial about the QGIS tool is available here.

+
+
+

Example

+

First we load all necessary packages.

+ +

Then we load the example dataset and prepare it into expected format of sf data.frame.

+ +

Let’s define variables necessary for KDE estimation, cell size of the resulting grid and band width of points.

+ +
+

Vector grid

+

Now we can prepare grid for KDE estimation. We prepare rectangular grid (hexagonal is the second option) with given cell size which is slightly bigger than convex hull of the data.

+ +

At this moment it is possible to calculate KDE using kde() function with specified settings.

+ +
## Using centroids instead of provided `grid` geometries to calculate KDE estimates.
+

The result can be visualized using tmap package.

+ +

+
+
+

Raster

+

Now we can prepare raster for KDE estimation. We prepare raster with given cell size which is slightly bigger than convex hull of the data.

+ +

At this moment it is possible to calculate KDE using kde() function with specified settings.

+ +

The result can be visualized using tmap package.

+ +

+
+
+ + + + + + + + + + + diff --git a/man/create_grid.Rd b/man/create_grid.Rd new file mode 100644 index 0000000..be3903e --- /dev/null +++ b/man/create_grid.Rd @@ -0,0 +1,57 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/create_grid.R +\name{create_grid_rectangular} +\alias{create_grid_rectangular} +\alias{create_grid_hexagonal} +\title{Create grid} +\usage{ +create_grid_rectangular( + geometry, + cell_size, + side_offset = 0, + only_inside = FALSE +) + +create_grid_hexagonal( + geometry, + cell_size, + side_offset = 0, + only_inside = FALSE +) +} +\arguments{ +\item{geometry}{\code{\link[sf]{sf}} \code{data.frame} containing geometry which should be cover by +the grid.} + +\item{cell_size}{\code{numeric} specifing the distance for equally spaced centers of polygons +(rectangular or hexagonal).} + +\item{side_offset}{\code{numeric} specifing the side offset, distance added to the convex hull +of input geometry to generate grid for KDE. Good estimate is usually the same value as band width of KDE.} + +\item{only_inside}{\code{logical} specifing if the grid cells should be generated only inside of the +geometry. Default value is \code{FALSE}.} +} +\value{ +\code{\link[sf]{sf}} \code{data.frame}. +} +\description{ +Create grid of equaly spaced rectangles or hexagons. The distance between centre points +in both x and y dimension is equal to \code{cell_size}. The function is effectively a wrapper around +\code{\link[sf]{st_make_grid}} with a little bit of preprocessing including generation of grid only inside +\code{\link[sf]{st_convex_hull}}. +} +\section{Functions}{ +\itemize{ +\item \code{create_grid_rectangular}: Create rectangular grid + +\item \code{create_grid_hexagonal}: Create hexagonal grid +}} + +\examples{ +library(sf) +nc <- st_read(system.file("shape/nc.shp", package="sf")) \%>\% st_transform(32031) +grid <- create_grid_hexagonal(nc, cell_size = 100000) +grid <- create_grid_rectangular(nc, cell_size = 100000, only_inside = TRUE) + +} diff --git a/man/create_raster.Rd b/man/create_raster.Rd new file mode 100644 index 0000000..0187979 --- /dev/null +++ b/man/create_raster.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/create_raster.R +\name{create_raster} +\alias{create_raster} +\title{Create raster} +\usage{ +create_raster(geometry, cell_size, side_offset = 0) +} +\arguments{ +\item{geometry}{\code{\link[sf]{sf}} \code{data.frame} containing geometry which should be cover by +the raster.} + +\item{cell_size}{\code{numeric} specifing the distance for equally spaced cells.} + +\item{side_offset}{\code{numeric} specifing the side offset, distance added to the convex hull +of input geometry to generate raster for KDE. Good estimate is usually the same value as band width of KDE.} +} +\value{ +\code{\link[raster]{Raster-class}} +} +\description{ +Create raster of equaly spaced cells. The distance between centre of cells +in both x and y dimension is equal to \code{cell_size}. +} +\examples{ +library(sf) +nc <- st_read(system.file("shape/nc.shp", package="sf")) \%>\% st_transform(32031) +raster <- create_raster(nc, cell_size = 100000) + +} diff --git a/man/kde.Rd b/man/kde.Rd new file mode 100644 index 0000000..b99d34d --- /dev/null +++ b/man/kde.Rd @@ -0,0 +1,65 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/kde.R +\name{kde} +\alias{kde} +\title{Kernel Density Estimation} +\usage{ +kde( + points, + band_width, + decay = 1, + kernel = "quartic", + scaled = FALSE, + grid, + cell_size +) +} +\arguments{ +\item{points}{\code{\link[sf]{sf}} \code{data.frame} containing only POINTS.} + +\item{band_width}{\code{numeric} specifing the band width for KDE.} + +\item{decay}{\code{numeric} specifing the decay parameter for \code{"triangular"} kernel. For +other kernels besides \code{"triangular"} the parameter is not used.} + +\item{kernel}{\code{character} specifing type of kernel to use. Available implemented kernels are +\code{"uniform", "quartic", "triweight", "epanechnikov", "triangular"}. Default is \code{"quartic"} and if +uknown kernel name is used it falls back to the default value.} + +\item{scaled}{\code{logical} specifing if the output values should be scaled. Default value is +\code{FALSE}.} + +\item{grid}{either \code{\link[sf]{sf}} \code{data.frame} (outcome of function +\code{\link{create_grid_rectangular}} or \code{\link{create_grid_hexagonal}}) or +\code{\link[raster]{Raster-class}} (outcome of function \code{\link{create_raster}}). +Does not have to be specified if \code{cell_size} is set.} + +\item{cell_size}{\code{numeric} specifing the distance for equal spaced points. Must be +higher than 0. Can be left out if \code{grid} is provided as \code{grid} is used instead. +The code used to generate grid is \code{\link{create_grid_rectangular}(points, cell_size, band_width)}.} +} +\value{ +either \code{\link[sf]{sf}} \code{data.frame} or \code{\link[raster]{Raster-class}} +depending on class of \code{grid} parameter. +} +\description{ +KDE for spatial data. The algorithm is heavily inspired by +\href{https://github.com/qgis/QGIS/blob/b3d2619976a69d7fb67b884492da491dfaba287c/src/analysis/raster/qgskde.cpp}{Heatmap tool} +in QGIS. The help for QGIS tools is provided \href{https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/qgis/interpolation.html#heatmap-kernel-density-estimation}{at the QGIS website}. +The a tutorial is provided \href{https://grindgis.com/software/heat-map-using-qgis}{here}. +} +\details{ +\code{grid} parameter specifies output of the function. KDE is calculated on the specified \code{grid}. +If grid is \code{\link[raster]{Raster-class}} then outcome is also \code{\link[raster]{Raster-class}}. +If grid is \code{\link[sf]{sf}} \code{data.frame} then outcome is also \code{\link[sf]{sf}} \code{data.frame}. +} +\examples{ +library(sf) +nc <- st_read(system.file("shape/nc.shp", package="sf")) \%>\% st_transform(32031) +grid <- create_grid_hexagonal(nc, cell_size = 100000) +points <- st_sample(nc, 500) \%>\% st_as_sf() +kde_estimate_grid <- kde(points, band_width = 150000, grid = grid) +raster <- create_raster(nc, cell_size = 100000) +kde_estimate_raster <- kde(points, band_width = 150000, grid = raster) + +} diff --git a/man/pipe.Rd b/man/pipe.Rd new file mode 100644 index 0000000..b7daf6a --- /dev/null +++ b/man/pipe.Rd @@ -0,0 +1,12 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils-pipe.R +\name{\%>\%} +\alias{\%>\%} +\title{Pipe operator} +\usage{ +lhs \%>\% rhs +} +\description{ +See \code{magrittr::\link[magrittr]{\%>\%}} for details. +} +\keyword{internal} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp new file mode 100644 index 0000000..eb3a03a --- /dev/null +++ b/src/RcppExports.cpp @@ -0,0 +1,33 @@ +// Generated by using Rcpp::compileAttributes() -> do not edit by hand +// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 + +#include + +using namespace Rcpp; + +// kde_estimate +NumericVector kde_estimate(NumericMatrix fishnet, NumericMatrix points, double bw, String kernel, bool scaled, double decay); +RcppExport SEXP _SpatialKDE_kde_estimate(SEXP fishnetSEXP, SEXP pointsSEXP, SEXP bwSEXP, SEXP kernelSEXP, SEXP scaledSEXP, SEXP decaySEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< NumericMatrix >::type fishnet(fishnetSEXP); + Rcpp::traits::input_parameter< NumericMatrix >::type points(pointsSEXP); + Rcpp::traits::input_parameter< double >::type bw(bwSEXP); + Rcpp::traits::input_parameter< String >::type kernel(kernelSEXP); + Rcpp::traits::input_parameter< bool >::type scaled(scaledSEXP); + Rcpp::traits::input_parameter< double >::type decay(decaySEXP); + rcpp_result_gen = Rcpp::wrap(kde_estimate(fishnet, points, bw, kernel, scaled, decay)); + return rcpp_result_gen; +END_RCPP +} + +static const R_CallMethodDef CallEntries[] = { + {"_SpatialKDE_kde_estimate", (DL_FUNC) &_SpatialKDE_kde_estimate, 6}, + {NULL, NULL, 0} +}; + +RcppExport void R_init_SpatialKDE(DllInfo *dll) { + R_registerRoutines(dll, NULL, CallEntries, NULL, NULL); + R_useDynamicSymbols(dll, FALSE); +} diff --git a/src/kde.cpp b/src/kde.cpp new file mode 100644 index 0000000..aa36006 --- /dev/null +++ b/src/kde.cpp @@ -0,0 +1,136 @@ +// source of the inspiration for the function +// https://github.com/qgis/QGIS/blob/b3d2619976a69d7fb67b884492da491dfaba287c/src/analysis/raster/qgskde.cpp +#include +#include + +using namespace Rcpp; + +// uniform kernel +double uniformKernel(double d, double bw, bool scaled){ + + if (scaled) { + double k = 2./(M_PI*bw); + return k * (0.5/bw); + } + else { + return 1.; + } +} + +//quarticKernel +double quarticKernel(double d, double bw, bool scaled){ + + if (scaled) { + double k = 116./(5.*M_PI*pow(bw, 2)); + return k*(15./16.)*pow(1.-pow(d/bw, 2),2); + } + else { + return pow(1.-pow(d/bw, 2), 2); + } +} + +//triweightKernel +double triweightKernel(double d, double bw, bool scaled){ + + if (scaled) { + double k = 128./(35.*M_PI*pow(bw, 2)); + return k * (35./32.)*pow(1.-pow(d/bw, 2), 3); + } + else { + return pow(1.-pow(d/bw, 2), 3); + } +} + +//epanechnikovKernel +double epanechnikovKernel(double d, double bw, bool scaled){ + + if (scaled) { + double k = 8./(3.*M_PI*pow(bw, 2)); + return k*(3./4.)*(1.-pow(d/bw, 2)); + } + else { + return (1.-pow(d/bw, 2)); + } +} + +double triangularKernel(double d, double bw, bool scaled, double decay){ + + if (scaled) { + if (decay >= 0){ + double k = 3./((1.+2.*decay)*M_PI*pow(bw,2)); + // Derived from Wand and Jones (1995), p. 175 (with addition of decay parameter) + return k*(1.-(1.-decay)*(d/bw)); + } + else { + // Non-standard or mathematically valid negative decay ("coolmap") + return (1.-(1.-decay)*(d/bw)); + } + } + else { + return(1.-(1.-decay)*(d/bw)); + } +} + +double kde_element(double d, double bw, String kernel, bool scaled, double decay){ + if (d <= bw) { + if (kernel == "uniform") { + return uniformKernel(d, bw, scaled); + } + else if(kernel == "quartic") { + return quarticKernel(d, bw, scaled); + } + else if(kernel == "triweight") { + return triweightKernel(d, bw, scaled); + } + else if(kernel == "epanechnikov") { + return epanechnikovKernel(d, bw, scaled); + } + else if(kernel == "triangular") { + return triangularKernel(d, bw, scaled, decay); + } + // default is uniform kernel + else{ + return uniformKernel(d, bw, scaled); + } + } else { + return 0; + } +} + +// [[Rcpp::export]] +NumericVector kde_estimate(NumericMatrix fishnet, + NumericMatrix points, + double bw, + String kernel, + bool scaled = false, + double decay = 1) { + + int nrow = fishnet.nrow(); + + NumericVector out(nrow); + + for (int i = 0; i < nrow; i++) { + + double d = 0; + + for (int j = 0; j < points.nrow(); j++) { + + NumericVector v1 = fishnet.row(i); + NumericVector v2 = points.row(j); + + NumericVector v3 = v1-v2; + + d += kde_element(sqrt(sum(pow(v3, 2.0))), + bw, + kernel, + scaled, + decay); + + } + + out(i) = d; + Rcpp::checkUserInterrupt(); + } + + return out; +} diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..7d974b6 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(SpatialKDE) + +test_check("SpatialKDE") diff --git a/tests/testthat/test-create_grids_raster.R b/tests/testthat/test-create_grids_raster.R new file mode 100644 index 0000000..b842cfb --- /dev/null +++ b/tests/testthat/test-create_grids_raster.R @@ -0,0 +1,50 @@ +context("Test grid and raster creation functions") +library(SpatialKDE) +library(sp) +library(sf) +library(raster) +library(dplyr) + +data(meuse) + +coordinates(meuse) <- ~x+y +proj4string(meuse) <- CRS("+init=epsg:28992") + +test_data <- meuse %>% + st_as_sf() %>% + select() + +test_that("create grid - correct", { + expect_s3_class(create_grid_rectangular(test_data, cell_size = 100), "sf") + expect_s3_class(create_grid_hexagonal(test_data, cell_size = 100), "sf") +}) + +test_that("create grid - wrong inputs", { + expect_error(create_grid_rectangular("a"), regex = "`geometry` must be of class `sf`.") + expect_error(create_grid_rectangular(test_data, cell_size = "-1"), + regex = "Currently it is `character` with value `-1`.") + expect_error(create_grid_rectangular(test_data, cell_size = -10), + regex = "Currently it is `numeric` with value `-10`.") + expect_error(create_grid_rectangular(test_data, cell_size = 1, side_offset = "-1"), + regex = "Currently it is `character` with value `-1`.") + expect_error(create_grid_rectangular(test_data, cell_size = 1, side_offset = -1), + regex = "Currently it is `numeric` with value `-1`.") + expect_error(create_grid_rectangular(test_data, cell_size = 1, side_offset = 1, only_inside = "a"), + regex = "Currently it is of type: `character`.") +}) + +test_that("create raster - correct", { + expect_s4_class(create_raster(test_data, cell_size = 100), "RasterLayer") +}) + +test_that("create raster - wrong inputs", { + expect_error(create_raster("a"), regex = "`geometry` must be of class `sf`.") + expect_error(create_raster(test_data, cell_size = "-1"), + regex = "Currently it is `character` with value `-1`.") + expect_error(create_raster(test_data, cell_size = -10), + regex = "Currently it is `numeric` with value `-10`.") + expect_error(create_raster(test_data, cell_size = 1, side_offset = "-1"), + regex = "Currently it is `character` with value `-1`.") + expect_error(create_raster(test_data, cell_size = 1, side_offset = -1), + regex = "Currently it is `numeric` with value `-1`.") +}) diff --git a/tests/testthat/test-kde.R b/tests/testthat/test-kde.R new file mode 100644 index 0000000..3e794f6 --- /dev/null +++ b/tests/testthat/test-kde.R @@ -0,0 +1,62 @@ +context("Test KDE") +library(SpatialKDE) +library(sp) +library(sf) +library(raster) +library(dplyr) + +data(meuse) + +coordinates(meuse) <- ~x+y +proj4string(meuse) <- CRS("+init=epsg:28992") + +test_data <- meuse %>% + st_as_sf() %>% + select() + +test_data_not_projected <- test_data %>% + st_transform(4326) + +test_data_polygons <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) %>% + st_transform(3969) + +test_grid <- create_grid_rectangular(test_data, cell_size = 100) + +test_grid_not_projected <- test_grid %>% st_transform(4326) + +test_raster <- create_raster(test_data, cell_size = 100) + +suppressWarnings(test_raster_not_projected <- projectRaster(test_raster, crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) + +test_that("kde - wrong inputs", { + expect_error(kde("a"), regexp = "Variable `points` must be of class `sf`, currently is has classes: character.") + expect_error(kde(test_data_not_projected), regexp = "Variable `points` layer must be projected.") + expect_error(kde(test_data_polygons), regexp = "Other types were found: MULTIPOLYGON.") + expect_error(kde(test_data), regexp = "Both variables `grid` and `cellsize` are not specified.") + expect_error(kde(test_data, cell_size = -1), regexp = "Currently it is `numeric` with value `-1`.") + expect_error(kde(test_data, cell_size = 100, band_width = -5), regexp = "Currently it is `numeric` with value `-5`.") + expect_warning(kde(test_data, cell_size = 100, band_width = 100, kernel = "wrong"), + regexp = "Unknown `kernel` used.") +}) + +test_that("kde - wrong inputs - grid input", { + expect_error(kde(test_data, cell_size = 100, band_width = 100, kernel = "quartic", + grid = test_grid_not_projected), + regexp = "Variable `grid` layer must be projected.") +}) + +test_that("kde - wrong inputs - raster input", { + expect_error(kde(test_data, cell_size = 100, band_width = 100, kernel = "quartic", + grid = test_raster_not_projected), + regexp = "Raster layer `grid` must be projected.") +}) + +test_that("results", { + expect_s3_class(kde(test_data, cell_size = 100, band_width = 100, kernel = "quartic", + grid = test_grid), + "sf") + + expect_s4_class(kde(test_data, cell_size = 100, band_width = 100, kernel = "quartic", + grid = test_raster), + "RasterLayer") +}) diff --git a/vignettes/SpatialKDE.Rmd b/vignettes/SpatialKDE.Rmd new file mode 100644 index 0000000..1a2c1c2 --- /dev/null +++ b/vignettes/SpatialKDE.Rmd @@ -0,0 +1,100 @@ +--- +title: "SpatialKDE quickstart" +author: "Jan Caha" +date: "`r Sys.time()`" +output: + rmarkdown::html_vignette +vignette: > + %\VignetteEngine{knitr::rmarkdown} + %\VignetteIndexEntry{SpatialKDE quickstart} + %\usepackage[utf8]{inputenc} +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE) +``` + +## Inspiration + +SpatialKDE implements kernel density estimation for spatial data with all the necessary settings, including selection of bandwidth, kernel type and underlying grid (cell size and shape). The algorithm is based on [Heatmap tool](https://docs.qgis.org/testing/en/docs/user_manual/processing_algs/qgis/interpolation.html#heatmap-kernel-density-estimation) from [QGIS](https://qgis.org/en/site/). The tool is adapted directly from [source code](https://github.com/qgis/QGIS/blob/b3d2619976a69d7fb67b884492da491dfaba287c/src/analysis/raster/qgskde.cpp). Example tutorial about the QGIS tool is [available here](https://grindgis.com/software/heat-map-using-qgis). + +## Example + +First we load all necessary packages. + +```{r, message=FALSE} +library(SpatialKDE) +library(sp) +library(sf) +library(dplyr) +library(tmap) +``` + +Then we load the example dataset and prepare it into expected format of `sf` `data.frame`. + +```{r} +data(meuse) + +meuse <- meuse %>% + st_as_sf(coords = c("x", "y"), dim = "XY") %>% + st_set_crs(28992) %>% + select() +``` + +Let's define variables necessary for KDE estimation, cell size of the resulting grid and band width of points. + +```{r} +cell_size <- 100 +band_width <- 150 +``` + +### Vector grid + +Now we can prepare grid for KDE estimation. We prepare rectangular grid (hexagonal is the second option) with given cell size which is slightly bigger than convex hull of the data. + +```{r} +grid_meuse <- meuse %>% + create_grid_rectangular(cell_size = cell_size, side_offset = band_width) +``` + +At this moment it is possible to calculate KDE using `kde()` function with specified settings. + +```{r} +kde <- meuse %>% + kde(band_width = band_width, kernel = "quartic", grid = grid_meuse) +``` + +The result can be visualized using `tmap` package. + +```{r} +tm_shape(kde) + + tm_polygons(col = "kde_value", palette = "viridis", title = "KDE Estimate") + +tm_shape(meuse) + + tm_bubbles(size = 0.1, col = "red") +``` + +### Raster + +Now we can prepare raster for KDE estimation. We prepare raster with given cell size which is slightly bigger than convex hull of the data. + +```{r} +raster_meuse <- meuse %>% + create_raster(cell_size = cell_size, side_offset = band_width) +``` + +At this moment it is possible to calculate KDE using `kde()` function with specified settings. + +```{r} +kde <- meuse %>% + kde(band_width = band_width, kernel = "triweight", grid = raster_meuse) +``` + +The result can be visualized using `tmap` package. + +```{r} +tm_shape(kde) + + tm_raster(palette = "viridis", title = "KDE Estimate") + +tm_shape(meuse) + + tm_bubbles(size = 0.1, col = "red") + +tm_layout(legend.outside = TRUE) +```