I'm trying to make IDW in R. In my tutorial, I found this approach which is easily understandable, expect of one point - why they used n = 6000 ?
How did they established the value to create a grid by spsample()??
library(maptools)
library(GISTools)
library(gstat)
require(deldir)
require(sp)
voronoipolygons = function(layer) {
crds = layer@coords
z = deldir(crds[,1], crds[,2])
w = tile.list(z)
polys = vector(mode='list', length=length(w))
for (i in seq(along=polys)) {
pcrds = cbind(w[[i]]$x, w[[i]]$y)
pcrds = rbind(pcrds, pcrds[1,])
polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
}
SP = SpatialPolygons(polys)
voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
y=crds[,2],
row.names=sapply(slot(SP, 'polygons'),
function(x) slot(x, 'ID'))))
return(voronoi)
}
library(gstat)
library(maptools)
data(fulmar)
fulmar.spdf <- SpatialPointsDataFrame(cbind(fulmar$x, fulmar$y), fulmar)
fulmar.spdf <- fulmar.spdf[fulmar.spdf$year == 1999, ]
fulmar.voro <-voronoipolygons(fulmar.spdf)
s.grid <- spsample(fulmar.voro, type = "regular", n = 6000) ### WHY ?? !!!
idw.est <- gstat::idw(fulmar ~ 1, fulmar.spdf, newdata=s.grid, idp = 1.0)
# extract the unique xand y locations in the grid
ux<-unique(coordinates(idw.est)[,1])
uy<-unique(coordinates(idw.est)[,2])
# extract the predicted values and format var1.pred into a matrix of gridded values
predmat <- matrix(idw.est$var1.pred, length(ux), length(uy))
When I tried to rerun the code for my data, I found that the length(ux) and length(uy) should be multiplication of length(idw.est) to correctly fill my predmat matrix otherwise I have a warning msg :
Warning message:
In matrix(idw.est$var1.pred, length(ux), length(uy)) :
data length [1003] is not a sub-multiple or multiple of the number of rows [49]
and my resulting interpolation is not correct.
Can you understand and explain me why did they used n = 6000?
Was it just by trial and error, or what trick can I use to set this number for any area?