[R] Maps in R

[Ott Toomet]

 | From: Nirmala Ravishankar <ravishan at fas.harvard.edu>
 | Date: Tue, 7 Jan 2003 23:35:19 -0500 (EST)
 | 
 | Is there a way to generate maps in R.  Specifically, I have calculated
 | estimates of intra-regional inequality for US states, and would like to
 | project that information onto a map.
 | 
 | Thanks,
 | Nirmala

There is a package called maps at

ftp://ftp.mcs.vuw.ac.nz/pub/statistics/map/

(maps_0.1-4.tar.gz worked well for me at linux, I am not quite sure
about windows).

There are not many different countries present but US is, both on
state and county level.  However, I have noticed problems with US maps
when experimenting with color codes.  

There seem not to a function to plot each state with a given color,
the main function map() draws different polygons with different color
but a state may contain many polygons.  I have written my own, I
include it below, it is made for Danish counties but should be easy to
change for e.g. US states.  

BTW, I have made a low-resolution map for Danish counties, if anybody
has interest.  

Best,

Ott

--------------------------------------------
This function plots Danish map with color codes as given by vector x.
x is a vector with names equal to county names.  The function finds
itself maximum and minimum value of x, sets the colors accordingly,
and draws a legend.

mapPlot.default <- function(x, main=NULL, ...) {
### plot a vector on the map using names of the vector
  mkPlot <- function(name, col) {
    map("danmark", name, col=col, fill=T, add=T)
   }
  mkCol <- function(val) {
    i <- as.integer(99*(val - x0)/(x1 - x0))
    i <- ifelse(i < 1, 1, i)
    i <- ifelse(i > 100, 100, i)
    col[i]
  }
  if(!any(search() == "package:maps")) {
    cat("Loading maps library...")
    library(maps, lib.loc="/home/otoomet/proge/R")
    cat("done\n")
  }
  x <- unlist(x)
  opar <- par(xpd=TRUE,
              las=0)
  on.exit(par(opar), add=TRUE)
  regionNames <- map("danmark", plot=FALSE, namesonly=TRUE)
  r <- map("danmark", ...)$range
  mtext(main, line=1, cex=1.5)
  col <- topo.colors(100)
  iRegions <- sapply(names(x), function(y) length(grep(y, regionNames)) > 0)
                                        # which regions exist on the map
  x <- x[iRegions]
  x0 <- range(x)[1]
  x1 <- range(x)[2]
  sapply(seq(along=x),
                   FUN=function(i) mkPlot(names(x[i]), mkCol(x[i])))
  lev <- zapsmall(pretty(x, n=10))
  legx <- r[2]
  legy <- r[4]
  legend(legx, legy, legend=as.character(lev), fill=mkCol(lev),
         bty="n")
}