[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)))

Obviously I'm missing something here. How does this generate  
contours? (It doesn't work for me).

Thanks for your patience.

More information about the R-sig-Geo mailing list