[GIS] How to do a Kernel Density with “population field” in r

arcgis-desktopkernel densitypopulationr

In ArcGIS, I can generate a kernel density map and weigh it by population field. For example, I can use city's population as weighting factor, so that more populated city has higher influence to the surrounding area.

I am trying to do the same in r. Right now I am using the package Splancs. However, I can only generate a kernel density map by treating every point the same. In other words, I cannot give more weights to certain points (e.g., a larger city).

Please let me know how to do it in r, and which package should I use.

Best Answer

I recent explored exactly this. Depending on what you want you could use the density funciton in the spatstat package. However, ppp.density returns an isotropic density, which is technically the intensity process, and there are no options for kernel type.

Here is a toy example using a raster to define extent and cell size. You can define the weights from the original sp points object. Please note: you really need to put thought into the sigma parameter (window size) and understand how it is specified. The default is "diggle" which, as seen in the example, converges on first order process. If you want 2nd order (local) spatial process you need to specify sigma as such. There are a number of bandwidth models available in spatstat to define sigma empirically. Please, use them.

I coerce back to a raster class object because honestly, it is just easier to deal with as far as talking to other sp class objects and writing to a raster format that a GIS can read. 

require(sp)
require(raster)
require(spatstat)

# create some data and plot
r <- raster(system.file("external/test.grd", package="raster"))
  pts <- sampleRandom(r, 1000, sp=TRUE) 

plot(r)
  plot(pts, pch=20, cex=0.50, add=TRUE) 

# function for coercing raster class object to spatstat im
raster.as.im <- function(im) { 
  as.im(as.matrix(im)[nrow(im):1,], xrange=bbox(im)[1,], 
        yrange=bbox(im)[2,]) 
    }   

# coerce raster to im
win <- raster.as.im(r)
  win <- as.mask(win, eps=res(r)[1])

# convert sp points to spatstat ppp object with raster defining window
x.ppp <- as.ppp(coordinates(pts), win) 

# Calculate density and coerce to raster object
den <- density(x.ppp, weights=pts@data[,1], diggle=TRUE)
  den <- raster(den)

plot(den)  
  plot(pts, pch=20, cex=0.50, add=TRUE)

This is just one approach and may not be exactly what you are after. If you would like a more generalized approach, akin to ArcGIS, let me know and I will expand my answer to include code for calling the density function in SAGA GIS (using RSAGA package), which allows for "Guassian" and "Quartic" kernels.

Related Question