krige{gstat} error with new data prediction
Herry, in k1<-krige(z ~ x+y, ~x+y, data=shp1, newdata=a1000, model=m2, nmax=40) shpl is alread a SpatialPointsDataFrame object, so you shouldn't specify it's coordinates, ~x+y. Rather, use k1<-krige(z ~ x+y, shp1, a1000, model=m2, nmax=40) I must admit that a more helpful error message would be useful here, as this starts to become a FAQ. -- Edzer
Alexander.Herr at csiro.au wrote:
Dear list,
I am using point data of an area as shapefile import (shp1). The
attribute data is log(AREA) as z reflecting areas size of the polygons.
Grid data (1000m pixel) as ascii grid imports (a1000) of the same
bounding box. Point and grid file stem from ArcInfo imported into R
using readShapePoints() and readAsciiGrid, this producing
SpatialPointsDataFrame and SpatialGridDataFrame.
I was hoping to predict the missing areas between the points with
kriging but to no avail. What am I doing wrong?
Code:
library(spatial)
library(maptools)
library(maps)
library(geoR)
library(gstat)
#v1<-variogram(z~x+y, loc=coordinates(shp1),data=shp1)
#XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXx
load("qldproperty.Rdata") #this has all data and v1
plot(v1,plot.numbers=T)
m1<-vgm(2.7,"Gau",200000,5.75)
plot(v1,plot.number=F, model=m1, ylim=c(4,9))
m2<-fit.variogram(v1,model=m1)
plot(v1,plot.number=F, model=m2, ylim=c(4,9))
#Xvalidation
shp1 at coords[,1]->x
shp1 at coords[,2]->y
shp1$z->z
as.data.frame(cbind(x,y,z))->dta1
kcv<-krige.cv(z~x+y, ~x+y, dta1, model=m2, nmax=40, nfold=50,
verbose=TRUE,maxdist=10000)
truehist(kcv$residual, xlim=c(-20,20), h=1, col="grey")
rownames(kcv[is.na(kcv$residual),])->exclude
kcv[!rownames(kcv) %in% exclude,]->tmp
#goodness of fit: Mean error (0), Root Mean square error (low), and mean
squared Dev.ratio. (1)
mean(tmp$residual);mean(tmp$residual**2);mean((tmp$residual)**2/tmp$var1
.var)
coordinates(tmp)<- ~x+y
bubble(tmp, "residual", main = "log(area): 50-fold validated krieging
residuals",fill=F)
k1<-krige(z ~ x+y, ~x+y, data=shp1, newdata=a1000, model=m2, nmax=40)
Error:
Error in function (classes, fdef, mtable) :
unable to find an inherited method for function "coordinates<-",
for signature "SpatialPointsDataFrame"
Running R2.4 on Mandrake linux 10.2 The data is at ftp://ftp.csiro.au/Herry/ file: qldproperty.Rdata Any help appreciated. Thanx Herry Dr Alexander Herr Spatial and statistical analyst CSIRO, Sustainable Ecosystems Davies Laboratory, University Drive, Douglas, QLD 4814 Private Mail Bag, Aitkenvale, QLD 4814 Phone/www (07) 4753 8510; 4753 8650(fax) Home: http://herry.ausbats.org.au Webadmin ABS: http://ausbats.org.au Sustainable Ecosystems: http://www.cse.csiro.au/ -------------------------------------------- -----Original Message----- From: r-sig-geo-bounces at stat.math.ethz.ch [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf Of r-sig-geo-request at stat.math.ethz.ch Sent: Tuesday, August 15, 2006 8:00 PM To: r-sig-geo at stat.math.ethz.ch Subject: R-sig-Geo Digest, Vol 36, Issue 12 Send R-sig-Geo mailing list submissions to r-sig-geo at stat.math.ethz.ch To subscribe or unsubscribe via the World Wide Web, visit
https://stat.ethz.ch/mailman/listinfo/r-sig-geo or, via email, send a message with subject or body 'help' to r-sig-geo-request at stat.math.ethz.ch You can reach the person managing the list at r-sig-geo-owner at stat.math.ethz.ch When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-Geo digest..." Today's Topics: 1. Linking lm residuls to specific polygon in shape files (David Martell) 2. Re: Linking lm residuls to specific polygon in shape files (Roger Bivand) 3. Re: Linking lm residuls to specific polygon in shape files (David Martell) ---------------------------------------------------------------------- Message: 1 Date: Mon, 14 Aug 2006 14:26:19 -0400 From: David Martell <martell at smokey.forestry.utoronto.ca> Subject: [R-sig-Geo] Linking lm residuls to specific polygon in shape files To: r-sig-geo at stat.math.ethz.ch Message-ID: <Pine.SOL.4.55.0608141416550.867 at smokey.forestry.utoronto.ca> Content-Type: TEXT/PLAIN; charset=US-ASCII I am using the lm.morantest procedure in spdep to test for spatial autocorrelation in my lm residuals. I am not sure where this procedure gets the residuals. I have prepared a GIS coverage that has a .shp file. I ran the lm model, produced the residuals, then put them them in the GIS coverage. When I run lm in my R program, it produces an lm object that also contains the residuals. Does lm.morantest get the residuals from the .shp file or the lm object. If it gets them from the latter, how does it determine which residual is assigned to which polygon in the .shp file. Thanks, Dave Martell _____________________________________________________________ David L. Martell, Professor Faculty of Forestry, University of Toronto 33 Willcocks Street, Toronto, Ontario, Canada, M5S 3B3 E-mail: martell at smokey.forestry.utoronto.ca URL: http://www.firelab.utoronto.ca Phone: (416) 978-6960 Fax: (416) 978-3834 Map: ES on campus map at http://oracle.osm.utoronto.ca/map/ ------------------------------ Message: 2 Date: Mon, 14 Aug 2006 22:12:55 +0200 (CEST) From: Roger Bivand <Roger.Bivand at nhh.no> Subject: Re: [R-sig-Geo] Linking lm residuls to specific polygon in shape files To: David Martell <martell at smokey.forestry.utoronto.ca> Cc: r-sig-geo at stat.math.ethz.ch Message-ID: <Pine.LNX.4.44.0608142156080.13097-100000 at reclus.nhh.no> Content-Type: TEXT/PLAIN; charset=US-ASCII On Mon, 14 Aug 2006, David Martell wrote: I am using the lm.morantest procedure in spdep to test for spatial autocorrelation in my lm residuals. I am not sure where this procedure gets the residuals. I have prepared a GIS coverage that has a .shp file. I ran the lm model, produced the residuals, then put them them in the GIS coverage. When I run lm in my R program, it produces an lm object that also contains the residuals. Does lm.morantest get the residuals from the .shp file or the lm object. If it gets them from the latter, how does it determine which residual is assigned to which polygon in the .shp file. lm.morantest() uses the whole "lm" object returned by lm(), so the residuals are by definition in the same order as the variables included in the model formula. If you need tight control, do something like this: xx <- readShapePoly(system.file("shapes/columbus.shp", package="maptools")[1], IDvar="NEIGNO") row.names(as(xx, "data.frame")) xx_nb <- poly2nb(as(xx, "SpatialPolygons")) attr(xx_nb, "region.id") lm_obj <- lm(CRIME ~ INC + HOVAL, data=xx) names(residuals(lm_obj)) lm.morantest(lm_obj, nb2listw(xx_nb), spChk=TRUE) where spChk=TRUE checks the equality of the names of the residuals of the "lm" object (taken from the data= object of the call to lm()) with the "region.id" attribute of the neighbour list. From the time the sapefile is read, it has nothing more to do with the SpatialPolygonsDataFrame object (here xx). Doing something like: xx$lm_obj_resids <- residuals(lm_obj) and then writePolyShape(xx, "newshape") includes them in a new shapefile. Does this clarify things? Roger Thanks, Dave Martell _____________________________________________________________ David L. Martell, Professor Faculty of Forestry, University of Toronto 33 Willcocks Street, Toronto, Ontario, Canada, M5S 3B3 E-mail: martell at smokey.forestry.utoronto.ca URL: http://www.firelab.utoronto.ca Phone: (416) 978-6960 Fax: (416) 978-3834 Map: ES on campus map at http://oracle.osm.utoronto.ca/map/ _______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no ------------------------------ Message: 3 Date: Mon, 14 Aug 2006 16:19:25 -0400 From: David Martell <martell at smokey.forestry.utoronto.ca> Subject: Re: [R-sig-Geo] Linking lm residuls to specific polygon in shape files To: Roger Bivand <Roger.Bivand at nhh.no> Cc: r-sig-geo at stat.math.ethz.ch Message-ID: <Pine.SOL.4.55.0608141617260.867 at smokey.forestry.utoronto.ca> Content-Type: TEXT/PLAIN; charset=US-ASCII Hello Roger, Yes, that's very helpful. Thanks very much, Dave Martell On Mon, 14 Aug 2006, Roger Bivand wrote: On Mon, 14 Aug 2006, David Martell wrote: I am using the lm.morantest procedure in spdep to test for spatial autocorrelation in my lm residuals. I am not sure where this procedure gets the residuals. I have prepared a GIS coverage that has a .shp file. I ran the lm model, produced the residuals, then put them them in the GIS coverage. When I run lm in my R program, it produces an lm object that also contains the residuals. Does lm.morantest get the residuals from the .shp file or the lm object. If it gets them from the latter, how does it determine which residual is assigned to which polygon in the .shp file. lm.morantest() uses the whole "lm" object returned by lm(), so the residuals are by definition in the same order as the variables included in the model formula. If you need tight control, do something like this: xx <- readShapePoly(system.file("shapes/columbus.shp", package="maptools")[1], IDvar="NEIGNO") row.names(as(xx, "data.frame")) xx_nb <- poly2nb(as(xx, "SpatialPolygons")) attr(xx_nb, "region.id") lm_obj <- lm(CRIME ~ INC + HOVAL, data=xx) names(residuals(lm_obj)) lm.morantest(lm_obj, nb2listw(xx_nb), spChk=TRUE) where spChk=TRUE checks the equality of the names of the residuals of the "lm" object (taken from the data= object of the call to lm()) with the "region.id" attribute of the neighbour list. From the time the sapefile is read, it has nothing more to do with the SpatialPolygonsDataFrame object (here xx). Doing something like: xx$lm_obj_resids <- residuals(lm_obj) and then writePolyShape(xx, "newshape") includes them in a new shapefile. Does this clarify things? Roger Thanks, Dave Martell _____________________________________________________________ David L. Martell, Professor Faculty of Forestry, University of Toronto 33 Willcocks Street, Toronto, Ontario, Canada, M5S 3B3 E-mail: martell at smokey.forestry.utoronto.ca URL: http://www.firelab.utoronto.ca Phone: (416) 978-6960 Fax: (416) 978-3834 Map: ES on campus map at http://oracle.osm.utoronto.ca/map/ _______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo -- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no _____________________________________________________________ David L. Martell, Professor Faculty of Forestry, University of Toronto 33 Willcocks Street, Toronto, Ontario, Canada, M5S 3B3 E-mail: martell at smokey.forestry.utoronto.ca URL: http://www.firelab.utoronto.ca Phone: (416) 978-6960 Fax: (416) 978-3834 Map: ES on campus map at http://oracle.osm.utoronto.ca/map/ ------------------------------ _______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo End of R-sig-Geo Digest, Vol 36, Issue 12 _______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo