Skip to content
Prev 6819 / 29559 Next

need to create US map with colors by state

On Thu, 5 Nov 2009, Greg Snow wrote:

            
Right. This is a data representation question, followed by a projection 
question (and maybe the offsetting of Hawaii and Alaska). An extended 
example (could someone maybe adapt as a movie?):

library(XML)
tf <- tempfile()
# grab some data
download.file("http://www.cdc.gov/flu/weekly/flureport.xml", tf)
tr0 <- xmlTreeParse(tf)
Chld <- xmlChildren(xmlRoot(tr0))
names(Chld)
# choose the most recent week
caption <- xmlAttrs(Chld[[56]])["subtitle"]
t56 <- xmlChildren(Chld[[56]])
names(t56)
# extract the states by abbreviation and activity level
abbrev <- character(length(t56))
label <- character(length(t56))
for (i in 1:length(t56)) {
   here <- xmlChildren(t56[[i]])
   abbrev[i] <- xmlValue(here[["abbrev"]])
   label[i] <- xmlValue(here[["label"]])
}
# get boundaries from the maps package
library(maps)
library(maptools)
sts <- map("state", fill=TRUE, plot=FALSE)
# convert to spatial polygons combining polygons with the same
# first name component
IDs <- sapply(strsplit(sts$names, ":"), function(x) x[1])
sp_sts <- map2SpatialPolygons(sts, IDs=IDs,
  proj4string=CRS("+proj=longlat +datum=WGS84"))
# match names and their abbreviations
data(state)
o <- match(row.names(sp_sts), tolower(state.name))
ABBs <- state.abb[o]
# add DC as we have a polygon but no abbreviation
ABBs[is.na(ABBs)] <- "DC"
# match the CDC and maps abbreviations
oo <- match(ABBs, abbrev)
# re-order the labels
label42 <- label[oo]
# change the SpatialPolygons names to their abbreviations
sp_sts1 <- spChFIDs(sp_sts, ABBs)
# make a SpatialPolygonsDataFrame
df <- data.frame(row.names=ABBs, label=factor(label42))
spdf <- SpatialPolygonsDataFrame(sp_sts1, df)
# and plot (could use better palette)
spplot(spdf, "label", main=caption)
# zooming in on DC to check if OK.
spplot(spdf, "label", xlim=c(-85, -75), ylim=c(35, 41), main=caption)

which match the display on:

http://www.cdc.gov/h1n1flu/updates/us/

for the states and entities for which there were boundaries in the source 
used. Other boundaries may be found as shapefiles on the US Census site.

To get to the projection used on the CDC webpage, use spTransform() in 
the rgdal package.

Access to US data continues to please, and lead to envy!

Roger