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()??


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], 
                                                         row.names=sapply(slot(SP, 'polygons'), 
                                                                          function(x) slot(x, 'ID'))))


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

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


Ignore this voronoi tessellation nonsense, it is just confusing the matter. The easiest way is to settle on a grain size (resolution), create an extent object using "raster::extent", create a raster object using the extent object and the defined resolution in the raster::raster function, assign incremental values and then coerce to a SpatialPointsDataFrame object using "as". This way you are not figuring out n but rather a supported grain. It would look something like this:


s = 360 # resolution
e <- raster::extent(172321.7, 203004.2, 72720.54, 110285.1)
r <- raster::raster(e, res=c(s,s))
  r[] <- 1:raster::ncell(r)

r <- as(r, "SpatialPointsDataFrame")

plot(r, pch=20, cex=0.5)

Note; you can pull your extent from any spatial (raster, sp) object eg., e <- extent(mypoints), it does not have to be explicitly defined with input coordinates.


The number of cells (n) used in a prediction grid for interpolation, is a consequence of defining its extent and resolution.

The extent is usually the area in which sample points were collected, so to avoid extrapolation; and the resolution can vary according to many things, for example: to what is being studied or due to the amount of sample points available within the area of interest.

Therefore, it could have been a different value from n = 6000, but I guess in that exercise they wanted to simulate a grid with greater number of slots than the available data (6000 > 729).
729 observations comes from the sample points in fulmar.spdf (nrow(fulmar.spdf)).

See the following post for a more didactic example about creating prediction grids:

Plotting map resulted from kriging in R

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.