[R-sig-Geo] Adding Contours to Map (revisited)
Michael Fuller
mmfuller at tiem.utk.edu
Fri Jul 21 19:50:47 CEST 2006
I'm trying to add filled contour lines of count data to a projected
map but I'm having difficulty getting the contours to match up
properly with the projection. A nearly identical problem was
described in a previous post (Nov 2005) by Jeanne Thibeault and I've
tried the approach recommended by Roger Bivand, but I've been
unsuccessful so far.
The raw data are the longitude-latitude coordinates (x,y) of 307
survey locations spread across the US. Each coordinate pair is
associated with a count value (z). I've generated a grid of
interpolated data using interp() {akima}
> dat.intrp <- interp(x,y,z)
But how can I display these as filled contours on a US map? Roger
suggested the following for Jeanne's interpolated data set called
"Sum.80.83.li" (also created using interp()):
<begin quote>
With a little work, you can use contourLines():
res <- contourLines(Sum.80.83.li)
contours_x <- unlist(sapply(res, function(x) c(x$x, NA)))
contours_y <- unlist(sapply(res, function(x) c(x$y, NA)))
contours <- list(x=contours_x, y=contours_y)
gets you contours you can give as an argument to mapproject(),
but you'll
need to handle label positions and values yourself based on res -
look at
str(res[[1]]) to get a view of what the raw contours look like.
plot(mapproject(contours, ...), type="l")
does the plotting.
<end quote>
I'm not sure what this transformation yields (I'm still a relative
newbie). It seems that 'contours_x' and 'contours_y' are vectors
whose values are 'NA', generated using unlist() and using sapply() to
replace the x and y values of 'res' with 'NA'?
For example:
junk <-c(1:5)
junk2 <-unlist(sapply(junk, function(x) c(x$x,NA)))
junk2
[1] NA NA NA NA NA
Obviously I'm missing something here. How does this generate
contours? (It doesn't work for me).
Thanks for your patience.
Mike
More information about the R-sig-Geo
mailing list