R : How to build heatmap with the leaflet package

I read a post about interactive maps with R using the `leaflet` package.

``````X=cbind(lng,lat)
kde2d <- bkde2D(X, bandwidth=c(bw.ucv(X[,1]),bw.ucv(X[,2])))

x=kde2d\$x1
y=kde2d\$x2
z=kde2d\$fhat
CL=contourLines(x , y , z)

m %>% addPolygons(CL[[5]]\$x,CL[[5]]\$y,fillColor = "red", stroke = FALSE)
``````

I am not familiar with the `bkde2D` function, so I'm wondering if this code could be generalized to any shapefiles?

What if each node has a specific weight, that we would like to represent on the heat map?

Are there other ways to create a heat map with `leaflet` map in R ?

• bke2d lets you do 2d binning (kernel density estimation) for a set of points (so lng/lat pairs work well). the ks package supports kernel smoothing for data from 1- to 6-dimensions. the akima package can do interpolation (useful when you need a regular grid). it might be worth reading up on the spatial task view for this before attempting to produce something that may not represent the data properly. – hrbrmstr Nov 3 '15 at 18:05
• ok, thanks for the link, I will definitely look this. Actually the bke2d function is not working that well with my data as it works in the example, and I can t figure why. – Felipe Nov 4 '15 at 17:02

Here's my approach for making a more generalized heat map in Leaflet using R. This approach uses `contourLines`, like the previously mentioned blog post, but I use `lapply` to iterate over all the results and convert them to general polygons. In the previous example it's up to the user to individually plot each polygon, so I would call this "more generalized" (at least this is the generalization I wanted when I read the blog post!).

``````## INITIALIZE
library("leaflet")
library("data.table")
library("sp")
library("rgdal")
# library("maptools")
library("KernSmooth")

infile <- "mvthefts.csv"

## Also, clean up variable names, and convert dates
if(!file.exists(infile)){
}
setnames(dat, tolower(colnames(dat)))
setnames(dat, gsub(" ", "_", colnames(dat)))
dat <- dat[!is.na(longitude)]
dat[ , date := as.IDate(date, "%m/%d/%Y")]

## MAKE CONTOUR LINES
## Note, bandwidth choice is based on MASS::bandwidth.nrd()
kde <- bkde2D(dat[ , list(longitude, latitude)],
bandwidth=c(.0045, .0068), gridsize = c(100,100))
CL <- contourLines(kde\$x1 , kde\$x2 , kde\$fhat)

## EXTRACT CONTOUR LINE LEVELS
LEVS <- as.factor(sapply(CL, `[[`, "level"))
NLEV <- length(levels(LEVS))

## CONVERT CONTOUR LINES TO POLYGONS
pgons <- lapply(1:length(CL), function(i)
Polygons(list(Polygon(cbind(CL[[i]]\$x, CL[[i]]\$y))), ID=i))
spgons = SpatialPolygons(pgons)

## Leaflet map with polygons
``````

Here's what you'll have at this point:

``````## Leaflet map with points and polygons
## Note, this shows some problems with the KDE, in my opinion...
## For example there seems to be a hot spot at the intersection of Mayfield and
## Fillmore, but it's not getting picked up.  Maybe a smaller bw is a good idea?

addCircles(lng = dat\$longitude, lat = dat\$latitude,
radius = .5, opacity = .2, col = "blue")
``````

And this is what the heat map with points would look like:

Here's an area that suggests to me that I need to tune some parameters or perhaps use a different kernel:

``````## Leaflet map with polygons, using Spatial Data Frame
## Initially I thought that the data frame structure was necessary
## This seems to give the same results, but maybe there are some
spgonsdf = SpatialPolygonsDataFrame(Sr = spgons,
data = data.frame(level = LEVS),
match.ID = TRUE)
• Thanks for this mini-tutorial. How did you select the `bandwidth` in `bkde2d()` ? – the_darkside Feb 21 '18 at 0:14
• @the_darkside great question. In reality I fiddle with it until I get something I like, I originally developed this map specifically to examine the bandwidth assumptions. In this case I actually used `MASS::bandwidth.nrd(dat\$latitude)` and `MASS::bandwidth.nrd(dat\$longitude)` as the starting points. See `?MASS::kde2d` documentation which links to `bandwith.nrd`. Also see `?KernSmooth::dpik` if you're interested for another approach. – geneorama Feb 21 '18 at 16:02
• if `gridsize = c(100,100)` does that mean there are a total of 10,000 cells? – the_darkside Feb 21 '18 at 17:25