r - How to plot a heat map on a spatial map

Question

I am new to spatial data analysis in R and would like to do something easy, I am still having difficulties... I have a big table with latitudes and longitudes

sample = structure(list(Longitude = c(-0.19117, -0.211708, -0.206458, 
-0.173862, -0.156618), Latitude = c(51.489096, 51.520075, 51.525301, 
51.482442, 51.495752), Location_Easting_OSGR = c(525680L, 524170L, 
524520L, 526900L, 528060L), Location_Northing_OSGR = c(178240L, 
181650L, 182240L, 177530L, 179040L)), .Names = c("Longitude", 
"Latitude", "Location_Easting_OSGR", "Location_Northing_OSGR"
), row.names = c(NA, -5L), class = c("data.table", "data.frame"
))

I got a map of the UK from GADM (level 2 of UK map).

I would like to be able to

1. plot points defined by longitude/latitude on the map
2. build a heat map that shows where the points are more concentrated...

Is it easy ? If not do you have some pointers (only UK please) Cheers

Answer 1

Is this what you had in mind?

Your sample was too small to demonstrate a heat map, so I created a bigger sample with artificial clusters at (long,lat) = (-1,52), (-2,54) and (-4.5,56). IMO the map would be more informative without the points.

Also, I downloaded the shapefile, not the .Rdata, and imported that. The reason is that you are much more likely to find shapefiles in other projects, and it is easy to import them into R.

setwd("< directory with all your files>")
library(rgdal)         # for readOGR(...)
library(ggplot2)
library(RColorBrewer)  # for brewer.pal(...)

sample <- data.frame(Longitude=c(-1+rnorm(50,0,.5),-2+rnorm(50,0,0.5),-4.5+rnorm(50,0,.5)),
                     Latitude =c(52+rnorm(50,0,.5),54+rnorm(50,0,0.5),56+rnorm(50,0,.5)))
UKmap  <- readOGR(dsn=".",layer="GBR_adm2")
map.df <- fortify(UKmap)

ggplot(sample, aes(x=Longitude, y=Latitude)) + 
  stat_density2d(aes(fill = ..level..), alpha=0.5, geom="polygon")+
  geom_point(colour="red")+
  geom_path(data=map.df,aes(x=long, y=lat,group=group), colour="grey50")+
  scale_fill_gradientn(colours=rev(brewer.pal(7,"Spectral")))+
  xlim(-10,+2.5) +
  coord_fixed()

Explanation:

This approach uses the ggplot package, which allows you to create layers and then render the map. The calls do the following:

ggplot -         establish `sample` as the default dataset and define (Longitude,Latitude) as (x,y)
stat_density2d - heat map layer; polygons with fill color based on relative frequency of points
geom_point -     the points
geom_path -      the map (boundaries of the admin regions)
scale_fill_gradientn - defines which colors to use for the fill
xlim -           x-axis limits
coord_fixed -    force aspect ratio = 1, so map is not distorted