[R-sig-Geo] krige{gstat} error with new data prediction

Edzer J. Pebesma e.pebesma at geo.uu.nl
Fri Nov 24 08:35:36 CET 2006


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
>




More information about the R-sig-Geo mailing list