# 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)
color = heat.colors(NLEV, NULL)[spgonsdf@data\$level])
• Scoured the interwebs trying to figure this out and this was by far the best example I found. Plugged it in to my code and it "just worked." Awesome. Thank you! – Jeff Allen Jan 31 '18 at 21:46
• Thanks! I've actually created a repo with several other map examples that might be useful for others github.com/geneorama/wnv_map_demo – geneorama Feb 9 '18 at 17:54
• 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

Building off of genorama's answer above, you can also convert the output of bkde2D into a raster rather than contour lines, using the fhat values as the raster cell values

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

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")]

## Create kernel density output
kde <- bkde2D(dat[ , list(longitude, latitude)],
bandwidth=c(.0045, .0068), gridsize = c(100,100))
# Create Raster from Kernel Density output
KernelDensityRaster <- raster(list(x=kde\$x1 ,y=kde\$x2 ,z = kde\$fhat))

#create pal function for coloring the raster
palRaster <- colorNumeric("Spectral", domain = KernelDensityRaster@data@values)

## Leaflet map with raster
colors = palRaster,
opacity = .8) %>%
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")

This is your output. Note that the low density values still show up as colored in the raster.

We can remove these low density cells with the following :

#set low density cells as NA so we can make them transparent with the colorNumeric function
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA

#create pal function for coloring the raster
palRaster <- colorNumeric("Spectral", domain = KernelDensityRaster@data@values, na.color = "transparent")

## Redraw the map
colors = palRaster,
opacity = .8) %>%
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")

Now any raster cell with a value of less than 1 is transparent.

If you want a binned raster, use the colorBin function rather than the colorNumeric function:

palRaster <- colorBin("Spectral", bins = 7, domain = KernelDensityRaster@data@values, na.color = "transparent")

## Leaflet map with raster
colors = palRaster,
opacity = .8) %>%
values = KernelDensityRaster@data@values,
title = "Kernel Density of Points")

To make it smoother, simply increase the gridsize in the bkde2D function. This increases the resolution of the generated raster. (I changed it to

gridsize = c(1000,1000)

Output:

• How can you convert the legend description “Kernel Density of Points” to something more intuitive, like “Thefts per square km”? I guess there is an equation linking the bandwidth, gridsize and projection, or perhaps even kdf\$fhat that describes the units. – fifthace Feb 14 '20 at 14:31
• Thanks! I really like this answer – geneorama Dec 14 '20 at 19:56
• Thank you very much for your answer. I got an error when running the code. The message is: "Error in get: object '.xts_chob' not found". Do you know why ? – Xiaoshi Jun 16 at 12:46
• See link for solution – Jacob Sanua Jun 17 at 19:11

An easy way of creating Leaflet heat maps in R is using the Leaflet.heat plugin. An excellent guide on how to use it can be found here. Hope you find it useful.