8

I would like to create a map of the USA (perhaps a heatmap) to show frequency of a certain characteristic among states. I am not sure of what package to use or if my data is in the proper form. My data is in the table tf

tf
 AB  AK  AL  AN  AR  AZ  CA  CO  CT  DC  DE  EN  FL  GA  HI  IA  ID  IL  IN  KS 
  1  21  31   1  12  56 316  53  31  16   7   1 335  63  11  42  29  73  40  2

For the most part, my abbreviations are US (aside from a few Canadian instances). What is the best suggested approach for graphically displaying this on a map?

Now how do I get granularity of less than 50 per color?

enter image description here

Brian Tompsett - 汤莱恩
  • 5,753
  • 72
  • 57
  • 129
megv
  • 1,421
  • 5
  • 24
  • 36
  • 1
    [Maptools](http://cran.r-project.org/web/packages/maptools/index.html) (`?maptools`) and [maps](http://cran.r-project.org/web/packages/maps/index.html) (`?maps`) will be a nice starting point. You may also find the [spatial task view](http://cran.r-project.org/web/views/Spatial.html) helpful. – Seb Dec 16 '11 at 17:37
  • 2
    possible duplicate of [State level unemployment in R](http://stackoverflow.com/questions/5385713/state-level-unemployment-in-r) – IRTFM Dec 16 '11 at 20:22

2 Answers2

8

two packages: maps, ggplot2. There is an excellent example at: ?map_data()

just to start with:

tf= structure(list(state = structure(1:14, .Label = c("AK", "AL", 
"AR", "AZ", "CA", "CO", "CT", "DE", "FL", "GA", "IA", "IL", "IN", 
"KS"), class = "factor"), num = c(21L, 31L, 12L, 56L, 316L, 53L, 
31L, 7L, 335L, 63L, 42L, 73L, 40L, 2L), region = structure(c(2L, 
1L, 4L, 3L, 5L, 6L, 7L, 8L, 9L, 10L, 13L, 11L, 12L, 14L), .Label = c("alabama", 
"alaska", "arizona", "arkansas", "california", "colorado", "connecticut", 
"delaware", "florida", "georgia", "illinois", "indiana", "iowa", 
"kansas"), class = "factor")), .Names = c("state", "num", "region"
), class = "data.frame", row.names = c(NA, -14L))

require(maps);require(ggplot2)

states <- map_data("state")
tfmerged <- merge(states, tf, sort = FALSE, by = "region")
tfmerged <- tfmerged[order(tfmerged$order), ]
qplot(long, lat, data = tfmerged, group = group, fill = num,
geom="polygon")

Then fill the rests of the states info.

mdsumner
  • 29,099
  • 6
  • 83
  • 91
aatrujillob
  • 4,738
  • 3
  • 19
  • 32
  • 1
    Would I just want to append an L to my numeric data then? – megv Dec 16 '11 at 18:24
  • @megv, no need to append an L, i just added that to show the data I used. What you might need to add is the name of the states so it matches with the states$region so the merge can work. – aatrujillob Dec 16 '11 at 18:38
  • Since I already have a table with state abbreviation and counts can can i just insert that directly into the structure. Sorry..I am a bit new to this. – megv Dec 16 '11 at 18:54
  • What you need is to make sure there is at least one key (a variable that matches in both datasets) when you are doing the merge. So if you want to keep the state abreviation I guess you would need to rename region or create a variable in states with that info. – aatrujillob Dec 16 '11 at 19:38
  • Also..what is the region pertain to in the structure? – megv Dec 16 '11 at 22:34
  • How would I change my color scale to include more bins..right now it every change of 50 – megv Dec 18 '11 at 19:44
  • Can you add an image of the result? – Roman Luštrik Dec 19 '11 at 13:42
  • I added the image..now how do I get it broken up to finer granularity than 50? – megv Dec 19 '11 at 21:05
  • This is complicated. You are actually drawing each and every state. when map package already has one build in. The simplest solution is to plot using maptools and just color the states – MySchizoBuddy Feb 09 '13 at 21:36
8

Another approach with spplot:

library(maps)
library(maptools)
library(sp)

First read the data and add a column with the names of the states:

txt <- "AB  AK  AL  AN  AR  AZ  CA  CO  CT  DC  DE  EN  FL  GA  HI  IA  ID  IL  IN  KS
    1  21  31   1  12  56 316  53  31  16   7   1 335  63  11  42  29  73  40  2"

dat <- stack(read.table(text = txt,  header = TRUE))
names(dat)[2] <-'state.abb'
dat$states <- tolower(state.name[match(dat$state.abb,  state.abb)])

Then you get the map and convert it to a SpatialPolygons:

mapUSA <- map('state',  fill = TRUE,  plot = FALSE)
nms <- sapply(strsplit(mapUSA$names,  ':'),  function(x)x[1])
USApolygons <- map2SpatialPolygons(mapUSA,  IDs = nms,  CRS('+proj=longlat'))

And now you add the information from your data:

idx <- match(unique(nms),  dat$states)
dat2 <- data.frame(value = dat$value[idx], state = unique(nms))
row.names(dat2) <- unique(nms)

USAsp <- SpatialPolygonsDataFrame(USApolygons,  data = dat2)

Finally you plot it:

spplot(USAsp['value'])

Added imageenter image description here

megv
  • 1,421
  • 5
  • 24
  • 36
Oscar Perpiñán
  • 4,491
  • 17
  • 28