[R-sig-Geo] contour lines on a plot

Virgilio Gomez Rubio Virgilio.Gomez at uclm.es
Fri Jan 23 20:42:13 CET 2009


Dear Maxime,

> Congratulations for the way you succeed in putting contourline at the 
> good place, although I'm not sure I understood the solution.
> I'm almost able to draw myplot... but not yet.
> The problem is that I can not use contourlines on a dataframe. I attache 
> my datafile to help you to see what it looks like.
> I tried to use contour instead of contourLines but it didn't worked.

The following linesw of code will get you some spatial interpolation using IDW (extension to ordinary kringing is very simple) and the contour lines you seek:

#
#
#

#Read data
library(sp)
library(gstat)
d<-read.table(file="Str65popK15coordK.txt", header=TRUE)
coordinates(d)<-~long+lat
proj4string(d)<-CRS("+proj=longlat")


#create à grid to perform kriging
x.range=c(-4.5,30)
y.range=c(40,55)
grd<-expand.grid(x=seq(from=x.range[1], to=x.range[2], by=1/6), 
   y=seq(from=y.range[1], to=y.range[2], by=1/6))
coordinates(grd)=~x+y
gridded(grd)=TRUE
proj4string(grd)<-CRS("+proj=longlat")

#IDW
g<-gstat(id="K1", formula=K1~1, data=d)
p<-predict(g,  newdata=grd)

#Get contour lines
library(maptools)
img<-as(p, "SpatialGridDataFrame")
img1<-as.image.SpatialGridDataFrame(img["K1.pred"])
CL<-contourLines(img1, levels=1:9/10)
CLSL<-ContourLines2SLDF(CL)

#Display results
pts<-list("sp.points", d,pch=20, col="black", cex=1)
lns<-list("sp.lines", CLSL )
spplot(p,zcol="K1.pred",col.regions=gray(0:100/100), cuts=40, 
   sp.layout=list(pts, lns), pretty=T)


Hope this helps,

Virgilio




More information about the R-sig-Geo mailing list