[R-sig-Geo] krige{gstat} error with new data prediction
Alexander.Herr at csiro.au
Alexander.Herr at csiro.au
Fri Nov 24 08:24:27 CET 2006
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
More information about the R-sig-Geo
mailing list