[R] help to make a map on R

[Brian.J.GREGOR@odot.state.or.us]

I've put together a few functions to help us make maps in our office.  The
first, import.poly, takes a ascii file export of a polygon coverage from
GeoMedia and imports it into a list that I've structured to use with the
polygon() function.  The second, plot.poly, makes maps using the polygon()
function. I also am working on functions to plot transportation networks
that I would be willing to share if anyone is interested.

import.poly <- function(x){
#  Open the text file and read into an object
   in.file <- file(x, "r")
   temp.file <- readLines(in.file)
   close(in.file)
#  Initialize a list to hold the results and vectors to hold all the
coordinates (for determining ranges)
#  poly <- lapply(1:length(temp.file), function(x) 0)
   poly <- list()
   all.x <- NULL
   all.y <- NULL
#  Each polygon is a character string in a vector.  Process each in turn.
   for(i in 1:length(temp.file)){
    temp.poly <- unlist(strsplit(temp.file[i], ",")) # create charater
vector from character string
    zone <- temp.poly[1] # extract the polygon identifier
    x <- as.numeric(temp.poly[seq(3, length(temp.poly), by=2)]) # extract a
vector of x coordinates
    y <- as.numeric(temp.poly[seq(2, length(temp.poly)-1, by=2)])  # extract
a vector of y coordinates
    poly[[i]] <- list("zone"=zone, "x"=x, "y"=y)  # make a list element of
the identifier and the xs and ys
    all.x <- c(all.x, x)  # concatenate the xs to a master list
    all.y <- c(all.y, y)  # concatenate the ys to a master list
   }
   x.range <- range(all.x)  # determine the range of all xs
   y.range <- range(all.y)  # determine the range of all ys
   result <- list("x.range"=x.range, "y.range"=y.range, "poly"=poly)  # put
the ranges and polygon info in a list
}

plot.poly <- function(geo, data, n, type="e", br=0, col="heat", leg.loc=0,
bnd=0, size=5, title=""){
#  The arguments are: geo = the geographic coverage
#  data = the data to be plotted by color.  Must be a vector with a names
attribute equal to the polygon names
#  n = the number of classes to depict
#  type = the types of class breaks: "e" = equal interval, "q" = quantile,
"c" = custom
#  br = a vector of custom breaks, only necessary if type="c"
#  col = name of an R color palette.  Supported types are "heat", "topo",
"terrain" and "rainbow"
#  leg.loc = a list (constructed using locator(1)) with the legend location.
This is optional.
#  bnd = a list (constructed using locator(2)) with the boundaries of the
area to be plotted.  This is optional.
#  size = the size of the maximum plot dimension in inches (default is 5
inches)
#  title = the title to be printed on the plot
#  Establish the plotting range if established with "bnd" argument or entire
coverage
   if(is.list(bnd)){xrange <- c(min(bnd$x), max(bnd$x)); yrange <-
c(min(bnd$y), max(bnd$y))}   
   else {xrange <- geo$x.range; yrange <- geo$y.range}
#  Scale the x and y dimenstions to be geographically correct   
   xlen <- abs(diff(xrange))*49  # 49 miles for each degree of longitude
   ylen <- abs(diff(yrange))*69  # 69 miles for each degree of latitude
   maxlen <- max(xlen, ylen)
   par(pin=size*c(xlen/maxlen, ylen/maxlen)) # set the plot size parameters
#  Plot a blank frame with title and axis labels
   plot(0, 0, xlim=xrange, ylim=yrange, type="n", xlab="longitude",
ylab="latitude", main=title)
#  Set up the color palette for plotting
   if(col=="heat") colors <- heat.colors(n)
   if(col=="rainbow") colors <- rainbow(n)
   if(col=="terrain") colors <- terrain.colors(n)
   if(col=="topo") colors <- topo.colors(n)
#  Set up the class breaks depending on the type parameter (custom, equal
interval, or quantile)
   if(type=="c") br <- br
   if(type=="e") br <- seq(min(data), max(data), length=n+1)
   if(type=="q") br <- quantile(data, probs=seq(0,1,1/n))
   cuts <- cut(data, br, labels=F, include.lowest=T)
   names(cuts) <- names(data) # name the vector of factors the same as the
data, used later in coloring
#  Go through a loop to plot each polygon
   for(i in 1:length(geo$poly)){
      color <- colors[cuts[geo$poly[[i]]$zone]] # choose the plot color
based on the class it is in
      polygon(geo$poly[[i]]$x, geo$poly[[i]]$y, col=color) # plot the
polygon with that color
   }
#  Plot the legend
   cuts <- cut(data, br, include.lowest=T) # do cuts() again to get
intervals for the cuts rather than integers 
   leg.txt <- as.character(levels(cuts)) # convert the intervals into text
to use in the legend
   if(leg.loc==0) leg.loc <- locator(1) # if the legend location has not
been specified wait for user to locate it
   legend(leg.loc$x, leg.loc$y, leg.txt, fill=colors, bg="white") # plot the
legend
}

Brian Gregor, P.E.
Transportation Planning Analysis Unit
Oregon Department of Transportation
Brian.J.GREGOR at odot.state.or.us
(503) 986-4120
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._