8

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

In this article, the author createD a heat map like this:

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 = leaflet() %>% addTiles() 
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
8
+50

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

inurl <- "https://data.cityofchicago.org/api/views/22s8-eq8h/rows.csv?accessType=DOWNLOAD"
infile <- "mvthefts.csv"

## LOAD DATA
## Also, clean up variable names, and convert dates
if(!file.exists(infile)){
    download.file(url = inurl, destfile = infile)
}
dat <- data.table::fread(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
leaflet(spgons) %>% addTiles() %>% 
    addPolygons(color = heat.colors(NLEV, NULL)[LEVS])

Here's what you'll have at this point: enter image description here

## 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?

leaflet(spgons) %>% addTiles() %>%
    addPolygons(color = heat.colors(NLEV, NULL)[LEVS]) %>%
    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:

enter image description here

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

enter image description here

## 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 
## advantages to using the data.frame, e.g. for adding more columns
spgonsdf = SpatialPolygonsDataFrame(Sr = spgons,
                                    data = data.frame(level = LEVS),
                                    match.ID = TRUE)
leaflet() %>% addTiles() %>%
    addPolygons(data = spgonsdf,
                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
  • 2
    @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
3

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.

2

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

inurl <- "https://data.cityofchicago.org/api/views/22s8-eq8h/rows.csv?accessType=DOWNLOAD"
infile <- "mvthefts.csv"

## LOAD DATA
## Also, clean up variable names, and convert dates
if(!file.exists(infile)){
  download.file(url = inurl, destfile = infile)
}
dat <- data.table::fread(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
leaflet() %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            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.

Raster Output

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
leaflet() %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            values = KernelDensityRaster@data@values, 
            title = "Kernel Density of Points")

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

Final Map

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
leaflet() %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            values = KernelDensityRaster@data@values, 
            title = "Kernel Density of Points")

Binned Raster Kernel Density

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:

Smoothed Raster

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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