[R-sig-Geo] need to create US map with colors by state

[Roger Bivand]

On Thu, 5 Nov 2009, Greg Snow wrote:

> The help page for state.vbm in the TeachingDemos package shows an 
> example of doing this.  The procedure is similar for other maps that can 
> be imported and plotted using the sp and maptools packages.  The tricky 
> part is to make sure that your color specifications line up properly 
> with the data and the polygons, this is most important when a single 
> state is made up of multiple polygons (one advantage of state.vbm is 
> that it is 1 polygon per state).

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

>
> Hope this helps,
>
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no