Skip to content

spplot: overlay colored points on grid

2 messages · Timothy W. Hilton, Edzer Pebesma

Original
Timothy W. Hilton
# 
R-sig-geo folks,

I'm trying to plot a SpatialPointsDataFrame over a SpatialGridDataFrame,
with both the points and the grid colored by the same color scale.  I
can get either one separately, but I haven't been able to combine them.
Code to reproduce my problem below.  The real data are results from
krige(), and I'd like to plot observations on top of the interpolated
field.

Ideally, I'd like to give the points a black outline, with a filled
background according to the value, like
plot(rnorm(10), pch=21, col="black", bg=cols)  ## cols from code below

Many thanks!
Tim

## set up a color scheme
cols <- c("#40004B", "#762A83", "#9970AB", "#C2A5CF", "#E7D4E8", "#D9F0D3", 
          "#A6DBA0", "#5AAE61", "#1B7837", "#00441B")
colbreaks <- seq(0, 1.0, by=0.1)

## make a 10x10 SpatialGridDataFrame w/ random data
grd <- SpatialGrid(grid = GridTopology(c(0,0), c(1,1), c(10,10)))
df <- SpatialGridDataFrame(grid=grd, data=data.frame(z=rnorm(nrow(coordinates(grd)))))

## make a SpatialPointsDataFrame with 15 points with random data
pts <- data.frame(x=round(runif(15) * 10),
                  y=round(runif(15) * 10),
                  z=runif(15))
coordinates(pts) <- ~x+y

## this plots the grid the way I want it plotted
spplot(df, col.regions=cols, at=colbreaks)

## this plots the points the way I'd like to
dev.new()
spplot(pts, zcol="z", col.regions=cols, cuts=colbreaks, key.space="right", pch=16, fill=TRUE)

## I can't figure out how to get the points colored by their value, like
## in the previous plot
dev.new()
spplot(df, col.regions=cols, at=colbreaks,
       sp.layout=list("sp.points", pts, pch=21, cuts=colbreaks),
       main="not quite right")
Edzer Pebesma
# 
Tim, maybe this works:

pts$col = findInterval(pts$z, colbreaks)
spplot(df, col.regions=cols, at=colbreaks,
       sp.layout=list(list("sp.points", pts, pch=21, col=1,cex=1.8),
       list("sp.points", pts, col = cols[pts$col], pch=16,cex=1.6)),
       main="maybe right")

plotting pts twice is a bit of a hack; you need to get the cex arguments
right by tweaking them somewhat.

Hth,
On 05/30/2011 11:00 PM, Timothy W. Hilton wrote: