Skip to content
Prev 425 / 29559 Next

kriging inside borders

On Fri, 20 May 2005, Chloe ARCHINARD wrote:

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

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