I am new on spatial kernel density estimation with r and need some suggestions. Saying, I like to estimate the density for some event occurring at a location, for example, the probability of occurrence of a disease in each state, or the probability of soybean yield at each county. And then use the spatial kernel estimation result to calculate the probability that the occurrence will be lower than some specified level (threshold). What R package has the functions for spatial kernel density estimation?

The section 8.4.2 in the book, "Geographically Weighted Regression: The Analysis of Spatially Varying Relationships", have a topic about geographically weighted Kernels.This method seems to be what I need. Do anyone know any R package with functions for geographically weighted kernels?

4 Answers 4


****Edit 04/03/2023 - There is a spatialKDE package (on CRAN), written in C++, that replicates the type of KDE available in GIS software such as QGIS and AcrPro. The KDE function available in spatialEco is now sf.kde and implements both a traditional unweighted/weighted density as well as the density implementation available in the ks package. Please read the ks documentation and associate literature to understand exactly what you are getting with these results. The eks package provides a spatial interface to ks via the sf package.

****Edit 11/08/2018 - There is now a function sp.kde in the spatialEco package that allows for weighted or unweighted kernel density estimates.

There is no simple implementation of a Kernal Density Estimate using weights in R. Most of the advice for KDE's are limited to spatial locations only. You can write a function to project results from the ks package to a grid, but this is not entirely straight forward. My best advice is to leverage existing implementations from a GIS. The best option I have found to date is calling SAGA GIS directly from R. This leverages the speed and efficiency of raster processing in a GIS without having to leave R.

Here is an example (require SAGA GIS v 2.1.0):


spath = "C:/Program Files/SAGA-GIS" 

# Working directory

# name of mask raster defining extent                          
rname <- paste(getwd(),"mask.img",sep="/")

cs = 120 # Cell size of mask raster

# Shapefile of points, with weights attributes (pop)     
pts <- paste(getwd(),"pts.shp",sep="/")
pop <- "mean_count"

# Output raster
( saga.env <- rsaga.env(path=spath, workspace=getwd()) )            

rsaga.geoprocessor("io_gdal", 0, list(GRIDS="tmpm.sgrd", FILES=rname), env=saga.env)
rsaga.geoprocessor("grid_gridding", 6, list(POINTS=pts, POPULATION=pop,  
                   RADIUS=s, KERNEL="quartic kernel", GRID_GRID="tmpm.sgrd", 
                   USER_SIZE=cs), env=saga.env)
rsaga.geoprocessor("io_gdal", 1, list(GRIDS="tmpkde.sgrd", FILE=outrast, FORMAT=6),

# Clean up temporary SAGA files                
file.remove(paste(getwd(), paste("tmpm",  c("sgrd", "sdat", "mgrd", "prj"), 
            sep="."), sep="/")) 
file.remove(paste(getwd(), paste("tmpkde",  c("sgrd", "sdat", "mgrd", "prj"), 
            sep="."), sep="/")) 

# Look at output in R           
r <- raster(outrast)
  • Really appreciate. Jeffrey. Jun 8, 2015 at 19:13
  • ks::kde explicitly accepts weights (its w argument).
    – whuber
    Sep 4, 2015 at 14:37
  • @whuber, I should have been clearer. I do state: "You can write a function to project results from the ks package to a grid, but this is not entirely straight forward." Which was to imply that ks::kde can take weights but there is not a spatial implementation of the packages kernel functions. Sep 5, 2015 at 18:54
  • Have you looked at predict.kde? Although slow, it is entirely straightforward. But perhaps I am misunderstanding, because I don't know what you mean by a "spatial implementation" of a function.
    – whuber
    Sep 7, 2015 at 17:01

I have also been looking for a proper way to perform a weighted bivariate kernel interpolation. The code below worked for me:

# Download an example dataset - those are tree logs in a 100x100m plot. I used the volume of log, as weight.

test <- read.csv("https://dl.dropboxusercontent.com/u/39606472/R_rep/test.csv")

# Evaluate effect of tree felt out and in the plot. By effect I mean a mixture between trunk volume and distance to tree on a regular grid. xmin/xmax = plot area, eval.points= where I want the effect to be evaluated

    kernel <-kde(x=test[,c("x1","y1")], xmin=c(-20,-20),xmax=c(120,120),eval.points = expand.grid(x=seq(5,95,10), y=seq(5,95,10)),w=res$vol)
        IDW <- data.frame(x=kernel[[2]]$x,y=kernel[[2]]$y,z=kernel[[3]])

enter image description here


I built something for this - see: https://stackoverflow.com/questions/20108870/implementing-a-different-kernel-for-2d-kernel-density-estimation-in-r#

From there -

I ended up modifying the kde2d function from the MASS library. Some significant revision was needed, as is shown below. That said, the code is very flexible, allowing an arbitrary 2-d kernel to be used. (rdist.earth() was used for the great circle distance, h is the chosen bandwidth, in this case, in km, and n is the number of grid points in each direction to be used. rdist.earth requires the "fields" library)

The function could be modified to perform calculations in more than 2d, but the grid gets large very fast in higher dimensions. (Not that it's small now.)

Comments and suggestions on elegance or performance are welcome!

kde2d_mod <- function (data, h, n = 200, lims = c(range(data$lat), range(data$lon))) {
#Data is a matrix: lon,lat for each source. (lon,lat to match rdist.earth format.)
print(Sys.time()) #for timing

nx <- dim(data)[1]
if (dim(data)[2] != 2) 
stop("data vectors have only lat-long data")
if (any(!is.finite(data))) 
stop("missing or infinite values in the data are not allowed")
if (any(!is.finite(lims))) 
stop("only finite values are allowed in 'lims'")
g<-grid(n,lims) #Function to create grid.

#The distance matrix gets large... Can we work around it? YES WE CAN!
#Allocate our output:

for (i in (1:sets)-1) {

   z[(i*10000+1):(min((i+1)*10000,dim(g)[1]))]<- apply( #Here is my kernel...

#Un-transpose the final data.

The key point here is that any kernel can be used in that inner loop; the downside is that this is evaluated at grid points, so a high-res grid is needed to run this; FFT would be great, but I didn't attempt it.

Grid Function:

grid<- function(n,lims) {
num <- rep(n, length.out = 2L)
gx <- seq.int(lims[1L], lims[2L], length.out = num[1L])
gy <- seq.int(lims[3L], lims[4L], length.out = num[2L])

v1<-matrix(v1, nrow=length(gy), ncol=length(gx))
v2<-t(matrix(v2, nrow=length(gx), ncol=length(gy)))

grid_out<-aperm(array(grid_out,dim=c(n,n,2)),c(3,2,1) ) #reshape

You can plot the values using image.plot, with the v1 and v2 matrices for your x,y points:

kde2d_mod_plot<-function(kde2d_mod_output,n,lims) ){
 num <- rep(n, length.out = 2L)
 gx <- seq.int(lims[1L], lims[2L], length.out = num[1L])
 gy <- seq.int(lims[3L], lims[4L], length.out = num[2L])

 v1<-matrix(v1, nrow=length(gy), ncol=length(gx))
 v2<-t(matrix(v2, nrow=length(gx), ncol=length(gy)))

 map('world', fill = FALSE,add=TRUE)

Take a look at density() function in spatstat package. Official site has a number of manuals and articles about this package (see Documentation). I would recommend to start with Analysing Spatial Point Patterns in R.

  • 1
    I do not believe that the is exactly what the OP is after. To quote Baddeley directly: "The result of density.ppp is not a probability density. It is an estimate of the intensity function of the point process that generated the point pattern data." Jun 8, 2015 at 18:58
  • I email Dr. Braddeley and he said the package did not have the function for such analysis. Thanks, SS_Rebelious. Jun 9, 2015 at 15:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.