[R-sig-Geo] krige{gstat} error with new data prediction
Edzer J. Pebesma
e.pebesma at geo.uu.nl
Tue Nov 28 09:37:23 CET 2006
Hi Herry, the problem is that your trend model specifies a linear trend
in x and y, the coordinate names of shp1, but that x and y cannot be
resolved for the prediction locations in a1000 because there they have
names coords.x1 and coords.x2 (probably chosen default by sp). Perhaps
you replace z~x+y with z~1 (you seem to have plenty of data), otherwise
make sure that in both data sets coordinates have the same names.
--
Edzer
Alexander.Herr at csiro.au wrote:
> Hi Edzer,
>
> the data are located at ftp://ftp.csiro.au/Herry, file
> qldproperty.Rdata.
>
> call is:
> k1<-krige(z ~ x+y, ~x+y, data=shp1, newdata=a1000, model=m2, nmax=40)
>
>
> SessionInfo:
> R version 2.4.0 (2006-10-03)
> i686-pc-linux-gnu
>
> locale:
> LC_CTYPE=en_US;LC_NUMERIC=C;LC_TIME=en_US;LC_COLLATE=en_US;LC_MONETARY=e
> n_US;LC_MESSAGES=en_US;LC_PAPER=en_US;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHON
> E=C;LC_MEASUREMENT=en_US;LC_IDENTIFICATION=C
>
> attached base packages:
> [1] "methods" "stats" "graphics" "grDevices" "utils"
> "datasets"
> [7] "base"
>
> other attached packages:
> geoR gstat MCMCpack MASS coda
> lattice
> "1.6-11" "0.9-34" "0.7-4" "7.2-29" "0.10-7"
> "0.14-13"
> mcmc RColorBrewer maps maptools sp
> foreign
> "0.5-1" "0.2-3" "2.0-32" "0.6-3" "0.9-4"
> "0.8-17"
> spatial
> "7.2-29"
>
>
> Thanx
> Herry
>
> Dr Alexander Herr - Herry
> 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: Edzer J. Pebesma [mailto:e.pebesma at geo.uu.nl]
> Sent: Monday, November 27, 2006 5:31 PM
> To: Herr, Alexander Herr - Herry (CSE, Townsville)
> Subject: Re: [R-sig-Geo] krige{gstat} error with new data prediction
>
> Alexander.Herr at csiro.au wrote:
>
>> Thanks Edzer,
>>
>> k1<-krige(z ~ x+y, shp1, a1000, model=m2, nmax=40) produces the
>> following:
>>
>> Error in eval(expr. envir. enclos): object "x" not found In addition:
>> Warning message:
>> grid has empty columns/rows in dimension 2 in: points2grid(points,
>> tolerance).
>>
>> The grid will have empty
>> I was hoping to predict values for within areas (currently assigned
>> bar some pixels stemming from spatial errors) and ignore those outside
>>
>
>
>> (ie the once that fill the grid boundary between study area and
>> bounding box). I presume the empty columns/rows refer to the NA
>> outside the study area?
>>
>> Any suggestions?
>>
>>
> Yes, provide me with data (off-list), such that I can reproduce your
> problem. Also send me the output of sessionInfo(), and make sure you use
> latest versions of everything.
> --
> Edzer
>
>> Thanks
>> 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: Edzer J. Pebesma [mailto:e.pebesma at geo.uu.nl]
>> Sent: Friday, November 24, 2006 5:36 PM
>> To: Herr, Alexander Herr - Herry (CSE, Townsville)
>> Cc: r-sig-geo at stat.math.ethz.ch
>> Subject: Re: [R-sig-Geo] 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$v
>>> a
>>> r1
>>> .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