Spatial Kernel Density Estimation with sf objects

[ratnanil]

To calculate kernel density of sf point features, I either need to convert sf -> sp -> ppp and use spatstat::density or use the function MASS::kde2d. I've created a wrapper around the latter to provide myself an easy access to the kernel density estimation with sf objects. Any chance something like this might be implemented in sf at one point?

st_kde <- function(points,cellsize, bandwith, extent = NULL){
  require(MASS)
  require(raster)
  require(sf)
  if(is.null(extent)){
    extent_vec <- st_bbox(points)[c(1,3,2,4)]
  } else{
    extent_vec <- st_bbox(extent)[c(1,3,2,4)]
  }
  
  n_y <- ceiling((extent_vec[4]-extent_vec[3])/cellsize)
  n_x <- ceiling((extent_vec[2]-extent_vec[1])/cellsize)
  
  extent_vec[2] <- extent_vec[1]+(n_x*cellsize)-cellsize
  extent_vec[4] <- extent_vec[3]+(n_y*cellsize)-cellsize
  
  coords <- st_coordinates(points)
  matrix <- kde2d(coords[,1],coords[,2],h = bandwith,n = c(n_x,n_y),lims = extent_vec)
  raster(matrix)
}


library(sf)
#> Linking to GEOS 3.6.2, GDAL 2.2.3, PROJ 4.9.3


points_sf <- st_as_sf(data.frame(x = rnorm(100), y = rnorm(100)), coords = c("x","y"))


points_kde <- st_kde(points_sf,0.05,1)
#> Loading required package: MASS
#> Loading required package: raster
#> Loading required package: sp
#> 
#> Attaching package: 'raster'
#> The following objects are masked from 'package:MASS':
#> 
#>     area, select

plot(points_kde)
plot(points_sf, add = TRUE)

^{Created on 2019-11-21 by the reprex package (v0.3.0)}

[edzer]

I'm not sure this kind of thing should go into sf; I am more in favor of using this kind of functionality from spatstat and write converters from sf to the classes used by spatstat. See also #1161.