R GIS layers not lining up (spplot/plot)

Question

I'm trying to make a choropleth in R using a polygon file, with a points layer of local care homes on top. When I plot each layer individually they're fine, however when I put them both together (using spplot for the polygons and plot for the points) the coordinates not longer line up, in that various points which should be in the area are now off the map (in the sea!). Is it possible that using an spplot/plot combo would cause the problem?

I'm confident there are no errors in the original files.

library(rgdal)
library(maps)
library(maptools)
library(RColorBrewer)

setwd(....)
shape <- readOGR(".","LSOA Boundaries")
setwd(".....")
values <- read.csv("lsoamappingvalues.csv")
plotme <- merge(shape,values, by.x="LSOA11CD",by.y="LSOA")

setwd(".....")
carehomes <- read.csv("care home points.csv")

x <- carehomes$LONG
y <- carehomes$LAT
homecords <- data.frame(x=as.numeric(x),y=as.numeric(y))
coordinates(homecords) <- c('x','y')

spplot(plotme, "POP", col.regions = colorRampPalette(brewer.pal(9, "Greens"))(100), col=NA) plot(homecords, add = TRUE, col='red', pch=19)

Answer 1

Unfortunately, you cannot use spplot with base plotting. This is because spplot is a wrapper for lattice, which is a high level plotting engine. I have always found lattice obtuse and believe that it is becoming obsolete. If you are going to invest time in learning a higher level plotting engine I would recommend ggplot2.

For you problem, there is no need to use spplot as, base low level plotting will work just fine. You can create a color vector that you can pass to the polygons or points using the "col" argument.

In this example I use an "ifelse" to create a colors vector based on ranges in a vector (copper) in the polygon data.

library(sp)
library(spatialEco)
library(RColorBrewer)

data(meuse)
coordinates(meuse) <- ~x+y
polys <- hexagons(meuse, 500) 
polys@data <- data.frame(polys@data, over(polys, meuse))
polys <- polys[!is.na(polys@data$copper),]

cols <- rev(brewer.pal(5, "Spectral"))
x <- polys@data$copper
colcode <- ifelse( x < 24, cols[1],
           ifelse( x >= 24 & x < 48, cols[2],
           ifelse( x >= 48 & x < 50, cols[3],            
           ifelse( x >= 50 & x < 69, cols[4],
           ifelse( x >= 69, cols[5], NA)))))

plot(polys, col=colcode)
  plot(meuse, pch=20, cex=0.75, add=TRUE) 
  box()
  legend("bottomright", legend=c("18-24","24-48","48-50","50-69","69-108"),
         fill=attr(colcode, "palette"), bty="n")
         title(main="Copper")

You could also use the "classIntervals" and "findColours" from the classInt package to create breaks and to create the color vector.

library(classInt)
( cuts <- classIntervals(polys@data$copper, style="fixed", 
               fixedBreaks=c(summary(polys@data$copper))) )
colcode <- findColours(cuts,  rev(brewer.pal(5, "Spectral")))   
plot(polys, col=colcode)

Answer 2

I really like using tmap for choropleth mapping in R. Here's the documentation. I think it's better for plotting maps.

Here's an example of a function I once wrote to plot a choropleth map.

function (function1_info, indicator, color_scheme, classes, classification, legend_title, map_title){

# The map plot will classify the data in quartiles, but the quartiles are based on the
# results from the other DataPrep function.

# load the R libraries needed

library(rgdal) #reads spatial data
library(tmap) # thematic mapping functions
library(RColorBrewer)


# get data from the shp_prep function into this function

data_prep <- function1_info

shapefile <- data_prep[[1]]

class_breaks <- data_prep[[2]]

Map <- tm_shape(shapefile) + 
    tm_polygons(indicator, palette = color_scheme, n = classes, style = classification, breaks = class_breaks, title = legend_title, cex = 0.9, legend.hist = TRUE, legend.hist.title = "Distribution", border.alpha = .5) + 
    tm_layout(title = map_title, title.size = .96, inner.margins = c(0, 0.1, 0, 0), outer.margins = 0.01, draw.frame = FALSE, title.position = c("right", "top")) + 
    tm_scale_bar()
Map