0

This is my first question here, so please be patient! Thanks in advance for taking your time.

I have the following problem, which I couldn't find an answer to: I am plotting income data for England in a map by travel-to-work area (TTWA). I'm using the classInt package to assign values into five classes by quantile, and then map these data using colors from the RColorBrewer package. However, when comparing the map with the raw data, the shading doesn't correspond to the right classes. For example, London has a value of 27168.83, but is shaded in the color corresponding to 18629-20166.

The images behind the links illustrate the problem. A snippet of the data frame, showing the earnings value for London. The map resultant from my code, with London highlighted

There must be a bug in the code that I just can't identify. Please see the code below. What am I doing wrong? 

# Clear workspace
rm(list = ls())
setwd("C:/Users/Carolin/OneDrive/PhD/Internal and international migration - paper 1/Data analysis/")

# Load packages
  library(sp)
  library(maptools)
  library(rgdal)
  library("classInt")
  library("RColorBrewer")
  library(rgeos)
  library(plyr)

# Read in shapefiles 
  TTWA <- readShapeSpatial("01_raw/Travel_to_Work_Areas_December_2011_Full_Clipped_Boundaries_in_United_Kingdom/Travel_to_Work_Areas_December_2011_Full_Clipped_Boundaries_in_United_Kingdom")
  TTWA <- TTWA[substr(TTWA$ttwa11cd, 1, 1) == "E", ]

# Load data from CSV
  data <- read.table(paste("02_intermediate/ttwa_workpl_tot_earnings_2016.csv", sep = ""), sep="\t", header = TRUE)

# Join attribute data to TTWA
  mapdata <- merge(data, TTWA@data, by.x="ttwa", by.y="ttwa11cd", all.y = TRUE)

# Set colours and breaks
  breaks <- classIntervals(mapdata$workpl_tot_earnings, n = 5, style = "quantile", dataPrecision = NULL)
  my_colours <- brewer.pal(5, "YlOrRd")

# Plot map
  plot(TTWA, col = my_colours[findInterval(mapdata$workpl_tot_earnings, breaks$brks, all.inside = FALSE)], axes = FALSE, border = FALSE)
  legend("topleft", legend = leglabs(round(breaks$brks, digits = 0)), fill = my_colours, bty = "o")

The shapefiles can be found here: http://geoportal.statistics.gov.uk/datasets/travel-to-work-areas-december-2011-full-clipped-boundaries-in-united-kingdom This is the structure of the ttwa_workpl_tot_earnings_2016.csv dataset: 

ttwa    workpl_tot_earnings  
E30000004   19066.27129483733  
E30000018   21016.34693877551  
E30000029   22202  
E30000039   18691  
E30000046   17585.72016081708  
E30000051   17873  
E30000054   18657.91722661888  
E30000061   18926.50491083942  
E30000064   17539.94306924291  
E30000070   20103  
E30000076   21062  
E30000093   20194.32707774799

Any help with this would be really appreciated.

Carolin Thol
  • 1
  • 1
  • can you include a sample data set that reproduces the error? See some suggestions [here](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) – De Novo Mar 08 '18 at 13:34
  • Thanks @DanHall, I've add the link to the shapefiles and an extract from the dataset. Does this help? – Carolin Thol Mar 08 '18 at 14:41
  • It's better if you follow the advice in [the link I shared in my first comment](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) and, for example, dput a minimal structure that will reproduce the error, instead of asking people to download and process a raw file. – De Novo Mar 08 '18 at 14:44
  • Sorry Dan, but I'm not getting an error: the code runs, but the results are obviously wrong. That's why I posted screenshots of my data frame and the map I'm getting, which do not correspond. I did read the page you suggested, but I don't think there's a way to represent shapefiles in this way. – Carolin Thol Mar 08 '18 at 14:56
  • by error i mean "thing that you didn't want to happen", e.g. values are wrong. – De Novo Mar 08 '18 at 15:02
  • If you could explain how I can post the shapefiles in another way, that would be very helpful. Otherwise I'm not sure how to reproduce the results. I've posted an extract of the data that I'm using. I'm using R x64 3.4.2 in RStudio on Windows 10. – Carolin Thol Mar 08 '18 at 15:07

0 Answers0