SpatialKDE

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.

Example

First we load all necessary packages.

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.

data(meuse)
str(meuse)

## 'data.frame':    155 obs. of  14 variables:
##  $ x      : num  181072 181025 181165 181298 181307 ...
##  $ y      : num  333611 333558 333537 333484 333330 ...
##  $ cadmium: num  11.7 8.6 6.5 2.6 2.8 3 3.2 2.8 2.4 1.6 ...
##  $ copper : num  85 81 68 81 48 61 31 29 37 24 ...
##  $ lead   : num  299 277 199 116 117 137 132 150 133 80 ...
##  $ zinc   : num  1022 1141 640 257 269 ...
##  $ elev   : num  7.91 6.98 7.8 7.66 7.48 ...
##  $ dist   : num  0.00136 0.01222 0.10303 0.19009 0.27709 ...
##  $ om     : num  13.6 14 13 8 8.7 7.8 9.2 9.5 10.6 6.3 ...
##  $ ffreq  : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
##  $ soil   : Factor w/ 3 levels "1","2","3": 1 1 1 2 2 2 2 1 1 2 ...
##  $ lime   : Factor w/ 2 levels "0","1": 2 2 2 1 1 1 1 1 1 1 ...
##  $ landuse: Factor w/ 15 levels "Aa","Ab","Ag",..: 4 4 4 11 4 11 4 2 2 15 ...
##  $ dist.m : num  50 30 150 270 380 470 240 120 240 420 ...

meuse1 <- 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.

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.

grid_meuse <- meuse1 %>% 
  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.

kde <- meuse1 %>% 
  kde(band_width = band_width, kernel = "quartic", grid = grid_meuse)

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

The result can be visualized using tmap package.

tm_shape(kde) + 
  tm_polygons(col = "kde_value", palette = "viridis", title = "KDE Estimate") +
tm_shape(meuse1) +
  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.

raster_meuse <- meuse1 %>% 
  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.

kde <- meuse1 %>% 
  kde(band_width = band_width, kernel = "triweight", grid = raster_meuse)

The result can be visualized using tmap package.

tm_shape(kde) + 
  tm_raster(palette = "viridis", title = "KDE Estimate") +
tm_shape(meuse1) +
  tm_bubbles(size = 0.1, col = "red") + 
tm_layout(legend.outside = TRUE)