[R-sig-Geo] krige{gstat} error with new data prediction
Alexander.Herr at csiro.au
Alexander.Herr at csiro.au
Mon Nov 27 03:19:07 CET 2006
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?
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$va
> 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