[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