[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