# R Spatial Kernel Density Estimation

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?

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")
require(ks)
``````

# 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[]\$x,y=kernel[]\$y,z=kernel[])
plot(test\$x1,test\$y1,cex=log(test\$vol),xlim=c(-20,120),ylim=c(-20,120))
points(test\$x1,test\$y1,cex=log(test\$vol))
`````` ****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):

``````require(RSAGA)
require(raster)

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

# Working directory

# name of mask raster defining extent

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
outrast=paste(getwd(),"pop_kde.img",sep="/")

# SET SAGA ENVIRONMENT
( saga.env <- rsaga.env(path=spath, workspace=getwd()) )

# CALCULATE KDE
s=18000 # SCALE (RADIUS) OF KDE
rsaga.geoprocessor("io_gdal", 0, list(GRIDS="tmpm.sgrd", FILES=rname), env=saga.env)
rsaga.geoprocessor("grid_gridding", 6, list(POINTS=pts, POPULATION=pop,
TARGET="tmpm.sgrd",USER_GRID=paste(getwd(),"tmpkde.sgrd",sep="/"),
USER_SIZE=cs), env=saga.env)
rsaga.geoprocessor("io_gdal", 1, list(GRIDS="tmpkde.sgrd", FILE=outrast, FORMAT=6),
env=saga.env)

# 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)
plot(r)
``````
• Really appreciate. Jeffrey. – user2037892 Jun 8 '15 at 19:13
• `ks::kde` explicitly accepts weights (its `w` argument). – whuber Sep 4 '15 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. – Jeffrey Evans Sep 5 '15 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 '15 at 17:01

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.

• 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." – Jeffrey Evans Jun 8 '15 at 18:58
• I email Dr. Braddeley and he said the package did not have the function for such analysis. Thanks, SS_Rebelious. – user2037892 Jun 9 '15 at 15:03

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)
if (dim(data) != 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'")
#Grid:
g<-grid(n,lims) #Function to create grid.

#The distance matrix gets large... Can we work around it? YES WE CAN!
sets<-ceiling(dim(g)/10000)
#Allocate our output:
z<-rep(as.double(0),dim(g))

for (i in (1:sets)-1) {
g_subset=g[(i*10000+1):(min((i+1)*10000,dim(g))),]
a_matrix<-rdist.earth(g_subset,data,miles=FALSE)

z[(i*10000+1):(min((i+1)*10000,dim(g)))]<- apply( #Here is my kernel...
a_matrix,1,FUN=function(X)
{sum(exp(-X^2/(2*(h^2))))/(2*pi*nx)}
)
rm(a_matrix)
}

print(Sys.time())
#Un-transpose the final data.
z<-t(matrix(z,n,n))
dim(z)<-c(n^2,1)
z<-as.vector(z)
return(z)
}
``````

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=rep(gy,length(gx))
v2=rep(gx,length(gy))
v1<-matrix(v1, nrow=length(gy), ncol=length(gx))
v2<-t(matrix(v2, nrow=length(gx), ncol=length(gy)))
grid_out<-c(unlist(v1),unlist(v2))

grid_out<-aperm(array(grid_out,dim=c(n,n,2)),c(3,2,1) ) #reshape
grid_out<-unlist(as.list(grid_out))
dim(grid_out)<-c(2,n^2)
grid_out<-t(grid_out)
return(grid_out)
}
``````

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=rep(gy,length(gx))
v2=rep(gx,length(gy))
v1<-matrix(v1, nrow=length(gy), ncol=length(gx))
v2<-t(matrix(v2, nrow=length(gx), ncol=length(gy)))

image.plot(v1,v2,matrix(kde2d_mod_output,n,n))