4

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?

Andre Silva
  • 10,259
  • 12
  • 54
  • 106
maycca
  • 3,376
  • 4
  • 30
  • 59

2 Answers2

2

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:

library(raster)
library(sp)

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)
  class(r)

r <- as(r, "SpatialPointsDataFrame")
  class(r)
  dim(r)

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.

Jeffrey Evans
  • 31,750
  • 2
  • 47
  • 94
1

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

Andre Silva
  • 10,259
  • 12
  • 54
  • 106