This short note shows how to create a regular grid onto which kriging or another prediction method can be applied.

We use the Meuse dataset as an example of the study area we want to predict over. We have some points in this area:

library(sp); library(gstat)
data(meuse)  # an example dataset supplied with the sp package
names(meuse)

##  [1] "x"       "y"       "cadmium" "copper"  "lead"    "zinc"    "elev"   
##  [8] "dist"    "om"      "ffreq"   "soil"    "lime"    "landuse" "dist.m"

coordinates(meuse) <- c("x", "y") # convert to a SpatialPointsDataFrame
summary(meuse)

## Object of class SpatialPointsDataFrame
## Coordinates:
##      min    max
## x 178605 181390
## y 329714 333611
## Is projected: NA 
## proj4string : [NA]
## Number of points: 155
## Data attributes:
##     cadmium           copper            lead            zinc       
##  Min.   : 0.200   Min.   : 14.00   Min.   : 37.0   Min.   : 113.0  
##  1st Qu.: 0.800   1st Qu.: 23.00   1st Qu.: 72.5   1st Qu.: 198.0  
##  Median : 2.100   Median : 31.00   Median :123.0   Median : 326.0  
##  Mean   : 3.246   Mean   : 40.32   Mean   :153.4   Mean   : 469.7  
##  3rd Qu.: 3.850   3rd Qu.: 49.50   3rd Qu.:207.0   3rd Qu.: 674.5  
##  Max.   :18.100   Max.   :128.00   Max.   :654.0   Max.   :1839.0  
##                                                                    
##       elev             dist               om         ffreq  soil   lime   
##  Min.   : 5.180   Min.   :0.00000   Min.   : 1.000   1:84   1:97   0:111  
##  1st Qu.: 7.546   1st Qu.:0.07569   1st Qu.: 5.300   2:48   2:46   1: 44  
##  Median : 8.180   Median :0.21184   Median : 6.900   3:23   3:12          
##  Mean   : 8.165   Mean   :0.24002   Mean   : 7.478                        
##  3rd Qu.: 8.955   3rd Qu.:0.36407   3rd Qu.: 9.000                        
##  Max.   :10.520   Max.   :0.88039   Max.   :17.000                        
##                                     NA's   :2                             
##     landuse       dist.m      
##  W      :50   Min.   :  10.0  
##  Ah     :39   1st Qu.:  80.0  
##  Am     :22   Median : 270.0  
##  Fw     :10   Mean   : 290.3  
##  Ab     : 8   3rd Qu.: 450.0  
##  (Other):25   Max.   :1000.0  
##  NA's   : 1

bbox(meuse)

##      min    max
## x 178605 181390
## y 329714 333611

We can define any extent for a grid; here we choose to cover the bounding box, and expand the grid to match the chosen resolution.

bbox(meuse)

##      min    max
## x 178605 181390
## y 329714 333611

res <- 120 # resolution
# round extremes to resolution
(x.min <- bbox(meuse)[1,1]%/%res*res)

## [1] 178560

(x.max <- (bbox(meuse)[1,2]+res)%/%res*res) # make sure it is outside the bbox

## [1] 181440

(y.min <- bbox(meuse)[2,1]%/%res*res)

## [1] 329640

(y.max <- (bbox(meuse)[2,2]+res)%/%res*res) # make sure it is outside the bbox

## [1] 333720

# grid of coordinates
grid <- expand.grid(x = seq(x.min, x.max, by=res),
              y = seq(y.min, y.max, by=res))
class(grid)

## [1] "data.frame"

coordinates(grid) <- c("x", "y")
class(grid)

## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"

gridded(grid) <- T; fullgrid(grid) <- T
class(grid)

## [1] "SpatialGrid"
## attr(,"package")
## [1] "sp"

Now we can interpolate over the grid, here using IDW just to show a result:

kres <- idw(zinc ~ 1, loc=meuse, newdata=grid)

## [inverse distance weighted interpolation]

spplot(kres, zcol="var1.pred")