[R-sig-Geo] kriging inside borders

[Roger Bivand]

On Fri, 20 May 2005, Roger Bivand wrote:

> On Fri, 20 May 2005, Chloe ARCHINARD wrote:
> 
> > My exact commands are :
> > library(maptools)
> > shapefile=read.shape('my.shp')
> > x=Map2poly(shapefile)
> > 
> > and the length is :
> > length(mappolys)
> > [1] 615
> 
> OK. So what I think you have is a grid of cells, possibly with cells 
> missing around the edges and maybe in the middle too. The splancs polygon 
> cannot have holes, by the way. I think that what you could do is to choose 
> the locations= values you need within the polygon grid first, and not use 
> borders=. I've tried to replicate your situation, and have made a 
> shapefile with rectangular grid cells with uneven edges and a hole. Then I 
> did:
> 
> library(sp)
> library(maptools)
> shapefile <- read.shape("<your file>.shp")
> x <- as.SpatialRings.Shapes(shapefile$Shapes, IDs=1:length(shapefile$Shapes))
> xp <- SpatialPoints(cbind(runif(50, 0, 3), runif(50, 0, 3)))
> 
> # here you do xp <- SpatialPoints(<your location points matrix>)
> 
> res <- overlay(x, xp)
> 
> # overlay() sets the value of res to the ID of the polygon the point falls 
> # in, or to NA if the point does not fall in any polygon.
> 
> plot(x, col="grey")
> points(xp[is.na(res)], col="red")
> points(xp[!is.na(res)], col="blue")
> 
> and you can choose the points for the prediction locations as:
> 
> <my prediction locations> <- xp[!is.na(res)]

Further to my own reply, the last line was not right:

> xx <- xp[!is.na(res)]
> class(xx)
[1] "SpatialPoints"
attr(,"package")
[1] "sp"
> class(coordinates(xx))
[1] "matrix"

so what you need is:

<my prediction locations> <- coordinates(xp[!is.na(res)])

which returns just the matrix of point coordinates without the other 
information.

This now works, using check.locations() called by krige.conv():

> check.locations(xx)
[1] "no"
> check.locations(coordinates(xx))
      coords.x1 coords.x2
 [1,] 2.8883985 1.5438226
 [2,] 2.1517145 1.4901741
 [3,] 0.8694240 1.9049170

...


> 
> Please check on the plot to see if overlay divides up the points correctly 
> first, just to be sure. If this works, it is will be a nice example of how 
> the coordination of spatial data functions in the "sp" package can make 
> things easier (if not, any reports on difficulties with "sp" are always 
> welcome).
> 
> Roger
> 
> >  But i'm not sure that it's the good polylist object i need because it's
> > a polygon grid in fact.
> 
> > Thank you for your help
> > Chloé
> > -----Message d'origine-----
> > De : Roger Bivand [mailto:Roger.Bivand at nhh.no] 
> > Envoyé : jeudi 19 mai 2005 19:07
> > À : Chloe ARCHINARD
> > Cc : r-sig-geo at stat.math.ethz.ch
> > Objet : Re: [R-sig-Geo] kriging inside borders
> > 
> > On Thu, 19 May 2005, Chloe ARCHINARD wrote:
> > 
> > 
> > > Hello, I'm using geoR package to compute kriging with the function
> > 
> > > krige.conv() and I used maptools (read.shape () + Map2poly () )  and
> > 
> > > splancs (borders) packages to include a object of class polylist for the
> > 
> > > attribute borders.
> > 
> > 
> > >  But I've got this message of error when 'borders' is in krige.conv () :  
> > 
> > > krige.conv: there are no prediction locations inside the borders.  My
> > 
> > > polylist is not good? The coordinates seam to be OK when I draw a map
> > 
> > > with this polylist object. And if I use 'borders' in the contour
> > 
> > > function : message error is :
> > 
> > 
> > > Error in as.data.frame.default(borders) : can't coerce class "polylist"
> > 
> > > into a data.frame
> > 
> > 
> > Firstly, it does help if the exact commands are given, rather than your 
> > 
> > choice of them. I think the main issue here is that a polylist object is a 
> > 
> > list of matrices with attributes, not a matrix. If your polylist object 
> > 
> > is:
> > 
> > 
> > x <- Map2poly(read.shape("my.shp"))
> > 
> > 
> > what is length(x)? If it is 1, you may be able to say x[[1]] to extract 
> > 
> > the first matrix in the list, but even then it may not work. If you can 
> > 
> > show us the commands you are giving, and the length of your polylist 
> > 
> > object, it should be possible to find something that works. 
> > 
> > 
> > 
> > > If someone could help me, thanks a lot!
> > 
> > 
> > > Chloé Archinard
> > 
> > > 
> > 
> > > 
> > 
> > 
> > -- 
> > 
> > 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
> > 
> > 
> > 
> > -- 
> > passerelle antivirus du campus CNRS de Montpellier
> > --
> > 
> > 
> > 
> > 
> 
> 

-- 
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